http://www.youtube.com/watch?v=HemL3kkgkVk
#include "def.h"
#include "math.h"
#include "float.h"
#include "2440addr.h"
#include "2440lib.h"
#include "2440slib.h"
#include "lcd.h"
#include <stdlib.h>
#include <math.h>
/* function prototypes FFT */
void fft(int N, double (*x)[2], double (*X)[2]);
void fft_rec(int N, int offset, int delta, double (*x)[2], double (*X)[2], double (*XX)[2]);
void ifft(int N, double (*x)[2], double (*X)[2]);
void FTT_PAINT (int nn);
void init_FFT (void);
#define TWO_PI (6.2831853071795864769252867665590057683943L)
#define NN 256
int N = NN ; /* number of points in FFT */
double (*x)[2]; /* pointer to time-domain samples */
double (*X)[2]; /* pointer to frequency-domain samples */
double dummy; /* scratch variable */
int iFFT = 0; /* generic index */
U32 X_antigou[256][2];
#define RGB(g,r,b) (unsigned int)( (r << 16) + (g << 8) + b )
void init_FFT (void)
{
x = malloc(2 * N * sizeof(double));
X = malloc(2 * N * sizeof(double));
Lcd_printf( 15 , 30 ,RGB( 0xff,0xff,0xff),RGB( 0x00,0x00,0x00),0,"0" );
Lcd_printf( 80 , 30 ,RGB( 0xff,0xff,0xff),RGB( 0x00,0x00,0x00),0,"1K" );
Lcd_printf( 148 , 30 ,RGB( 0xff,0xff,0xff),RGB( 0x00,0x00,0x00),0,"2K" );
Lcd_printf( 216 , 30 ,RGB( 0xff,0xff,0xff),RGB( 0x00,0x00,0x00),0,"3K" );
Lcd_printf( 284 , 30 ,RGB( 0xff,0xff,0xff),RGB( 0x00,0x00,0x00),0,"4K" );
Lcd_printf( 352 , 30 ,RGB( 0xff,0xff,0xff),RGB( 0x00,0x00,0x00),0,"5K" );
Lcd_printf( 420 , 30 ,RGB( 0xff,0xff,0xff),RGB( 0x00,0x00,0x00),0,"6K" );
Lcd_printf( 488 , 30 ,RGB( 0xff,0xff,0xff),RGB( 0x00,0x00,0x00),0,"7K" );
Lcd_printf( 556 , 30 ,RGB( 0xff,0xff,0xff),RGB( 0x00,0x00,0x00),0,"8K" );
Lcd_printf( 624 , 30 ,RGB( 0xff,0xff,0xff),RGB( 0x00,0x00,0x00),0,"9K" );
Lcd_printf( 692 , 30 ,RGB( 0xff,0xff,0xff),RGB( 0x00,0x00,0x00),0,"10K" );
Lcd_printf( 760 , 30 ,RGB( 0xff,0xff,0xff),RGB( 0x00,0x00,0x00),0,"11K" );
}
void FTT_PAINT (int nn)
{
int m;
int mmm;
U32 X1;
fft(nn, x, X);
mmm=0;
for (m = 1; m <= (nn/2)-1 ; m++ )
{ X1 = sqrt( pow(X[m][0],2) + pow(X[m][1],2) );
if((X1) >= 415 )
{X1 = 415;}
mmm = mmm+6;
Glib_Line( 10+mmm , 476 , 10+mmm , 476-X_antigou[m][0] ,RGB(0X00,0X00,0X00));
Glib_Line( 10+mmm , 476 , 10+mmm , 476-X1 ,RGB(0XFA,0X4F,0X02));
Glib_Line( 10+mmm+1 , 476 , 10+mmm+1 , 476-X_antigou[m][0] ,RGB(0X00,0X00,0X00));
Glib_Line( 10+mmm+1 , 476 , 10+mmm+1 , 476-X1 ,RGB(0XFA,0X4F,0X02));
Glib_Line( 10+mmm+2 , 476 , 10+mmm+2 , 476-X_antigou[m][0] ,RGB(0X00,0X00,0X00));
Glib_Line( 10+mmm+2 , 476 , 10+mmm+2 , 476-X1 ,RGB(0XFA,0X4F,0X02));
Glib_Line( 10+mmm+3 , 476 , 10+mmm+3 , 476-X_antigou[m][0] ,RGB(0X00,0X00,0X00));
Glib_Line( 10+mmm+3 , 476 , 10+mmm+3 , 476-X1 ,RGB(0XFA,0X4F,0X02));
Glib_Line( 10+mmm+4 , 476 , 10+mmm+4 , 476-X_antigou[m][0] ,RGB(0X00,0X00,0X00));
Glib_Line( 10+mmm+4 , 476 , 10+mmm+4 , 476-X1 ,RGB(0XFA,0X4F,0X02));
Glib_Line( 10+mmm+5 , 476 , 10+mmm+5 , 476-X_antigou[m][0] ,RGB(0X00,0X00,0X00));
X_antigou[m][1] = X1;
}
for (m = 1; m <= (nn/2)-1 ; m++ )
{ X_antigou[m][0] = X_antigou[m][1]; }
}
/*----------------------------------------------------------------------------
fft.c - fast Fourier transform and its inverse (both recursively)
Copyright (C) 2004, Jerome R. Breitenbach. All rights reserved.
The author gives permission to anyone to freely copy, distribute, and use
this file, under the following conditions:
- No changes are made.
- No direct commercial advantage is obtained.
- No liability is attributed to the author for any damages incurred.
----------------------------------------------------------------------------*/
/* FFT */
void fft(int N, double (*x)[2], double (*X)[2])
{
double (*XX)[2] = malloc(2 * N * sizeof(double));
fft_rec(N, 0, 1 , x, X, XX);
free(XX);
}
/* FFT recursion */
void fft_rec(int N, int offset, int delta, double (*x)[2], double (*X)[2], double (*XX)[2])
{
int N2 = N/2; /* half the number of points in FFT */
int k; /* generic index */
double cs, sn; /* cosine and sine */
int k00, k01, k10, k11; /* indices for butterflies */
double tmp0, tmp1; /* temporary storage */
if(N != 2) /* Perform recursive step. */
{
fft_rec(N2, offset, 2*delta, x, XX, X);
fft_rec(N2, offset+delta, 2*delta, x, XX, X);
for(k=0; k<N2; k++)
{
k00 = offset + k*delta; k01 = k00 + N2*delta;
k10 = offset + 2*k*delta; k11 = k10 + delta;
cs = cos(TWO_PI*k/(double)N); sn = sin(TWO_PI*k/(double)N);
tmp0 = cs * XX[k11][0] + sn * XX[k11][1];
tmp1 = cs * XX[k11][1] - sn * XX[k11][0];
X[k01][0] = XX[k10][0] - tmp0;
X[k01][1] = XX[k10][1] - tmp1;
X[k00][0] = XX[k10][0] + tmp0;
X[k00][1] = XX[k10][1] + tmp1;
}
}
else
{
k00 = offset; k01 = k00 + delta;
X[k01][0] = x[k00][0] - x[k01][0];
X[k01][1] = x[k00][1] - x[k01][1];
X[k00][0] = x[k00][0] + x[k01][0];
X[k00][1] = x[k00][1] + x[k01][1];
}
}
/* IFFT */
void ifft(int N, double (*x)[2], double (*X)[2])
{
int N2 = N/2; /* half the number of points in IFFT */
int i; /* generic index */
double tmp0, tmp1; /* temporary storage */
/* Calculate IFFT via reciprocity property of DFT. */
fft(N, X, x);
x[0][0] = x[0][0]/N; x[0][1] = x[0][1]/N;
x[N2][0] = x[N2][0]/N; x[N2][1] = x[N2][1]/N;
for(i=1; i<N2; i++)
{
tmp0 = x[i][0]/N; tmp1 = x[i][1]/N;
x[i][0] = x[N-i][0]/N; x[i][1] = x[N-i][1]/N;
x[N-i][0] = tmp0; x[N-i][1] = tmp1;
}
}
