svcore  1.9
FFTapi.cpp
Go to the documentation of this file.
00001 /* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
00002 
00003 /*
00004     Sonic Visualiser
00005     An audio file viewer and annotation editor.
00006     Centre for Digital Music, Queen Mary, University of London.
00007     This file copyright 2006 Chris Cannam and QMUL.
00008     FFT code from Don Cross's public domain FFT implementation.
00009     
00010     This program is free software; you can redistribute it and/or
00011     modify it under the terms of the GNU General Public License as
00012     published by the Free Software Foundation; either version 2 of the
00013     License, or (at your option) any later version.  See the file
00014     COPYING included with this distribution for more information.
00015 */
00016 
00017 #include "FFTapi.h"
00018 
00019 #ifndef HAVE_FFTW3F
00020 
00021 #include <cmath>
00022 #include <iostream>
00023 
00024 void
00025 fft(unsigned int n, bool inverse, double *ri, double *ii, double *ro, double *io)
00026 {
00027     if (!ri || !ro || !io) return;
00028 
00029     unsigned int bits;
00030     unsigned int i, j, k, m;
00031     unsigned int blockSize, blockEnd;
00032 
00033     double tr, ti;
00034 
00035     if (n < 2) return;
00036     if (n & (n-1)) return;
00037 
00038     double angle = 2.0 * M_PI;
00039     if (inverse) angle = -angle;
00040 
00041     for (i = 0; ; ++i) {
00042         if (n & (1 << i)) {
00043             bits = i;
00044             break;
00045         }
00046     }
00047 
00048     int *table = new int[n];
00049 
00050     for (i = 0; i < n; ++i) {
00051         
00052         m = i;
00053         
00054         for (j = k = 0; j < bits; ++j) {
00055             k = (k << 1) | (m & 1);
00056             m >>= 1;
00057         }
00058         
00059         table[i] = k;
00060     }
00061 
00062     if (ii) {
00063         for (i = 0; i < n; ++i) {
00064             ro[table[i]] = ri[i];
00065             io[table[i]] = ii[i];
00066         }
00067     } else {
00068         for (i = 0; i < n; ++i) {
00069             ro[table[i]] = ri[i];
00070             io[table[i]] = 0.0;
00071         }
00072     }
00073 
00074     blockEnd = 1;
00075 
00076     for (blockSize = 2; blockSize <= n; blockSize <<= 1) {
00077 
00078         double delta = angle / (double)blockSize;
00079         double sm2 = -sin(-2 * delta);
00080         double sm1 = -sin(-delta);
00081         double cm2 = cos(-2 * delta);
00082         double cm1 = cos(-delta);
00083         double w = 2 * cm1;
00084         double ar[3], ai[3];
00085 
00086         for (i = 0; i < n; i += blockSize) {
00087 
00088             ar[2] = cm2;
00089             ar[1] = cm1;
00090 
00091             ai[2] = sm2;
00092             ai[1] = sm1;
00093 
00094             for (j = i, m = 0; m < blockEnd; j++, m++) {
00095 
00096                 ar[0] = w * ar[1] - ar[2];
00097                 ar[2] = ar[1];
00098                 ar[1] = ar[0];
00099 
00100                 ai[0] = w * ai[1] - ai[2];
00101                 ai[2] = ai[1];
00102                 ai[1] = ai[0];
00103 
00104                 k = j + blockEnd;
00105                 tr = ar[0] * ro[k] - ai[0] * io[k];
00106                 ti = ar[0] * io[k] + ai[0] * ro[k];
00107 
00108                 ro[k] = ro[j] - tr;
00109                 io[k] = io[j] - ti;
00110 
00111                 ro[j] += tr;
00112                 io[j] += ti;
00113             }
00114         }
00115 
00116         blockEnd = blockSize;
00117     }
00118 
00119 /* fftw doesn't normalise, so nor will we
00120 
00121     if (inverse) {
00122 
00123         double denom = (double)n;
00124 
00125         for (i = 0; i < n; i++) {
00126             ro[i] /= denom;
00127             io[i] /= denom;
00128         }
00129     }
00130 */
00131     delete[] table;
00132 }
00133 
00134 struct fftf_plan_ {
00135     int size;
00136     int inverse;
00137     float *real;
00138     fftf_complex *cplx;
00139 };
00140 
00141 fftf_plan
00142 fftf_plan_dft_r2c_1d(int n, float *in, fftf_complex *out, unsigned)
00143 {
00144     if (n < 2) return 0;
00145     if (n & (n-1)) return 0;
00146     
00147     fftf_plan_ *plan = new fftf_plan_;
00148     plan->size = n;
00149     plan->inverse = 0;
00150     plan->real = in;
00151     plan->cplx = out;
00152     return plan;
00153 }
00154 
00155 fftf_plan
00156 fftf_plan_dft_c2r_1d(int n, fftf_complex *in, float *out, unsigned)
00157 {
00158     if (n < 2) return 0;
00159     if (n & (n-1)) return 0;
00160     
00161     fftf_plan_ *plan = new fftf_plan_;
00162     plan->size = n;
00163     plan->inverse = 1;
00164     plan->real = out;
00165     plan->cplx = in;
00166     return plan;
00167 }
00168 
00169 void
00170 fftf_destroy_plan(fftf_plan p)
00171 {
00172     delete p;
00173 }
00174 
00175 void
00176 fftf_execute(const fftf_plan p)
00177 {
00178     float *real = p->real;
00179     fftf_complex *cplx = p->cplx;
00180     int n = p->size;
00181     int forward = !p->inverse;
00182 
00183     double *ri = new double[n];
00184     double *ro = new double[n];
00185     double *io = new double[n];
00186 
00187     double *ii = 0;
00188     if (!forward) ii = new double[n];
00189 
00190     if (forward) {
00191         for (int i = 0; i < n; ++i) {
00192             ri[i] = real[i];
00193         }
00194     } else {
00195         for (int i = 0; i < n/2+1; ++i) {
00196             ri[i] = cplx[i][0];
00197             ii[i] = cplx[i][1];
00198             if (i > 0) {
00199                 ri[n-i] = ri[i];
00200                 ii[n-i] = -ii[i];
00201             }
00202         }
00203     }
00204 
00205     fft(n, !forward, ri, ii, ro, io);
00206 
00207     if (forward) {
00208         for (int i = 0; i < n/2+1; ++i) {
00209             cplx[i][0] = ro[i];
00210             cplx[i][1] = io[i];
00211         }
00212     } else {
00213         for (int i = 0; i < n; ++i) {
00214             real[i] = ro[i];
00215         }
00216     }
00217 
00218     delete[] ri;
00219     delete[] ro;
00220     delete[] io;
00221     if (ii) delete[] ii;
00222 }
00223 
00224 #endif