svcore
1.9
|
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