From 438cfc61ee3fae8875dbbf24a44fa7d75969f2b3 Mon Sep 17 00:00:00 2001
From: Ludovic Pouzenc <ludovic@pouzenc.fr>
Date: Wed, 6 Jun 2012 21:42:40 +0000
Subject: Ajout Makefile, résolution crash malloc (tableau statique à la
place). Données audio 2 channels : apparament avec un offset style [i*nchan]
on obtient bcp de valeurs à 0... FFT : valeurs en sortie toutes pétées, faut
débugguer. Conso CPU : correcte
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit
git-svn-id: file:///var/svn/2012-violon-leds/trunk@13 6be1fa4d-33ac-4c33-becc-79fcb3794bb6
---
tests/test5/fft.c | 153 ++++++++++++++++++++++++++++++++++++++++++++++++++++++
1 file changed, 153 insertions(+)
create mode 100644 tests/test5/fft.c
(limited to 'tests/test5/fft.c')
diff --git a/tests/test5/fft.c b/tests/test5/fft.c
new file mode 100644
index 0000000..902d7f0
--- /dev/null
+++ b/tests/test5/fft.c
@@ -0,0 +1,153 @@
+#include "fft.h"
+#include <math.h>
+#include <stdlib.h>
+
+#define MaxFastBits 16
+
+int **gFFTBitTable = NULL;
+
+void FFT(int NumSamples,
+ gboolean InverseTransform,
+ float *RealIn, float *ImagIn, float *RealOut, float *ImagOut)
+{
+ int NumBits; /* Number of bits needed to store indices */
+ int i, j, k, n;
+ int BlockSize, BlockEnd;
+
+ double angle_numerator = 2.0 * M_PI;
+ double tr, ti; /* temp real, temp imaginary */
+/*
+ if (!IsPowerOfTwo(NumSamples)) {
+ fprintf(stderr, "%d is not a power of two\n", NumSamples);
+ exit(1);
+ }
+*/
+ if (!gFFTBitTable)
+ InitFFT();
+
+ if (!InverseTransform)
+ angle_numerator = -angle_numerator;
+
+ NumBits = NumberOfBitsNeeded(NumSamples);
+
+ /*
+ ** Do simultaneous data copy and bit-reversal ordering into outputs...
+ */
+ for (i = 0; i < NumSamples; i++) {
+ j = FastReverseBits(i, NumBits);
+ RealOut[j] = RealIn[i];
+ ImagOut[j] = (ImagIn == NULL) ? 0.0 : ImagIn[i];
+ }
+
+ /*
+ ** Do the FFT itself...
+ */
+
+ BlockEnd = 1;
+ for (BlockSize = 2; BlockSize <= NumSamples; BlockSize <<= 1) {
+
+ double delta_angle = angle_numerator / (double) BlockSize;
+
+ double sm2 = sin(-2 * delta_angle);
+ double sm1 = sin(-delta_angle);
+ double cm2 = cos(-2 * delta_angle);
+ double cm1 = cos(-delta_angle);
+ double w = 2 * cm1;
+ double ar0, ar1, ar2, ai0, ai1, ai2;
+
+ for (i = 0; i < NumSamples; i += BlockSize) {
+ ar2 = cm2;
+ ar1 = cm1;
+
+ ai2 = sm2;
+ ai1 = sm1;
+
+ for (j = i, n = 0; n < BlockEnd; j++, n++) {
+ ar0 = w * ar1 - ar2;
+ ar2 = ar1;
+ ar1 = ar0;
+
+ ai0 = w * ai1 - ai2;
+ ai2 = ai1;
+ ai1 = ai0;
+
+ k = j + BlockEnd;
+ tr = ar0 * RealOut[k] - ai0 * ImagOut[k];
+ ti = ar0 * ImagOut[k] + ai0 * RealOut[k];
+
+ RealOut[k] = RealOut[j] - tr;
+ ImagOut[k] = ImagOut[j] - ti;
+
+ RealOut[j] += tr;
+ ImagOut[j] += ti;
+ }
+ }
+ BlockEnd = BlockSize;
+ }
+
+ /*
+ ** Need to normalize if inverse transform...
+ */
+
+ if (InverseTransform) {
+ float denom = (float) NumSamples;
+
+ for (i = 0; i < NumSamples; i++) {
+ RealOut[i] /= denom;
+ ImagOut[i] /= denom;
+ }
+ }
+}
+
+void InitFFT()
+{
+ gFFTBitTable = malloc(MaxFastBits*sizeof(int));
+
+ int len = 2;
+ int b, i;
+ for (b=1; b<=MaxFastBits; b++) {
+ gFFTBitTable[b-1]=malloc(len*sizeof(int));
+
+ for (i=0; i<len; i++)
+ gFFTBitTable[b-1][i] = ReverseBits(i, b);
+
+ len <<= 1;
+ }
+}
+
+int NumberOfBitsNeeded(int PowerOfTwo)
+{
+ int i;
+
+/*
+ if (PowerOfTwo < 2) {
+ fprintf(stderr, "Error: FFT called with size %d\n", PowerOfTwo);
+ exit(1);
+ }
+*/
+ for (i = 0;; i++)
+ if (PowerOfTwo & (1 << i))
+ return i;
+}
+
+inline int FastReverseBits(int i, int NumBits)
+{
+ if (NumBits <= MaxFastBits)
+ return gFFTBitTable[NumBits - 1][i];
+ else
+ return ReverseBits(i, NumBits);
+}
+
+int ReverseBits(int index, int NumBits)
+{
+ int i, rev;
+
+ for (i = rev = 0; i < NumBits; i++) {
+ rev = (rev << 1) | (index & 1);
+ index >>= 1;
+ }
+
+ return rev;
+}
+
+
--
cgit v1.2.3