Marsyas
0.6.0-alpha
|
00001 /* 00002 ** Copyright (C) 1998-2004 George Tzanetakis <gtzan@cs.uvic.ca> 00003 ** 00004 ** This program is free software; you can redistribute it and/or modify 00005 ** it under the terms of the GNU General Public License as published by 00006 ** the Free Software Foundation; either version 2 of the License, or 00007 ** (at your option) any later version. 00008 ** 00009 ** This program is distributed in the hope that it will be useful, 00010 ** but WITHOUT ANY WARRANTY; without even the implied warranty of 00011 ** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 00012 ** GNU General Public License for more details. 00013 ** 00014 ** You should have received a copy of the GNU General Public License 00015 ** along with this program; if not, write to the Free Software 00016 ** Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. 00017 */ 00018 00019 #include <marsyas/fft.h> 00020 00021 using namespace Marsyas; 00022 00023 //typedef struct { mrs_real re ; mrs_real im ; } complex ; //lmartins: not needed and can create name clashing! 00024 00025 /* 00026 * bitreverse places real array x containing N/2 complex values 00027 * into bit-reversed order 00028 */ 00029 00030 void 00031 fft::bitreverse(mrs_real x[], int N ) 00032 { 00033 mrs_real rtemp, itemp ; 00034 int i, j, m ; 00035 for ( i = j = 0 ; i < N ; i += 2, j += m ) { 00036 if ( j > i ) { 00037 rtemp = x[j] ; itemp = x[j+1] ; /* complex exchange */ 00038 x[j] = x[i] ; x[j+1] = x[i+1] ; 00039 x[i] = rtemp ; x[i+1] = itemp ; 00040 } 00041 for ( m = N>>1 ; m >= 2 && j >= m ; m >>= 1 ) 00042 j -= m ; 00043 } 00044 } 00045 00046 00047 00048 /* 00049 * COMPUTES A 2*N SIZED FOURIER TRANSFORM 00050 * 00051 * If forward is true, rfft replaces 2*N real data points in x with 00052 * N complex values representing the positive frequency half of their 00053 * Fourier spectrum, with x[1] replaced with the real part of the Nyquist 00054 * frequency value. If forward is false, rfft expects x to contain a 00055 * positive frequency spectrum arranged as before, and replaces it with 00056 * 2*N real values. N MUST be a power of 2. 00057 */ 00058 void 00059 fft::rfft( mrs_real x[], int N, int forward ) 00060 { 00061 mrs_real c1, c2, h1r, h1i, h2r, h2i, wr, wi, wpr, wpi, temp, theta ; 00062 mrs_real xr, xi ; 00063 int i, i1, i2, i3, i4, N2p1 ; 00064 theta = PI/N ; 00065 wr = 1. ; 00066 wi = 0. ; 00067 c1 = 0.5 ; 00068 if ( forward ) 00069 { 00070 c2 = -0.5 ; 00071 cfft( x, N, forward ) ; 00072 xr = x[0] ; 00073 xi = x[1] ; 00074 } 00075 else 00076 { 00077 c2 = 0.5 ; 00078 theta = -theta ; 00079 xr = x[1] ; 00080 xi = 0. ; 00081 x[1] = 0. ; 00082 } 00083 wpr = (mrs_real)(-2.*pow( sin( 0.5*theta ), 2. )); 00084 wpi = sin( theta ) ; 00085 N2p1 = (N<<1) + 1 ; 00086 for ( i = 0 ; i <= N>>1 ; ++i ) 00087 { 00088 i1 = i<<1 ; 00089 i2 = i1 + 1 ; 00090 i3 = N2p1 - i2 ; 00091 i4 = i3 + 1 ; 00092 if ( i == 0 ) { 00093 h1r = c1*(x[i1] + xr ) ; 00094 h1i = c1*(x[i2] - xi ) ; 00095 h2r = -c2*(x[i2] + xi ) ; 00096 h2i = c2*(x[i1] - xr ) ; 00097 x[i1] = h1r + wr*h2r - wi*h2i ; 00098 x[i2] = h1i + wr*h2i + wi*h2r ; 00099 xr = h1r - wr*h2r + wi*h2i ; 00100 xi = -h1i + wr*h2i + wi*h2r ; 00101 } 00102 else { 00103 h1r = c1*(x[i1] + x[i3] ) ; 00104 h1i = c1*(x[i2] - x[i4] ) ; 00105 h2r = -c2*(x[i2] + x[i4] ) ; 00106 h2i = c2*(x[i1] - x[i3] ) ; 00107 x[i1] = h1r + wr*h2r - wi*h2i ; 00108 x[i2] = h1i + wr*h2i + wi*h2r ; 00109 x[i3] = h1r - wr*h2r + wi*h2i ; 00110 x[i4] = -h1i + wr*h2i + wi*h2r ; 00111 } 00112 wr = (temp = wr)*wpr - wi*wpi + wr ; 00113 wi = wi*wpr + temp*wpi + wi ; 00114 } 00115 if ( forward ) 00116 x[1] = xr ; 00117 else 00118 cfft( x, N, forward ) ; 00119 } 00120 00121 00122 /* 00123 * cfft replaces real array x containing NC complex values 00124 * (2*NC real values alternating real, imagininary, etc.) 00125 * by its Fourier transform if forward is true, or by its 00126 * inverse Fourier transform if forward is false, using a 00127 * recursive Fast Fourier transform method due to Danielson 00128 * and Lanczos. NC MUST be a power of 2. 00129 */ 00130 00131 void 00132 fft::cfft(mrs_real x[], int NC, int forward ) 00133 { 00134 mrs_real wr, wi, wpr, wpi, theta, scale ; 00135 int mmax, ND, m, i, j, delta ; 00136 ND = NC<<1 ; 00137 bitreverse( x, ND ) ; 00138 for ( mmax = 2 ; mmax < ND ; mmax = delta ) { 00139 delta = mmax<<1 ; 00140 theta = TWOPI/( forward? mmax : -mmax ) ; 00141 wpr = (mrs_real)(-2.*pow( sin( 0.5*theta ), 2. )); 00142 wpi = sin( theta ) ; 00143 wr = 1. ; 00144 wi = 0. ; 00145 for ( m = 0 ; m < mmax ; m += 2 ) { 00146 register mrs_real rtemp, itemp ; 00147 for ( i = m ; i < ND ; i += delta ) 00148 { 00149 j = i + mmax ; 00150 rtemp = wr*x[j] - wi*x[j+1] ; 00151 itemp = wr*x[j+1] + wi*x[j] ; 00152 x[j] = x[i] - rtemp ; 00153 x[j+1] = x[i+1] - itemp ; 00154 x[i] += rtemp ; 00155 x[i+1] += itemp ; 00156 } 00157 wr = (rtemp = wr)*wpr - wi*wpi + wr ; 00158 wi = wi*wpr + rtemp*wpi + wi ; 00159 } 00160 } 00161 /* 00162 * scale output 00163 */ 00164 00165 scale = forward ? (mrs_real)1.0/ND : (mrs_real)2.0 ; 00166 { register mrs_real *xi=x, *xe=x+ND ; 00167 while ( xi < xe ) 00168 *xi++ *= scale ; 00169 } 00170 00171 } 00172 00173 00174 00175 00176 00177 00178 00179 00180