realfft.c

数字信号处理中的实序列的快速傅里叶变换的实现

#include "stdafx.h" 
#include "math.h" 
void rfft(int n, double x[]); 

int main(int argc, char* argv[]) 
{ 
int i,n; 
double x[64]; 
n=64; 
for(i=0;i <10;i++) x[i]=0; 
for(i=0;i <n;i++) x[i]=exp(-(i-10)/15.0)*sin(6.2831853*(i-10)/16.0); 
rfft(n,x); 
printf("\n DISCRETE FOURIER TRANSFORM\n"); 
printf(" %10.7f",x[0]); 
printf("\n"); 
printf(" %10.7f + J%10.7f\n",x[1],x[n-1]); 
for(i=2;i <n/2;i+=2) 
{ 
printf(" %10.7f + J%10.7f\n",x[i],x[n-i]); 
printf(" %10.7f + J%10.7f\n",x[i+1],x[n-i-1]); 
//printf("\n"); 
} 
printf(" %10.7f",x[n/2]); 
printf("\n"); 
printf(" %10.7f + J%10.7f\n",x[n/2-1],-x[n/2+1]); 
for(i=2;i <n/2;i+=2) 
{ 
printf(" %10.7f + J%10.7f\n",x[n/2-i],-x[n/2+i]); 
printf(" %10.7f + J%10.7f\n",x[n/2-i-1],-x[n/2+i+1]); 
//printf("\n"); 
} 
return 0; 
} 

void rfft(int n, double x[]) //实序列快速傅里叶变换 
{ 
int i,j,k,m,i1,i2,i3,i4,n1,n2,n4; 
double a,e,cc,ss,xt,t1,t2; 
for(j=1,i=1;i <16;i++) 
{ 
m=i; 
j= 2*j; 
if(j==n) break; 
} 
n1=n-1; 
for(j=0,i=0;i <n1;i++) 
{ 
if(i <j) 
{ 
xt=x[j]; 
x[j]=x[i]; 
x[i]=xt; 
} 
k=n/2; 
while (k <(j+1)) { 
j=j-k; 
k=k/2; 
} 
j=j+k; 
} 
for(i=0;i <n;i+=2) 
{ 
xt=x[i]; 
x[i]=xt+x[i+1]; 
x[i+1]=xt-x[i+1]; 
} 
n2=1; 
for(k=2;k <=m;k++) 
{ 
n4=n2; 
n2=2*n4; 
n1=2*n2; 
e=6.28318530718/n1; 
for(i=0;i <n;i+=n1) 
{ 
xt=x[i]; 
x[i]=xt+x[i+n2]; 
x[i+n2]=xt-x[i+n2]; 
x[i+n2+n4]=-x[i+n2+n4];// 
a=e; 
for(j=1;j <(n4-1);j++) 
{ 
i1=i+j; 
i2=i-j+n2; 
i3=i+j+n2; 
i4=i-j+n1; 
cc=cos(a); 
ss=sin(a); 
a=a+e; 
t1=cc*x[i3]+ss*x[i4]; 
t2=ss*x[i3]-cc*x[i4]; 
x[i4]=x[i2]-t2; 
x[i3]=-x[i2]-t2;// 
x[i2]=x[i1]-t1; 
x[i1]=x[i1]+t1; 
} 
} 
} 
}