fft的C语言源代码
fft c语言程序
fft c语言程序以下是一个基于C语言的FFT（快速傅里叶变换）程序示例：```c#include <stdio.h>#include <math.h>#define PI 3.14159265void FFT(double complex *x, int n ) {if ( n <= 1 ) return;double complex even ( n /2 );double complex odd ( n /2 );for ( int i =0 ; i < n / 2 ; i++) {even (i) = x[2 * i];odd (i) = x[2 * i + 1 ];FFT(even , n /2 );FFT(odd , n /2 );}for ( int i =0 ; i < n / 2 ; i++) {double complex t = cexp(-I *2 *PI * i / n ) * odd (i); x[i] = even (i) + t ;x[i + n /2 ] = even (i) - t ;}}int main() {double complex x[] = {1 , 2 , 3 , 4 , 5 , 6 , 7 , 8 };int n = sizeof(x) / sizeof(double complex );FFT(x, n );for ( int i =0 ; i < n ; i++) {printf( "%f + %fi\n" , creal(x[i]), cimag(x[i]));}return 0 ;}```这个程序使用了递归的方法来实现快速傅里叶变换。
它接受一个复数数组`x`作为输入，并计算其傅里叶变换。
在主函数中，定义了一个复数数组`x`，并调用`FFT`函数对其进行变换。
最后，输出变换后的结果。
fft傅立叶变换C语言
fft傅立叶变换C语言#include "stdio.h"#include <math.h> //調用源（頭）文件struct compx /*定義一個複數結構*/{float real;float imag;};struct compx s[ 257 ]; //FFT輸入輸出均從是s[1]開始存入struct compx EE(struct compx,struct compx); //定義複數相乘結構void FFT(struct compx xin,int N); /*定義FFT函數*/struct compx EE(struct compx a1,struct compx b2) //兩複數相乘的程序{struct compx b3; //b3保存兩複數間的結果b3.real=a1.real*b2.real-a1.imag*b2.imag; //兩複數間的運算b3.imag=a1.real*b2.imag+a1.imag*b2.real;return(b3); /*返回結果*/}void FFT(struct compx xin,int N) /*FFT函數體*/{int f,m,nv2,nm1,i,k,j=1,l; /*定義變量*/struct compx v,w,t; /*定義結構變量*/nv2=N/2; /*最高位值的權值*/f=N; /*f為中間變量*/for(m=1;(f=f/2)!=1;m++){;} /*求級數m*/nm1=N-1; /*nm1為數組長度*/for(i=1;i<=nm1;i++) /*倒序*/{if(i<j) {t=xin[ j ];xin[j]=xin[ i ];xin[ i ] =t;} /*i<j則換位*/k=nv2; /*k為倒序中相應位置的權值*/while(k<j) {j=j-k;k=k/2;} /*k<j時最高為變為0*/j=j+k; /* j為數組中的位數，是一個十進制數*/}{int le,lei,ip; //變量初始化,le為序列長度float pi;for(l=1;l<=m;l++) /*l控制級數*/{le=pow(2,l); /*le等於2的l次方*/lei=le/2; /*蝶形兩節點間的距離*/pi=3.14159265;v.real=1.0; // 此次的v運於複數的初始化v.imag=0.0; w.real=cos(pi/lei); /*旋轉因子*/w.imag=-sin(pi/lei);for(j=1;j<=lei;j++) //外循環控制蝶行運算的級數{for(i=j;i<=N;i=i+le) //內循環控制每級間的運算次數{ip=i+lei; /*蝶形運算的下一個節點*/t=EE(xin[ ip ],v); /*第一個旋轉因子*/xin[ ip ].real=xin[ i ].real-t.real; /*蝶形計算*/xin[ ip ].imag=xin[ i ].imag-t.imag;xin[ i ].real=xin[ i ].real+t.real;xin[ i ].imag=xin[ i ].imag+t.imag;}v=EE(v,w); //調用EE複數相乘程序，結果給下次的循環}}}}main() /*定義主函數*/{int N,i; //變量初始化，N為總點數，i為每點數printf("shu ru N de ge shu N="); /*提示輸入*/scanf("%d",&N); /*輸入N*/ for(i=1;i<=N;i++) /*輸入*/{printf("di %d ge shu real=",i);getchar();scanf("%f",&s[ i ].real);getchar();printf("\n");printf("di %d ge shu imag=",i);scanf("%f",&s[ i ].imag);printf("\n");}FFT(s,N); /*調用FFt*/for(i=1;i<=N;i++) /*輸出*/{printf("%f",s[ i ].real);printf(" + ");printf("%f",s[ i ].imag);printf("j");printf("\n");}}。
c语言快速傅里叶变换处理数据
c语言快速傅里叶变换处理数据快速傅里叶变换(FFT)是一种常见的信号处理算法，可以在O(nlogn)的时间内处理一个长度为n的序列。
在c语言中，可以使用外部库来进行FFT计算。
下面是一个使用fftw库进行FFT处理的示例代码：```c#include <fftw3.h> // 使用fftw库进行FFTint main(){int n = 8;double signal[] = {0.707, 1.0, 0.707, 0.0, -0.707, -1.0, -0.707, 0.0}; // 输入信号fftw_complex* in = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * n);fftw_complex* out = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * n);for (int i = 0; i < n; i++) {in[i][0] = signal[i]; // 实部in[i][1] = 0; // 虚部设为0}fftw_plan p = fftw_plan_dft_1d(n, in, out,FFTW_FORWARD, FFTW_ESTIMATE);fftw_execute(p);for (int i = 0; i < n; i++) {printf("%f + %fj\n", out[i][0], out[i][1]); // 输出FFT结果}fftw_destroy_plan(p);fftw_free(in);fftw_free(out);return 0;}```此代码中，输入的信号长度为8，信号值存储在signal数组中。
首先需要分配内存来存储FFT所需的输入和输出数据。
然后，使用fftw_plan_dft_1d函数创建一个FFT计算计划，并使用fftw_execute 函数执行FFT计算。
C语言实现FFT(快速傅里叶变换)
#include <iom128.h>#include <intrinsics.h>/*********************************************************************快速福利叶变换C函数函数简介：此函数是通用的快速傅里叶变换C语言函数，移植性强，以下部分不依赖硬件。
此函数采用联合体的形式表示一个复数，输入为自然顺序的复数（输入实数是可令复数虚部为0），输出为经过FFT变换的自然顺序的复数使用说明：使用此函数只需更改宏定义FFT_N的值即可实现点数的改变，FFT_N的应该为2的N次方，不满足此条件时应在后面补0函数调用：FFT(s);时间：2010-2-20版本：Ver1.0参考文献：**********************************************************************/#include<math.h>#define PI 3.1415926535897932384626433832795028841971 //定义圆周率值#define FFT_N 128 //定义福利叶变换的点数struct compx {float real,imag;}; //定义一个复数结构struct compx s[FFT_N]; //FFT输入和输出：从S[1]开始存放，根据大小自己定义/*******************************************************************函数原型：struct compx EE(struct compx b1,struct compx b2)函数功能：对两个复数进行乘法运算输入参数：两个以联合体定义的复数a,b输出参数：a和b的乘积，以联合体的形式输出*******************************************************************/struct compx EE(struct compx a,struct compx b){struct compx c;c.real=a.real*b.real-a.imag*b.imag;c.imag=a.real*b.imag+a.imag*b.real;return(c);}/*****************************************************************函数原型：void FFT(struct compx *xin,int N)函数功能：对输入的复数组进行快速傅里叶变换（FFT）输入参数：*xin复数结构体组的首地址指针，struct型*****************************************************************/void FFT(struct compx *xin){int f,m,nv2,nm1,i,k,l,j=0;struct compx u,w,t;nv2=FFT_N/2; //变址运算，即把自然顺序变成倒位序，采用雷德算法nm1=FFT_N-1;for(i=0;i<nm1;i++){if(i<j) //如果i<j,即进行变址{t=xin[j];xin[j]=xin[i];xin[i]=t;}k=nv2; //求j的下一个倒位序while(k<=j) //如果k<=j,表示j的最高位为1{j=j-k; //把最高位变成0k=k/2; //k/2，比较次高位，依次类推，逐个比较，直到某个位为0 }j=j+k; //把0改为1}{int le,lei,ip; //FFT运算核，使用蝶形运算完成FFT运算f=FFT_N;for(l=1;(f=f/2)!=1;l++) //计算l的值，即计算蝶形级数;for(m=1;m<=l;m++) // 控制蝶形结级数{ //m表示第m级蝶形，l为蝶形级总数l=log（2）Nle=2<<(m-1); //le蝶形结距离，即第m级蝶形的蝶形结相距le点lei=le/2; //同一蝶形结中参加运算的两点的距离u.real=1.0; //u为蝶形结运算系数，初始值为1u.imag=0.0;w.real=cos(PI/lei); //w为系数商，即当前系数与前一个系数的商w.imag=-sin(PI/lei);for(j=0;j<=lei-1;j++) //控制计算不同种蝶形结，即计算系数不同的蝶形结{for(i=j;i<=FFT_N-1;i=i+le) //控制同一蝶形结运算，即计算系数相同蝶形结{ip=i+lei; //i，ip分别表示参加蝶形运算的两个节点t=EE(xin[ip],u); //蝶形运算，详见公式xin[ip].real=xin[i].real-t.real;xin[ip].imag=xin[i].imag-t.imag;xin[i].real=xin[i].real+t.real;xin[i].imag=xin[i].imag+t.imag;}u=EE(u,w); //改变系数，进行下一个蝶形运算}}}}/************************************************************函数原型：void main()函数功能：测试FFT变换，演示函数使用方法输入参数：无输出参数：无************************************************************/void main(){int i;for(i=0;i<FFT_N;i++) //给结构体赋值{s[i].real=sin(2*3.141592653589793*i/FFT_N); //实部为正弦波FFT_N点采样，赋值为1 s[i].imag=0; //虚部为0}FFT(s); //进行快速福利叶变换for(i=0;i<FFT_N;i++) //求变换后结果的模值，存入复数的实部部分s[i].real=sqrt(s[i].real*s[i].real+s[i].imag*s[i].imag);while(1);}#include <iom128.h>#include <intrinsics.h>/*********************************************************************快速福利叶变换C程序包函数简介：此程序包是通用的快速傅里叶变换C语言函数，移植性强，以下部分不依赖硬件。
FFT及IFFTC语言实现
FFT及IFFTC语言实现FFT（快速傅里叶变换）和IFFT（快速傅里叶逆变换）是傅里叶变换在计算机科学和信号处理中的高效实现方法。
它们广泛应用于图像处理、音频处理、通信系统等领域。
下面将详细介绍FFT和IFFT的C语言实现。
首先，让我们了解一下DFT（离散傅里叶变换）。
DFT将一个离散的时间域序列转换为离散的频域序列，其定义如下：其中，N表示序列的长度，x(n)是输入序列，X(k)是输出序列。
FFT是DFT的一种高效实现方法，它利用了序列的对称性质，将操作的复杂度从O(N^2)降低到O(NlogN)。
IFFT则是FFT的逆过程，可以将频域序列恢复为时间域序列。
以下是FFT的C语言实现代码：```c#include <stdio.h>#include <math.h>typedef structdouble real;double imag;result.real = a.real * b.real - a.imag * b.imag;result.imag = a.real * b.imag + a.imag * b.real;return result;result.real = a.real + b.real; result.imag = a.imag + b.imag; return result;result.real = a.real - b.real; result.imag = a.imag - b.imag; return result;if (N <= 1)return;}for (int i = 0; i < N / 2; i++) even[i] = x[2 * i];odd[i] = x[2 * i + 1];}fft(even, N / 2);fft(odd, N / 2);for (int k = 0; k < N / 2; k++)x[k] = add(even[k], t);x[k + N / 2] = subtract(even[k], t); }//完整的IFFT实现代码略，与FFT相似int mai//输入序列x(n)int N = sizeof(x) / sizeof(x[0]);//调用FFTfft(x, N);//输出频域序列X(k)for (int k = 0; k < N; k++)printf("X(%d) = %.2f + %.2fi\n", k, x[k].real, x[k].imag);}//调用IFFT//...return 0;```IFFT的实现与FFT类似，只需将其中几处计算公式略作变换即可。
FFT算法的C语言编程
#include <stdio.h>#include <math.h>#include <stdlib.h>#define PI 3.1415926/*This program is used to conduct FFT or IFFT transform based on 2 discomposition data can be loaded by typing on the keyboard or from file.*/typedefstruct{double real;doubleimag;}complex;long NN;long t2t(long n, int M){longnn=0;for(inti=0;i<M;i++){nn=nn+n%2*pow(2,M-1-i);n=n/2;}returnnn;}intdefM(long n){int M=0;while(pow(2,M)<n){M++;}return M;}voiddatain(complex x[],long N){int flag;longi;char s[20];FILE *fpi;doublea,b;printf("Load data from keybord(input 0) or file(input 1)\n");scanf("%d",&flag);if(flag==0){printf("Please enter the data in the form of a,b--a+j*b ending with'*'\n");for(i=0;getchar()!='*';i++)scanf("%lf%lf",&x[i].real,&x[i].imag);while(i<N){x[i].real=0;x[i].imag=0;i++;}}elseif(flag==1){printf("Please input the name of the file with its File type suffix-(eg:data.txt)\n");gets(s);gets(s);//To avoid the quick response to enterfpi=fopen(s,"r");if(!fpi){printf("Open file error");}i=0;while(!feof(fpi)){fscanf(fpi,"%lf",&a);fscanf(fpi,"%lf",&b);x[i].real=a;x[i].imag=b;i++;}fclose(fpi);while(i<N){x[i].real=0;x[i].imag=0;i++;}}}voiddataout(complex x[],long N){longi;for(i=0;i<N;i++)printf("%.3lf+j*%.3lf\n",x[i].real,x[i].imag);}void reverse(complex x[],complex xi[],long n,int m) {longi;long temp;for(i=0;i<n;i++){temp=t2t(i,m);xi[i].real=x[temp].real;xi[i].imag=x[temp].imag;}}complexcom_add(complex x1,complex x2){complex y;y.real=x1.real+x2.real;y.imag=x1.imag+x2.imag;return y;}complexcom_sub(complex x1,complex x2){complex y;y.real=x1.real-x2.real;y.imag=x1.imag-x2.imag;return y;}complexcom_mul(complex x1,complex x2){complex y;y.real=x1.real*x2.real-x1.imag*x2.imag;y.imag=x1.real*x2.imag+x2.real*x1.imag;return y;}complexcom_div(complex x1,complex x2){complex y;double down=x2.real*x2.real+x2.imag*x2.imag;y.real=(x1.real*x2.real+x1.imag*x2.imag)/down;y.imag=(x2.real*x1.imag-x1.real*x2.imag)/down;return y;}voidfft(complex x[],long NN, int M){intm,i,j;longNw,Ni;complexWm,temp;for(m=0;m<M;m++){Nw=pow(2,m+1);Ni=pow(2,m);for(i=0;i<Ni;i++){for(j=i;j<NN;j=j+Nw){Wm.real=cos(2*i*PI/Nw);Wm.imag=-sin(2*i*PI/Nw);temp=com_mul(x[j+Ni],Wm);x[j+Ni]=com_sub(x[j],temp);x[j]=com_add(x[j],temp);}}}}voidifft(complex x[],long NN, int M){intm,i,j;longNw,Ni;complexWm,temp;for(m=0;m<M;m++){Nw=pow(2,m+1);Ni=pow(2,m);for(i=0;i<Ni;i++){for(j=i;j<NN;j=j+Nw){Wm.real=cos(2*i*PI/Nw);Wm.imag=-sin(2*i*PI/Nw);temp=com_mul(x[j+Ni],Wm);x[j+Ni]=com_sub(x[j],temp);x[j]=com_add(x[j],temp);}}}for(i=0;i<NN;i++){x[i].real=x[i].real/NN;x[i].imag=-x[i].imag/NN;}}int main(){longN,i;intM,flag;doublea,b;printf("Please enter the length of the sequence\n");scanf("%ld",&N);M=defM(N);NN=pow(2,M);complexxy[NN];complex xyi[NN];datain(xy,N);for(i=N;i<NN;i++){xy[i].real=0;xy[i].imag=0;}reverse(xy,xyi,NN,M);printf("Please choose FFT(1) or IFFT(-1)\n");scanf("%d",&flag);if(flag==1)fft(xyi,NN,M);elseif(flag==-1)ifft(xyi,NN,M);printf("The initial sequence is:\n");dataout(xy,N);printf("The FFT results are as follows:\n"); dataout(xyi,N);return 0;}。
最简洁的FFT代码(C 实现)
int D,L,J,K; double P;
for(L=1;L<=M;L++)
//M 级碟形运算
{
D=pow(2,(L-1));
for(J=0;J<=D-1;J++)
//每一级有 D 个旋转因子
{
P=J*pow(2,M-L);
double real=cos(2*pi*P/N);
double imag=-sin(2*pi*P/N);
#include<iostream.h> #include<math.h> const int M=3; const int N=8; int change(int e) {
//M N 用常变量 //计算倒序数
int b[M]; int sum=0; for(int i=0;i<M;i++) {
b[M-i-1]=e%2; e=e/2; sum=(sum+b[M-i-1])*2; } return sum/2; } int main() { cout<<"N="<<N<<" "<<"M="<<M<<endl; int a[N]; double A[N][2],B[N][2]; cout<<"原数组为："<<endl; for(int j=0;j<N;j++) { A[j][0]=(j+1)*2; A[j][1]=0; cout<<"A["<<j<<"][0]="<<A[j][0]<<" cout<<"A["<<j<<"][1]="<<A[j][1]<<" if((j+1)%2==0)
FFT算法C语言程序代码(可打印修改)
//wLen——FFT 的长度 Struct complexData{
//定义一个复数结构
float re;
float im;
};
Void gen_w_r2(struct complexData *twiFactor,int wLen)
{
int iFactor;
float stepFactor;
stepFactor=2.0*pi/wLen;
short n2,ie,ia,i,j,k,m; float rtemp,itemp,c,s;
n2=n; ie=1; for(k=n;k>1;k>>=1) //loop c
{ n2>>=1;ia=0; for(j=0;j<ie;j++) //loop b { c=w[2*j]; s=w[2*j+1]; for(i=0;i<n2;i++) //loop a { m=ia+n2; rtemp=c*x[2*m]+s*x[2*m+1]; itemp= c*x[2*m+1]-s*x[2*m]; x[2*m]=x[2*ia]-rtemp; x[2*m+1]=x[2*ia+1]-itemp; x[2*ia]=x[2*ia]+rtemp; x[2*ia+1]=x[2*ia+1]+itemp; ia++; } ia+=n2; } ie<<=1; }
}
}
2、在运行 FFT 之前，对输入序列进行倒序变换，代码如下：
//bitRevData——指向位变换序列的指针
//revLen——FFT 长度 Void bit_rev(struct complexData *bitRevData,int revLen)
FFT最详细的源代码和解释
我自己的一些详细标注，有利于深入了解FFT，后面附加几位网友对FFT的理解及源代码，让广大朋友更迅速的掌握FFT#include "DSP281x_Device.h" // DSP281x Headerfile Include File，添加所有头文件#include "DSP281x_Examples.h" // DSP281x Examples Include File，条件编译而已#include "f2812a.h" //一些变量的宏定义而已#include"math.h"#define PI 3.1415926 //前变后常#define SAMPLENUMBER 128//#define SAMPLENUMBER 512void InitForFFT();void MakeWave();//void FFT(float dataR[SAMPLENUMBER],float dataI[SAMPLENUMBER]);int INPUT[SAMPLENUMBER],DATA[SAMPLENUMBER];float fWaveR[SAMPLENUMBER],fWaveI[SAMPLENUMBER],w[SAMPLENUMBER];float sin_tab[SAMPLENUMBER],cos_tab[SAMPLENUMBER];//逐级计算FFT，一级一级递推void FFT(float dataR[SAMPLENUMBER],float dataI[SAMPLENUMBER]){int x0,x1,x2,x3,x4,x5,x6,xx;int i,j,k,b,p,L;float TR,TI,temp;/********** following code invert sequence ************///倒序for ( i=0;i<SAMPLENUMBER;i++ )//就是码位倒置嘛，二进制各个位独立出来再反向{ //128七位二进制表示，/号代表右移嘛x0=x1=x2=x3=x4=x5=x6=0;x0=i&0x01; x1=(i/2)&0x01; x2=(i/4)&0x01; x3=(i/8)&0x01;x4=(i/16)&0x01; x5=(i/32)&0x01; x6=(i/64)&0x01;xx=x0*64+x1*32+x2*16+x3*8+x4*4+x5*2+x6;//最低位，最高位反过来dataI[xx]=dataR[i];}for ( i=0;i<SAMPLENUMBER;i++ ){dataR[i]=dataI[i]; dataI[i]=0; //对应过来}/************** following code FFT *******************/for ( L=1;L<=7;L++ ){ /* for(1) */b=1; i=L-1;/* b的意义非常重大,b表示当前层不同旋转因子的个数*/while ( i>0 ){b=b*2; i--;} /* b= 2^(L-1) */for ( j=0;j<=b-1;j++ ) /* for (2) */{p=1; i=7-L;while ( i>0 ) /* p=pow(2,7-L)*j; */{p=p*2; i--;}p=p*j;for ( k=j;k<128;k=k+2*b ) /* for (3) */{TR=dataR[k]; TI=dataI[k]; temp=dataR[k+b];dataR[k]=dataR[k]+dataR[k+b]*cos_tab[p]+dataI[k+b]*sin_tab[p];dataI[k]=dataI[k]-dataR[k+b]*sin_tab[p]+dataI[k+b]*cos_tab[p];dataR[k+b]=TR-dataR[k+b]*cos_tab[p]-dataI[k+b]*sin_tab[p];dataI[k+b]=TI+temp*sin_tab[p]-dataI[k+b]*cos_tab[p]; //递推嘛，防止立马调用结果，} /* END for (3) */ //引入一个中间变量存原始值，} /* END for (2) */ //防止上一步对下一步的影响} /* END for (1) */for ( i=0;i<SAMPLENUMBER/2;i++ ) //对称性，前半部分即可{w[i]=sqrt(dataR[i]*dataR[i]+dataI[i]*dataI[i]);}} /* END FFT */main(){int i;InitForFFT();MakeWave();for ( i=0;i<SAMPLENUMBER;i++ ){fWaveR[i]=INPUT[i];fWaveI[i]=0.0f;w[i]=0.0f;}FFT(fWaveR,fWaveI);//输入波形进行FFT变换,此处引入起始实参即可递推下去for ( i=0;i<SAMPLENUMBER;i++ ){DATA[i]=w[i];//变换后的波形转换到输出接口}while ( 1 );// break point}//旋转因子事先初始化好，方便调用void InitForFFT(){int i;for ( i=0;i<SAMPLENUMBER;i++ ){sin_tab[i]=sin(PI*2*i/SAMPLENUMBER);cos_tab[i]=cos(PI*2*i/SAMPLENUMBER);//旋转因子事先初始化好，方便调用}}//利用这个，确实能产生各种各样的谐波void MakeWave() //利用这个，确实能产生各种各样的谐波{int i;for ( i=0;i<SAMPLENUMBER;i++ )//1024是相应的幅值嘛，只要弄出个基波，以之为标准即可{INPUT[i]=sin(PI*2*i/SAMPLENUMBER)*1024+sin(PI*2*i/SAMPLENUMBER*3)*1024/3+sin(PI*2*i/SAMPLENUMBER*5)*1024/5+sin(PI*2*i/SAMPLENUMBER*7)*1024/7+sin(PI*2*i/SAMPLENUMBER*9)*1024/9;//INPUT[i]=sin(PI*2*i/SAMPLENUMBER*3)*1024;}}X[0] = x[0]*(1-0) + x[1]*(1-0) = x[0] + 1*x[1]; X[1] = x[0]*(1-0) + x[1]*(-1-0) = x[0] - 1*x[1];这就是单个2点蝶形算法.FFT实现流程图分析(N=8, 以8点信号为例)FFT implementation of an 8-point DFT as two 4-point DFTs and four 2-pointDFTs8点FFT流程图(Layer表示层, gr表示当前层的颗粒)下面以LayerI为例.LayerI部分, 具有4个颗粒, 每个颗粒2个输入(注意2个输入的来源, 由时域信号友情提供, 感谢感谢☺)我们将输入x[k]分为两部分x_r[k], x_i[k]. 具有实部和虚部, 时域信号本没有虚部的, 因此可以让x_i[k]为0. 那么为什么还要画蛇添足分为实部和虚部呢? 这是因为LayerII, LayerIII的输入是复数, 为了编码统一而强行分的.当然你编码时可以判断当前层是否为1来决定是否分. 但是我想每个人最后都会倾向分的.旋转因子 tw = cos(2*PI*k/N)-j*sin(2*PI*k/N); 也可以分为实部和虚部, 令其为tw_r, tw_i;则tw = tw_r - j*tw_i;X[k] = (x_r[k] + j*x_i[k]) + (tw_r–j*tw_i) * (x_r[k+N/2]+j*x_i[k+N/2])则X_R[k] = x_r[k] + tw_r*x_r[k+N/2] + tw_i*x_i[k+N/2];X_I[k] = x_i[k] - tw_i*x_r[k+N/2] + tw_r*x_i[k+N/2];LayerII部分, 具有2个颗粒, 每个颗粒4个输入(注意4个输入的来源, 由LayerI友情提供, 感谢感谢☺)LayerIII部分, 具有1个颗粒, 每个颗粒8个输入(注意8个输入的来源, 由LayerII友情提供, 感谢感谢☺)LayerI, LayerII, LayerIII从左往右, 蝶形信号运算流非常明显!假令输入为x[k], x[k+N/2], 输出为X[k], X[k+N/2]. x[k]分解为x_r[k], x_i[k]部分则该蝶形运算为X[k]= (x_r[k]-j*x_i[k]) +(x_r[k+N/2]-j*x_i[k+N/2])*(cos(2*PI*k/N)-j*sin(2*PI*k/N));再令cos(2*PI*k/N)为tw1, sin(2*PI*k/N)为tw2 则X[k] = (x_r[k]-j*x_i[k]) + (x_r[k+N/2]-j*x_i[k+N/2])*(tw1-j*tw2);X_R[k] = x_r[k] + x_r[k+N/2]*tw1 - x_i[k+N/2]*tw2;X_I[K] = x_i[k]x_r[k] = x_r[k] + x_r[k+b]*tw1 + x_i[k+b]*tw2;x_i[k] = x_i[k] - x_r[k+b]*tw2 + x_i[k+b]*tw1;譬如8点输入x[8]1. 先分割成2部分: x[0], x[2], x[4], x[6] 和 x[1], x[3], x[5], x[7]2. 信号x[0], x[2], x[4], x[6]再分割成x[0], x[4] 和 x[2], x[6]信号x[1], x[3], x[5], x[7]再分割成x[1], x[5] 和 x[3], x[7]3. 无法分割了, 已经分割成2点了☺.如上图:在LayerI的时候, 我们是对2点进行DFT.( 一共4次DFT )输入为 x[0]&x[4]; x[2]&x[6]; x[1]&x[5]; x[3]&x[7]输出为 y[0],y[1]; Y[2],y[3]; Y[4],y[5]; Y[6],y[7];流程:I. 希望将输入直接转换为x[0], x[4], x[2], x[6], x[1], x[5], x[3], x[7]的顺序II. 对转换顺序后的信号进行4次DFT步骤I代码实现/*** 反转算法.这个算法效率比较低!先用起来在说, 之后需要进行优化.*/static void bitrev( void ){int p=1, q, i;int bit_rev[ N ];float xx_r[ N ];bit_rev[ 0 ] = 0;while( p < N ){for(q=0; q<p; q++){bit_rev[ q ] = bit_rev[ q ] * 2;bit_rev[ q + p ] = bit_rev[ q ] + 1;}p *= 2;}for(i=0; i<N; i++) xx_r[ i ] = x_r[ i ];for(i=0; i<N; i++) x_r[i] = xx_r[ bit_rev[i] ];}// ------------------------ 此刻序列x重排完毕------------------------步骤II代码实现int j;float TR; // 临时变量float tw1; // 旋转因子/* 两点DFT */for(k=0; k<N; k+=2){// 两点DFT简化告诉我们tw1=1TR = x_r[k]; // TR就是A, x_r[k+b]就是B.x_r[k] = TR + tw1*x_r[k+b];x_r[k+b] = TR - tw1*x_r[k+b];}在LayerII的时候, 我们希望得到z, 就需要对y进行DFT.y[0],y[2]; y[1],y[3]; y[4],y[6]; y[5],y[7];z[0], z[1]; z[2],z[3]; z[4],z[5]; z[6],z[7];在LayerIII的时候, 我们希望得到v, 就需要对z进行DFT.z[0],z[4]; z[1],z[5]; z[2],z[6]; z[3],z[7];v[0],v[1]; v[2],v[3]; v[4],v[5]; v[6],v[7];准备令输入为x[s], x[s+N/2], 输出为y[s], y[s+N/2]这个N绝对不是上面的8, 这个N是当前颗粒的输入样本总量对于LayerI而言N是2; 对于LayerII而言N是4; 对于LayerIII而言N是8复数乘法:(a+j*b) * (c+j*d)实部 = a*c – bd;虚部 = ad + bc;旋转因子:实现(C描述)#include <stdio.h>#include <math.h>#include <stdlib.h>//#include "complex.h"//--------------------------------------------------------------------------#define N 8 //64#define M 3 //6 //2^m=N#define PI 3.1415926//--------------------------------------------------------------------------float twiddle[N/2] = {1.0, 0.707, 0.0, -0.707};float x_r[N] = {1, 1, 1, 1, 0, 0, 0, 0};float x_i[N]; //N=8/*float twiddle[N/2] = {1, 0.9951, 0.9808, 0.9570, 0.9239, 0.8820, 0.8317, 0.7733,0.7075,0.6349, 0.5561, 0.4721, 0.3835, 0.2912, 0.1961, 0.0991 ,0.0000,-0.0991,-0.1961,-0.2912,-0.3835,-0.472 1,-0.5561,-0.6349,-0.7075,-0.7733,0.8317,-0.8820,-0.9239,-0.9570,-0.9808,-0.9951}; //N=64float x_r[N]={1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,};float x_i[N];*/FILE *fp;// ----------------------------------- func-----------------------------------/*** 初始化输出虚部*/static void fft_init( void ){int i;for(i=0; i<N; i++) x_i[i] = 0.0;}/*** 反转算法.将时域信号重新排序.* 这个算法有改进的空间*/static void bitrev( void ){int p=1, q, i;int bit_rev[ N ]; //float xx_r[ N ]; //bit_rev[ 0 ] = 0;while( p < N ){for(q=0; q<p; q++){bit_rev[ q ] = bit_rev[ q ] * 2;bit_rev[ q + p ] = bit_rev[ q ] + 1;}p *= 2;}for(i=0; i<N; i++) xx_r[ i ] = x_r[ i ];for(i=0; i<N; i++) x_r[i] = xx_r[ bit_rev[i] ]; }/* ------------ add by sshc625 ------------ */static void bitrev2( void ){return ;}/* */void display( void ){printf("/n/n");int i;for(i=0; i<N; i++)printf("%f/t%f/n", x_r[i], x_i[i]);}/****/void fft1( void ){ fp = fopen("log1.txt", "a+");int L, i, b, j, p, k, tx1, tx2;float TR, TI, temp; // 临时变量float tw1, tw2;/* 深M. 对层进行循环. L为当前层, 总层数为M. */for(L=1; L<=M; L++){fprintf(fp,"----------Layer=%d----------/n", L);/* b的意义非常重大,b表示当前层的颗粒具有的输入样本点数 */b = 1;i = L - 1;while(i > 0){b *= 2;i--;}// -------------- 是否外层对颗粒循环, 内层对样本点循环逻辑性更强一些呢! --------------/** outter对参与DFT的样本点进行循环* L=1, 循环了1次(4个颗粒, 每个颗粒2个样本点)* L=2, 循环了2次(2个颗粒, 每个颗粒4个样本点)* L=3, 循环了4次(1个颗粒, 每个颗粒8个样本点)*/for(j=0; j<b; j++){/* 求旋转因子tw1 */p = 1;i = M - L; // M是为总层数, L为当前层.while(i > 0){p = p*2;i--;}p = p * j;tx1 = p % N;tx2 = tx1 + 3*N/4;tx2 = tx2 % N;// tw1是cos部分, 实部; tw2是sin部分, 虚数部分.tw1 = ( tx1>=N/2)? -twiddle[tx1-N/2] : twiddle[ tx1 ]; tw2 = ( tx2>=N/2)? -twiddle[tx2-(N/2)] : twiddle[tx2];/** inner对颗粒进行循环* L=1, 循环了4次(4个颗粒, 每个颗粒2个输入)* L=2, 循环了2次(2个颗粒, 每个颗粒4个输入)* L=3, 循环了1次(1个颗粒, 每个颗粒8个输入)*/for(k=j; k<N; k=k+2*b){TR = x_r[k]; // TR就是A, x_r[k+b]就是B.TI = x_i[k];temp = x_r[k+b];/** 如果复习一下 (a+j*b)(c+j*d)两个复数相乘后的实部虚部分别是什么* 就能理解为什么会如下运算了, 只有在L=1时候输入才是实数, 之后层的* 输入都是复数, 为了让所有的层的输入都是复数, 我们只好让L=1时候的* 输入虚部为0* x_i[k+b]*tw2是两个虚数相乘*/fprintf(fp, "tw1=%f, tw2=%f/n", tw1, tw2);x_r[k] = TR + x_r[k+b]*tw1 + x_i[k+b]*tw2;x_i[k] = TI - x_r[k+b]*tw2 + x_i[k+b]*tw1;x_r[k+b] = TR - x_r[k+b]*tw1 - x_i[k+b]*tw2;x_i[k+b] = TI + temp*tw2 - x_i[k+b]*tw1;fprintf(fp, "k=%d, x_r[k]=%f, x_i[k]=%f/n", k, x_r[k],x_i[k]);fprintf(fp, "k=%d, x_r[k]=%f, x_i[k]=%f/n", k+b, x_r[k+b], x_i[k+b]);} //} //} //}/*** ------------ add by sshc625 ------------* 该实现的流程为* for( Layer )* for( Granule )* for( Sample )*****/void fft2( void ){ fp = fopen("log2.txt", "a+");int cur_layer, gr_num, i, k, p;float tmp_real, tmp_imag, temp; // 临时变量, 记录实部float tw1, tw2;// 旋转因子,tw1为旋转因子的实部cos部分, tw2为旋转因子的虚部sin部分.int step; // 步进int sample_num; // 颗粒的样本总数(各层不同, 因为各层颗粒的输入不同)/* 对层循环 */for(cur_layer=1; cur_layer<=M; cur_layer++){/* 求当前层拥有多少个颗粒(gr_num) */gr_num = 1;i = M - cur_layer;while(i > 0){i--;gr_num *= 2;}/* 每个颗粒的输入样本数N' */sample_num = (int)pow(2, cur_layer);/* 步进. 步进是N'/2 */step = sample_num/2;/* */k = 0;/* 对颗粒进行循环 */for(i=0; i<gr_num; i++){/** 对样本点进行循环, 注意上限和步进*/for(p=0; p<sample_num/2; p++){// 旋转因子, 需要优化...tw1 = cos(2*PI*p/pow(2, cur_layer)); tw2 = -sin(2*PI*p/pow(2, cur_layer));tmp_real = x_r[k+p];tmp_imag = x_i[k+p];temp = x_r[k+p+step];/*(tw1+jtw2)(x_r[k]+jx_i[k])** real : tw1*x_r[k] - tw2*x_i[k]* imag : tw1*x_i[k] + tw2*x_r[k]* 我想不抽象出一个* typedef struct {* double real; // 实部* double imag; // 虚部* } complex; 以及针对complex的操作* 来简化复数运算是否是因为效率上的考虑!*//* 蝶形算法 */x_r[k+p] = tmp_real + ( tw1*x_r[k+p+step] -tw2*x_i[k+p+step] );x_i[k+p] = tmp_imag + ( tw2*x_r[k+p+step] +tw1*x_i[k+p+step] );/* X[k] = A(k)+WB(k)* X[k+N/2] = A(k)-WB(k) 的性质可以优化这里*/// 旋转因子, 需要优化...tw1 = cos(2*PI*(p+step)/pow(2, cur_layer));tw2 = -sin(2*PI*(p+step)/pow(2, cur_layer));x_r[k+p+step] = tmp_real + ( tw1*temp -tw2*x_i[k+p+step] );x_i[k+p+step] = tmp_imag + ( tw2*temp +tw1*x_i[k+p+step] );printf("k=%d, x_r[k]=%f, x_i[k]=%f/n", k+p, x_r[k+p], x_i[k+p]);printf("k=%d, x_r[k]=%f, x_i[k]=%f/n", k+p+step,x_r[k+p+step], x_i[k+p+step]);}/* 开跳!:) */k += 2*step;}}}/** 后记:* 究竟是颗粒在外层循环还是样本输入在外层, 好象也差不多, 复杂度完全一样.* 但以我资质愚钝花费了不少时间才弄明白这数十行代码.* 从中我发现一个于我非常有帮助的教训, 很久以前我写过一部分算法, 其中绝大多数都是递归.* 将数据量减少, 减少再减少, 用归纳的方式来找出数据量加大代码的规律* 比如FFT* 1. 先写死LayerI的代码; 然后再把LayerI的输出作为LayerII的输入, 又写死代码; ......* 大约3层就可以统计出规律来. 这和递归也是一样, 先写死一两层, 自然就出来了! * 2. 有的功能可以写伪代码, 不急于求出结果, 降低复杂性, 把逻辑结果定出来后再添加. * 比如旋转因子就可以写死, 就写1.0. 流程出来后再写旋转因子.* 寥寥数语, 我可真是流了不少汗! Happy!*/void dft( void ){int i, n, k, tx1, tx2;float tw1,tw2;float xx_r[N],xx_i[N];/** clear any data in Real and Imaginary result arrays prior to DFT */for(k=0; k<=N-1; k++)xx_r[k] = xx_i[k] = x_i[k] = 0.0;// caculate the DFTfor(k=0; k<=(N-1); k++){for(n=0; n<=(N-1); n++){tx1 = (n*k);tx2 = tx1+(3*N)/4;tx1 = tx1%(N);tx2 = tx2%(N);if(tx1 >= (N/2))tw1 = -twiddle[tx1-(N/2)];elsetw1 = twiddle[tx1];if(tx2 >= (N/2))tw2 = -twiddle[tx2-(N/2)];elsetw2 = twiddle[tx2];xx_r[k] = xx_r[k]+x_r[n]*tw1;xx_i[k] = xx_i[k]+x_r[n]*tw2;}xx_i[k] = -xx_i[k];}// displayfor(i=0; i<N; i++)printf("%f/t%f/n", xx_r[i], xx_i[i]);}//---------------------------------------------------------------------------int main( void ){fft_init( );bitrev( );// bitrev2( );//fft1( );fft2( );display( );system( "pause" );// dft();return 1;}#include <iom128.h>#include <intrinsics.h>/*********************************************************************快速福利叶变换C函数函数简介：此函数是通用的快速傅里叶变换C语言函数，移植性强，以下部分不依赖硬件。
C语言、Matlab实现FFT几种编程实例..
/*******************************************************************
函数原型：struct compx EE(struct compx b1,struct compx b2)
函数功能：对两个复数进行乘法运算
输入参数：两个以联合体定义的复数a,b
输出参数：a和b的乘积，以联合体的形式输出
*******************************************************************/
struct compx EE(struct compx a,struct compx b)
while(1);
}
%////二、
%/////////////////////////////////////////////////////////////////////////////////////////////////////////////
%///////////////////////////////////////////////////////////////////////////////////////////////////////////
{
double real;
double img;
}complex;
void fft(); /*快速傅里叶变换*/
void ifft(); /*快速傅里叶逆变换*/
void initW();
void change();
void add(complex ,complex ,complex *); /*复数加法*/
快速傅里叶变换(FFT)源代码
疯狂代码 / ĵ:http://VC/Article20803.htm才发现原来这堆变换说白了只是些数字游戏也没用到啥 复变知识最后用C实现了下FFT也算告段落代码如下:
{ tmp = pData1[i]; pData1[i] = pData1[rev]; pData1[rev] = tmp; } }
// output result data for(i=0; i<n; i) cout<<pData1[i].re<<"\t"; cout<<endl; for(i=0; i<n; i) cout<<pData1[i].im<<"\t"; cout<<endl;
# <iostream> # <fstream> # <math.h>
using std;
const double PI = 3.14159265358979323846;
n; // 数据个数 = 2logn次方 logn;
/// 复数结构体 struct stCompNum { double re; double im; };
cnt <<= 1;
stCompNum* pTmp = NULL; pTmp = pData1; pData1 = pData2; pData2 = pTmp; }
// resort for(i=0; i<n; i) { rev = reverseBits(i, logn); stCompNum tmp; (rev > i)
0000000000000000
2009-2-12 3:42:50 疯狂代码 /
FFT-C快速傅里叶变换超级详细的原代码
快速傅立叶变换（FFT）的C++实现收藏标准的离散傅立叶DFT 变换形式如：y k=Σj=0n-1a jωn-kj = A (ωn-k).（ωn k为复数1 的第k 个n 次方根，且定义多项式A (x)=Σj=0n-1a j x j）而离散傅立叶逆变换IDFT (Inverse DFT)形式如：a j=(Σk=0n-1y kωn kj)/n .yk=Σj=0n-1 ajωn-kj = A (ωn-k).（ωnk 为复数1 的第k 个n 次方根，且定义多项式 A (x) = Σj=0n-1 ajxj ）而离散傅立叶逆变换IDFT (Inverse DFT)形式如：aj=(Σk=0n-1 ykωnkj)/n .以下不同颜色内容为引用并加以修正：快速傅立叶变换（Fast Fourier Transform，FFT）是离散傅立叶变换（Discrete Fourier transform，DFT）的快速算法，它是根据离散傅立叶变换的奇、偶、虚、实等特性，对离散傅立叶变换的算法进行改进获得的。
它对傅立叶变换的理论并没有新的发现，但是对于在计算机系统或者说数字系统中应用离散傅立叶变换，可以说是进了一大步。
设Xn 为N 项的复数序列，由DFT 变换，任一Xi 的计算都需要N 次复数乘法和N -1 次复数加法，而一次复数乘法等于四次实数乘法和两次实数加法，一次复数加法等于两次实数加法，即使把一次复数乘法和一次复数加法定义成一次“运算”（四次实数乘法和四次实数加法），那么求出N 项复数序列的Xi ，即N 点DFT 变换大约就需要N2 次运算。
当N =1024 点甚至更多的时候，需要N2 = 1048576 次运算，在FFT 中，利用ωn 的周期性和对称性，把一个N 项序列（设N 为偶数），分为两个N / 2 项的子序列，每个N / 2点DFT 变换需要(N / 2)2 次运算，再用N 次运算把两个N / 2点的DFT 变换组合成一个N 点的DFT 变换。
c语言编写fft
#include<stdio.h>#include<math.h>#define PI 3.14159//逆序子程序，xr[]和xi[]分别输入序列的实部和虚部，n为序列长度（FFT点数）void ss(float xr[],float xi[],int n){int i=0,j,s;float a,bj;for(j=1;j<n;j++){for(s=n/2;s<=i;s=s/2){i=i-s;}i=i+s;if(i>j){a=xr[i];bj=xi[i];xi[i]=xi[j];xr[i]=xr[j];xr[j]=a;xi[j]=bj;}}}//------------------------------------------------------------------------//基2时分FFT子程序，xr[]和xi[]分别输入序列的实部和虚部,调用后返回为FFT的实部//和虚部，n为序列长度（FFT点数），m=1为FFT，m=-1为IFFT，PI=3.14159void fft(float xr[],float xi[],int m,int n){int loop1,loop2,loop3,le,le1,ip,i=1;float wr,wi,tr,ti,ur1,ur,ui,ls;ls=log(n)/log(2); //n=2^lsss(xr,xi,n); //逆序for(loop1=1;loop1<=ls;loop1++) //分级计算FFT，共ls级蝶形运算{ur=1.0;ui=0.0; //每一级的W（2^loop 0）le=pow(2,loop1); //le=2^loop每个蝶群中相同次序运算蝶的左前翅地址间隔le1=le/2; //le1:每个蝶群中所含的运算蝶个数wr=cos(PI*m/le1);wi=-sin(PI*m/le1);for(loop2=0;loop2<le1;loop2++) //计算每个蝶群中的第loop2个运算蝶{for(loop3=loop2;loop3<n;loop3=loop3+le) //由第一个蝶群开始，计算蝶群中//的第loop2个运算蝶{ip=loop3+le1;tr=ur*xr[ip]-ui*xi[ip];ti=ur*xi[ip]+ui*xr[ip];xr[ip]=xr[loop3]-tr;xi[ip]=xi[loop3]-ti;xr[loop3]=tr+xr[loop3];xi[loop3]=ti+xi[loop3];}ur1=wr*ur-ui*wi;ui=wr*ui+ur*wi;ur=ur1;}}if(m==-1){for(i=0;i<n;i++){xr[i]=xr[i]/n;xi[i]=xi[i]/n;}}}//-------------------------------------------------------------------。
fft c语言代码
以下是使用C语言实现快速傅里叶变换（FFT）的示例代码：c复制代码#include<stdio.h>#include<math.h>#define PI 3.14void fft(double* x, double* out, int n) {if (n == 1) {out[0] = x[0];return;}double w_n = 2 * PI / n;double w_m = exp(-1j * w_n / 2);double w = 1.0;for (int k = 0; k < n / 2; k++) {double t = w * x[2 * k];out[k] = t + x[2 * k + 1];out[k + n / 2] = t - x[2 * k + 1];w *= w_m;}fft(x, out, n / 2);fft(&x[n / 2], &out[n / 2], n / 2);}int main() {int n = 8; // 输入序列长度为8double x[] = {1, 1, 1, 1, 0, 0, 0, 0}; // 输入序列double out[n]; // 输出序列fft(x, out, n); // 进行FFT变换for (int i = 0; i < n; i++) {printf("%f ", out[i]); // 输出FFT变换结果}printf("\n");return0;}以上代码中，fft函数实现了快速傅里叶变换，main函数用于测试。
在fft函数中，首先判断序列长度是否为1，如果是则直接返回序列本身；否则，将序列分为奇数位和偶数位两个子序列，分别进行FFT变换，然后将两个子序列的FFT结果进行合并。
在循环中，使用w表示旋转因子，w_n表示旋转角度，w_m表示旋转因子的乘积因子。
C语言写的FFT代码
tw1 = cos(2*PI*(p+step)/pow(2, cur_layer));
tw2 = -sin(2*PI*(p+step)/pow(2, cur_layer));
x_r[k+p+step] = tmp_real + ( tw1*temp - tw2*x_i[k+p+step] );
/**
*初始化输出虚部
*/
static voidfft_ir(i=0; i<N; i++) x_i[i] = 0.0;
}
/**
*反转算法.将时域信号重新排序.
*这个算法有改进的空间
*/
static voidbitrev(void)
{
intp=1, q, i;
intbit_rev[ N ];
printf("k=%d, x_r[k]=%f, x_i[k]=%f\n", k+p, x_r[k+p], x_i[k+p]);
printf("k=%d, x_r[k]=%f, x_i[k]=%f\n", k+p+step, x_r[k+p+step], x_i[k+p+step]);
}
/*开跳!:) */
x_r[k+p] = tmp_real + ( tw1*x_r[k+p+step] - tw2*x_i[k+p+step] );
x_i[k+p] = tmp_imag + ( tw2*x_r[k+p+step] + tw1*x_i[k+p+step] );
FFT算法的c++实现
FFT算法源代码（C++）#include <iostream>#include <iomanip>#include "math.h"using namespace std;double pi = 3.1415926535897932;class complex{public://无参构造函数complex(){re=0;im=0;}//有参构造函数complex(double real,double imag){re=real;im=imag;}//加法complex operator + (complex& c){return complex( re + c.re , im + c.im );}//减法complex operator - (complex& c){return complex( re - c.re , im - c.im );}//乘法complex operator * (complex& c){return complex( (re * c.re)-(im * c.im) , (re * c.im)+(im * c.re) );}//除法complex operator / (complex& c){return complex( ( re*c.re + im*c.im )/( c.re*c.re + c.im*c.im ),((im * c.re)-(re * c.im))/((c.re*c.re)+(c.im*c.im)) );}//除2void half(){re=re/2;im=im/2;}//输出到ostreamvoid insert(ostream& out){if(re>=0) cout<<" ";out<<re<< ((im>=0)?'+':'-') <<'j'<< ((im>=0)?im:(0-im));}//显示void show(){cout<<re<< ((im>=0)?'+':'-') <<'j'<< ((im>=0)?im:(0-im)) ;}//设值void setValue(double real,double imag){if(real>0.00000001||real< -0.00000001) re=real;else re=0;if(imag>0.00000001||imag< -0.00000001) im=imag;else im=0;}friend complex c_pow(complex c,int y);private:double re,im;};//流输出ostream& operator << (ostream& out, complex c){c.insert(out);return out;}//乘方complex c_pow(complex c,int y){complex returnValue = c ;for(int i = 1 ; i < y ; i++ ){returnValue = returnValue * c ;}return returnValue;}//获得Wnvoid getW(complex w[],int len){for(int i=0 ; i<len ; i++ ){w[i].setValue( cos(0 - pi*2*(i)/len) , sin(0 - pi*2*(i)/len) );}}//获得mWnvoid getMW(complex w[],int len){for(int i=0 ; i<len ; i++ ){w[i].setValue( cos(pi*2*(i)/len) , sin(pi*2*(i)/len) );}}//倒序重排void resort(complex x[],int len){//初始化int half=len/2;int p1=0,p2;int k;complex temp;//重排for( p2=0 ; p2+1<len ; p2++ ){if(p2<p1 && p1+1<len){temp=x[p1];x[p1]=x[p2];x[p2]=temp;}k=half;while(k<p1+1 && p1+1<len){p1=p1-k;k=k/2;}p1=p1+k;}}//基2时间FFTvoid FFT2t(complex x[],complex y[],int len){//初始化int m=0,n,k,j,q;int tlen=len/2;while(tlen!=0){m++;tlen=tlen/2;}resort(x,len); //倒序重排for(j=0;j<len;j++){y[j]=x[j];}complex *xt = new complex[len];complex *w = new complex[len];getW(w,len);//开始迭代运算q=len;for(int i=1;i<=m;i++){q=q/2;n=len/q/2;for(j=0;j<len;j++){xt[j]=y[j];}if(i!=1){ for(j=0;j<q;j++){for(k=0;k<n;k++){xt[n+j*2*n+k]=xt[n+j*2*n+k]*w[k*q];}}for(j=0;j<n;j++){for(k=0;k<q;k++){y[j+k*2*n]=xt[j+k*2*n]+xt[j+n+k*2*n];}}for(;j<2*n;j++){for(k=0;k<q;k++){y[j+k*2*n]=xt[j-n+k*2*n]-xt[j+k*2*n];}}}delete[]xt,w;}//基2频率FFT和IFFTvoid FFT(complex x[],complex y[],int len,bool IFFT = false) {//初始化int m=0,n,k,j,q;int tlen=len/2;while(tlen!=0){m++;tlen=tlen/2;}for(j=0;j<len;j++){y[j]=x[j];}complex *xt = new complex[len];complex *w = new complex[len];if(IFFT==true) //若是IFFT则获取Wn序列{getMW(w,len);}else //若是FFT则获取mWn序列{getW(w,len);}//开始迭代运算for(int i=1;i<=m;i++){n=n/2;q=len/n/2;for(j=0;j<len;j++){xt[j]=y[j];}for(j=0;j<n;j++){for(k=0;k<q;k++){y[j+k*2*n]=xt[j+k*2*n]+xt[j+n+k*2*n];}}for(;j<2*n;j++){for(k=0;k<q;k++){y[j+k*2*n]=xt[j-n+k*2*n]-xt[j+k*2*n];}}if(i!=m){ for(j=0;j<q;j++){for(k=0;k<n;k++){y[n+j*2*n+k]=y[n+j*2*n+k]*w[k*q];}}}if(IFFT==true)//若是IFFT则除以2{for(j=0;j<len;j++){y[j].half();}}}resort(y,len); //倒序重排delete[]xt,w;}void main(){//初始化complex q(0.9,0.3);complex x[32],y[32],z[32],w[32],v[32],w1;w1.setValue(1,0);getW(w,32); //获取Wn序列x[0].setValue(1,0);for(int i=1;i<32;i++)//定原序列的值{x[i]=x[i-1]*q;}//输出原序列cout<<"原序列："<<endl;for(i=0;i<32;i+=2){cout<<setiosflags(ios::fixed)<<setprecision(6)<<x[i]<<"\t"<<x[i+1]<<endl;}for(i=0;i<32;i++) //进行公式计算{complex a=w1-(x[31]*q),b=w1-q*w[i];v[i]=a/b;}//输出公式计算值cout<<"公式计算值："<<endl;for(i=0;i<32;i+=2){cout<<v[i]<<"\t"<<v[i+1]<<endl;}FFT2t(x,y,32); //进行基2时间变换//输出基2时间变换后序列cout<<"基2时间变换后序列："<<endl;for(i=0;i<32;i+=2){cout<<y[i]<<"\t"<<y[i+1]<<endl;}FFT(y,z,32,1); //进行反变换//输出反变换后的序列cout<<"反变换后的序列："<<endl;for(i=0;i<32;i+=2)cout<<z[i]<<"\t"<<z[i+1]<<endl;}}运行结果。