Question

请教大虾一个fft的问题

Answer 1

把你的程序发过来看下，时序FFT算法还是挺简单的啊。
我这有两个，一个是在网上搜的，一个是我自己写的。
复数这个类自己写一个就可以了。

// FFT_t.cpp : 定义控制台应用程序的入口点。
//
#include "stdafx.h"
#include "complex.h"
#include "cmath"
#include "iostream"
#include "stdio.h"
#include "fstream"
using namespace std;
#include "stdio.h"
#define NUM 128

int changeadd(int k,int m)//倒序方法一，直接算出其值
{
//m是二进制数的位数
//这个地主要注意的是当m很大时，K1可能会超出范围
int k1=0;
for(int i=0;i<m;i++)
{
if(k&(int)pow(2.0,i))
k1+=pow(2.0,m-i-1);
}
return k1;
}
void FFT(complex *in,complex * out,int M)
{
ofstream file("fft.dat");
int N=pow(2.0,M);
out[0]=in[0];
out[N-1]=in[N-1];
complex t;
for(int i=1;i<N-1;i++)
{
if(i<changeadd(i,M))
{
t=in[changeadd(i,M)];
out[changeadd(i,M)]=in[i];
out[i]=t;
}
}
file<<"in"<<endl;
for(int i=0;i<N;i++)
file<<in[i];
//cout<<"out"<<endl;
//for(int i=0;i<N;i++)
// cout<<out[i];
for(int L=1;L<=M;L++)//L第几层蝶形，
{
int B=pow(2.0,L-1);//B同一个蝶形的两个点间的距离

for(int J=0;J<=B-1;J++)
{
int P=pow(2.0,M-L)*J;

complex temp(cos(2*PI*P/N),-sin(2*PI*P/N));
file<<temp;
for(int K=J;K<者族此=N-1;K+=pow(2.0,L))
{
//if(L==1&&J==0&&K==2){__asm int 3}
out[K]=out[K]+out[K+B]*temp;
out[K+B]=out[K]-out[K+B]*temp*2;
}
}
//观察蝶形的具体样子
//file<<"第"<<L<<"层"<<endl;
//for(int i=0;i<N;i++)
//{
// file<<out[i];
//}
//file<<endl;
}
file.close();
}
//////////////////////////////////////////////////////
////在网上找到的一种算法
void _FFT(complex Input[],int Length,int isign)
{
//isign代表是正还是逆变换 -1表示正变换，1表示逆变换首迅
ofstream _file("_fft.dat");
int l,i,m,mr=0;
complex t;
float tm,pisign=isign*PI;
for(m=1;m<Length;m++)/穗缺/倒序
{
l=Length>>1;
while(mr+l>=Length) l>>=1;
mr=mr%l+l;
if(mr>m)
{
t=Input[m];
Input[m]=Input[mr];
Input[mr]=t;
}
}
_file<<"in"<<endl;
for(int i=0;i<Length;i++)
_file<<Input[i];
/////蝶形计算
l=1;
while(l<Length)
{
for(m=0;m<l;m++)
{
tm=pisign/l*m;
_file<<complex(cos(tm),sin(tm));
for(i=m;i<Length;i+=l<<1)
{
//if(l==1&&m==0&&i==2){__asm int 3}
t=Input[i+l]*complex(cos(tm),sin(tm));
Input[i+l]=Input[i]-t;
Input[i]+=t;
}
}
l<<=1;
// _file<<"第"<<l<<"层"<<endl;
// for(int i=0;i<Length;i++)
// {
//
// _file<<Input[i];
// }
// _file<<endl;
}
if(isign==1) for(l=0;l<Length;l++) Input[l]/=Length;
_file.close();
}
//////////////////////////////////////
int _tmain(int argc, _TCHAR* argv[])
{
complex in[NUM];
complex _in[NUM];
double t;
int x=log(NUM*1.0)/log(2.0);
for(int i=0;i<NUM;i++)
{
t=sin(0.01*2*PI*i)+sin(0.02*2*PI*i)+sin(0.4*2*PI*i);
_in[i]=in[i]=complex(t,0);
}
FFT(in,in,x);
_FFT(_in,NUM,-1);
//for(int i=0;i<NUM;i++)
// cout<<in[i];
//int x;
//while(1)
//{
// cin>>x;
// if(x==0)break;
// //changeadd;
// printf("%x\n",x);
//printf("%x\n",changeadd(2,7));
//}
return 0;
}
里边的调试代码我也没有删除，能帮助你理解