源：用于ARM上的FFT与IFFT源代码（C语言，不依赖特定平台）
代码在2011年全国电子大赛结束后（2011年9月3日）发布，多个版本，注释详细。
/******************************************************************************* ** 程序名称：快速傅里叶变换(FFT) ** 程序描述：本程序实现快速傅里叶变换 ** 程序作者：宋元瑞 ** 最后修改：2011年4月5日 *******************************************************************************/ #include <stdio.h> #include <math.h> #define PI 3.141592653589 //圆周率，12位小数 #define N 8 //傅里叶变换的点数 #define M 3 //蝶形运算的级数，N = 2^M typedef double ElemType; //原始数据序列的数据类型,可以在这里设置 typedef struct //定义复数结构体 { ElemType real,imag; }complex; complex data[N]; //定义存储单元，原始数据与负数结果均使用之 ElemType result[N]; //存储FFT后复数结果的模 //变址 void ChangeSeat(complex *DataInput) { int nextValue,nextM,i,k,j=0; complex temp; nextValue=N/2; //变址运算，即把自然顺序变成倒位序，采用雷德算法 nextM=N-1; for (i=0;i<nextM;i++) { if (i<j) //如果i<j,即进行变址 { temp=DataInput[j]; DataInput[j]=DataInput[i]; DataInput[i]=temp; } k=nextValue; //求j的下一个倒位序 while (k<=j) //如果k<=j,表示j的最高位为1 { j=j-k; //把最高位变成0 k=k/2; //k/2，比较次高位，依次类推，逐个比较，直到某个位为0 } j=j+k; //把0改为1 } } /* //变址 void ChangeSeat(complex *DataInput) { complex Temp[N]; int i,n,New_seat; for(i=0; i<N; i++) { Temp[i].real = DataInput[i].real; Temp[i].imag = DataInput[i].imag; } for(i=0; i<N; i++) { New_seat = 0; for(n=0;n<M;n++) { New_seat = New_seat | (((i>>n) & 0x01) << (M-n-1)); } DataInput[New_seat].real = Temp[i].real; DataInput[New_seat].imag = Temp[i].imag; } } */ //复数乘法 complex XX_complex(complex a, complex b) { complex temp; temp.real = a.real * b.real-a.imag*b.imag; temp.imag = b.imag*a.real + a.imag*b.real; return temp; } //FFT void FFT(void) { int L=0,B=0,J=0,K=0; int step=0; ElemType P=0,T=0; complex W,Temp_XX; //ElemType TempResult[N]; ChangeSeat(data); for(L=1; L<=M; L++) { B = 1<<(L-1);//B=2^(L-1) for(J=0; J<=B-1; J++) { P = (1<<(M-L))*J;//P=2^(M-L) *J step = 1<<L;//2^L for(K=J; K<=N-1; K=K+step) { W.real = cos(2*PI*P/N); W.imag = -sin(2*PI*P/N); Temp_XX = XX_complex(data[K+B],W); data[K+B].real = data[K].real - Temp_XX.real; data[K+B].imag = data[K].imag - Temp_XX.imag; data[K].real = data[K].real + Temp_XX.real; data[K].imag = data[K].imag + Temp_XX.imag; } } } } void IFFT(void) { int L=0,B=0,J=0,K=0; int step=0; ElemType P=0,T=0; complex W,Temp_XX; //ElemType TempResult[N]; ChangeSeat(data); for(L=1; L<=M; L++) { B = 1<<(L-1);//B=2^(L-1) for(J=0; J<=B-1; J++) { P = (1<<(M-L))*J;//P=2^(M-L) *J step = 1<<L;//2^L for(K=J; K<=N-1; K=K+step) { W.real = cos(2*PI*P/N); W.imag = sin(2*PI*P/N);//逆运算，这里跟FFT符号相反 Temp_XX = XX_complex(data[K+B],W); data[K+B].real = data[K].real - Temp_XX.real; data[K+B].imag = data[K].imag - Temp_XX.imag; data[K].real = data[K].real + Temp_XX.real; data[K].imag = data[K].imag + Temp_XX.imag; } } } } int main(int argc, char *argv[]) { int i = 0; for(i=0; i<N; i++)//制造输入序列 { data[i].real = sin(2*PI*i/N); printf("%lf ",data[i]); } printf(" "); FFT();//进行FFT计算 printf(" "); for(i=0; i<N; i++) {printf("%lf ",sqrt(data[i].real*data[i].real+data[i].imag*data[i].imag));} IFFT();//进行FFT计算 printf(" "); for(i=0; i<N; i++) {printf("%lf ",data[i].real/N);} printf(" "); /*for(i=0; i<N; i++) {printf("%lf ",data[i].imag/N);} printf(" ");*/ /*for(i=0; i<N; i++) {printf("%lf ",sqrt(data[i].real*data[i].real+data[i].imag*data[i].imag)/N);}*/ return 0; }
/******************************************************************************* ** 程序名称：快速傅里叶变换(FFT) ** 程序描述：本程序实现快速傅里叶变换 ** 性能提升：修正了IFFT的bug,使用宏定义改变N大小 ** 程序版本：V6.5 ** 程序作者：宋元瑞 ** 最后修改：2011年5月16日 *******************************************************************************/ #include <stdio.h> #include <math.h> #define PI 3.14159265358979323846264338327950288419716939937510 //圆周率，50位小数 #define PI2 6.28318530717958647692528676655900576839433879875021 #define N 1024 //傅里叶变换的点数 #define M 10 //蝶形运算的级数，N = 2^M #define Npart2 512 //创建正弦函数表时取PI的1/2 #define Npart4 256 //创建正弦函数表时取PI的1/4 #define FFT_RESULT(x) (sqrt(data[x].real*data[x].real+data[x].imag*data[x].imag)) #define IFFT_RESULT(x) (data[x].real/N) typedef float ElemType; //原始数据序列的数据类型,可以在这里设置 typedef struct //定义复数结构体 { ElemType real,imag; }complex; complex data[N]; //定义存储单元，原始数据与负数结果均使用之 ElemType SIN_TABLE[Npart4+1]; //产生模拟原始数据输入 void InputData(void) { int i; for(i=0; i<N; i++)//制造输入序列 { data[i].real = sin(2*PI*i/N); //正弦波 data[i].imag = 0; } } //创建正弦函数表 void CREATE_SIN_TABLE(void) { int i=0; for(i=0; i<=Npart4; i++) { SIN_TABLE[i] = sin(PI*i/Npart2);//SIN_TABLE[i] = sin(PI2*i/N); } } ElemType Sin_find(ElemType x) { int i = (int)(N*x); i = i>>1; if(i>Npart4)//注意：i已经转化为0~N之间的整数了！ {//不会超过N/2 i = Npart2 - i;//i = i - 2*(i-Npart4); } return SIN_TABLE[i]; } ElemType Cos_find(ElemType x) { int i = (int)(N*x); i = i>>1; if(i<Npart4)//注意：i已经转化为0~N之间的整数了！ {//不会超过N/2 //i = Npart4 - i; return SIN_TABLE[Npart4 - i]; } else//i>Npart4 && i<N/2 { //i = i - Npart4; return -SIN_TABLE[i - Npart4]; } } //变址 void ChangeSeat(complex *DataInput) { int nextValue,nextM,i,k,j=0; complex temp; nextValue=N/2; //变址运算，即把自然顺序变成倒位序，采用雷德算法 nextM=N-1; for (i=0;i<nextM;i++) { if (i<j) //如果i<j,即进行变址 { temp=DataInput[j]; DataInput[j]=DataInput[i]; DataInput[i]=temp; } k=nextValue; //求j的下一个倒位序 while (k<=j) //如果k<=j,表示j的最高位为1 { j=j-k; //把最高位变成0 k=k/2; //k/2，比较次高位，依次类推，逐个比较，直到某个位为0 } j=j+k; //把0改为1 } } //复数乘法 /*complex XX_complex(complex a, complex b) { complex temp; temp.real = a.real * b.real-a.imag*b.imag; temp.imag = b.imag*a.real + a.imag*b.real; return temp; }*/ //FFT运算函数 void FFT(void) { int L=0,B=0,J=0,K=0; int step=0, KB=0; //ElemType P=0; ElemType angle; complex W,Temp_XX; ChangeSeat(data);//变址 //CREATE_SIN_TABLE(); for(L=1; L<=M; L++) { step = 1<<L;//2^L B = step>>1;//B=2^(L-1) for(J=0; J<B; J++) { //P = (1<<(M-L))*J;//P=2^(M-L) *J angle = (double)J/B; //这里还可以优化 W.imag = -Sin_find(angle); //用C++该函数课声明为inline W.real = Cos_find(angle); //用C++该函数课声明为inline //W.real = cos(angle*PI); //W.imag = -sin(angle*PI); for(K=J; K<N; K=K+step) { KB = K + B; //Temp_XX = XX_complex(data[KB],W); //用下面两行直接计算复数乘法，省去函数调用开销 Temp_XX.real = data[KB].real * W.real-data[KB].imag*W.imag; Temp_XX.imag = W.imag*data[KB].real + data[KB].imag*W.real; data[KB].real = data[K].real - Temp_XX.real; data[KB].imag = data[K].imag - Temp_XX.imag; data[K].real = data[K].real + Temp_XX.real; data[K].imag = data[K].imag + Temp_XX.imag; } } } } //IFFT运算函数 void IFFT(void) { int L=0,B=0,J=0,K=0; int step=0, KB=0; //ElemType P=0; ElemType angle; complex W,Temp_XX; ChangeSeat(data);//变址 //CREATE_SIN_TABLE(); for(L=1; L<=M; L++) { step = 1<<L;//2^L B = step>>1;//B=2^(L-1) for(J=0; J<B; J++) { //P = (1<<(M-L))*J;//P=2^(M-L) *J angle = (double)J/B; //这里还可以优化 W.imag = Sin_find(angle); //用C++该函数课声明为inline W.real = Cos_find(angle); //用C++该函数课声明为inline //W.real = cos(angle*PI); //W.imag = -sin(angle*PI); for(K=J; K<N; K=K+step) { KB = K + B; //Temp_XX = XX_complex(data[KB],W); //用下面两行直接计算复数乘法，省去函数调用开销 Temp_XX.real = data[KB].real * W.real-data[KB].imag*W.imag; Temp_XX.imag = W.imag*data[KB].real + data[KB].imag*W.real; data[KB].real = data[K].real - Temp_XX.real; data[KB].imag = data[K].imag - Temp_XX.imag; data[K].real = data[K].real + Temp_XX.real; data[K].imag = data[K].imag + Temp_XX.imag; } } } } //主函数 int main(int argc, char *argv[]) { int i = 0; CREATE_SIN_TABLE(); //创建正弦函数表 ,这句只需在程序开始时执行一次 InputData(); //输入原始数据 ，此处用公式模拟；实际应用时为AD采样数据 FFT(); //进行 FFT计算 for(i=0; i<N; i++) {printf("%f ",FFT_RESULT(i));}/**/ IFFT();//进行 IFFT计算 for(i=0; i<N; i++) {printf("%f ",IFFT_RESULT(i));}/**/ return 0; }
/******************************************************************************* ** 程序名称：快速傅里叶变换(FFT) ** 程序描述：本程序实现快速傅里叶变换及其逆变换 ** 性能提升：修正了FFT的bug,变量重新命名,并将 N_FFT改为动态值 ** 程序版本：V6.6 ** 程序作者：宋元瑞 ** 最后修改：2011年5月16日 *******************************************************************************/ #include <stdio.h> #include <stdlib.h> #include <math.h> //#define OUTPRINT printf //#define DEL /##/ #define RESULT(x) sqrt(data_of_N_FFT[x].real*data_of_N_FFT[x].real+data_of_N_FFT[x].imag*data_of_N_FFT[x].imag) #define PI 3.14159265358979323846264338327950288419716939937510 //圆周率，50位小数 #define PI2 6.28318530717958647692528676655900576839433879875021 int N_FFT=0; //傅里叶变换的点数 int M_of_N_FFT=0; //蝶形运算的级数，N = 2^M int Npart2_of_N_FFT=0; //创建正弦函数表时取PI的1/2 int Npart4_of_N_FFT=0; //创建正弦函数表时取PI的1/4 typedef float ElemType; //原始数据序列的数据类型,可以在这里设置 typedef struct //定义复数结构体 { ElemType real,imag; }complex_of_N_FFT,*ptr_complex_of_N_FFT; ptr_complex_of_N_FFT data_of_N_FFT=NULL;//开辟存储单元，原始数据与负数结果均使用之 ElemType* SIN_TABLE_of_N_FFT=NULL; //产生模拟原始数据输入 void InputData(void) { int i; for (i=0; i<N_FFT; i++)//制造输入序列 { data_of_N_FFT[i].real = sin(2*PI*i/N_FFT); //正弦波 data_of_N_FFT[i].imag = 0; printf("%f ",data_of_N_FFT[i].real); } } //创建正弦函数表 void CREATE_SIN_TABLE(void) { int i=0; for (i=0; i<=Npart4_of_N_FFT; i++) { SIN_TABLE_of_N_FFT[i] = sin(PI*i/Npart2_of_N_FFT);//SIN_TABLE[i] = sin(PI2*i/N); } } ElemType Sin_find(ElemType x) { int i = (int)(N_FFT*x); i = i>>1; if (i>Npart4_of_N_FFT)//注意：i已经转化为0~N之间的整数了！ { //不会超过N/2 i = Npart2_of_N_FFT - i;//i = i - 2*(i-Npart4); } return SIN_TABLE_of_N_FFT[i]; } ElemType Cos_find(ElemType x) { int i = (int)(N_FFT*x); i = i>>1; if (i<Npart4_of_N_FFT)//注意：i已经转化为0~N之间的整数了！ { //不会超过N/2 //i = Npart4 - i; return SIN_TABLE_of_N_FFT[Npart4_of_N_FFT - i]; } else//i>Npart4 && i<N/2 { //i = i - Npart4; return -SIN_TABLE_of_N_FFT[i - Npart4_of_N_FFT]; } } //变址 void ChangeSeat(complex_of_N_FFT *DataInput) { int nextValue,nextM,i,k,j=0; complex_of_N_FFT temp; nextValue=N_FFT/2; //变址运算，即把自然顺序变成倒位序，采用雷德算法 nextM=N_FFT-1; for (i=0;i<nextM;i++) { if (i<j) //如果i<j,即进行变址 { temp=DataInput[j]; DataInput[j]=DataInput[i]; DataInput[i]=temp; } k=nextValue; //求j的下一个倒位序 while (k<=j) //如果k<=j,表示j的最高位为1 { j=j-k; //把最高位变成0 k=k/2; //k/2，比较次高位，依次类推，逐个比较，直到某个位为0 } j=j+k; //把0改为1 } } //复数乘法 /*complex XX_complex(complex a, complex b) { complex temp; temp.real = a.real * b.real-a.imag*b.imag; temp.imag = b.imag*a.real + a.imag*b.real; return temp; }*/ //FFT运算函数 void FFT(void) { int L=0,B=0,J=0,K=0; int step=0, KB=0; //ElemType P=0; ElemType angle; complex_of_N_FFT W,Temp_XX; ChangeSeat(data_of_N_FFT);//变址 //CREATE_SIN_TABLE(); for (L=1; L<=M_of_N_FFT; L++) { step = 1<<L;//2^L B = step>>1;//B=2^(L-1) for (J=0; J<B; J++) { //P = (1<<(M-L))*J;//P=2^(M-L) *J angle = (double)J/B; //这里还可以优化 W.imag = -Sin_find(angle); //用C++该函数课声明为inline W.real = Cos_find(angle); //用C++该函数课声明为inline //W.real = cos(angle*PI); //W.imag = -sin(angle*PI); for (K=J; K<N_FFT; K=K+step) { KB = K + B; //Temp_XX = XX_complex(data[KB],W); //用下面两行直接计算复数乘法，省去函数调用开销 Temp_XX.real = data_of_N_FFT[KB].real * W.real-data_of_N_FFT[KB].imag*W.imag; Temp_XX.imag = W.imag*data_of_N_FFT[KB].real + data_of_N_FFT[KB].imag*W.real; data_of_N_FFT[KB].real = data_of_N_FFT[K].real - Temp_XX.real; data_of_N_FFT[KB].imag = data_of_N_FFT[K].imag - Temp_XX.imag; data_of_N_FFT[K].real = data_of_N_FFT[K].real + Temp_XX.real; data_of_N_FFT[K].imag = data_of_N_FFT[K].imag + Temp_XX.imag; } } } } //IFFT运算函数 void IFFT(void) { int L=0,B=0,J=0,K=0; int step=0, KB=0; //ElemType P=0; ElemType angle; complex_of_N_FFT W,Temp_XX; ChangeSeat(data_of_N_FFT);//变址 //CREATE_SIN_TABLE(); for (L=1; L<=M_of_N_FFT; L++) { step = 1<<L;//2^L B = step>>1;//B=2^(L-1) for (J=0; J<B; J++) { //P = (1<<(M-L))*J;//P=2^(M-L) *J angle = (double)J/B; //这里还可以优化 W.imag = Sin_find(angle); //用C++该函数课声明为inline W.real = Cos_find(angle); //用C++该函数课声明为inline //W.real = cos(angle*PI); //W.imag = -sin(angle*PI); for (K=J; K<N_FFT; K=K+step) { KB = K + B; //Temp_XX = XX_complex(data[KB],W); //用下面两行直接计算复数乘法，省去函数调用开销 Temp_XX.real = data_of_N_FFT[KB].real * W.real-data_of_N_FFT[KB].imag*W.imag; Temp_XX.imag = W.imag*data_of_N_FFT[KB].real + data_of_N_FFT[KB].imag*W.real; data_of_N_FFT[KB].real = data_of_N_FFT[K].real - Temp_XX.real; data_of_N_FFT[KB].imag = data_of_N_FFT[K].imag - Temp_XX.imag; data_of_N_FFT[K].real = data_of_N_FFT[K].real + Temp_XX.real; data_of_N_FFT[K].imag = data_of_N_FFT[K].imag + Temp_XX.imag; } } } } //初始化FFT程序 //N_FFT是 FFT的点数，必须是2的次方 void Init_FFT(int N_of_FFT) { int i=0; int temp_N_FFT=1; N_FFT = N_of_FFT; //傅里叶变换的点数 ，必须是 2的次方 M_of_N_FFT = 0; //蝶形运算的级数，N = 2^M for (i=0; temp_N_FFT<N_FFT; i++) { temp_N_FFT = 2*temp_N_FFT; M_of_N_FFT++; } printf(" %d ",M_of_N_FFT); Npart2_of_N_FFT = N_FFT/2; //创建正弦函数表时取PI的1/2 Npart4_of_N_FFT = N_FFT/4; //创建正弦函数表时取PI的1/4 //ptr_complex_of_N_FFT data_of_N_FFT=NULL;//开辟存储单元，原始数据与负数结果均使用之 data_of_N_FFT = (ptr_complex_of_N_FFT)malloc(N_FFT * sizeof(complex_of_N_FFT)); //ptr_complex_of_N_FFT SIN_TABLE_of_N_FFT=NULL; SIN_TABLE_of_N_FFT = (ElemType *)malloc((Npart4_of_N_FFT+1) * sizeof(ElemType)); CREATE_SIN_TABLE(); //创建正弦函数表 } //结束 FFT运算，释放相关内存 void Close_FFT(void) { if (data_of_N_FFT != NULL) { free(data_of_N_FFT); //释放内存 data_of_N_FFT = NULL; } if (SIN_TABLE_of_N_FFT != NULL) { free(SIN_TABLE_of_N_FFT); //释放内存 SIN_TABLE_of_N_FFT = NULL; } } //主函数 int main(int argc, char *argv[]) { int i = 0; Init_FFT(1024); //初始化各项参数，此函数只需执行一次 InputData(); //输入原始数据 FFT(); //进行 FFT计算 printf(" "); for (i=0; i<N_FFT; i++) { printf("%f ",RESULT(i)); } IFFT();//进行 IFFT计算 printf(" "); for (i=0; i<N_FFT; i++) { printf("%f ",data_of_N_FFT[i].real/N_FFT); } Close_FFT(); //结束 FFT运算，释放相关内存 return 0; }
/******************************************************************************* ** 模块名称：快速傅里叶变换(FFT) ** 模块描述：本程序实现快速傅里叶变换 ** 性能提升：已达到网上同类程序最高速度 ** 模块版本：V6.0 ** 模块作者：宋元瑞 ** 最后修改：2011年5月6日 *******************************************************************************/ /******************************************************************************* ** 程序说明： FFT运算输入参数 source_Data(i) 是一个N大小的数组（注意是小括号） FFT运算输出结果 result_Data(i) 是一个N大小的数组（注意是小括号） ** 调用举例： int main(int argc, char *argv[]) { int i = 0; ///以下为FFT运算的调用方式 FFT_prepare(); //为FFT做好准备，此函数只需程序开始时执行一次即可，切勿写在循环中 while(1) { for(i=0;i<N_FFT;i++) //输入原始数据 {source_Data(i) = sin(2*PI*i/N_FFT);}//注意inputData后面是小括号 FFT(); //进行FFT计算 //输出结果：XXX = result_Data(i); } return 0; } *******************************************************************************/ #ifndef SYRFFT_H //#pragma once //#include <stdio.h> #include <math.h> #define FFT_prepare() CREATE_SIN_TABLE(); //创建正弦函数表 #define source_Data(i) data_FFT[i].imag //接收输入数据的数组，大小为N #define result_Data(i) sqrt(data_FFT[i].real*data_FFT[i].real+data_FFT[i].imag*data_FFT[i].imag)//FFT结果 #define PI 3.14159265358979323846264338327950288419716939937510 //圆周率，50位小数 #define PI2 6.28318530717958647692528676655900576839433879875021 #define N_FFT 1024 //傅里叶变换的点数 #define M_of_N_FFT 10 //蝶形运算的级数，N = 2^M #define Npart2_of_N_FFT 512 //创建正弦函数表时取PI的1/2 #define Npart4_of_N_FFT 256 //创建正弦函数表时取PI的1/4 typedef float ElemType; //原始数据序列的数据类型,可以在这里设置 typedef struct //定义复数结构体 { ElemType real,imag; }complex_of_N_FFT; complex_of_N_FFT data_FFT[N_FFT]; //存放要进行FFT运输的数据，运算结果也存放这里 ElemType SIN_TABLE[Npart4_of_N_FFT+1]; //产生模拟原始数据输入 /* void InputData(void) { int i; for(i=0; i<N_FFT; i++) //制造输入序列 { source_Data(i) = sin(2*PI*i/N_FFT); //模拟输入正弦波 //data_FFT[i].imag = sin(2*PI*i/N); //正弦波 //printf("%f ",data_FFT[i].imag); } }*/ //创建正弦函数表 void CREATE_SIN_TABLE(void) { int i=0; for(i=0; i<=Npart4_of_N_FFT; i++) { SIN_TABLE[i] = sin(PI*i/Npart2_of_N_FFT);//SIN_TABLE[i] = sin(PI2*i/N); } } ElemType Sin_find(ElemType x) { int i = (int)(N_FFT*x); i = i>>1; if(i>Npart4_of_N_FFT)//注意：i已经转化为0~N之间的整数了！ {//不会超过N/2 i = Npart2_of_N_FFT - i;//i = i - 2*(i-Npart4); } return SIN_TABLE[i]; } ElemType Cos_find(ElemType x) { int i = (int)(N_FFT*x); i = i>>1; if(i<Npart4_of_N_FFT)//注意：i已经转化为0~N之间的整数了！ {//不会超过N/2 //i = Npart4 - i; return SIN_TABLE[Npart4_of_N_FFT - i]; } else//i>Npart4 && i<N/2 { //i = i - Npart4; return -SIN_TABLE[i - Npart4_of_N_FFT]; } } //变址 void ChangeSeat(complex_of_N_FFT *DataInput) { //变址前data_FFT[i].imag存储了输入数据，变址后data_FFT[i].real存储了输入数据 int i,n,New_seat; for(i=0; i<N_FFT; i++) { New_seat = 0; for(n=0;n<M_of_N_FFT;n++) { New_seat = New_seat | (((i>>n) & 0x01) << (M_of_N_FFT-n-1)); } DataInput[New_seat].real = DataInput[i].imag; } for(i=0; i<N_FFT; i++) { DataInput[i].imag = 0; } } //复数乘法 /* complex_of_N_FFT XX_complex(complex_of_N_FFT a, complex_of_N_FFT b) { complex_of_N_FFT temp; temp.real = a.real * b.real-a.imag*b.imag; temp.imag = b.imag*a.real + a.imag*b.real; return temp; }*/ //FFT运算函数 void FFT(void) { int L=0,B=0,J=0,K=0; int step=0, KB=0; //ElemType P=0; ElemType angle; complex_of_N_FFT W,Temp_XX; ChangeSeat(data_FFT); //变址 for(L=1; L<=M_of_N_FFT; L++) { step = 1<<L;//2^L B = step>>1;//B=2^(L-1) for(J=0; J<B; J++) { //P = (1<<(M-L))*J;//P=2^(M-L) *J angle = (double)J/B; //这里还可以优化 W.imag = -Sin_find(angle); //用C++该函数课声明为inline W.real = Cos_find(angle); //用C++该函数课声明为inline //W.real = cos(angle*PI); //W.imag = -sin(angle*PI); for(K=J; K<N_FFT; K=K+step) { KB = K + B; //Temp_XX = XX_complex(data_FFT[KB],W); //用下面两行直接计算复数乘法，省去函数调用开销 Temp_XX.real = data_FFT[KB].real * W.real-data_FFT[KB].imag*W.imag; Temp_XX.imag = W.imag*data_FFT[KB].real + data_FFT[KB].imag*W.real; data_FFT[KB].real = data_FFT[K].real - Temp_XX.real; data_FFT[KB].imag = data_FFT[K].imag - Temp_XX.imag; data_FFT[K].real = data_FFT[K].real + Temp_XX.real; data_FFT[K].imag = data_FFT[K].imag + Temp_XX.imag; } } } } #define SYRFFT_H #endif
/******************************************************************************* ** 程序名称：快速傅里叶变换(FFT) ** 程序描述：本程序实现快速傅里叶变换 ** 程序版本：V6.0 ** 程序作者：宋元瑞 ** 最后修改：2011年5月6日 *******************************************************************************/ #include <stdio.h> #include "syrFFT_H.h" //包含FFT运算头文件 //主函数 int main(int argc, char *argv[]) { int i = 0; //以下3句为FFT运算的调用函数 FFT_prepare(); //为FFT做好准备，此函数只需程序开始时执行一次即可，切勿写在循环中 //while(1) { for(i=0;i<N_FFT;i++)//模拟输入 {source_Data(i) = sin(2*PI*i/N_FFT);}//注意inputData后面是小括号 FFT(); for(i=0; i<N_FFT; i++)//输出结果 {printf("%lf ",result_Data(i));} } return 0; }
#ifndef FFT_H #include <stdlib.h> #include <math.h> #define RESULT(x) sqrt(data_of_N_FFT[x].real*data_of_N_FFT[x].real+data_of_N_FFT[x].imag*data_of_N_FFT[x].imag) #define PI 3.14159265358979323846264338327950288419716939937510 //圆周率，50位小数 #define PI2 6.28318530717958647692528676655900576839433879875021 int N_FFT=0; //傅里叶变换的点数 int M_of_N_FFT=0; //蝶形运算的级数，N = 2^M int Npart2_of_N_FFT=0; //创建正弦函数表时取PI的1/2 int Npart4_of_N_FFT=0; //创建正弦函数表时取PI的1/4 typedef float ElemType; //原始数据序列的数据类型,可以在这里设置 typedef struct //定义复数结构体 { ElemType real,imag; }complex_of_N_FFT,*ptr_complex_of_N_FFT; ptr_complex_of_N_FFT data_of_N_FFT=NULL;//开辟存储单元，原始数据与负数结果均使用之 ElemType* SIN_TABLE_of_N_FFT=NULL; //创建正弦函数表 void CREATE_SIN_TABLE(void) { int i=0; for(i=0; i<=Npart4_of_N_FFT; i++) { SIN_TABLE_of_N_FFT[i] = sin(PI*i/Npart2_of_N_FFT);//SIN_TABLE[i] = sin(PI2*i/N); } } ElemType Sin_find(ElemType x) { int i = (int)(N_FFT*x); i = i>>1; if(i>Npart4_of_N_FFT)//注意：i已经转化为0~N之间的整数了！ {//不会超过N/2 i = Npart2_of_N_FFT - i;//i = i - 2*(i-Npart4); } return SIN_TABLE_of_N_FFT[i]; } ElemType Cos_find(ElemType x) { int i = (int)(N_FFT*x); i = i>>1; if(i<Npart4_of_N_FFT)//注意：i已经转化为0~N之间的整数了！ {//不会超过N/2 //i = Npart4 - i; return SIN_TABLE_of_N_FFT[Npart4_of_N_FFT - i]; } else//i>Npart4 && i<N/2 { //i = i - Npart4; return -SIN_TABLE_of_N_FFT[i - Npart4_of_N_FFT]; } } //变址 void ChangeSeat(complex_of_N_FFT *DataInput) { int nextValue,nextM,i,k,j=0; complex_of_N_FFT temp; nextValue=N_FFT/2; //变址运算，即把自然顺序变成倒位序，采用雷德算法 nextM=N_FFT-1; for (i=0;i<nextM;i++) { if (i<j) //如果i<j,即进行变址 { temp=DataInput[j]; DataInput[j]=DataInput[i]; DataInput[i]=temp; } k=nextValue; //求j的下一个倒位序 while (k<=j) //如果k<=j,表示j的最高位为1 { j=j-k; //把最高位变成0 k=k/2; //k/2，比较次高位，依次类推，逐个比较，直到某个位为0 } j=j+k; //把0改为1 } } //复数乘法 /*complex XX_complex(complex a, complex b) { complex temp; temp.real = a.real * b.real-a.imag*b.imag; temp.imag = b.imag*a.real + a.imag*b.real; return temp; }*/ //FFT运算函数 void FFT(void) { int L=0,B=0,J=0,K=0; int step=0, KB=0; //ElemType P=0; ElemType angle; complex_of_N_FFT W,Temp_XX; ChangeSeat(data_of_N_FFT);//变址 //CREATE_SIN_TABLE(); for(L=1; L<=M_of_N_FFT; L++) { step = 1<<L;//2^L B = step>>1;//B=2^(L-1) for(J=0; J<B; J++) { //P = (1<<(M-L))*J;//P=2^(M-L) *J angle = (double)J/B; //这里还可以优化 W.imag = -Sin_find(angle); //用C++该函数课声明为inline W.real = Cos_find(angle); //用C++该函数课声明为inline //W.real = cos(angle*PI); //W.imag = -sin(angle*PI); for(K=J; K<N_FFT; K=K+step) { KB = K + B; //Temp_XX = XX_complex(data[KB],W); //用下面两行直接计算复数乘法，省去函数调用开销 Temp_XX.real = data_of_N_FFT[KB].real * W.real-data_of_N_FFT[KB].imag*W.imag; Temp_XX.imag = W.imag*data_of_N_FFT[KB].real + data_of_N_FFT[KB].imag*W.real; data_of_N_FFT[KB].real = data_of_N_FFT[K].real - Temp_XX.real; data_of_N_FFT[KB].imag = data_of_N_FFT[K].imag - Temp_XX.imag; data_of_N_FFT[K].real = data_of_N_FFT[K].real + Temp_XX.real; data_of_N_FFT[K].imag = data_of_N_FFT[K].imag + Temp_XX.imag; } } } } //IFFT运算函数 void IFFT(void) { int L=0,B=0,J=0,K=0; int step=0, KB=0; //ElemType P=0; ElemType angle; complex_of_N_FFT W,Temp_XX; ChangeSeat(data_of_N_FFT);//变址 //CREATE_SIN_TABLE(); for(L=1; L<=M_of_N_FFT; L++) { step = 1<<L;//2^L B = step>>1;//B=2^(L-1) for(J=0; J<B; J++) { //P = (1<<(M-L))*J;//P=2^(M-L) *J angle = (double)J/B; //这里还可以优化 W.imag = Sin_find(angle); //用C++该函数课声明为inline W.real = Cos_find(angle); //用C++该函数课声明为inline //W.real = cos(angle*PI); //W.imag = -sin(angle*PI); for(K=J; K<N_FFT; K=K+step) { KB = K + B; //Temp_XX = XX_complex(data[KB],W); //用下面两行直接计算复数乘法，省去函数调用开销 Temp_XX.real = data_of_N_FFT[KB].real * W.real-data_of_N_FFT[KB].imag*W.imag; Temp_XX.imag = W.imag*data_of_N_FFT[KB].real + data_of_N_FFT[KB].imag*W.real; data_of_N_FFT[KB].real = data_of_N_FFT[K].real - Temp_XX.real; data_of_N_FFT[KB].imag = data_of_N_FFT[K].imag - Temp_XX.imag; data_of_N_FFT[K].real = data_of_N_FFT[K].real + Temp_XX.real; data_of_N_FFT[K].imag = data_of_N_FFT[K].imag + Temp_XX.imag; } } } } //初始化FFT程序 //N_FFT是 FFT的点数，必须是2的次方 void Init_FFT(int N_of_FFT) { int i=0; int temp_N_FFT=1; N_FFT = N_of_FFT; //傅里叶变换的点数 ，必须是 2的次方 M_of_N_FFT = 0; //蝶形运算的级数，N = 2^M for(i=0; temp_N_FFT<N_FFT; i++) { temp_N_FFT = 2*temp_N_FFT; M_of_N_FFT++; } //printf(" %d ",M_of_N_FFT); Npart2_of_N_FFT = N_FFT/2; //创建正弦函数表时取PI的1/2 Npart4_of_N_FFT = N_FFT/4; //创建正弦函数表时取PI的1/4 //ptr_complex_of_N_FFT data_of_N_FFT=NULL;//开辟存储单元，原始数据与负数结果均使用之 data_of_N_FFT = (ptr_complex_of_N_FFT)malloc(N_FFT * sizeof(complex_of_N_FFT)); //ptr_complex_of_N_FFT SIN_TABLE_of_N_FFT=NULL; SIN_TABLE_of_N_FFT = (ElemType *)malloc((Npart4_of_N_FFT+1) * sizeof(ElemType)); CREATE_SIN_TABLE(); //创建正弦函数表 } //结束 FFT运算，释放相关内存 void Close_FFT(void) { if(data_of_N_FFT != NULL) { free(data_of_N_FFT); //释放内存 data_of_N_FFT = NULL; } if(SIN_TABLE_of_N_FFT != NULL) { free(SIN_TABLE_of_N_FFT); //释放内存 SIN_TABLE_of_N_FFT = NULL; } } #define FFT_H #endif
/******************************************************************************* ** 程序名称：快速傅里叶变换(FFT) ** 程序描述：本程序实现快速傅里叶变换及其逆变换 ** 性能提升：修正了FFT的bug ** 程序版本：V6.6 ** 程序作者：宋元瑞 ** 最后修改：2011年5月16日 *******************************************************************************/ #include "jxust_fft6_6.h" //产生模拟原始数据输入 ,在实际应用时替换为AD采样数据 void InputData(void) { int i; for(i=0; i<N_FFT; i++)//制造输入序列 { data_of_N_FFT[i].real = sin(2*PI*i/N_FFT); //输入采样数据 data_of_N_FFT[i].imag = 0; } } //主函数 ，示例如何调用FFT运算 int main(int argc, char *argv[]) { int i = 0; Init_FFT(1024); //①初始化各项参数，此函数只需执行一次 ; 参数为FFT的点数，必须为2的次方 InputData(); //②输入原始数据 ,此处在实际应用时替换为AD采样数据 FFT(); //③进行 FFT计算 //for(i=0; i<N_FFT; i++)//④输出FFT频谱结果 sqrt(data_of_N_FFT[i].real*data_of_N_FFT[i].real+data_of_N_FFT[i].imag*data_of_N_FFT[i].imag) //{printf("%f ",RESULT(i));} IFFT();//进行 IFFT计算 //for(i=0; i<N_FFT; i++)//(可选步骤)⑤输出 IFFT结果 ，滤波时会用到；暂时不用 //{printf("%f ",data_of_N_FFT[i].real/N_FFT);} Close_FFT(); //结束 FFT运算，释放相关内存 ;此函数在彻底结束FFT运算后调用， return 0; }
#ifndef SYRFFT_6_55_H #include <math.h> #define FFT_RESULT(x) (sqrt(data_of_N_FFT[x].real*data_of_N_FFT[x].real+data_of_N_FFT[x].imag*data_of_N_FFT[x].imag)) #define IFFT_RESULT(x) (data_of_N_FFT[x].real/N_FFT) #define PI 3.14159265358979323846264338327950288419716939937510 //圆周率，50位小数 #define PI2 6.28318530717958647692528676655900576839433879875021 #define N_FFT 1024 //傅里叶变换的点数 #define M_of_N_FFT 10 //蝶形运算的级数，N = 2^M #define Npart2_of_N_FFT 512 //创建正弦函数表时取PI的1/2 #define Npart4_of_N_FFT 256 //创建正弦函数表时取PI的1/4 typedef float ElemType; //原始数据序列的数据类型,可以在这里设置 typedef struct //定义复数结构体 { ElemType real,imag; }complex_of_N_FFT,*ptr_complex_of_N_FFT; complex_of_N_FFT data_of_N_FFT[N_FFT]; //定义存储单元，原始数据与负数结果均使用之 ElemType SIN_TABLE_of_N_FFT[Npart4_of_N_FFT+1]; //创建正弦函数表 void CREATE_SIN_TABLE(void) { int i=0; for(i=0; i<=Npart4_of_N_FFT; i++) { SIN_TABLE_of_N_FFT[i] = sin(PI*i/Npart2_of_N_FFT);//SIN_TABLE[i] = sin(PI2*i/N); } } ElemType Sin_find(ElemType x) { int i = (int)(N_FFT*x);//注意：i已经转化为0~N之间的整数了！ i = i>>1;// i = i / 2; if(i>Npart4_of_N_FFT) {//根据FFT相关公式，sin()参数不会超过PI， 即i不会超过N/2 i = Npart2_of_N_FFT - i;//i = i - 2*(i-Npart4); } return SIN_TABLE_of_N_FFT[i]; } ElemType Cos_find(ElemType x) { int i = (int)(N_FFT*x);//注意：i已经转化为0~N之间的整数了！ i = i>>1; if(i < Npart4_of_N_FFT ) { //不会超过N/2 //i = Npart4 - i; return SIN_TABLE_of_N_FFT[Npart4_of_N_FFT - i]; } else //i > Npart4 && i < N/2 { //i = i - Npart4; return -SIN_TABLE_of_N_FFT[i - Npart4_of_N_FFT]; } } //变址 void ChangeSeat(complex_of_N_FFT *DataInput) { int nextValue,nextM,i,k,j=0; complex_of_N_FFT temp; nextValue=N_FFT/2; //变址运算，即把自然顺序变成倒位序，采用雷德算法 nextM=N_FFT-1; for (i=0;i<nextM;i++) { if (i<j) //如果i<j,即进行变址 { temp=DataInput[j]; DataInput[j]=DataInput[i]; DataInput[i]=temp; } k=nextValue; //求j的下一个倒位序 while (k<=j) //如果k<=j,表示j的最高位为1 { j=j-k; //把最高位变成0 k=k/2; //k/2，比较次高位，依次类推，逐个比较，直到某个位为0 } j=j+k; //把0改为1 } } //复数乘法 /*complex XX_complex(complex a, complex b) { complex temp; temp.real = a.real * b.real-a.imag*b.imag; temp.imag = b.imag*a.real + a.imag*b.real; return temp; }*/ //FFT运算函数 void FFT(void) { int L=0,B=0,J=0,K=0; int step=0, KB=0; //ElemType P=0; ElemType angle; complex_of_N_FFT W,Temp_XX; ChangeSeat(data_of_N_FFT);//变址 //CREATE_SIN_TABLE(); for(L=1; L<=M_of_N_FFT; L++) { step = 1<<L;//2^L B = step>>1;//B=2^(L-1) for(J=0; J<B; J++) { //P = (1<<(M-L))*J;//P=2^(M-L) *J angle = (double)J/B; //这里还可以优化 W.imag = -Sin_find(angle); //用C++该函数课声明为inline W.real = Cos_find(angle); //用C++该函数课声明为inline //W.real = cos(angle*PI); //W.imag = -sin(angle*PI); for(K=J; K<N_FFT; K=K+step) { KB = K + B; //Temp_XX = XX_complex(data[KB],W); //用下面两行直接计算复数乘法，省去函数调用开销 Temp_XX.real = data_of_N_FFT[KB].real * W.real-data_of_N_FFT[KB].imag*W.imag; Temp_XX.imag = W.imag*data_of_N_FFT[KB].real + data_of_N_FFT[KB].imag*W.real; data_of_N_FFT[KB].real = data_of_N_FFT[K].real - Temp_XX.real; data_of_N_FFT[KB].imag = data_of_N_FFT[K].imag - Temp_XX.imag; data_of_N_FFT[K].real = data_of_N_FFT[K].real + Temp_XX.real; data_of_N_FFT[K].imag = data_of_N_FFT[K].imag + Temp_XX.imag; } } } } //IFFT运算函数 void IFFT(void) { int L=0,B=0,J=0,K=0; int step=0, KB=0; //ElemType P=0; ElemType angle; complex_of_N_FFT W,Temp_XX; ChangeSeat(data_of_N_FFT);//变址 //CREATE_SIN_TABLE(); for(L=1; L<=M_of_N_FFT; L++) { step = 1<<L;//2^L B = step>>1;//B=2^(L-1) for(J=0; J<B; J++) { //P = (1<<(M-L))*J;//P=2^(M-L) *J angle = (double)J/B; //这里还可以优化 W.imag = Sin_find(angle); //用C++该函数课声明为inline W.real = Cos_find(angle); //用C++该函数课声明为inline //W.real = cos(angle*PI); //W.imag = -sin(angle*PI); for(K=J; K<N_FFT; K=K+step) { KB = K + B; //Temp_XX = XX_complex(data[KB],W); //用下面两行直接计算复数乘法，省去函数调用开销 Temp_XX.real = data_of_N_FFT[KB].real * W.real-data_of_N_FFT[KB].imag*W.imag; Temp_XX.imag = W.imag*data_of_N_FFT[KB].real + data_of_N_FFT[KB].imag*W.real; data_of_N_FFT[KB].real = data_of_N_FFT[K].real - Temp_XX.real; data_of_N_FFT[KB].imag = data_of_N_FFT[K].imag - Temp_XX.imag; data_of_N_FFT[K].real = data_of_N_FFT[K].real + Temp_XX.real; data_of_N_FFT[K].imag = data_of_N_FFT[K].imag + Temp_XX.imag; } } } } //初始化FFT程序 void Init_FFT() { CREATE_SIN_TABLE(); //创建正弦函数表 } //结束 FFT运算，释放相关内存 void Close_FFT(void) { } #define SYRFFT_6_55_H #endif
/******************************************************************************* ** 程序名称：快速傅里叶变换(FFT) ** 程序描述：本程序实现快速傅里叶变换 ** 性能提升：修正了IFFT的bug,变量重新命名,使用宏定义改变N大小 ** 程序版本：V6.55 ** 程序作者：宋元瑞 ** 最后修改：2011年5月22日 *******************************************************************************/ #include "syrFFT_6_55.h" //产生模拟原始数据输入 ,在实际应用时替换为AD采样数据 void InputData(void) { int i; for(i=0; i<N_FFT; i++)//制造输入序列 { data_of_N_FFT[i].real = sin(2*PI*i/N_FFT); //输入采样数据 data_of_N_FFT[i].imag = 0; } } //主函数 ，示例如何调用FFT运算 int main(int argc, char *argv[]) { int i = 0; Init_FFT(); //①初始化各项参数，此函数只需执行一次 ; 修改FFT的点数去头文件的宏定义处修改 InputData(); //②输入原始数据 ,此处在实际应用时替换为AD采样数据 FFT(); //③进行 FFT计算 //for(i=0; i<N_FFT; i++)//④输出FFT频谱结果 sqrt(data_of_N_FFT[i].real*data_of_N_FFT[i].real+data_of_N_FFT[i].imag*data_of_N_FFT[i].imag) //{printf("%f ",FFT_RESULT(i));} IFFT();//进行 IFFT计算 //for(i=0; i<N_FFT; i++)//(可选步骤)⑤输出 IFFT结果 ，滤波时会用到；暂时不用 //{printf("%f ",IFFT_RESULT(i));} Close_FFT(); //结束 FFT运算，此版本此句无用，可不写 return 0; }
回复更精彩，请到原出处关注。