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
| #include <stdio.h>
#include <stdlib.h>
#include <math.h>
void FFT(int number);
#define SWAP(a,b) tempr=(a);(a)=(b);(b)=tempr
double tableau[2048]; // tableau pour 1024 points (les imagineurs =0)
void main()
{
FILE *fichier=NULL;
long lsb, msb, i;
long data, amplitude;
fichier=fopen("1Khz.data", "rb");
for (i=0; i<2048; i++)
{data=fgetc(fichier);
lsb=data; //codage little indian
data=fgetc(fichier);
msb=data;
data=msb*256+lsb; // la donnee
//conversion signé nonsigné
data=data-32768;
if (data<0) data =data+65536;
tableau[i]=data; // ranger dans le tableau
}
FFT(1024); //on fait la FFT!
i=0;
while (i<1024) // la moitié de la tableau (fréquences positives)
{
amplitude =sqrt(tableau[i]*tableau[i]+tableau[i+1]*tableau[i+1]);
printf("amplitude de fréquence No %ld est: %ld\n", i, amplitude);
i=i+2;
}
fclose(fichier);
system("PAUSE");
}
void FFT( int number)
{
long n,mmax,m,j,istep,i;
long wtemp,wr,wpr,wpi,wi,theta,tempr,tempi;
n=number* 2;
j=0;
for (i=0;i<n/2;i+=2) {
if (j > i) {
SWAP(tableau[j],tableau[i]);
SWAP(tableau[j+1],tableau[i+1]);
if((j/2)<(n/4)){
SWAP(tableau[(n-(i+2))],tableau[(n-(j+2))]);
SWAP(tableau[(n-(i+2))+1],tableau[(n-(j+2))+1]);
}
}
m=n/2;
while (m >= 2 && j >= m) {
j -= m;
m = m/2;
}
j += m;
}
mmax=2;
while (n > mmax)
{
istep = mmax<< 1;
theta=(2*3.141592653589793/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*tableau[j-1]-wi*tableau[j];
tempi=wr*tableau[j]+wi*tableau[j-1];
tableau[j-1]=tableau[i-1]-tempr;
tableau[j]=tableau[i]-tempi;
tableau[i-1] += tempr;
tableau[i] += tempi;
}
wr=(wtemp=wr)*wpr-wi*wpi+wr;
wi=wi*wpr+wtemp*wpi+wi;
}
mmax=istep;
}
} |
Partager