快速傅里叶变换程序设计

1 设计任务描述1.1 设计题目快速傅里叶变换程序设计设计要求1.2.1 设计目的1）理解FFT的算法以及利用DSP实现的方法。

2）能熟练的调试程序并能观察其结果。

3）熟悉TMS320C54x系列DSP芯片的软件设计方法。

基本要求1）研究FFT原理以及利用DSP实现的方法。

2）编写FFT程序。

3）调试程序，观察结果。

2 设计思路2.1 FFT 算法原理若给定由N 个信号样本{x （0），x （1），…，x （N -1）}组成的信号序列x （n ），DFT 可用式2-1给出：1()()N nkNn X k x n W-==∑ k =0,1，…，N -1（2-1）式2-1中，nk N W 称为旋转因子或蝶形因子，nkN W =2/j nk Neπ-。

从中可以看出：当信号样本为复数时，计算单个()X k 需经过N 次复数乘法和N -1次复数加法运算，相当于4N 次实数乘法和2（2N -1）次实数加法。

完成全部N 点DFT 共需2N 次复数乘法和N （N -1）复数加法运算。

可见，随着N 不断增加，整个DFT运算量是相当庞大的，而FFT 算法通过对计算过程的深入分析，利用旋转因子nkNW 具有的周期性与对称性，实现了降低运算复杂度的目的。

当序列长度N 为偶数时，信号序列x （n ）可被分解为奇、偶两个子序列，相应的N 点DFT 被分解为两个N /2点的DFT ：()()()kN X k G k W H k =+ k =0,1, …，N /2-1（2-2）(/2)()()k N X N k G k W H k +=- k =0,1, …，N /2-1（2-3）式（2-2）和（2-3）中，G(k)和()H k 分别表示x （n ）分解后得到的N /2点偶序列点奇序列的DFT 。

式（2-2）和式（2-3）表明，只要求出G(k)和()H k ，x （n ）前N /2点和后N /2点的DFT 就得到了，整个序列的DFT 也就得到了。

#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语言函数，移植性强，以下部分不依赖硬件。

快速傅里叶变换程序设计快速傅里叶变换（Fast Fourier Transform，简称FFT）是一种用于将一个函数在时域中的离散序列转换为频域中的连续谱的数学算法。

它在信号处理、图像处理、通信等领域有着广泛的应用。

本文将介绍快速傅里叶变换的原理，并给出一个简单的基于C语言的FFT程序设计。

1.快速傅里叶变换原理快速傅里叶变换通过分而治之的策略，将一个长度为N的序列分成两个长度为N/2的序列，并通过递归地计算这两个序列的FFT得到最终的结果。

FFT的核心思想是将离散序列的DFT（离散傅里叶变换）分解为两个子问题的DFT，并利用旋转因子的性质实现计算的重复利用，从而大大减少了计算量。

2.FFT的具体实现步骤（1）将输入序列分为偶数下标和奇数下标序列，分别求它们的DFT；（2）通过旋转因子将奇数下标序列转换为频域中的虚部；（3）将两个子问题的FFT结果合并得到最终的结果。

3.基于C语言的FFT程序设计下面是一个简单的基于C语言的FFT程序设计，用于计算一个长度为N的离散序列的FFT。

```c#include <stdio.h>#include <math.h>typedef structdouble real;double imag;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++)double theta = -2 * PI * k / N;t.real = cos(theta) * odd[k].real + sin(theta) * odd[k].imag; t.imag = cos(theta) * odd[k].imag - sin(theta) * 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;}int maiint N;printf("请输入序列的长度：");scanf("%d", &N);for (int i = 0; i < N; i++)printf("请输入第%d个元素的实部和虚部：", i+1);scanf("%lf %lf", &x[i].real, &x[i].imag);}fft(x, N);printf("FFT结果为：\n");for (int i = 0; i < N; i++)printf("第%d个元素：%.2f+%.2fi\n", i+1, x[i].real,x[i].imag);}return 0;```在实际使用中，我们可以根据需要对该程序进行优化，如增加输入提示、处理边界情况、改进算法等，从而更好地满足实际需求。

快速傅立叶变换（FFT ）的实现一、实验目的1.了解FFT 的原理及算法；2.了解DSP 中FFT 的设计及编程方法；3.熟悉FFT 的调试方法；二、实验原理FFT 是一种高效实现离散付立叶变换的算法，把信号从时域变换到频域，在频域分析处理信息。

对于长度为N 的有限长序列x (n )，它的离散傅里叶变换为：(2/)j N nk N W e π-=，称为旋转因子，或蝶形因子。

在x (n )为复数序列的情况下，计算X (k )：对某个k 值，需要N 次复数乘法、(N -1)次复数加法；对所有N 个k 值，需要2N 次复数乘法和N (N -1)次复数加法。

对于N 相当大时（如1024）来说，直接计算它的DFT 所作的计算量是很大的，FFT 的基本思想在于： 利用2()j nk N N W e π-=的周期性即：k N k N N W W +=对称性：/2k k N N N W W +=-将原有的N 点序列分成两个较短的序列，这些序列的DFT 可以很简单的组合起来得到原序列的DFT 。

按时间抽取的FFT ——DIT FFT 信号流图如图5.1所示：图5.1 时间抽取的FFT —DIT FFT 信号流图FFT 算法主要分为以下四步。

第一步 输入数据的组合和位倒序∑=-=10)()(N n nk N W n x k X把输入序列作位倒序是为了在整个运算最后的输出中得到的序列是自然顺序。

第二步 实现N 点复数FFT第一级蝶形运算；第二级蝶形运算；第三级至log2N 级蝶形运算；FFT 运算中的旋转因子N W 是一个复数，可表示：为了实现旋转因子N W 的运算，在存储空间分别建立正弦表和余弦表，每个表对应从0度到180度，采用循环寻址来对正弦表和余弦表进行寻址。

第三步 功率谱的计算X (k )是由实部()R X k 和虚部()I X k 组成的复数：()()()R I X k X k jX k =+；计算功率谱时只需将FFT 变换好的数据，按照实部()R X k 和虚部()I X k 求它们的平方和，然后对平方和进行开平方运算。