原创 c语言快速傅里叶变换

#include<stdio.h>
#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;i<length;i++)     //实现倒叙排列
{
     k="i";
  j=0;
  for(m=0;m<jishu;m++)          
  {
   j="j"*2+k%2;
   k/=2;
  }
     if(i<j)
  {
     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;j<kind;j++)  //每种蝶形运算的起始点为j,共有kind种
    {
     for(i=j;i<length;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";i<8;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;i<8;i++)             //输出结果
      printf("X(%d)=  %lf  +   %lf i\n",i+1,x.GetR(),x.GetI());
}