又苦看了两天书,倒是出来了
作了三个序列的实验,均作的8点的:
x1=2*i+j0 // 0<=i<=7
x2=3*i=j0 // 0<=i<=7
x3=2*i=j3*i // 0<=i<=7
fft(x1,8)的结果
56.0000
-8.0000+19.3137j
-8.0000+8.0000j
-8.0000+3.3137j
-8.0000
-8.0000-3.3137j
-8.0000-8.0000j
-8.0000-19.3137j
fft(x2,8)的结果
84.0000
-12.0000+28.9706j
-12.0000+12.0000j
-12.0000+4.9706j
-12.0000
-12.0000-4.9706j
-12.0000-12.0000j
-12.0000-28.9706j
fft(x3,8)的结果经分离后的两个序列
56.0000
-8.0000+19.3137j
-8.0000+8.0000j
-8.0000+3.3137j
-8.0000
-8.0000-3.3137j
-8.0000-8.0000j
-8.0000-19.3137j
//******************
84.0000
-12.0000-28.9706j
-12.0000-12.0000j
-12.0000-4.9706j
-12.0000
-12.0000+4.9706j
-12.0000+12.0000j
-12.0000+28.9706j
仿真用的程序:
声明:FFT主程序是在网上找的经修改后的程序,作者的相关信息均在哈
/*时间抽选基2FFT及IFFT算法C语言实现*/
/*Author :Junyi Sun*/
/*Copyright 2004-2005*/
/*Mail:ccnusjy@yahoo.com.cn*/
#include <iostream.h>
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
/*定义复数类型*/
typedef struct
{
float real;
float img;
}complex;
void fft(complex *dt,int n); /*快速傅里叶变换*/
// void ifft(); /*快速傅里叶反变换*/
void initW(complex *W,int n); /*初始化变换核*/
void depart(complex *x1,complex *x2,int n);
void add(complex ,complex ,complex *); /*复数加法*/
void mul(complex ,complex ,complex *); /*复数乘法*/
void sub(complex ,complex ,complex *); /*复数减法*/
void output(complex *dt,int n); /*输出结果*/
double PI; /*圆周率*/
double PI2;// 2*PI
int main()
{
int i,method;
int size_x=0; /*输入序列的大小，在本程序中仅限2的次幂*/
double a;
complex *x,*y;
system("cls");
//PI=atan(1)*4;
PI=3.141592653589793;
PI2=6.283185307179586;
cout<<"print N=
";
cin>>size_x;
a=PI2/size_x;
x=new complex;
y=new complex;
if(x==NULL || y==NULL)
exit(0);
cout<<"x.real=2*i x.img=0
";
for(i=0;i<size_x;i++)
{
x.real=2*i;
x.img=0;
}
fft(x,size_x);
output(x,size_x);
cout<<"any key NO. ";
cin>>i;
cout<<"x.real=3*i x.img=0
";
for(i=0;i<size_x;i++)
{
x.real=3*i;
x.img=0;
}
fft(x,size_x);
output(x,size_x);
cout<<"any key NO. ";
cin>>i;
cout<<"2*i+j3*i
";
for(i=0;i<size_x;i++)
{
x.real=2*i;
x.img=3*i;
}
fft(x,size_x);
depart(x,y,size_x); //数据分离
output(x,size_x);
output(y,size_x);
delete[] y;
delete[] x;
return 0;
}
/*快速傅里叶变换*/
void fft(complex *dt,int n)
{
int i=0,j=0,k=0,r=0;
complex up,down,product,*W;
double t;
/*变址计算，将x(n)码位倒置*/
//***********************
for(i=0;i<n;i++)
{
k=i;j=0;
t=(log(n)/log(2));
while( (t--)>0 )
{
j=j<<1;
j|=(k & 1);
k=k>>1;
}
if(i<j)
{
product=dt;
dt=dt[j];
dt[j]=product;
}
}
//*********************
W=new complex[n];
if(W==NULL)
exit(0);
initW(W,n);
//**********************************
i=j=k=0;
t=(log(n)/log(2));
for(i=0;i<t ;i++)
{ //*一级蝶形运算
r=1<<i;//取 r
for(j=0;j<n;j+= 2*r )
{ //*一组蝶形运算
for(k=0;k<r;k++)
{ //*一个蝶形运算
mul(dt[j+k+r],W[n*k/2/r],&product);
add(dt[j+k],product,&up);
sub(dt[j+k],product,&down);
dt[j+k]=up;
dt[j+k+r]=down;
}
}
}
delete[] W;
}
void depart(complex *x1,complex *x2,int n)
{
float temp1_data_r,temp1_data_i,temp2_data_r,temp2_data_i;
int i;
complex *s_x1;
s_x1=new complex[n];
if(s_x1==NULL)
exit(0);
for(i=0;i<n;i++)
s_x1=x1;
for(i=1;i<n;i++)
{
temp1_data_r=(s_x1.real+s_x1[n-i].real)*0.5;
temp2_data_i=(s_x1.real-s_x1[n-i].real)*0.5;
temp2_data_r=(s_x1.img+s_x1[n-i].img)*0.5;
temp1_data_i=(s_x1.img-s_x1[n-i].img)*0.5;
x1.real=temp1_data_r;
x1.img=temp1_data_i;
x2.real=temp2_data_r;
x2.img=temp2_data_i;
}
x1[0].real=s_x1[0].real;
x1[0].img=0;
x2[0].real=s_x1[0].img;
x2[0].img=0;
delete[] s_x1;
}
/*初始化变换核,cos,sin表*/
void initW(complex *W,int n)
{
/*输入序列,变换核*/
int i;
double a;
a=PI2/n;
for(i=0;i<n;i++)
{
W.real=cos(a*i);
W.img=-1*sin(a*i);
}
}
/*输出傅里叶变换的结果*/
void output(complex *dt,int n)
{
int i;
printf("The result are as follows
");
for(i=0;i<n;i++)
{
printf("%.4f",dt.real);
if(dt.img>=0.0001)
printf("+%.4fj
",dt.img);
else if(fabs(dt.img)<0.0001)
printf("
");
else
printf("%.4fj
",dt.img);
}
}
void add(complex a,complex b,complex *c)
{
c->real=a.real+b.real;
c->img=a.img+b.img;
}
void mul(complex a,complex b,complex *c)
{
c->real=a.real*b.real - a.img*b.img;
c->img=a.real*b.img + a.img*b.real;
}
void sub(complex a,complex b,complex *c)
{
c->real=a.real-b.real;
c->img=a.img-b.img;
}
//可能FFT程序中还有一些BUG,请指教一下
|
楼主
标题置顶
标题高亮
