C语言编写FFT程序

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语言#include &quot;stdio.h&quot;#include &lt;math.h&gt; //調用源（頭）文件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&lt;=nm1;i++) /*倒序*/{if(i&lt;j) {t=xin[ j ];xin[j]=xin[ i ];xin[ i ] =t;} /*i&lt;j則換位*/k=nv2; /*k為倒序中相應位置的權值*/while(k&lt;j) {j=j-k;k=k/2;} /*k&lt;j時最高為變為0*/j=j+k; /* j為數組中的位數，是一個十進制數*/}{int le,lei,ip; //變量初始化,le為序列長度float pi;for(l=1;l&lt;=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&lt;=lei;j++) //外循環控制蝶行運算的級數{for(i=j;i&lt;=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(&quot;shu ru N de ge shu N=&quot;); /*提示輸入*/scanf(&quot;%d&quot;,&amp;N); /*輸入N*/ for(i=1;i&lt;=N;i++) /*輸入*/{printf(&quot;di %d ge shu real=&quot;,i);getchar();scanf(&quot;%f&quot;,&amp;s[ i ].real);getchar();printf(&quot;\n&quot;);printf(&quot;di %d ge shu imag=&quot;,i);scanf(&quot;%f&quot;,&amp;s[ i ].imag);printf(&quot;\n&quot;);}FFT(s,N); /*調用FFt*/for(i=1;i&lt;=N;i++) /*輸出*/{printf(&quot;%f&quot;,s[ i ].real);printf(&quot; + &quot;);printf(&quot;%f&quot;,s[ i ].imag);printf(&quot;j&quot;);printf(&quot;\n&quot;);}}。

#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>#defineFFT_N128//定义傅立叶变换的点数structcompx{floatreal,imag;};//定义一个复数结构structcompxs[FFT_N];//FFT输入和输出：从S[1]开始存放，根据大小自己定义/*******************************************************************函数原型：structcompxEE(structcompxb1,structcompxb2)函数功能：对两个复数进行乘法运算输入参数：两个以联合体定义的复数a,b输出参数：a和b的乘积，以联合体的形式输出*******************************************************************/ structcompxEE(structcompxa,structcompxb){structcompxc;c.real=a.real*b.real-a.imag*b.imag;c.imag=a.real*b.imag+a.imag*b.real;return(c);}/*****************************************************************函数原型：voidFFT(structcompx*xin,intN)函数功能：对输入的复数组进行快速傅里叶变换（FFT）输入参数：*xin复数结构体组的首地址指针，struct型*****************************************************************/ voidFFT(structcompx*xin){intf,m,nv2,nm1,i,k,l,j=0;structcompxu,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}{intle,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);//改变系数，进行下一个蝶形运算}}}}/************************************************************函数原型：voidmain()函数功能：测试FFT变换，演示函数使用方法输入参数：无输出参数：无************************************************************/ voidmain(){inti;for(i=0;i<FFT_N;i++)//给结构体赋值{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（快速傅里叶变换）和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类似，只需将其中几处计算公式略作变换即可。

#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;}。

#include "myapp.h"#include "ICETEK-VC5509-EDU.h"#include "scancode.h"#include <math.h>#define PI 3.1415926#define SAMPLENUMBER 128voidInitForFFT();void MakeWave();int INPUT[SAMPLENUMBER],DATA[SAMPLENUMBER];floatfWaveR[SAMPLENUMBER],fWaveI[SAMPLENUMBER],w[SAMPLENUMBER ];float sin_tab[SAMPLENUMBER],cos_tab[SAMPLENUMBER];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);for ( i=0;i<SAMPLENUMBER;i++ ){DATA[i]=w[i];}while ( 1 ); // break point}void FFT(float dataR[SAMPLENUMBER],float dataI[SAMPLENUMBER]) {int x0,x1,x2,x3,x4,x5,x6,xx;inti,j,k,b,p,L;float TR,TI,temp;/********** following code invert sequence ************/for ( i=0;i<SAMPLENUMBER;i++ ){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;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 */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++ ){INPUT[i]=sin(PI*2*i/SAMPLENUMBER*3)*1024;}}#include"myapp.h"voidCLK_init(){ioport unsigned int *clkmd;clkmd=(unsigned int *)0x1c00;*clkmd =0x2033; // 0x2033;//0x2413;// 200MHz=0x2513}voidSetDSPPLL(unsigned intuPLL){ioport unsigned int *clkmd;clkmd=(unsigned int *)0x1c00;*clkmd =uPLL;}voidTMCR_reset( void ){ioport unsigned int *TMCR_MGS3=(unsigned int *)0x07FE; ioport unsigned int *TMCR_MM =(unsigned int *)0x07FF; *TMCR_MGS3 =0x510;*TMCR_MM =0x000;}。

