1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125
| #include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define PI 3.1415926535897932385
#define SWAP(a,b) tempr=(a);(a)=(b);(b)=tempr
void FFT(float data[],unsigned long N,int isign);
int main()
{
int i, j=5; // j pour le signal carre
float sample, amplitude, A[2048];
/* géneration d'un signal sinusoidal 1000hz echantilloné à 44100*/
for (i=0; i<1024; i++)
{ sample= 32000*sin((i/44.1)*2*PI); // 44100:1000=44.1 (44.1 echantillons par periode)
A[2*i]=sample;
A[2*i+1]=0;
}
/* ******* ou un signal rectangle (plusieurs harmoniques)
for( i = 0 ; i < 1024 ; i++)
{
if( i % 44 == 0 ) // 44 echantillons par période
{
if ( j == 5 ) // amplitude = 5 volt
j = -5;
else
j = 5;
}
A[2*i] = j;
A[2*i+1] =0;
} ******* */
FFT(A,1024,1);
/* Affiche la TF du signal d'entrée */
i = 0 ;
while ( i < 256) //les 256 premières fréquences
{
amplitude = sqrt(A[i]*A[i]+A[i+1]*A[i+1]);
printf(" l'amplitude à l'indice %ld est: %f \n",i, amplitude);
i=i+2;
}
system("PAUSE");
return 0;
}
void FFT(float data[], unsigned long N, int isign)
{
unsigned long n,mmax,m,j,istep,i;
double wtemp,wr,wpr,wpi,wi,theta,tempr,tempi;
n=N * 2;
j=0;
for (i=0;i<n/2;i+=2) {
if (j > i) {
SWAP(data[j],data[i]);
SWAP(data[j+1],data[i+1]);
if((j/2)<(n/4)){
SWAP(data[(n-(i+2))],data[(n-(j+2))]);
SWAP(data[(n-(i+2))+1],data[(n-(j+2))+1]);
}
}
m=n/2;
while (m >= 2 && j >= m) {
j -= m;
m = m/2;
}
j += m;
}
mmax=2;
//external loop
while (n > mmax)
{
istep = mmax<< 1;
theta=isign*(6.28318530717959/mmax);
wtemp=sin(0.5*theta);
wpr = -2.0*wtemp*wtemp;
wpi=sin(theta);
wr=1.0;
wi=0.0;
//internal loops
for (m=1;m<mmax;m+=2) {
for (i= m;i<=n;i+=istep) {
j=i+mmax;
tempr=wr*data[j-1]-wi*data[j];
tempi=wr*data[j]+wi*data[j-1];
data[j-1]=data[i-1]-tempr;
data[j]=data[i]-tempi;
data[i-1] += tempr;
data[i] += tempi;
}
wr=(wtemp=wr)*wpr-wi*wpi+wr;
wi=wi*wpr+wtemp*wpi+wi;
}
mmax=istep;
}
} |