C＋＋编写的FFT算法函数（两种方法）
读研的时候用C++写的两种方法的FFT算法，FFT和IFFT是浙江大学出版社《数值分析》书的方法，记不得名字了。FFT2和IFFT2是一般信号处理书上讲的蝶形算法，都经过调试与matlab校验过。需要说明的是函数中用到的MatrixC是我自己用C++构造的一个复数矩阵类，Complex也是自己构造的复数类，各自封装了一些方法以便简化主程序，提高可读性。void FFT(MatrixC &fk,MatrixC &Fn,int N){
const double pi=3.1415926;
int k,n,P,q;
P=(int)(log10(N)/log10(2));
MatrixC temp(fk);
for(q=1;q<=P;q++)
{
for(k=0;k<(int)(pow(2,(P-q)));k++)
{
for(n=0;n<(int)(pow(2,(q-1)));n++)
{
Fn((int)(k*pow(2,q)+n),0)=temp((int)(k*pow(2,q-1)+n),0)
+temp((int)(k*pow(2,q-1)+n+pow(2,P-1)),0);
Fn((int)(k*pow(2,q)+n+pow(2,q-1)),0)=(temp((int)(k*pow(2,q-1)+n),0)
-temp((int)(k*pow(2,q-1)+n+pow(2,P-1)),0))
*Complex(cos(2*pi*k*pow(2,q-1)/N),-sin(2*pi*k*pow(2,q-1)/N));
}
}
temp=Fn;
}
}
void IFFT(MatrixC &fk,MatrixC &Fn,int N){
const double pi=3.1415926;
int k,n,P,q;
P=(int)(log10(N)/log10(2));
MatrixC temp(Fn);
for(q=1;q<=P;q++)
{
for(k=0;k<(int)(pow(2,(P-q)));k++)
{
for(n=0;n<(int)(pow(2,(q-1)));n++)
{
fk((int)(k*pow(2,q)+n),0)=temp((int)(k*pow(2,q-1)+n),0)
+temp((int)(k*pow(2,q-1)+n+pow(2,P-1)),0);
fk((int)(k*pow(2,q)+n+pow(2,q-1)),0)=(temp((int)(k*pow(2,q-1)+n),0)
-temp((int)(k*pow(2,q-1)+n+pow(2,P-1)),0))
*Complex(cos(2*pi*k*pow(2,q-1)/N),sin(2*pi*k*pow(2,q-1)/N));
}
}
temp=fk;
}
fk=fk*(1.0/N);
}
void FFT2(MatrixC &fk,MatrixC &Fn,int N){
MatrixC ifk(fk);
int n=(int)(log10(N)/log10(2));
int i;
//inverse the order of the index of fk
int order,inverse,temp;
for(order=0;order<N;order++){
inverse=0;
for(i=0;i<n;i++){
temp=(order>>i)&1;
temp=temp<<(n-1-i);
inverse=inverse|temp;
}
ifk(order,0)=fk(inverse,0);
}
//FFT algorithm like butterfly
Fn=ifk;
int j,k,kk,d,w;
const double pi=3.1415926;
Complex W,Fn_temp;
for(i=1;i<=n;i++){
w=(int)pow(2,i);
d=w/2;
for(j=0;j<d;j++){
k=j;
W=Complex(cos(-2*pi/pow(2,i)*j),sin(-2*pi/pow(2,i)*j));
while(k<N){
kk=k+d;
Fn_temp=Fn(k,0);
Fn(k,0)=Fn(k,0)+W*Fn(kk,0);
Fn(kk,0)=Fn_temp-W*Fn(kk,0);
k+=w;
}
}
}
}
void IFFT2(MatrixC &fk,MatrixC &Fn,int N){
MatrixC iFn(Fn);
int n=(int)(log10(N)/log10(2));
int i;
//inverse the order of the index of fk
int order,inverse,temp;
for(order=0;order<N;order++){
inverse=0;
for(i=0;i<n;i++){
temp=(order>>i)&1;
temp=temp<<(n-1-i);
inverse=inverse|temp;
}
iFn(order,0)=Fn(inverse,0);
}
//FFT algorithm like butterfly
fk=iFn;
int j,k,kk,d,w;
const double pi=3.1415926;
Complex W,fk_temp;
for(i=1;i<=n;i++){
w=(int)pow(2,i);
d=w/2;
for(j=0;j<d;j++){
k=j;
W=Complex(cos(-2*pi/pow(2,i)*j),-sin(-2*pi/pow(2,i)*j));
while(k<N){
kk=k+d;
fk_temp=fk(k,0);
fk(k,0)=fk(k,0)+W*fk(kk,0);
fk(kk,0)=fk_temp-W*fk(kk,0);
k+=w;
}
}
}
fk=fk*(1.0/N);
}
[ 本帖最后由 花如月 于 2007-12-17 18:13 编辑 ] 其实现在数学工具强大，我当时自己编的原因：一是对算法本身感兴趣，二是考验下自己的C++应用能力。收获是：我现在对FFT的理论和物理意义都比较了解了。呵呵。
有点难以看懂，我现在都看不懂了，当时注释写太少了（不好的编程习惯），不好意思！
[ 本帖最后由 hyl2323 于 2007-12-17 18:06 编辑 ] 学习了，多谢分享！！
页:
[1]