rfft c语言程序rfft C语言程序是一种用于实现快速傅里叶变换（FFT）的C语言程序。

FFT是一种数学算法，用于将一个时域信号转换为频域信号，广泛应用于信号处理、图像处理以及其他领域。

在C语言中，我们可以使用库函数来实现rfft程序。

一个常用的库是FFTW （Fastest Fourier Transform in the West），它提供了高效的FFT算法和相应的C语言接口。

以下是一个简单示例，演示如何使用rfft函数来进行信号的快速傅里叶变换：```c#include <stdio.h>#include <fftw3.h>#define N 1024int main() {// 创建输入和输出数组double in[N], out[N];// 创建FFTW执行计划fftw_plan plan;plan = fftw_plan_r2r_1d(N, in, out, FFTW_R2HC, FFTW_ESTIMATE);// 填充输入数组（示例中为简单的正弦波）for (int i = 0; i < N; i++) {in[i] = sin(2 * M_PI * i / N);}// 执行傅里叶变换fftw_execute(plan);// 输出结果for (int i = 0; i < N / 2 + 1; i++) {printf("频率：%d，幅度：%f\n", i, out[i]);}// 释放内存和计划fftw_destroy_plan(plan);fftw_cleanup();return 0;}```这段代码使用了FFTW库中的`fftw_plan_r2r_1d`函数来创建执行计划，指定了输入和输出数组的大小。

然后，通过填充输入数组（可以根据需求修改）并调用`fftw_execute`函数执行傅里叶变换。

最后，输出结果为频率和幅度。

DIT-基2FFT的浮点C语言程序：1、生成旋转因子，复数结构,旋转因子Wn=exp(—j*2*pi/N）//twiFactor——指向旋转因子矩阵的指针//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;for(iFactor=0;iFactor〈(wLen>〉1）；iFactor++）{twiFactor[iFactor］。

re=cos（stepFactor＊iFactor)；twiFactor[iFactor].im=sin(stepFactor＊iFactor); //W［n］=exp(j＊2＊pi*n/N），n=0，1，…，(N/2—1)｝｝2、在运行FFT之前，对输入序列进行倒序变换,代码如下：//bitRevData——指向位变换序列的指针//revLen—-FFT长度Void bit_rev(struct complexData ＊bitRevData，int revLen)｛struct complexData tempRev；int iRev，jRev，kRev，halfLen;halfLen=revLen>〉1；jRev=0;for(iRev=0；iRev〈(revLen—1）;iRev++）{If（iRev<jRev）{tempRev=bitRevData[jRev］；bitRevData[jRev]=bitRevData[iRev];bitRevData［iRev］= tempRev；}kRev=halfLen；while(kRev〈=jRev）｛jRev=jRev—kRev;kRev=kRev>>1;｝｝｝3、FFT计算.有3个循环体，分别为a内循环，b中间循环，c外循环,内循环实现蝶形结计算，循环a和b完成所有的蝶形结运算，而循环c则表示完成FFT算法所需要的级数。

#include <iom128.h>#include <intrinsics.h>#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>#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>**********************************************************************/#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/4+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 create_sin_tab(float *sin_t)函数功能：创建一个正弦采样表，采样点数与福利叶变换点数相同输入参数：*sin_t存放正弦表的数组指针输出参数：无******************************************************************/void create_sin_tab(float *sin_t){int i;for(i=0;i<=FFT_N/4;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/4)a=SIN_TAB[n];else if(n>FFT_N/4&&n<FFT_N/2){n-=FFT_N/4;a=SIN_TAB[FFT_N/4-n];}else if(n>=FFT_N/2&&n<3*FFT_N/4){n-=FFT_N/2;a=-SIN_TAB[n];}else if(n>=3*FFT_N/4&&n<3*FFT_N){n=FFT_N-n;a=-SIN_TAB[n];}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）Nle=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); }。