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