arotineoffft.txt

从网上下载的一个和matlab实现快速傅立叶变换相关的程序

下面给出FFT和IFFT的算法：

struct FFT_paras 
{ 
int N; // 输入序列长度 

float *pRin; //输入的实部，返回模 

float *pIin; //输入的虚部，返回幅角 

float *pRout; //返回FFT或IFFT结果的实部 

float *pIout; //返回FFT或IFFT结果的虚部 

bool do_IFFT; // FFT或IFFT的标志，TRUE 求 IFFT；FALSE 求 FFT 

bool get_all; // 是否输出模和幅角，TRUE 输出模和幅角；FALSE不输出 
}; 


void kbfft(FFT_paras para1) 
{ 
int it,m,is,i,j,nv,lo; 
float p,q,s,vr,vi,poddr,poddi; 

bool do_IFFT=para1.do_IFFT; 
bool get_all=para1.get_all; 
float* pr=para1.pRin ; 
float* pi=para1.pIin ; 
float* fr=para1.pRout ; 
float* fi=para1.pIout ; 

int n=para1.N; 
int k=int(log10(n)/log10(2)); 

for (it=0;it<=n-1;it ) 
{ 
m=it; is=0; 
for (i=0;i<=k-1;i ) 
{j=m/2;is=2*is (m-2*j); m=j;} 
fr[it]=pr[is]; 
fi[it]=pi[is]; 
} 

pr[0]=float(1.0); 
pi[0]=float(0.0); 
p=float(6.283185306/(1.0*n)); 
pr[1]=float(cos(p)); 
pi[1]=float(-sin(p)); 

if(do_IFFT) pi[1]=-1*pi[1];//如果是求逆快速傅立叶变换 

for(i=2;i<=n-1;i ) 
{ 
p=pr[i-1]*pr[1]; 
q=pi[i-1]*pi[1]; 
s=(pr[i-1] pi[i-1])*(pr[1] pi[1]); 
pr[i]=p-q; 
pi[i]=s-p-q; 
} 

for (it=0;it<=n-2;it=it 2) 
{ 
vr=fr[it]; 
vi=fi[it]; 
fr[it]=vr fr[it 1]; 
fi[it]=vi fi[it 1]; 
fr[it 1]=vr-fr[it 1]; 
fi[it 1]=vi-fi[it 1]; 
} 

m=n/2; nv=2; 
for(lo=k-2;lo>=0;lo--) 
{ 
m=m/2; nv=2*nv; 
for(it=0;it<=(m-1)*nv;it=it nv) 
for(j=0;j<=(nv/2)-1;j ) 
{ 
p=pr[m*j]*fr[it j nv/2]; 
q=pi[m*j]*fi[it j nv/2]; 
s=pr[m*j] pi[m*j]; 
s=s*(fr[it j nv/2] fi[it j nv/2]); 
poddr=p-q; 
poddi=s-p-q; 
fr[it j nv/2]=fr[it j]-poddr; 
fi[it j nv/2]=fi[it j]-poddi; 
fr[it j]=fr[it j] poddr; 
fi[it j]=fi[it j] poddi; 
} 
} 

if (do_IFFT) //如果是求逆快速傅立叶变换 
for(i=0;i<=n-1;i ) 
{ 
fr[i]=float(fr[i]/n); 
fi[i]=float(fi[i]/n); 
} 

if(get_all) 
for(i=0;i<=n-1;i ) 
{ 
pr[i]=float(sqrt(fr[i]*fr[i] fi[i]*fi[i])); 

if(fabs(fr[i])<0.000001*fabs(fi[i])) 
{ 
if ((fi[i]*fr[i])>0) 
pi[i]=float(3.1415926/2); 
else 
pi[i]=float(-3.1415926/2); 
} 
else 
pi[i]=float(atan(fi[i]/fr[i]));//幅角单位弧度 
} 

}