diff options
Diffstat (limited to 'tests/test5/fft.c')
-rw-r--r-- | tests/test5/fft.c | 80 |
1 files changed, 79 insertions, 1 deletions
diff --git a/tests/test5/fft.c b/tests/test5/fft.c index 902d7f0..a9957a7 100644 --- a/tests/test5/fft.c +++ b/tests/test5/fft.c @@ -2,12 +2,18 @@ #include <math.h> #include <stdlib.h> +/* Everything here comes from Audacity 1.3.13 + (orignally in C++ and with more genericity and functionnality) + Original Author : Dominic Mazzoni + Licenced under GPL 2.0 (see LICENCE) +*/ + #define MaxFastBits 16 int **gFFTBitTable = NULL; void FFT(int NumSamples, - gboolean InverseTransform, + int InverseTransform, float *RealIn, float *ImagIn, float *RealOut, float *ImagOut) { int NumBits; /* Number of bits needed to store indices */ @@ -151,3 +157,75 @@ int ReverseBits(int index, int NumBits) } +/* + * PowerSpectrum + * + * This function computes the same as RealFFT, above, but + * adds the squares of the real and imaginary part of each + * coefficient, extracting the power and throwing away the + * phase. + * + * For speed, it does not call RealFFT, but duplicates some + * of its code. + */ + +void PowerSpectrum(float *In, float *Out) +{ + int i; + + float theta = M_PI / 128; + + float tmpReal[128]; + float tmpImag[128]; + float RealOut[128]; + float ImagOut[128]; + + for (i = 0; i < 128; i++) { + tmpReal[i] = In[2 * i]; + tmpImag[i] = In[2 * i + 1]; //FIXME : WTFFFF ? + tmpImag[i]=0; + } + + FFT(128, 0, tmpReal, tmpImag, RealOut, ImagOut); + + float wtemp = sin(0.5 * theta); + + float wpr = -2.0 * wtemp * wtemp; + float wpi = -1.0 * sin(theta); + float wr = 1.0 + wpr; + float wi = wpi; + + int i3; + + float h1r, h1i, h2r, h2i, rt, it; + for (i = 1; i < 128 / 2; i++) { + + i3 = 128 - i; + + h1r = 0.5 * (RealOut[i] + RealOut[i3]); + h1i = 0.5 * (ImagOut[i] - ImagOut[i3]); + h2r = 0.5 * (ImagOut[i] + ImagOut[i3]); + h2i = -0.5 * (RealOut[i] - RealOut[i3]); + + rt = h1r + wr * h2r - wi * h2i; + it = h1i + wr * h2i + wi * h2r; + + Out[i] = rt * rt + it * it; + + rt = h1r - wr * h2r + wi * h2i; + it = -h1i + wr * h2i + wi * h2r; + + Out[i3] = rt * rt + it * it; + + wr = (wtemp = wr) * wpr - wi * wpi + wr; + wi = wi * wpr + wtemp * wpi + wi; + } + + rt = (h1r = RealOut[0]) + ImagOut[0]; + it = h1r - ImagOut[0]; + Out[0] = rt * rt + it * it; + rt = RealOut[128 / 2]; + it = ImagOut[128 / 2]; + Out[128 / 2] = rt * rt + it * it; +} + |