Chris@191: /* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
Chris@191: 
Chris@191: /*
Chris@191:     QM DSP Library
Chris@191: 
Chris@191:     Centre for Digital Music, Queen Mary, University of London.
Chris@191:     This file by Chris Cannam.
Chris@191: 
Chris@191:     This program is free software; you can redistribute it and/or
Chris@191:     modify it under the terms of the GNU General Public License as
Chris@191:     published by the Free Software Foundation; either version 2 of the
Chris@191:     License, or (at your option) any later version.  See the file
Chris@191:     COPYING included with this distribution for more information.
Chris@191: */
Chris@191: 
Chris@191: #include "DCT.h"
Chris@191: 
Chris@191: #include <cmath>
Chris@191: 
Chris@191: DCT::DCT(int n) :
Chris@191:     m_n(n),
Chris@191:     m_scaled(n, 0.0),
Chris@191:     m_time2(n*4, 0.0),
Chris@191:     m_freq2r(n*4, 0.0),
Chris@191:     m_freq2i(n*4, 0.0),
Chris@191:     m_fft(n*4)
Chris@191: {
Chris@191:     m_scale = m_n * sqrt(2.0 / m_n);
Chris@191: }
Chris@191: 
Chris@191: DCT::~DCT()
Chris@191: {
Chris@191: }
Chris@191: 
Chris@191: void
Chris@191: DCT::forward(const double *in, double *out)
Chris@191: {
Chris@191:     for (int i = 0; i < m_n; ++i) {
Chris@191: 	m_time2[i*2 + 1] = in[i];
Chris@191: 	m_time2[m_n*4 - i*2 - 1] = in[i];
Chris@191:     }
Chris@191: 
Chris@191:     m_fft.forward(m_time2.data(), m_freq2r.data(), m_freq2i.data());
Chris@191: 
Chris@191:     for (int i = 0; i < m_n; ++i) {
Chris@191: 	out[i] = m_freq2r[i];
Chris@191:     }
Chris@191: }
Chris@191: 
Chris@191: void
Chris@191: DCT::forwardUnitary(const double *in, double *out)
Chris@191: {
Chris@191:     forward(in, out);
Chris@191:     for (int i = 0; i < m_n; ++i) {
Chris@191: 	out[i] /= m_scale;
Chris@191:     }
Chris@191:     out[0] /= sqrt(2.0);
Chris@191: }
Chris@191: 
Chris@191: void
Chris@191: DCT::inverse(const double *in, double *out)
Chris@191: {
Chris@191:     for (int i = 0; i < m_n; ++i) {
Chris@191: 	m_freq2r[i] = in[i];
Chris@191:     }
Chris@191:     for (int i = 0; i < m_n; ++i) {
Chris@191: 	m_freq2r[m_n*2 - i] = -in[i];
Chris@191:     }
Chris@191:     m_freq2r[m_n] = 0.0;
Chris@191: 
Chris@191:     for (int i = 0; i <= m_n*2; ++i) {
Chris@191: 	m_freq2i[i] = 0.0;
Chris@191:     }
Chris@191:     
Chris@191:     m_fft.inverse(m_freq2r.data(), m_freq2i.data(), m_time2.data());
Chris@191: 
Chris@191:     for (int i = 0; i < m_n; ++i) {
Chris@191: 	out[i] = m_time2[i*2 + 1];
Chris@191:     }
Chris@191: }
Chris@191: 
Chris@191: void
Chris@191: DCT::inverseUnitary(const double *in, double *out)
Chris@191: {
Chris@191:     for (int i = 0; i < m_n; ++i) {
Chris@191: 	m_scaled[i] = in[i] * m_scale;
Chris@191:     }
Chris@191:     m_scaled[0] *= sqrt(2.0);
Chris@191:     inverse(m_scaled.data(), out);
Chris@191: }
Chris@191: