基2fft算法.txt

DIT-FFT算法子程序 FFT长度为2m

简介：本程序包括三个部分，头文件fft.h和主程序fft.c；另外还有一段程序是用fft计算实际题目,计算两个有限序列的线性卷积.

一、头文件fft.h
#ifndef FFT_DEF
#define FFT_DEF 1
#include<math.h>
#include<malloc.h>
typedef struct{
float real;
float imag;
}COMPLEX;
extern void ifft(COMPLEX *x,int m);
extern void fft(COMPLEX *x,int m);

/*快速傅里叶变换fft*/
void fft(COMPLEX *x,int m)
{
static COMPLEX *w;
static int mstore=0;
static int n=1;
COMPLEX u,temp,tm;
COMPLEX *xi,*xip,*xj,*wptr;
int i,j,k,l,le,windex;
double arg,w_real,w_imag,wrecur_real,wrecur_imag,wtemp_real;
if(m!=mstore)
{ if(mstore!=0) free(w);
mstore=m;
if(m==0) return;
n=1<<m; le=n/2;
w=(COMPLEX*)malloc((le-1)*sizeof(COMPLEX));
if(!w){ printf("
Unable to allocte complex w array 
"); exit(1); }
arg=4.0*atan(1.0)/le; xj=w;
wrecur_real=w_real=cos(arg); wrecur_imag=w_imag=-sin(arg ); 
for(j=1;j<le;j++)
{ xj->real=(float)wrecur_real; xj->imag=(float)wrecur_imag;
xj++;
wtemp_real=wrecur_real*w_real-wrecur_imag*w_imag;
wrecur_imag=wrecur_real*w_imag+wrecur_imag*w_real;
wrecur_real=wtemp_real;
}
}
le=n; windex=1;
for(l=0;l<m;l++)
{ le=le/2;
for(i=0;i<n;i=i+2*le)
{ xi=x+i; xip=xi+le;
temp.real=(xi->real+xip->real); temp.imag=(xi->imag+xip->imag);
xip->real=(xi->real-xip->real); xip->imag=(xi->imag-xip->imag);
*xi=temp;
}
wptr=w+windex-1;
for(j=1;j<le;j++)
{ u=*wptr;
for(i=j;i<n;i=i+2*le)
{ xi=x+i; xip=xi+le;
temp.real=(xi->real+xip->real); temp.imag=(xi->imag+xip->imag);
tm.real=xi->real-xip->real; tm.imag=xi->imag-xip->imag;
xip->real=(tm.real*u.real-tm.imag*u.imag);
xip->imag=(tm. real*u.imag+tm.imag*u.real);
*xi=temp;
}
wptr=wptr+windex;
}
windex=2*windex;
}
for(i=0;i<n;++i)
{ j=0;
for(k=0;k<m;++k)
j=(j<<1)|(1&(i>>k));
if(i<j)
{ xi=x+i; xj=x+j; temp=*xj; *xj=*xi; *xi=temp; }
}
}

/*快速傅里叶反变换ifft*/
void ifft(COMPLEX *x,int m)
{ 
static COMPLEX *w;
static int mstore=0;
static int n=1;
COMPLEX u,temp,tm;
COMPLEX *xi,*xip,*xj,*wptr;
int i,j,k,l,le,windex;
double arg,w_real,w_imag,wrecur_real,wrecur_imag,wtemp_real;
if(m!=mstore)
{ if(mstore!=0) free(w);
mstore=m;
if(m==0) return;
n=1<<m; le=n/2;
w=(COMPLEX*)malloc((le-1)*sizeof(COMPLEX));
if(!w){ printf("
Unable to allocte complex w array 
"); exit(1); }
arg=4.0*atan(1.0)/le; xj=w;
wrecur_real=w_real=cos(arg); wrecur_imag=w_imag=-sin(arg ); 
for(j=1;j<le;j++)
{ xj->real=(float)wrecur_real; xj->imag=(float)wrecur_imag;
xj++;
wtemp_real=wrecur_real*w_real-wrecur_imag*w_imag;
wrecur_imag=wrecur_real*w_imag+wrecur_imag*w_real;
wrecur_real=wtemp_real;
}
}
le=n; windex=1;
for(l=0;l<m;l++)
{ le=le/2;
for(i=0;i<n;i=i+2*le)
{ xi=x+i; xip=xi+le;
temp.real=(xi->real+xip->real); temp.imag=(xi->imag+xip->imag);
xip->real=(xi->real-xip->real); xip->imag=(xi->imag-xip->imag);
*xi=temp;
}
wptr=w+windex-1;
for(j=1;j<le;j++)
{ u=*wptr;
for(i=j;i<n;i=i+2*le)
{ xi=x+i; xip=xi+le;
temp.real=(xi->real+xip->real); temp.imag=(xi->imag+xip->imag);
tm.real=xi->real-xip->real; tm.imag=xi->imag-xip->imag;
xip->real=(tm.real*u.real-tm.imag*u.imag);
xip->imag=(tm.real*u.imag+tm.imag*u.real);
*xi=temp;
}
wptr=wptr+windex;
}
windex=2*windex;
}
for(i=0;i<n;++i)
{ j=0;
for(k=0;k<m;++k)
j=(j<<1)|(1&(i>>k));
if(i<j)
{ xi=x+i; xj=x+j; temp=*xj; *xj=*xi; *xi=temp; }
}
scale=(float)(1.0/n);
for(i=0;i<n;i++)
{ x[i].real=scale*x[i].real; x[i].imag=scale*x[i].imag; }
}
#endif

二、主程序fft.c
#include"FFT.h"
#include"stdio.h"
void main()
{
COMPLEX x[16],y[16];
int i;
printf("Input 16 point for fft(like this x+yj ):
");
for(i=0;i<16;i++)
{ scanf("%f+%fj",&x[i].real,&x[i].imag); }
printf("Input 16 point for ifft(like this x+yj ):
");
for(i=0;i<16;i++)
{ scanf("%f+%fj",&y[i].real,&y[i].imag); }

fft(x,4);
printf("The FFT result is:
");
for(i=0;i<16;i++)
{
printf("%f",x[i].real);
if(x[i].imag>=0)printf("+");
printf("%fj ",x[i].imag);
if((i+1)%3==0) printf("
");
}
printf("
");
ifft(y,4);
printf("The IFFT result is:
");
for(i=0;i<16;i++)
{
printf("%f",y[i].real);
if(y[i].imag>=0)printf("+");
printf("%fj ",y[i].imag);
if((i+1)%3==0) printf("
");
}
}
测试数据：fft: 0，1，2，3，4，5，6，7，8，9，10，11，12，13，14，15
ifft:16，15，14，13，12，11，10，9，8，7，6，5，4，3，2，1
运行结果：.
The FFT result is:
120.000000+0.000000j -8.000000+40.218716j -8.000000+19.313709j
-8.000000+11.972845j -8.000000+8.000000j -8.000000+5.345429j
-8.000000+3.313708j -8.000000+1.591297j -8.000000+0.000000j
-8.000000-1.591297j -8.000000-3.313708j -8.000000-5.345428j
-8.000000-8.000000j -8.000000-11.972846j -8.000000-19.313709j
-8.000000-40.218716j

The IFFT result is:
8.500000+0.000000j 0.500000+2.513670j 0.500000+1.207107j
0.500000+0.748303j 0.500000+0.500000j 0.500000+0.334089j
0.500000+0.207107j 0.500000+0.099456j 0.500000+0.000000j
0.500000-0.099456j 0.500000-0.207107j 0.500000-0.334089j
0.500000-0.500000j 0.500000-0.748303j 0.500000-1.207107j
0.500000-2.513670j

三、用fft程序实现x(n)=1, 0<=n<=15 0, 其他,h(n) =a^n, 0<=n<=10 0, 其他的线性卷积。
#include"FFT.h"
#include"stdio.h"
main()
{ COMPLEX x[32],h[32],r[32];
float a;
int i;
printf("Input the value of a :
");
scanf("%f",&a);
for(i=0;i<32;i++)
{
if(i<16) x[i].real=1;
else x[i].real=0;
if(i<11) {
if(i==0) h[i].real=1;
else h[i].real=h[i-1].real*a;
}
else h[i].real=0;
h[i].imag=0; x[i].imag=0;
}
fft(x,5); fft(h,5);
for(i=0;i<32;i++)
{
r[i].real=x[i].real*h[i].real-x[i].imag*h[i].imag;
r[i].imag=x[i].real*h[i].imag+x[i].imag*h[i].real;
}
ifft(r,5);
printf("The result is:
");
for(i=0;i<32;i++)
{
printf("%f",r[i].real);
if(r[i].imag>=0)printf("+");
printf("%fj ",r[i].imag);
if((i+1)%3==0) printf("
");
}
}
运行结果如下：
Input the value of a :
.0.8
The result is:
1.000000+0.000000j 1.800000+0.000000j 2.440000+0.000000j
2.952000-0.000000j 3.361600+0.000000j 3.689280-0.000000j
3.951424+0.000000j 4.161139-0.000000j 4.328911-0.000000j
4.463129-0.000000j 4.570503-0.000000j 4.570503-0.000000j
4.570503-0.000000j 4.570503-0.000000j 4.570503-0.000000j
4.570503-0.000000j 3.570503+0.000000j 2.770503+0.000000j
2.130503+0.000000j 1.618503+0.000000j 1.208904+0.000000j
0.881223+0.000000j 0.619079+0.000000j 0.409364+0.000000j
0.241592+0.000000j 0.107374+0.000000j 0.000000+0.000000j
0.000000+0.000000j 0.000000+0.000000j 0.000000+0.000000j
0.000000+0.000000j 0.000000+0.00000