1 /* 2 * Copyright (C) 2010 The Android Open Source Project 3 * 4 * Licensed under the Apache License, Version 2.0 (the "License"); you may not 5 * use this file except in compliance with the License. You may obtain a copy of 6 * the License at 7 * 8 * http://www.apache.org/licenses/LICENSE-2.0 9 * 10 * Unless required by applicable law or agreed to in writing, software 11 * distributed under the License is distributed on an "AS IS" BASIS, WITHOUT 12 * WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the 13 * License for the specific language governing permissions and limitations under 14 * the License. 15 */ 16 17 #include "Fft.h" 18 19 #include <stdlib.h> 20 #include <math.h> 21 22 Fft::Fft(void) : mCosine(0), mSine(0), mFftSize(0), mFftTableSize(0) { } 23 24 Fft::~Fft(void) { 25 fftCleanup(); 26 } 27 28 /* Construct a FFT table suitable to perform a DFT of size 2^power. */ 29 void Fft::fftInit(int power) { 30 fftCleanup(); 31 fftMakeTable(power); 32 } 33 34 void Fft::fftCleanup() { 35 delete [] mSine; 36 delete [] mCosine; 37 mSine = NULL; 38 mCosine = NULL; 39 mFftTableSize = 0; 40 mFftSize = 0; 41 } 42 43 /* z <- (10 * log10(x^2 + y^2)) for n elements */ 44 int Fft::fftLogMag(float *x, float *y, float *z, int n) { 45 float *xp, *yp, *zp, t1, t2, ssq; 46 47 if(x && y && z && n) { 48 for(xp=x+n, yp=y+n, zp=z+n; zp > z;) { 49 t1 = *--xp; 50 t2 = *--yp; 51 ssq = (t1*t1)+(t2*t2); 52 *--zp = (ssq > 0.0)? 10.0 * log10((double)ssq) : -200.0; 53 } 54 return 1; //true 55 } else { 56 return 0; // false/fail 57 } 58 } 59 60 int Fft::fftMakeTable(int pow2) { 61 int lmx, lm; 62 float *c, *s; 63 double scl, arg; 64 65 mFftSize = 1 << pow2; 66 mFftTableSize = lmx = mFftSize/2; 67 mSine = new float[lmx]; 68 mCosine = new float[lmx]; 69 scl = (M_PI*2.0)/mFftSize; 70 for (s=mSine, c=mCosine, lm=0; lm<lmx; lm++ ) { 71 arg = scl * lm; 72 *s++ = sin(arg); 73 *c++ = cos(arg); 74 } 75 mBase = (mFftTableSize * 2)/mFftSize; 76 mPower2 = pow2; 77 return(mFftTableSize); 78 } 79 80 81 /* Compute the discrete Fourier transform of the 2**l complex sequence 82 * in x (real) and y (imaginary). The DFT is computed in place and the 83 * Fourier coefficients are returned in x and y. 84 */ 85 void Fft::fft( float *x, float *y ) { 86 float c, s, t1, t2; 87 int j1, j2, li, lix, i; 88 int lmx, lo, lixnp, lm, j, nv2, k=mBase, im, jm, l = mPower2; 89 90 for (lmx=mFftSize, lo=0; lo < l; lo++, k *= 2) { 91 lix = lmx; 92 lmx /= 2; 93 lixnp = mFftSize - lix; 94 for (i=0, lm=0; lm<lmx; lm++, i += k ) { 95 c = mCosine[i]; 96 s = mSine[i]; 97 for ( li = lixnp+lm, j1 = lm, j2 = lm+lmx; j1<=li; 98 j1+=lix, j2+=lix ) { 99 t1 = x[j1] - x[j2]; 100 t2 = y[j1] - y[j2]; 101 x[j1] += x[j2]; 102 y[j1] += y[j2]; 103 x[j2] = (c * t1) + (s * t2); 104 y[j2] = (c * t2) - (s * t1); 105 } 106 } 107 } 108 109 /* Now perform the bit reversal. */ 110 j = 1; 111 nv2 = mFftSize/2; 112 for ( i=1; i < mFftSize; i++ ) { 113 if ( j < i ) { 114 jm = j-1; 115 im = i-1; 116 t1 = x[jm]; 117 t2 = y[jm]; 118 x[jm] = x[im]; 119 y[jm] = y[im]; 120 x[im] = t1; 121 y[im] = t2; 122 } 123 k = nv2; 124 while ( j > k ) { 125 j -= k; 126 k /= 2; 127 } 128 j += k; 129 } 130 } 131 132 /* Compute the discrete inverse Fourier transform of the 2**l complex 133 * sequence in x (real) and y (imaginary). The DFT is computed in 134 * place and the Fourier coefficients are returned in x and y. Note 135 * that this DOES NOT scale the result by the inverse FFT size. 136 */ 137 void Fft::ifft(float *x, float *y ) { 138 float c, s, t1, t2; 139 int j1, j2, li, lix, i; 140 int lmx, lo, lixnp, lm, j, nv2, k=mBase, im, jm, l = mPower2; 141 142 for (lmx=mFftSize, lo=0; lo < l; lo++, k *= 2) { 143 lix = lmx; 144 lmx /= 2; 145 lixnp = mFftSize - lix; 146 for (i=0, lm=0; lm<lmx; lm++, i += k ) { 147 c = mCosine[i]; 148 s = - mSine[i]; 149 for ( li = lixnp+lm, j1 = lm, j2 = lm+lmx; j1<=li; 150 j1+=lix, j2+=lix ) { 151 t1 = x[j1] - x[j2]; 152 t2 = y[j1] - y[j2]; 153 x[j1] += x[j2]; 154 y[j1] += y[j2]; 155 x[j2] = (c * t1) + (s * t2); 156 y[j2] = (c * t2) - (s * t1); 157 } 158 } 159 } 160 161 /* Now perform the bit reversal. */ 162 j = 1; 163 nv2 = mFftSize/2; 164 for ( i=1; i < mFftSize; i++ ) { 165 if ( j < i ) { 166 jm = j-1; 167 im = i-1; 168 t1 = x[jm]; 169 t2 = y[jm]; 170 x[jm] = x[im]; 171 y[jm] = y[im]; 172 x[im] = t1; 173 y[im] = t2; 174 } 175 k = nv2; 176 while ( j > k ) { 177 j -= k; 178 k /= 2; 179 } 180 j += k; 181 } 182 } 183 184 int Fft::fftGetSize(void) { return mFftSize; } 185 186 int Fft::fftGetPower2(void) { return mPower2; } 187
