Chris@226: /* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
Chris@226: 
Chris@226: /*
Chris@226:     Sonic Visualiser
Chris@226:     An audio file viewer and annotation editor.
Chris@226:     Centre for Digital Music, Queen Mary, University of London.
Chris@226:     This file copyright 2006 Chris Cannam and QMUL.
Chris@226:     
Chris@226:     This program is free software; you can redistribute it and/or
Chris@226:     modify it under the terms of the GNU General Public License as
Chris@226:     published by the Free Software Foundation; either version 2 of the
Chris@226:     License, or (at your option) any later version.  See the file
Chris@226:     COPYING included with this distribution for more information.
Chris@226: */
Chris@226: 
Chris@226: #ifndef _FFT_API_H_
Chris@226: #define _FFT_API_H_
Chris@226: 
Chris@226: #ifdef HAVE_FFTW3F
Chris@226: 
Chris@226: #include <fftw3.h>
Chris@226: 
Chris@226: #define fftf_complex fftwf_complex
Chris@226: #define fftf_malloc fftwf_malloc
Chris@226: #define fftf_free fftwf_free
Chris@226: #define fftf_plan fftwf_plan
Chris@226: #define fftf_plan_dft_r2c_1d fftwf_plan_dft_r2c_1d
Chris@226: #define fftf_plan_dft_c2r_1d fftwf_plan_dft_c2r_1d
Chris@226: #define fftf_execute fftwf_execute
Chris@226: #define fftf_destroy_plan fftwf_destroy_plan
Chris@226: 
Chris@226: #else
Chris@226: 
Chris@226: // Provide a fallback FFT implementation if FFTW3f is not available.
Chris@226: 
Chris@226: typedef float fftf_complex[2];
Chris@226: #define fftf_malloc malloc
Chris@226: #define fftf_free free
Chris@226: 
Chris@226: struct fftf_plan_;
Chris@226: typedef fftf_plan_ *fftf_plan;
Chris@226: 
Chris@226: fftf_plan fftf_plan_dft_r2c_1d(int n, float *in, fftf_complex *out, unsigned);
Chris@226: fftf_plan fftf_plan_dft_c2r_1d(int n, fftf_complex *in, float *out, unsigned);
Chris@226: void fftf_execute(const fftf_plan p);
Chris@226: void fftf_destroy_plan(fftf_plan p);
Chris@226: 
Chris@226: #define FFTW_ESTIMATE 0
Chris@293: #define FFTW_MEASURE 0
Chris@226: 
Chris@226: #endif
Chris@226: 
Chris@1091: #include <vector>
Chris@1091: #include <complex>
Chris@1131: #include <mutex>
Chris@1091: 
Chris@1091: class FFTForward // with fft shift but not window
Chris@1091: {
Chris@1131:     static std::mutex m_mutex;
Chris@1131:     
Chris@1091: public:
Chris@1091:     FFTForward(int size) :
Chris@1131:         m_size(size)
Chris@1131:     {
Chris@1131:         std::lock_guard<std::mutex> lock(m_mutex);
Chris@1131:         m_input = (float *)fftf_malloc(size * sizeof(float));
Chris@1131:         m_output = (fftf_complex *)fftf_malloc((size/2 + 1) * sizeof(fftf_complex));
Chris@1131:         m_plan = fftf_plan_dft_r2c_1d(size, m_input, m_output, FFTW_ESTIMATE);
Chris@1131:     }
Chris@1091: 
Chris@1091:     ~FFTForward() {
Chris@1131:         std::lock_guard<std::mutex> lock(m_mutex);
Chris@1091:         fftf_destroy_plan(m_plan);
Chris@1091:         fftf_free(m_input);
Chris@1091:         fftf_free(m_output);
Chris@1091:     }
Chris@1091: 
Chris@1091:     std::vector<std::complex<float> > process(std::vector<float> in) const {
Chris@1091:         const int hs = m_size/2;
Chris@1091:         for (int i = 0; i < hs; ++i) {
Chris@1091:             m_input[i] = in[i + hs];
Chris@1095:         }
Chris@1095:         for (int i = 0; i < hs; ++i) {
Chris@1091:             m_input[i + hs] = in[i];
Chris@1091:         }
Chris@1091:         fftf_execute(m_plan);
Chris@1091:         std::vector<std::complex<float> > result;
Chris@1091:         result.reserve(hs + 1);
Chris@1091:         for (int i = 0; i <= hs; ++i) {
Chris@1093:             result.push_back({ m_output[i][0], m_output[i][1] });
Chris@1091:         }
Chris@1091:         return result;
Chris@1091:     }
Chris@1091: 
Chris@1091: private:
Chris@1091:     int m_size;
Chris@1091:     float *m_input;
Chris@1091:     fftf_complex *m_output;
Chris@1091:     fftf_plan m_plan;
Chris@1091: };
Chris@1091: 
Chris@226: #endif
Chris@226: