#includestdio.h
我们提供的服务有：成都网站设计、成都网站制作、微信公众号开发、网站优化、网站认证、额济纳ssl等。为上千余家企事业单位解决了网站和推广的问题。提供周到的售前咨询和贴心的售后服务，是有科学管理、有技术的额济纳网站制作公司
#include math.h
class complex //定义一个类,实现复数的所有操作
{
double Real,Image; //实部与虚部
public:
complex(double r="0",double i="0"){Real=r;Image=i;}
double GetR(){return Real;} //取出实部
double GetI(){return Image;} //取出虚部
complex operator + (complex ); //复数加法
complex operator - (complex ); //复数减法
complex operator * (complex ); //复数乘法
void operator =(complex ); //复数 赋值
};
complex complex::operator + (complex c) //复数加法
{
complex t;
t.Real=Real+c.Real;
t.Image=Image+c.Image;
return t;
}
complex complex::operator - (complex c) //复数减法
{
complex t;
t.Real=Real-c.Real;
t.Image=Image-c.Image;
return t;
}
complex complex::operator * (complex c) //复数乘法
{
complex t;
t.Real=Real*c.Real-Image*c.Image;
t.Image=Real*c.Image+Image*c.Real;
return t;
}
void complex::operator = (complex c) //复数 赋值
{
Real=c.Real;
Image=c.Image;
}
void fft(complex a[],int length,int jishu) //实现fft的函数
{
const double PI="3".141592653589793;
complex u,Wn,t;
int i,j,k,m,kind,distance,other;
double tmp;
for(i=0;ilength;i++) //实现倒叙排列
{
k="i";
j=0;
for(m=0;mjishu;m++)
{
j="j"*2+k%2;
k/=2;
}
if(ij)
{
t="a";
a=a[j];
a[j]=t;
}
}
for(m=1;m=jishu;m++) //第m级蝶形运算,总级数为jishu
{
kind = (int)pow(2,m-1); //第m级有2^(m-1)种蝶形运算
distance = 2*kind; //同种蝶形结相邻距离为2^m
u=complex(1,0); //旋转因子初始值为 1
tmp=PI/kind;
Wn=complex(cos(tmp),-sin(tmp));//旋转因子Wn
for(j=0;jkind;j++) //每种蝶形运算的起始点为j,共有kind种
{
for(i=j;ilength;i+=distance) //同种蝶形运算
{
other=i+kind;//蝶形运算的两个因子对应单元下标的距离为2^(m-1)
t=a[other]*u; // 蝶形运算的乘积项
a[other]=a-t; //蝶形运算
a=a+t; //蝶形运算
}
u="u"*Wn; //修改旋转因子,多乘一个基本DFT因子WN
}
}
}
void main(void)
{
double a,b;
complex x[8]; //此程序以8点序列测试
printf("8点序列:\n");
for(int i="0";i8;i++) //初始化并输出原始序列
{
x=complex(i,i+1);
printf("x(%d) = %lf + %lf i\n",i+1,x.GetR(),x.GetI());
}
fft(x,8,3); //调用fft函数
printf("fft变换的结果为:\n");
for(i=0;i8;i++) //输出结果
printf("X(%d)= %lf + %lf i\n",i+1,x.GetR(),x.GetI());
}
我自己做了一个程序，和你想象的不一样，
首先用函数产生一个序列f[n]，然后调用FFT:
void __stdcall FFT(
long N, // Serial Length, N 0 for DFT, N 0 for iDFT - inversed Discrete Fourier Transform
double * inputReal, double * inputImaginary, // inputs
double * AmplitudeFrequences, double * PhaseFrequences) // outputs
比如 FFT(n, input, 0, FA, FP)
然后用printf把f[n]，FA[n]和FP[n]打印出来，生成一个文本文件，这个文件可以直接粘到EXCEL里面去，然后用EXCEL生成图表就一目了然了，非常清楚精准
另一种方法更直接，就是我把FFT做成了一个动态链接库wfft.dll，然后你打开EXCEL，在第一列产生一个自动生成的函数值，比如cos(2pi*w) + cos(16pi *w)，然后用宏调用我这个动态链接库，就在另外两列自动生成了幅频和相频数列，选择这两个序列就可以自动生成曲线和图表了
需要的话我可以把我做的样例FFT/EXCEL发给你
float ar[1024],ai[1024];/* 原始数据实部，虚部 */
float a[2050];
void fft(int nn) /* nn数据长度 */
{
int n1,n2,i,j,k,l,m,s,l1;
float t1,t2,x,y;
float w1,w2,u1,u2,z;
float fsin[10]={0.000000,1.000000,0.707107,0.3826834,0.1950903,0.09801713,0.04906767,0.02454123,0.01227154,0.00613588,};
float fcos[10]={-1.000000,0.000000,0.7071068,0.9238796,0.9807853,0.99518472,0.99879545,0.9996988,0.9999247,0.9999812,};
switch(nn)
{
case 1024: s=10; break;
case 512: s=9; break;
case 256: s=8; break;
}
n1=nn/2; n2=nn-1;
j=1;
for(i=1;i=nn;i++)
{
a[2*i]=ar[i-1];
a[2*i+1]=ai[i-1];
}
for(l=1;ln2;l++)
{
if(lj)
{
t1=a[2*j];
t2=a[2*j+1];
a[2*j]=a[2*l];
a[2*j+1]=a[2*l+1];
a[2*l]=t1;
a[2*l+1]=t2;
}
k=n1;
while (kj)
{
j=j-k;
k=k/2;
}
j=j+k;
}
for(i=1;i=s;i++)
{
u1=1;
u2=0;
m=(1i);
k=m1;
w1=fcos[i-1];
w2=-fsin[i-1];
for(j=1;j=k;j++)
{
for(l=j;lnn;l=l+m)
{
l1=l+k;
t1=a[2*l1]*u1-a[2*l1+1]*u2;
t2=a[2*l1]*u2+a[2*l1+1]*u1;
a[2*l1]=a[2*l]-t1;
a[2*l1+1]=a[2*l+1]-t2;
a[2*l]=a[2*l]+t1;
a[2*l+1]=a[2*l+1]+t2;
}
z=u1*w1-u2*w2;
u2=u1*w2+u2*w1;
u1=z;
}
}
for(i=1;i=nn/2;i++)
{
ar[i]=4*a[2*i+2]/nn; /* 实部 */
ai[i]=-4*a[2*i+3]/nn; /* 虚部 */
a[i]=4*sqrt(ar[i]*ar[i]+ai[i]*ai[i]); /* 幅值 */
}
}