c@120: /* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
c@120: /*
c@120:     Constant-Q library
c@120:     Copyright (c) 2013-2014 Queen Mary, University of London
c@120: 
c@120:     Permission is hereby granted, free of charge, to any person
c@120:     obtaining a copy of this software and associated documentation
c@120:     files (the "Software"), to deal in the Software without
c@120:     restriction, including without limitation the rights to use, copy,
c@120:     modify, merge, publish, distribute, sublicense, and/or sell copies
c@120:     of the Software, and to permit persons to whom the Software is
c@120:     furnished to do so, subject to the following conditions:
c@120: 
c@120:     The above copyright notice and this permission notice shall be
c@120:     included in all copies or substantial portions of the Software.
c@120: 
c@120:     THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
c@120:     EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
c@120:     MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
c@120:     NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY
c@120:     CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF
c@120:     CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
c@120:     WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
c@120: 
c@120:     Except as contained in this notice, the names of the Centre for
c@120:     Digital Music; Queen Mary, University of London; and Chris Cannam
c@120:     shall not be used in advertising or otherwise to promote the sale,
c@120:     use or other dealings in this Software without prior written
c@120:     authorization.
c@120: */
c@120: 
c@120: #include "FFT.h"
c@120: 
c@122: #include "MathUtilities.h"
c@120: 
c@120: #include "kiss_fft.h"
c@120: #include "kiss_fftr.h"
c@120: 
c@120: #include <cmath>
c@120: 
c@120: #include <iostream>
c@120: 
c@120: #include <stdexcept>
c@120: 
c@120: class FFT::D
c@120: {
c@120: public:
c@120:     D(int n) : m_n(n) {
c@120:         m_planf = kiss_fft_alloc(m_n, 0, NULL, NULL);
c@120:         m_plani = kiss_fft_alloc(m_n, 1, NULL, NULL);
c@120:         m_kin = new kiss_fft_cpx[m_n];
c@120:         m_kout = new kiss_fft_cpx[m_n];
c@120:     }
c@120: 
c@120:     ~D() {
c@120:         kiss_fft_free(m_planf);
c@120:         kiss_fft_free(m_plani);
c@120:         delete[] m_kin;
c@120:         delete[] m_kout;
c@120:     }
c@120: 
c@120:     void process(bool inverse,
c@120:                  const double *ri,
c@120:                  const double *ii,
c@120:                  double *ro,
c@120:                  double *io) {
c@120: 
c@120:         for (int i = 0; i < m_n; ++i) {
c@120:             m_kin[i].r = ri[i];
c@120:             m_kin[i].i = (ii ? ii[i] : 0.0);
c@120:         }
c@120: 
c@120:         if (!inverse) {
c@120: 
c@120:             kiss_fft(m_planf, m_kin, m_kout);
c@120: 
c@120:             for (int i = 0; i < m_n; ++i) {
c@120:                 ro[i] = m_kout[i].r;
c@120:                 io[i] = m_kout[i].i;
c@120:             }
c@120: 
c@120:         } else {
c@120: 
c@120:             kiss_fft(m_plani, m_kin, m_kout);
c@120: 
c@120:             double scale = 1.0 / m_n;
c@120: 
c@120:             for (int i = 0; i < m_n; ++i) {
c@120:                 ro[i] = m_kout[i].r * scale;
c@120:                 io[i] = m_kout[i].i * scale;
c@120:             }
c@120:         }
c@120:     }
c@120:     
c@120: private:
c@120:     int m_n;
c@120:     kiss_fft_cfg m_planf;
c@120:     kiss_fft_cfg m_plani;
c@120:     kiss_fft_cpx *m_kin;
c@120:     kiss_fft_cpx *m_kout;
c@120: };        
c@120: 
c@120: FFT::FFT(int n) :
c@120:     m_d(new D(n))
c@120: {
c@120: }
c@120: 
c@120: FFT::~FFT()
c@120: {
c@120:     delete m_d;
c@120: }
c@120: 
c@120: void
c@120: FFT::process(bool inverse,
c@120:              const double *p_lpRealIn, const double *p_lpImagIn,
c@120:              double *p_lpRealOut, double *p_lpImagOut)
c@120: {
c@120:     m_d->process(inverse,
c@120:                  p_lpRealIn, p_lpImagIn,
c@120:                  p_lpRealOut, p_lpImagOut);
c@120: }
c@120:     
c@120: class FFTReal::D
c@120: {
c@120: public:
c@120:     D(int n) : m_n(n) {
c@120:         if (n % 2) {
c@120:             throw std::invalid_argument
c@120:                 ("nsamples must be even in FFTReal constructor");
c@120:         }
c@120:         m_planf = kiss_fftr_alloc(m_n, 0, NULL, NULL);
c@120:         m_plani = kiss_fftr_alloc(m_n, 1, NULL, NULL);
c@120:         m_c = new kiss_fft_cpx[m_n];
c@120:     }
c@120: 
c@120:     ~D() {
c@120:         kiss_fftr_free(m_planf);
c@120:         kiss_fftr_free(m_plani);
c@120:         delete[] m_c;
c@120:     }
c@120: 
c@120:     void forward(const double *ri, double *ro, double *io) {
c@120: 
c@120:         kiss_fftr(m_planf, ri, m_c);
c@120: 
c@120:         for (int i = 0; i <= m_n/2; ++i) {
c@120:             ro[i] = m_c[i].r;
c@120:             io[i] = m_c[i].i;
c@120:         }
c@120: 
c@120:         for (int i = 0; i + 1 < m_n/2; ++i) {
c@120:             ro[m_n - i - 1] =  ro[i + 1];
c@120:             io[m_n - i - 1] = -io[i + 1];
c@120:         }
c@120:     }
c@120: 
c@120:     void forwardMagnitude(const double *ri, double *mo) {
c@120: 
c@120:         double *io = new double[m_n];
c@120: 
c@120:         forward(ri, mo, io);
c@120: 
c@120:         for (int i = 0; i < m_n; ++i) {
c@120:             mo[i] = sqrt(mo[i] * mo[i] + io[i] * io[i]);
c@120:         }
c@120: 
c@120:         delete[] io;
c@120:     }
c@120: 
c@120:     void inverse(const double *ri, const double *ii, double *ro) {
c@120: 
c@120:         // kiss_fftr.h says
c@120:         // "input freqdata has nfft/2+1 complex points"
c@120: 
c@120:         for (int i = 0; i < m_n/2 + 1; ++i) {
c@120:             m_c[i].r = ri[i];
c@120:             m_c[i].i = ii[i];
c@120:         }
c@120:         
c@120:         kiss_fftri(m_plani, m_c, ro);
c@120: 
c@120:         double scale = 1.0 / m_n;
c@120: 
c@120:         for (int i = 0; i < m_n; ++i) {
c@120:             ro[i] *= scale;
c@120:         }
c@120:     }
c@120: 
c@120: private:
c@120:     int m_n;
c@120:     kiss_fftr_cfg m_planf;
c@120:     kiss_fftr_cfg m_plani;
c@120:     kiss_fft_cpx *m_c;
c@120: };
c@120: 
c@120: FFTReal::FFTReal(int n) :
c@120:     m_d(new D(n)) 
c@120: {
c@120: }
c@120: 
c@120: FFTReal::~FFTReal()
c@120: {
c@120:     delete m_d;
c@120: }
c@120: 
c@120: void
c@120: FFTReal::forward(const double *ri, double *ro, double *io)
c@120: {
c@120:     m_d->forward(ri, ro, io);
c@120: }
c@120: 
c@120: void
c@120: FFTReal::forwardMagnitude(const double *ri, double *mo)
c@120: {
c@120:     m_d->forwardMagnitude(ri, mo);
c@120: }
c@120: 
c@120: void
c@120: FFTReal::inverse(const double *ri, const double *ii, double *ro)
c@120: {
c@120:     m_d->inverse(ri, ii, ro);
c@120: }
c@120: 
c@120: 
c@120: