Digital signal Processing C language (3)------FFT

Fast Fourier transform

Fft.c

#include"math.h"voidFftDouble*x,Double*y,intNintSign ) {    inti,j,k,l,m,n1,n2; DoubleC,c1,e,s,s1,t,tr,ti;  for(j=1, i=1;i< -; i++) {m=i; J=2*J; if(j==n) Break; } N1=n-1;  for(j=0, i=0; i<n1;i++)    {        if(i<j) {TR=X[j]; TI=Y[j]; X[J]=X[i]; Y[J]=Y[i]; X[i]=TR; Y[i]=ti; } k=n/2;  while(k< (j+1) ) {J=j-K; K=k/2; } J=j+K; } N1=1;  for(l=1; l<=m;l++) {N1=2*N1; N2=n1/2; E=3.14159265359/N2; C=1.0; S=0.0; C1=cos (e); S1=-sign*sin (e);  for(j=0; j<n2;j++)        {             for(i=j;i<n;i+=N1) {k=i+N2; TR=c*x[k]-s*Y[k]; TI=c*y[k]+s*X[k]; X[K]=x[i]-TR; Y[K]=y[i]-ti; X[i]=x[i]+TR; Y[i]=y[i]+ti; } t=C; C=c*c1-s*S1; S=t*s1+s*C1; }    }    if(sign==-1)    {         for(i=0; i<n;i++) {X[i]/=N; Y[i]/=N; }    }}

#include <QCoreApplication>#include"fft.c"intMainintargcChar*argv[])    {qcoreapplication D (argc, argv); intI,j,n; Doublea1,a2,c,c1,c2,d1,d2,q1,q2,w,w1,w2; Doublex[ +],y[ +],a[ +],b[ +]; N= +; A1=0.9; A2=0.3; x[0]=1.0; y[0]=0.0;  for(i=1; i<n;i++) {X[i]=a1*x[i-1]-a2*y[i-1]; Y[i]=a2*x[i-1]+a1*y[i-1]; } printf ("\ input\n");  for(i=0; i<n/2; i++)    {         for(j=0;j<2; j + +) {printf ("%10.7f+j%10.7f", x[2*i+j],y[2*i+J]); } printf ("\ n"); } Q1=x[n-1]; Q2=y[n-1];  for(i=0; i<n;i++) {W=6.28348530718/n*i; W1=cos (w); W2=-sin (w); C1=1.0-a1*w1+a2*W2; C2=a1*w2+a2*W1; C=c1*c1+c2*C2; D1=1.0-a1*q1+a2*Q2; D2=a1*q2+a2*Q1; A[i]= (C1*D1+C2*D2)/C; B[i]= (C2*D1-C1*D2)/C; } printf ("\ Theoretical dft\n");  for(i=0; i<n/2; i++)    {         for(j=0;j<2; j + +) {printf ("%10.7f+j%10.7f", a[2*i+j],b[2*i+J]); } printf ("\ n"); } FFT (X,y,n,1); printf ("\ndft\n");  for(i=0; i<n/2; i++)    {         for(j=0;j<2; j + +) {printf ("%10.7f+j%10.7f", x[2*i+j],y[2*i+J]); } printf ("\ n"); } FFT (X,y,n,-1); printf ("\nidft\n");  for(i=0; i<n/2; i++)    {         for(j=0;j<2; j + +) {printf ("%10.7f+j%10.7f", x[2*i+j],y[2*i+J]); } printf ("\ n"); }    returnd.exec ();}

Digital signal Processing C language (3)------FFT