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;
+}
+