三种C语言实现的FFT例程
- 格式:doc
- 大小:70.50 KB
- 文档页数:13
傅里叶变换c程序傅里叶变换是一种用于分析和处理信号的重要数学工具,广泛应用于信号处理、图像处理、通信系统等领域。
在数学上,傅里叶变换可以将一个时域上的连续或离散信号转换为频域上的连续或离散信号,从而提供了一种将信号从时域转换为频域表示的方法。
在计算机科学中,傅里叶变换也有很大的应用。
通过使用傅里叶变换算法,可以对信号进行频谱分析,提取信号的频率成分,并在频域上进行滤波、去噪等处理。
傅里叶变换的计算可以使用多种方法,包括连续傅里叶变换(FFT)和离散傅里叶变换(DFT)。
为了实现傅里叶变换的计算,可以使用C语言编写相应的程序。
下面是一个简单的傅里叶变换C程序的示例:```c#include <stdio.h>#include <stdlib.h>#include <math.h>#define PI 3.14159265typedef struct {double real;double imag;} Complex;void fft(Complex* x, int N) {if (N <= 1) {}Complex* even = (Complex*) malloc(N/2 * sizeof(Complex)); Complex* odd = (Complex*) malloc(N/2 * sizeof(Complex)); 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++) {Complex t;double angle = -2 * PI * k / N;t.real = cos(angle) * odd[k].real + sin(angle) * odd[k].imag; t.imag = cos(angle) * odd[k].imag - sin(angle) * odd[k].real; x[k].real = even[k].real + t.real;x[k].imag = even[k].imag + t.imag;x[k + N/2].real = even[k].real - t.real;x[k + N/2].imag = even[k].imag - t.imag;}free(even);}int main(void) {int N = 8;Complex* x = (Complex*) malloc(N * sizeof(Complex));// 初始化输入信号for (int i = 0; i < N; i++) {x[i].real = i;x[i].imag = 0;}fft(x, N);printf("傅里叶变换结果:\n");for (int i = 0; i < N; i++) {printf("X[%d] = %.2f + %.2fi\n", i, x[i].real, x[i].imag);}free(x);return 0;}```该程序演示了如何通过C代码实现傅里叶变换。
傅里叶在c语言中的使用傅里叶变换是一种在信号处理、图像处理等领域具有重要应用的算法。
它可以帮助我们将复杂的信号分解成一系列简单的正弦和余弦函数的叠加,从而更好地分析信号的频率特性。
在C语言中,傅里叶变换有着广泛的应用,下面我们将介绍傅里叶变换的基本概念、实现方法以及应用场景。
一、傅里叶变换的基本概念傅里叶变换是一种将时间域或空间域的信号转换为频域信号的数学方法。
其基本原理是将一个复杂的信号分解成无数个简单的正弦和余弦函数的叠加,这种叠加称为傅里叶级数。
在实际应用中,为了减少计算量,通常只取级数的前几项进行计算。
傅里叶变换的数学表达式如下:X(f) = ∫(-∞,∞) x(t) * e^(-jωt) dt其中,x(t)表示时域信号,X(f)表示频域信号,ω表示角频率,j表示虚数单位。
二、傅里叶变换在C语言中的实现方法1.离散傅里叶变换(DFT)离散傅里叶变换是傅里叶变换的一种离散形式,适用于离散信号的处理。
在C语言中,可以使用以下步骤实现DFT:(1)预处理:对输入信号进行窗函数处理,以减少频谱泄漏和旁瓣干扰。
(2)计算:按照DFT的计算公式,对输入信号的每个样本进行傅里叶变换。
(3)后处理:对变换结果进行幅度谱和相位谱的计算,并进行归一化处理。
2.快速傅里叶变换(FFT)快速傅里叶变换是一种高效计算离散傅里叶变换的方法,其时间复杂度为O(n log n)。
在C语言中,可以使用以下步骤实现FFT:(1)初始化:根据输入信号的长度,构建FFT递归函数。
(2)基2递归:按照FFT递归函数,对输入信号进行分组,并计算每组的傅里叶变换结果。
(3)合并:将每组的傅里叶变换结果合并,得到最终的傅里叶变换结果。
三、傅里叶变换的应用场景傅里叶变换在信号处理、图像处理等领域具有广泛的应用,如音频信号分析、图像滤波、模式识别等。
通过傅里叶变换,我们可以更好地分析信号的频率特性,从而为后续的处理和分析提供便利。
四、C语言实现傅里叶变换的实战案例以下是一个简单的C语言实现离散傅里叶变换的示例:```c#include <stdio.h>#include <stdlib.h>#include <math.h>void fft(float *in, float *out, int n) {// 基2递归实现FFT}int main() {int n = 8; // 采样点数float x[] = {1, 1, 1, 0, 0, 0, 0, 0}; // 输入信号float x_fft[n]; // 傅里叶变换结果fft(x, x_fft, n);// 输出傅里叶变换结果for (int i = 0; i < n; i++) {printf("x[%d] = %f", i, x_fft[i]);}return 0;}```五、总结与展望本文介绍了傅里叶变换在C语言中的基本概念、实现方法和应用场景。
在C语言中实现DFT(离散傅里叶变换)算法可以使用库函数fft(Fast Fourier Transform,快速傅里叶变换),也可以自己实现。
以下是一个简单的DFT算法的C语言实现:
void idft(double *x, int n)
{
double complex *y;
int i;
for (i = 0; i < n; i++) {
y[i] = creal(x[i]);
y[i] *= pow(-1.0, i * 2) * x[i];
}
for (i = 0; i < n; i++) {
x[i] = creal(y[i]);
y[i] = complex_conj(y[i]);
}
}
这个算法的基本思路是:首先将输入向量x中的所有元素求实部和虚部,然后将虚部乘以-1.0的i*2次方,再将实部和虚部相加,得到DFT输出向量y中的每个元素。
最后将y中的实部和虚部交换,得到原始向量x的反转DFT输出。
在这个算法中,我们使用了C语言中的complex类型来表示复数,使用creal()函数来获取复数的实部,使用complex_conj()函数来计算复数的共轭。
需要注意的是,这个算法只适用于n为偶数的情况,因为DFT输出向量中的元素需要和输入向量中的元素一一对应。
如果n为奇数,则需要对输入向量进行填充或删除,以使其满足偶数长度。
一、概述傅里叶变换是信号处理和数据压缩中常用的数学工具,它可以将时域信号转换为频域信号,从而便于分析和处理。
而快速傅里叶变换(FFT)则是一种高效的计算傅里叶变换的方法,可以大大提高计算效率,广泛应用于信号处理、图像处理、通信系统等领域。
二、傅里叶变换原理傅里叶变换的基本思想是将一个时域信号分解为不同频率的正弦和余弦函数的叠加,从而得到该信号的频谱图。
具体来说,对于一个连续信号x(t),它的傅里叶变换X(ω)定义为:X(ω) = ∫[0,∞]x(t)e^(-jωt)dt其中,ω为频率变量,X(ω)表示在频率ω处的信号能量。
而对于离散信号x[n],它的傅里叶变换X[k]则定义为:X[k] = ∑[n=0,N-1]x[n]e^(-j2πkn/N)其中,N为信号的采样点数,k为频率域的序号。
上述公式称为离散傅里叶变换(DFT),计算复杂度为O(N^2)。
而快速傅里叶变换则通过巧妙的算法设计,将计算复杂度降低到O(NlogN)。
三、快速傅里叶变换算法概述快速傅里叶变换的算法最早由Cooley和Tukey在1965年提出,它的基本思想是将一个长度为N的DFT分解为两个长度为N/2的DFT的组合,通过递归地分解和合并,最终实现对整个信号的快速计算。
下面我们来介绍一种常用的快速傅里叶变换算法:递归式分治法。
四、递归式分治法递归式分治法是一种高效的计算DFT的方法,它的基本思想是将长度为N的DFT分解为两个长度为N/2的DFT,并通过递归地调用自身,最终实现对整个信号的傅里叶变换。
具体来说,假设有一个长度为N的信号x[n],对其进行快速傅里叶变换的过程可以分为如下几个步骤:1. 将长度为N的信号x[n]分为长度为N/2的偶数序号和奇数序号的两个子信号x_even[n]和x_odd[n];2. 对子信号x_even[n]和x_odd[n]分别进行快速傅里叶变换,得到它们的频域表示X_even[k]和X_odd[k];3. 结合X_even[k]和X_odd[k],计算原信号的频域表示X[k]。
C语言实现FFTFFT(快速傅里叶变换)是一种常用的算法,用于计算离散傅里叶变换(DFT)的快速算法。
它在信号处理、图像处理、通信等方面有广泛的应用。
C语言提供了一组标准库函数来支持FFT算法的实现。
下面以C语言为例,展示如何实现FFT算法。
1.理解DFT首先,我们需要理解离散傅里叶变换(DFT)的概念。
DFT将时域离散信号转换为频域离散信号,它的计算公式如下:其中N是信号的长度,k表示频域的频率,n表示时域的时间。
离散信号经过DFT变换后,可以得到相应频率的幅度和相位信息。
2. Cooley-Tukey算法FFT算法采用了Cooley-Tukey算法的思想,它的主要思路是将DFT问题递归地分解为更小的DFT问题。
这样可以减少计算量,提高计算效率。
Cooley-Tukey算法具体如下:-如果信号的长度N是2的整数次幂,将原始信号分为偶数索引和奇数索引的两部分。
-对分离出的偶数索引部分和奇数索引部分分别进行DFT变换。
-最后将两个结果进行合并。
下面是一个简单的例子,展示了如何使用C语言实现FFT算法:```c#include <stdio.h>#include <math.h>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] = even[k] + t;x[k + N/2] = even[k] - t;}int maiint N = 8; // 信号长度//输入信号x[0]=1;x[1]=1;x[2]=1;x[3]=1;x[4]=0;x[5]=0;x[6]=0;x[7]=0;fft(x, N);//打印频域的幅度信息for (int k = 0; k < N; k++)printf("Amplitude at frequency %d: %f\n", k, cabs(x[k]));}return 0;```在上面的示例中,我们首先定义了一个fft函数,用于实现FFT算法。
傅里叶变换是一种在数学和工程中广泛使用的技术,它可以将一个信号分解为不同频率的成分。
在C语言中,我们可以使用库函数来实现傅里叶变换。
下面是一个简单的C语言程序,使用库函数进行傅里叶变换:c复制代码#include <stdio.h>#include <math.h>#include <complex.h>#include <fftw3.h>int main() {int N = 1024; // 采样点数fftw_complex *in, *out; // 输入和输出数组fftw_plan p; // 傅里叶变换计划// 创建输入和输出数组in = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N);out = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N);// 生成一个简单的正弦波信号作为输入信号for (int i = 0; i < N; i++) {in[i][0] = sin(2 * M_PI * i / N);in[i][1] = 0;}// 创建傅里叶变换计划p = fftw_plan_dft_1d(N, in, out, FFTW_FORWARD, FFTW_ESTIMATE);// 执行傅里叶变换fftw_execute(p);// 输出结果for (int i = 0; i < N; i++) {printf("%f + %fi\n", out[i][0], out[i][1]);}// 释放内存和关闭计划fftw_destroy_plan(p);fftw_free(in);fftw_free(out);return 0;}在这个程序中,我们首先定义了一个简单的正弦波信号作为输入信号。
/*此算法为FFT相位差算法的C语言实现,经验证可用(网上有很多此类程序无法运行);前半部分为FFT算法的子程序,也可用作其它用途*/#include<stdio.h>#include<stdlib.h>#include<math.h>typedef struct{double real,imag;}complex;complex cexp(complex a){complex z;z.real=exp(a.real)*cos(a.imag);z.imag=exp(a.real)*sin(a.imag);return(z);}void mcmpfft(complex x[],int n,int isign){complex t,z,ce;double pisign;int mr,m,l,j,i,nn;for(i=1;i<=16;i++)/*n must be power of2*/{nn=(int)pow(2,i);if(n==nn)break;}if(i>16){printf("N is not a power of2!\n");return;}z.real=0.0;pisign=4*isign*atan(1.);/*pisign的值为+180度或-180度*/mr=0;for(m=1;m<n;m++){l=n;while(mr+l>=n)l=l/2;mr=mr%l+l;if(mr<=m)continue;t.real=x[m].real;t.imag=x[m].imag;x[m].real=x[mr].real;x[m].imag=x[mr].imag;x[mr].real=t.real;x[mr].imag=t.imag;}l=1;while(1){if(l>=n){if(isign==-1)/*isign=-1For Forward Transform*/return;for(j=0;j<n;j++)/*Inverse Transform*/{x[j].real=x[j].real/n;x[j].imag=x[j].imag/n;}return;}for(m=0;m<l;m++)/*完成当前级所有蝶形运算*/{for(i=m;i<n;i=i+2*l)/*完成当前级相同W因子的所有蝶形运算*/{z.imag=m*pisign/l;ce=cexp(z);t.real=x[i+l].real*ce.real-x[i+l].imag*ce.imag;t.imag=x[i+l].real*ce.imag+x[i+l].imag*ce.real;x[i+l].real=x[i].real-t.real;/*原位运算*/x[i+l].imag=x[i].imag-t.imag;x[i].real=x[i].real+t.real;x[i].imag=x[i].imag+t.imag;}}l=2*l;/*确定下一级蝶形运算中W因子个数*/}}void main(){int i,N=1024,k;float pi=3.14159265,x[10000],A[10000],fi=3,f0=89,fs=450,max=0,o,oo,fi1,f00;complex s[2000];for(i=0;i<N;i++){x[i]=1*sin(2*pi*f0*i/fs+fi*pi/180);s[i].real=x[i];s[i].imag=0;}mcmpfft(s,N,-1);for(i=0;i<N;i++){A[i]=sqrt(s[i].real*s[i].real+s[i].imag*s[i].imag);if(A[i]>max){max=A[i];k=i;}}max=max*2/N;fi1=atan(s[k].imag/s[k].real)*180/pi+90;if(s[k].real>s[k].imag){o=(s[k+1].real-s[k-1].real)/(2*s[k].real-s[k-1].real-s[k+1].real);} else{o=(s[k+1].imag-s[k-1].imag)/(2*s[k].imag-s[k-1].imag-s[k+1].imag);}f00=(k-o)*fs/N;oo=(f00-f0)/f0*100;printf("fi=%f\nA=%f\nk0=%d\nf0=%f\n%f%",fi1,max,k,f00,oo);getch();}。
fft c语言程序以下是一个简单的FFT C语言程序的示例:```#include <stdio.h>#include <complex.h>#include <math.h>// 定义复数类型typedef double complex cplx;#define N 8void fft(cplx buf[], cplx out[], int n, int step){if (step < n){fft(out, buf, n, step * 2);fft(out + step, buf + step, n, step * 2);for (int i = 0; i < n; i += 2 * step){cplx t = cexp(-I * M_PI * i / n) * out[i + step]; buf[i / 2] = out[i] + t;buf[(i + n) / 2] = out[i] - t;}}}int main(){// 输入离散频域信号cplx buf[N] = {1, 2, 3, 4, 4, 3, 2, 1};// 创建输出数组cplx out[N];// 进行FFT变换fft(buf, out, N, 1);// 打印输出结果for (int i = 0; i < N; i++){printf("out[%d]: %.2f + %.2fi\n", i, creal(out[i]),cimag(out[i]));}return 0;}```该程序实现了一个基于蝶形算法的快速傅里叶变换算法。
它接受一个长度为N的复数数组作为输入,并将计算结果存储在另一个复数数组中。
在主函数中,我们定义了一个长度为N 的离散频域信号,并通过调用fft函数进行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语言函数,移植性强,以下部分不依赖硬件。
此程序包采用联合体的形式表示一个复数,输入为自然顺序的复数(输入实数是可令复数虚部为0),输出为经过FFT变换的自然顺序的复数.此程序包可在初始化时调用create_sin_tab()函数创建正弦函数表,以后的可采用查表法计算耗时较多的sin和cos运算,加快可计算速度使用说明:使用此函数只需更改宏定义FFT_N的值即可实现点数的改变,FFT_N的应该为2的N次方,不满足此条件时应在后面补0。
若使用查表法计算sin值和cos值,应在调用FFT函数前调用create_sin_tab()函数创建正弦表函数调用:FFT(s);时间:2010-2-20版本:Ver1.1参考文献:**********************************************************************/#include<math.h>#define FFT_N 128 //定义福利叶变换的点数#define PI 3.1415926535897932384626433832795028841971 //定义圆周率值struct compx {float real,imag;}; //定义一个复数结构struct compx s[FFT_N]; //FFT输入和输出:从S[0]开始存放,根据大小自己定义float SIN_TAB[FFT_N/2]; //定义正弦表的存放空间/*******************************************************************函数原型: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 create_sin_tab(float *sin_t)函数功能:创建一个正弦采样表,采样点数与福利叶变换点数相同输入参数:*sin_t存放正弦表的数组指针输出参数:无******************************************************************/ void create_sin_tab(float *sin_t){int i;for(i=0;i<FFT_N/2;i++)sin_t[i]=sin(2*PI*i/FFT_N);}/****************************************************************** 函数原型:void sin_tab(float pi)函数功能:采用查表的方法计算一个数的正弦值输入参数:pi 所要计算正弦值弧度值,范围0--2*PI,不满足时需要转换输出参数:输入值pi的正弦值******************************************************************/ float sin_tab(float pi){int n;float a;n=(int)(pi*FFT_N/2/PI);if(n>=0&&n<FFT_N/2)a=SIN_TAB[n];else if(n>=FFT_N/2&&n<FFT_N)a=-SIN_TAB[n-FFT_N/2];return a;}/****************************************************************** 函数原型:void cos_tab(float pi)函数功能:采用查表的方法计算一个数的余弦值输入参数:pi 所要计算余弦值弧度值,范围0--2*PI,不满足时需要转换输出参数:输入值pi的余弦值******************************************************************/ float cos_tab(float pi){float a,pi2;pi2=pi+PI/2;if(pi2>2*PI)pi2-=2*PI;a=sin_tab(pi2);return a;}/*****************************************************************函数原型: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)N le=2<<(m-1); //le蝶形结距离,即第m级蝶形的蝶形结相距le点lei=le/2; //同一蝶形结中参加运算的两点的距离u.real=1.0; //u为蝶形结运算系数,初始值为1u.imag=0.0;//w.real=cos(PI/lei); //不适用查表法计算sin值和cos值// w.imag=-sin(PI/lei);w.real=cos_tab(PI/lei); //w为系数商,即当前系数与前一个系数的商w.imag=-sin_tab(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;create_sin_tab(SIN_TAB);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语言函数,移植性强,以下部分不依赖硬件。