cannam@0: /* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
cannam@0: 
cannam@0: /*
cannam@0:     QM DSP Library
cannam@0: 
cannam@0:     Centre for Digital Music, Queen Mary, University of London.
cannam@0:     This file is based on Don Cross's public domain FFT implementation.
cannam@0: */
cannam@0: 
cannam@0: #include "FFT.h"
cannam@55: 
cannam@55: #include "maths/MathUtilities.h"
cannam@55: 
Chris@129: #include "kiss_fft.h"
Chris@129: #include "kiss_fftr.h"
Chris@129: 
cannam@0: #include <cmath>
cannam@0: 
cannam@55: #include <iostream>
cannam@55: 
Chris@129: #include <stdexcept>
Chris@129: 
Chris@129: class FFT::D
Chris@129: {
Chris@129: public:
Chris@129:     D(int n) : m_n(n) {
Chris@129:         m_planf = kiss_fft_alloc(m_n, 0, NULL, NULL);
Chris@129:         m_plani = kiss_fft_alloc(m_n, 1, NULL, NULL);
Chris@129:         m_kin = new kiss_fft_cpx[m_n];
Chris@129:         m_kout = new kiss_fft_cpx[m_n];
Chris@129:     }
Chris@129: 
Chris@129:     ~D() {
Chris@129:         kiss_fft_free(m_planf);
Chris@129:         kiss_fft_free(m_plani);
Chris@129:         delete[] m_kin;
Chris@129:         delete[] m_kout;
Chris@129:     }
Chris@129: 
Chris@129:     void process(bool inverse,
Chris@129:                  const double *ri,
Chris@129:                  const double *ii,
Chris@129:                  double *ro,
Chris@129:                  double *io) {
Chris@129: 
Chris@129:         for (int i = 0; i < m_n; ++i) {
Chris@129:             m_kin[i].r = ri[i];
Chris@129:             m_kin[i].i = (ii ? ii[i] : 0.0);
Chris@129:         }
Chris@129: 
Chris@129:         if (!inverse) {
Chris@129: 
Chris@129:             kiss_fft(m_planf, m_kin, m_kout);
Chris@129: 
Chris@129:             for (int i = 0; i < m_n; ++i) {
Chris@129:                 ro[i] = m_kout[i].r;
Chris@129:                 io[i] = m_kout[i].i;
Chris@129:             }
Chris@129: 
Chris@129:         } else {
Chris@129: 
Chris@129:             kiss_fft(m_plani, m_kin, m_kout);
Chris@129: 
Chris@129:             double scale = 1.0 / m_n;
Chris@129: 
Chris@129:             for (int i = 0; i < m_n; ++i) {
Chris@129:                 ro[i] = m_kout[i].r * scale;
Chris@129:                 io[i] = m_kout[i].i * scale;
Chris@129:             }
Chris@129:         }
Chris@129:     }
Chris@129:     
Chris@129: private:
Chris@129:     int m_n;
Chris@129:     kiss_fft_cfg m_planf;
Chris@129:     kiss_fft_cfg m_plani;
Chris@129:     kiss_fft_cpx *m_kin;
Chris@129:     kiss_fft_cpx *m_kout;
Chris@129: };        
Chris@129: 
Chris@114: FFT::FFT(int n) :
Chris@129:     m_d(new D(n))
cannam@0: {
cannam@0: }
cannam@0: 
cannam@0: FFT::~FFT()
cannam@0: {
Chris@129:     delete m_d;
cannam@0: }
cannam@0: 
Chris@129: void
Chris@129: FFT::process(bool inverse,
Chris@129:              const double *p_lpRealIn, const double *p_lpImagIn,
Chris@129:              double *p_lpRealOut, double *p_lpImagOut)
Chris@129: {
Chris@129:     m_d->process(inverse,
Chris@129:                  p_lpRealIn, p_lpImagIn,
Chris@129:                  p_lpRealOut, p_lpImagOut);
Chris@129: }
Chris@129:     
Chris@129: class FFTReal::D
Chris@129: {
Chris@129: public:
Chris@129:     D(int n) : m_n(n) {
Chris@129:         if (n % 2) {
Chris@129:             throw std::invalid_argument
Chris@129:                 ("nsamples must be even in FFTReal constructor");
Chris@129:         }
Chris@129:         m_planf = kiss_fftr_alloc(m_n, 0, NULL, NULL);
Chris@129:         m_plani = kiss_fftr_alloc(m_n, 1, NULL, NULL);
Chris@129:         m_c = new kiss_fft_cpx[m_n];
Chris@129:     }
Chris@129: 
Chris@129:     ~D() {
Chris@129:         kiss_fftr_free(m_planf);
Chris@129:         kiss_fftr_free(m_plani);
Chris@129:         delete[] m_c;
Chris@129:     }
Chris@129: 
Chris@129:     void forward(const double *ri, double *ro, double *io) {
Chris@129: 
Chris@129:         kiss_fftr(m_planf, ri, m_c);
Chris@129: 
Chris@129:         for (int i = 0; i <= m_n/2; ++i) {
Chris@129:             ro[i] = m_c[i].r;
Chris@129:             io[i] = m_c[i].i;
Chris@129:         }
Chris@129: 
Chris@129:         for (int i = 0; i + 1 < m_n/2; ++i) {
Chris@129:             ro[m_n - i - 1] =  ro[i + 1];
Chris@129:             io[m_n - i - 1] = -io[i + 1];
Chris@129:         }
Chris@129:     }
Chris@129: 
Chris@129:     void inverse(const double *ri, const double *ii, double *ro) {
Chris@129: 
Chris@129:         for (int i = 0; i < m_n; ++i) {
Chris@129:             m_c[i].r = ri[i];
Chris@129:             m_c[i].i = ii[i];
Chris@129:         }
Chris@129:         
Chris@129:         kiss_fftri(m_plani, m_c, ro);
Chris@129: 
Chris@129:         double scale = 1.0 / m_n;
Chris@129: 
Chris@129:         for (int i = 0; i < m_n; ++i) {
Chris@129:             ro[i] *= scale;
Chris@129:         }
Chris@129:     }
Chris@129: 
Chris@129: private:
Chris@129:     int m_n;
Chris@129:     kiss_fftr_cfg m_planf;
Chris@129:     kiss_fftr_cfg m_plani;
Chris@129:     kiss_fft_cpx *m_c;
Chris@129: };
Chris@129: 
Chris@114: FFTReal::FFTReal(int n) :
Chris@129:     m_d(new D(n)) 
cannam@0: {
cannam@64: }
cannam@0: 
cannam@64: FFTReal::~FFTReal()
cannam@64: {
Chris@129:     delete m_d;
cannam@64: }
cannam@64: 
cannam@64: void
Chris@129: FFTReal::forward(const double *ri, double *ro, double *io)
cannam@64: {
Chris@129:     m_d->forward(ri, ro, io);
cannam@64: }
cannam@64: 
Chris@114: void
Chris@129: FFTReal::inverse(const double *ri, const double *ii, double *ro)
Chris@114: {
Chris@129:     m_d->inverse(ri, ii, ro);
Chris@114: }
Chris@114: 
cannam@64: 
Chris@129: