实验题目：实验1 FFT 算法的应用
姓 名： 学 号： 上课时间： FFT 算法的应用
1． 实验目的：
离散傅氏变换（DFT ）的目的是把信号由时域变换到频域，从而可以在频域分析处理信息，得到的结果再由逆DFT 变换到时域。FFT 是DFT 的一种快速算法。在数字信号处理系统中，FFT 作为一个非常重要的工具经常使用，甚至成为DSP 运算能力的一个考核因素。
本实验通过使用MATLAB 函数中的FFT 命令计算离散时间信号的频谱，以加深对离散信号的DFT 的理解及其FFT 算法的运用。
2． 实验要求：
对实验内容中给定的序列求给定点数N 的FFT 和IFFT ，利用MATLAB 编程完成计算，绘出相应图形。并与理论计算相比较，说明实验结果的原因。
3． 实验原理：
一．数字滤波器设计：
（一）基—2按时间抽取FFT 算法
对于有限长离散数字信号{x[n]}，0 ≤ n ≤ N-1，其离散谱{x[k]}可以由离
()1
,...,1,0][)2(
1
-==--=∑N k e
n x k X nk N
j N n π
散付氏变换（DFT ）求得。DFT 的定义为
可以方便的把它改写为如下形式： 不难看出，W N 是周期性的，且周期为N ，即
W N 的周期性是DFT 的关键性质之一。为了强调起见，常用表达式W N 取代W 以便明确其周期是N 。
由DFT 的定义可以看出，在x[n]为复数序列的情况下，完全直接运算N 点DFT 需要（N-1）2次复数乘法和N （N-1）次加法。因此，对于一些相当大的N 值（如1024）来说，直接计算它的DFT 所作的计算量是很大的。FFT 的基本思想在于，将原有的N 点序列序列分成两个较短的序列，这些序列的DFT 可以很简单的组合起来得到原序列的DFT 。例如，若N 为偶数，将原有的N 点序列分成两个（N/2）点序列，那么计算N 点DFT 将只需要约[(N/2)2 ·2]=N 2/2次复数乘法。即比直接计算少作一半乘法。因子（N/2）2表示直接计算（N/2）点DFT 所需要的乘法次数，而乘数2代表必须完成两个DFT 。上述处理方法可以反复使用，即（N/2）点的DFT 计算也可以化成两个（N/4）点的DFT （假定N/2为偶数），从而又少作一半的乘法。这样一级一级的划分下去一直到最后就划分成两点的FFT 运算的情况。比如，一个N = 8点的FFT 运算按照这种方法来计算FFT 可以用下面的流程图来表示：
x(0)
x(1)
x(2)
x(3)x(4)
x(5)
x(6)
x(7)
X(7)
X(6)
X(5)
X(4)X(3)X(2)X(1)X(0)
()1
,...,1,0][1
0-==∑-=N k W n x k X nk N
N n ...
2,1,0,))((±±==++l m W W nk
N
lN k mN n N
关于蝶形结运算的具体原理及其推导可以参照讲义，在此就不再赘述。按频率抽取的FFT的原理也可查阅相关资料，这里就不再推导了。
二．使用到的MATLAB命令：
函数fft(x)可以计算R点序列的R点DFT值；而fft(x,N)则计算R点序列的N点DFT，若R>N，则直接截取R点DFT的前N点，若R 4．实验内容： (1)计算一个实数序列[]1,0256 =≤≤的1024点FFT，注意使用将此序列组 x n n 合成一复数序列后再计算的方法。 (1)实验分析： 求某实信号y(n)的复谱,可认为是将实信号加上数值为零的虚部变成复信号(x(n)+j0),再用FFT求其离散付里叶变换。这种作法很不经济,因为把实序列变成复序列,存储器要增加一倍,且计算机运行时,即使虚部为零,也要进行涉及虚部的运算,浪费了运算量。合理的解决方法是利用复数据FFT 对实数据进行有效计算。 用一个N点的FFT运算获得一个2N点实序列的DFT 设x(n)是2N点的实序列,现人为地将x(n)分为偶数组x1(n)和奇数组x2(n) x1(n)=x(2n) n=0,1,…,N-1 x2(n)=x(2n+1) n=0,1,…,N-1 然后将x1(n)及x2(n)组成一个复序列y(n)=x1(n)+jx2(k) 通过N点FFT运算可得到 Y(k)=X1(k)+jX2(k) 根据前面的讨论，得到 为求2N点x（n）所对应的X（k），需求出X（k）与X1（k），X2（k）的关系 而 所以 X(k)=X1(k)+nk W X2(k)。 N 这样，由x1(n)及x2(n)组成复序列，经FFT运算求得Y(k)后，再利用共轭对称性求出X1(k),X2(k)，最后利用上式求出X(k),从而达到了用一个N点的FFT计算一个2N点实序列DFT的目的。 (2)matlab源程序： function [ xk ] = gongshijun2( ) x1=[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,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,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,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,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,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,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,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]; x2=[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,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,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,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,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,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,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]; yn=x1+1i*x2; yk=fft(yn,512); x3=x1(end:-1:1); x4=x2(end:-1:1); yn=x3+1i*x4; ynk=fft(yn,512); x1k=0.5*(yk+conj(ynk)); x2k=-(1i)*0.5*(yk-conj(ynk)); xk=x1k+x2k*exp(1i*360/1024*512); m=0:length(xk)-1; plot(m,xk,'b-'); end (3)实验结果截图： (4)实验图截图： (2) 分别计算两个实数序列 5() c o s , 128 16 x n n n π =≤≤和5()sin ,012816 x n n n π =≤≤的128点FFT ，注意使用将此二序列组合成一复数序列后再计算的方法。 (1)实验分析： 一个N 点FFT 同时计算两个N 点实序列的DFT 设x1(n),x2(n)是彼此独立的两个N 点实序列,且X1(k)=DFT[x1(n)]，X2(k)=DFT[x2(n)] 可通过一次FFT 运算同时获得X1(k),X2(k)。算法如下： 首先将x1(n),x2(n)分别当作一复序列的实部及虚部,令 x(n)=x1(n)+jx2(n) 通过FFT 运算可获得x(n)的DFT 值 X(k)=DFT[x1(n)]+jDFT[x2(n)]=X1(k)+jX2(k) 利用离散付里叶变换的共轭对称性 有了x(n)的FFT 运算结果X(k),由上式即可得到X1(k),X2(k)的值。 则，由题目的设 x1(n)=5()cos ,012816x n n n π=≤≤,x2(n)=5()sin ,012816 x n n n π=≤≤ 组合成复数序列x(n)=1280,16 5sin 165cos ≤≤+n n j n π π，根据实验原理就可以求得相应的结果。 (2)matlab 源程序： function [gg,aa] = gongshijun1() n=0:128; xn=cos(56.25*n)+1i*sin(56.25*n); xk=fft(xn,128); xnk=cos(56.25*(128-n))+1i*sin(56.25*(128-n)); xnkk=fft(xnk,128); gg=0.5*(xk+conj(xnkk)); aa=-(1i)*0.5*(xk-conj(xnkk)); m=0:length(gg)-1; plot(m,gg,'k:',m,aa,'b-'); end (3)实验结果截图： （4）实验图截图： (3) 利用DFT 的方式计算下面两序列的线性卷积： g[n]={3, 4, -2, 0, 1, -4},h[n]={1, -3, 0, 4, -2, 3} (1)实验分析： 用FFT 法也就是用圆周卷积来代替线性卷积，为了不产生混跌，其必要条件是使用g(n),h(n)都补零值点，补到至少N=M+L-1(M=6，L=6)，即 116,05 0),({)(≤≤≤≤=n n n g n g ,11 6,05 0),({)(≤≤≤≤=n n n h n h 然后计算圆周卷积 )()()(n h n x n y ?= 这时，y(n)就能代表线性卷积的结果。 用FFT 计算y(n)值得步骤如下： (1)求H(k)=DFT[h(n)],N 点； (2)求G(k)=DFT[g(n)],N 点； (3)计算Y(k)=H(k)G(k)； (4)求y(n)=IDFT[Y(k)],N点。 (2)matlab源程序： function gg = gongshijun( ) g= [3,4,-2,0,1,-4,0,0,0,0,0,0]; h= [1,-3,0,4,-2,3,0,0,0,0,0,0]; Gk= fft(g,11); Hk= fft(h,11); Yk= Gk.*Hk; gg= ifft(Yk,11); n=0:length(gg)-1; plot(n,gg,'k.'); (2)实验结果和实验图截图： 班级：学号：姓名 实验二FFT算法的MATLAB实现 （一）实验目的： （1）掌握用matlab进行FFT在数字信号处理中的高效率应用。 （2）学习用FFT对连续信号和时域离散信号进行谱分析。 （二）实验内容及运行结果： 题1：若x(n)=cos(nπ/6)是一个N=12的有限序列，利用MATLAB计算它的DFT 并进行IDFT变换同时将原图与IDFT变换后的图形进行对比。当求解IFFT变换中，采样点数少于12时，会产生什么问题。 程序代码： N=12; n=0:11; Xn=cos(n*pi/6); k=0:11; nk=n'*k; WN=exp(-j*2*pi/N) WNnk=WN.^nk XK=Xn*WNnk; figure(1) stem(Xn) figure(2) stem(abs(XK)) 运行结果： IFFT变换中，当采样点数少于12时图像如下图显示： 分析：由图像可以看出，当采样点数小于12时，x(n)的频谱不变，周期为6，而XK 的频谱图发生改变。 题2：对以下序列进行谱分析 132()()103()8470x n R n n n x n n n =+≤≤?? =-≤≤??? 其他n 选择FFT 的变换区间N 为8和16点两种情况进行频谱分析，分别打印其幅频特 性曲线并进行对比、分析和讨论。 ㈠ 程序代码： x=ones(1,3);nx=0:2; x1k8=fft(x,8); F=(0:length(x1k8)-1)'*2/length(x1k8); %进行对应的频率转换 stem(f,abs(x1k8));%8点FFT title('8点FFTx_1(n)'); xlabel('w/pi'); ylabel('幅度'); N=8时： 按时间抽取的基2FFT 算法分析及MATLAB 实现 一、DIT-FFT 算法的基本原理 基2FFT 算法的基本思想是把原始的N 点序列依次分解成一系列短序列，充分利用旋转因子的周期性和对称性，分别求出这些短序列对应的DFT ，再进行适当的组合，得到原N 点序列的DFT ，最终达到减少运算次数，提高运算速度的目的。 按时间抽取的基2FFT 算法，先是将N 点输入序列x(n)在时域按奇偶次序分解成2个N/2点序列x1(n)和x2(n)，再分别进行DFT 运算，求出与之对应的X1(k)和X2(k)，然后利用图1所示的运算流程进行蝶形运算，得到原N 点序列的DFT 。只要N 是2的整数次幂，这种分解就可一直进行下去，直到其DFT 就是本身的1点时域序列。 图1 DIT-FFT 蝶形运算流图 二、DIT-FFT 算法的运算规律及编程思想 1.原位计算 [ 对N=M 2点的FFT 共进行M 级运算，每级由N/2个蝶形运算组成。在同一级中，每个蝶的输入数据只对本蝶有用，且输出节点与输入节点在同一水平线上，这就意味着每算完一个蝶后，所得数据可立即存入原输入数据所占用的数组元素(存储单元)，经过M 级运算后，原来存放输入序列数据的N 个存储单元中可依次存放X(k)的N 个值，这种原位(址)计算的方法可节省大量内存。 2.旋转因子的变化规律 N 点DIT ―FFT 运算流图中，每个蝶形都要乘以旋转因子p W N ，p 称为旋转因子的指数。例如N ＝8 ＝3 2 时各级的旋转因子： 第一级：L=1， 有1个旋转因子：p W N =J /4W N =J 2L W J=0 第二级：L=2，有2个旋转因子：p W N =J /2W N =J 2L W J=0,1 第三级：L=3，有4个旋转因子：p W N =J W N =J 2L W J=0,1,2,3 对于N ＝M 2的一般情况，第L 级共有1-L 2个不同的旋转因子： p W N =J 2L W J=0,1,2,… ,1-L 2 －1 L 2=M 2×M -L 2= N ·M -L 2 故： 按照上面两式可以确定第L 级运算的旋转因子 \ 3、同一级中，同一旋转因子对应蝶形数目 第L 级FFT 运算中，同一旋转因子用在L -M 2 个蝶形中； 4、同一级中，蝶形运算使用相同旋转因子之间相隔的“距离” 第L 级中，蝶距：D=L 2； 5、同一蝶形运算两输入数据的距离 在输入倒序，输出原序的FFT 变换中，第L 级的每一个蝶形的2个输入数据相距：B=1 -L 2。 6、码位颠倒 输入序列x(n)经过M 级时域奇、偶抽选后，输出序列X(k)的顺序和输入序列的顺序关系为倒位关系。 ' 1 引言 CCS(Code Composer Studio)是TI公司的DSP集成开发环境。它提供了环境配置、源文件编辑、程序调试、跟踪和分析等工具,帮助用户在一个软件环境下完成编辑、编译链接、调试和数据分析等工作。与TI提供的早期软件开发工具相比,利用CCS能够加快软件开发进程,提高工作效率。CCS一般工作在两种模式下:软件仿真器和与硬件开发板相结合的在线编程。前者可以脱离DSP芯片,在PC机上模拟DSP指令集与工作机制,主要用于前期算法实现和调试。后者实时运行在DSP芯片上,可以在线编制和调试应用程序。 2 C语言和汇编语言的混合编程 TMS320 C5000系列的软件设计通常有三种方法: (1) 用C语言开发; (2) 用汇编语言开发; (3) C和汇编的混合开发。 其中用C语言开发具有兼容性和可移植的优点,有利于缩短开发周期和减少开发难度,但是在运算量较大的情况下,C代码的效率还是无法和手工编写的汇编代码的效率相比,比如FFT运算,用汇编语言开发的效率高,程序执行速度快,而且可以合理利用芯片的硬件资源,但是开发难度较大,开发周期长,而且可读性和可移植性差。C和汇编的混合编程则可以充分利用前两者的优点,以达到最佳利用DSP资源的目的。但是,采用C和汇编语言混合编程必须遵循相关函数调用规则和寄存器调用规则,否则会给程序的开发带来意想不到的问题。 2.1 C语言和汇编语言混合编程的四种方法 (1) 独立编写汇编程序和C程序,分开编译或汇编成各自的目标代码模块,再用链接器将二者链接起来。这种方法比较灵活,但是设计者必须自己维护各汇编模块的入口和出口代码,自己计算传递的参数在堆栈中的偏移量,工作量较大,但是能做到对程序的绝对控制。 (2) 在C程序中使用汇编程序中定义的变量和常数。 (3) 在C程序中内嵌汇编语句。这种方法可以实现C语言无法实现的一些硬件控制功能,如修改中断控制寄存器。 (4) 将C语言编译生成相应的汇编代码,手工修改和优化C编译器生成的汇编代码。采用这种方法可以控制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 变换。这样变换以后，总的运算次数就变成N + 2 * (N / 2)2 = N + N2 / 2。继续上面的例子，N =1024 时，总的运算次数就变成了525312 次，节省了大约50% 的运算量。而如果我们将这种“一分为二”的思想不断进行下去，直到分成两两一组的DFT 运算单元，那么N 点的DFT 变换就只需要N * log2N 次的运算，N = 1024 点时，运算量仅有10240 次，是先前的直接算法的1% ，点数越多，运算量的节约就越大，这就是FFT 的优越性。 FFT 的实现可以自顶而下，采用递归，但是对于硬件实现成本高，对于软件实现都不够高效，改用迭代较好，自底而上地解决问题。感觉和归并排序的迭代版很类似，不过先要采用“位反转置换”的方法把Xi 放到合适的位置，设i 和j 互为s = log2N 位二进制的回文数，假设s = 3, i = (110)2 = 6, j = (011)2 = 3, 如果i ≠j , 那么Xi 和Xj 应该互换位置。（关于这个回文数的生成，是很有趣而且是很基本的操作，想当初偶初学C++ 的时候就有这样的习题。）当“位反转置换”完成后，先将每一个Xi 看作是独立的多项式，然后两个两个地将它们合并成一个多项式（每个多项式有2 项），合并实际上是“蝶形运算”（Butterfly Operation, 参考《算法导论》吧^_^），继续合并（第二次的每个多项式有4 项），直到只剩 A1=str2double(get(handles.edit8,'String')); A2=str2double(get(handles.edit9,'String')); F1=str2double(get(handles.edit10,'String')); F2=str2double(get(handles.edit11,'String')); Fs=str2double(get(handles.edit12,'String')); N=str2double(get(handles.edit13,'String')); t=[0:1/Fs:(N-1)/Fs]; x=A1*sin(2*pi*F1*t)+A2*sin(2*pi*F2*t); %信号x的离散值 axes(handles.axes1) %在axes1中作原始信号图 plot(x); grid on m=nextpow2(x);N=2^m; % 求x的长度对应的2的最低幂次m if length(x) 实验2 FFT 算法实现 2.1 实验目的 1、 加深对快速傅里叶变换的理解。 2、 掌握FFT 算法及其程序的编写。 3、 掌握算法性能评测的方法。 2.2 实验原理 一个连续信号)(t x a 的频谱可以用它的傅立叶变换表示为 dt e t x j X t j a a Ω-+∞ ∞-?= Ω)()( （2-1） 如果对该信号进行理想采样，可以得到采样序列 )()(nT x n x a = （2-2） 同样可以对该序列进行z 变换，其中T 为采样周期 ∑+∞∞--=n z n x z X )()( （2-3） 当ωj e z =的时候，我们就得到了序列的傅立叶变换 ∑+∞∞-=n j j e n x e X ωω)()( （2-4） 其中ω称为数字频率，它和模拟域频率的关系为 s f T /Ω=Ω=ω （2-5） 式中的s f 是采样频率。上式说明数字频率是模拟频率对采样率s f 的归一化。同模拟域的情况相似，数字频率代表了序列值变化的速率，而序列的傅立叶变换称为序列的频谱。序列的傅立叶变换和对应的采样信号频谱具有下式的对应关系。 ∑+∞∞--=)2(1)(T m j X T e X a j πωω （2-6） 即序列的频谱是采样信号频谱的周期延拓。从式（2-6）可以看出，只要分析采样序列的频谱，就可以得到相应的连续信号的频谱。注意：这里的信号必须是带限信号，采样也必须满 足Nyquist 定理。 在各种信号序列中，有限长序列在数字信号处理中占有很重要的地位。无限长的序列也往往可以用有限长序列来逼近。对于有限长的序列我们可以使用离散傅立叶变换（DFT ），这一变换可以很好地反应序列的频域特性，并且容易利用快速算法在计算机上实现当序列的长度是N 时，我们定义离散傅立叶变换为： ∑-===10)()]([)(N n kn N W n x n x DFT k X （2-7） 其中N j N e W π 2-=，它的反变换定义为： ∑-=-==10)(1)]([)(N k kn N W k X N k X IDFT n x （2-8） 根据式（2-3）和（2-7）令k N W z -=，则有 ∑-====-10)]([)(|)(N n nk N W z n x DFT W n x z X k N （2-9） 可以得到k N j k N e W z z X k X π2|)()(===-,k N W -是z 平面单位圆上幅角为k N πω2=的点，就是将单位圆进行N 等分以后第k 个点。所以，X(k)是z 变换在单位圆上的等距采样，或者说是序列傅立叶变换的等距采样。时域采样在满足Nyquist 定理时，就不会发生频谱混淆；同样地，在频率域进行采样的时候，只要采样间隔足够小，也不会发生时域序列的混淆。 DFT 是对序列傅立叶变换的等距采样，因此可以用于序列的频谱分析。在运用DFT 进行频谱分析的时候可能有三种误差，分析如下： （1）混淆现象 从式（2-6）中可以看出，序列的频谱是采样信号频谱的周期延拓，周期是2π/T ，因此当采样速率不满足Nyquist 定理，即采样频率T f s /1=小于两倍的信号（这里指的是实信号）频率时，经过采样就会发生频谱混淆。这导致采样后的信号序列频谱不能真实地反映原信号的频谱。所以，在利用DFT 分析连续信号频谱的时候，必须注意这一问题。避免混淆现象的唯一方法是保证采样的速率足够高，使频谱交叠的现象不出现。这就告诉我们，在确定信号的采样频率之前，需要对频谱的性质有所了解。在一般的情况下，为了保证高于折叠频率的分量不会出现，在采样之前，先用低通模拟滤波器对信号进行滤波。 （2）泄漏现象 实际中的信号序列往往很长，甚至是无限长序列。为了方便，我们往往用截短的序列来近似它们。这样可以使用较短的DFT 来对信号进行频谱分析。这种截短等价于给原信号序列乘以一个矩形窗函数。而矩形窗函数的频谱不是有限带宽的，从而它和原信号的频谱进行卷积以后会扩展原信号的频谱。值得一提的是，泄漏是不能和混淆完全分离开的，因为泄露导致频谱的扩展，从而造成混淆。为了减小泄漏的影响，可以选择适当的窗函数使频谱的扩散减到最小。 （3）栅栏效应 因为DFT 是对单位圆上z 变换的均匀采样，所以它不可能将频谱视为一个连续函数。 STM32F103 12-15元左右 本文将以一个实例来介绍如何使用STM32提供的DSP库函数进行FFT。 1.FFT运算效率 使用STM32官方提供的DSP库进行FFT，虽然在使用上有些不灵活（因为它是基4的FFT，所以FFT的点数必须是4^n），但其执行效率确实非常高效，看图1所示的FFT运算效率测试数据便可见一斑。该数据来自STM32 DSP库使用文档。 图1 FFT运算效率测试数据 由图1可见，在STM32F10x系列处理器上，如果使用72M的系统主频，进行64点的FFT运算，仅仅需要0.078ms而已。如果是进行1024点的FFT运算，也才需要2.138ms。 2.如何使用STM32提供的DSP库函数 2.1下载STM32的DSP库 大家可以从网上搜索下载得到STM32的DSP库，这里提供一个下载的地址：https://https://www.doczj.com/doc/e510264500.html,/public/STe2ecommunities/mcu/Lists/cortex_ mx_stm32/DispForm.aspx?ID=30831&RootFolder=%2fpublic%2fST e2ecommunities%2fmcu%2fLists%2fcortex%5fmx%5fstm32%2fST M32F10x%20DSP%20library%2c%20where%20is%20it 2.2添加DSP库到自己的工程项目中 下载得到STM32的DSP库之后，就可以将其添加到自己的工程项目中了。 其中，inc文件夹下的stm32_dsp.h和table_fft.h两个文件是必须添加的。stm32_dsp.h是STM32的DSP库的头文件。 src文件夹下的文件可以有选择的添加（用到那个添加那个即可）。因为我只用到了256点的FFT，所以这里我只添加了cr4_fft_256_stm32.s文件。添加完成后的项目框架如图2所示。 图2 项目框架 2.3模拟采样数据 根据采样定理，采样频率必须是被采样信号最高频率的2倍。这里，我要采集的是音频信号，音频信号的频率范围是20Hz到20KHz，所以我使用的采用频率是44800Hz。那么在进行256点FFT时，将得到44800Hz / 256 = 175Hz的频率分辨率。 为了验证FFT运算结果的正确性，这里我模拟了一组采样数据，并将该采样数据存放到了long类型的lBufInArray数组中，且该数组中每个元素的高16 位存储采样数据的实部，低16位存储采样数据的虚部(总是为0)。 为什么要这样做呢？是因为后面要调用STM32的DSP库函数，需要传入的参数规定了必须是这样的数据格式。 下面是具体的实现代码： 1 /****************************************************************** 2函数名称:InitBufInArray() 3函数功能:模拟采样数据，采样数据中包含3种频率正弦波(350Hz，8400Hz，18725Hz) 4参数说明: 5备注:在lBufInArray数组中，每个数据的高16位存储采样数据的实部， 6低16位存储采样数据的虚部(总是为0) 7作者:博客园依旧淡然（https://www.doczj.com/doc/e510264500.html,/menlsh/） 8 *******************************************************************/ 128点FFT 算法设计方法 X(k)=127 n 0x n =∑（）w nk 128, X （k ）=63n 0x n =∑（2）w 1282nk +63n 0x n =∑（2+1） w 128(2n+1)k =63n 0x n =∑（2）w 64nk +(63n 0x n =∑（2+1） w 64nk )w 1281k =H(k)+G(k) w 1281k , H(k)= 63 n 0h n =∑（）w 64nk ,h(n)=x(2n) G(k)= 63 n 0g n =∑（） w 64nk ,g(n)=x(2n+1) H(k) = 63n 0h n =∑（） w 64nk =7a 0=∑(7b 0h a b =∑（+8）)w 64(a+8b)k H(k)=H(k 0+8k 1)= 7a 0=∑(7b 0 h a b =∑（+8）)w 64(a+8b)(k0+8k1) =7a 0=∑(7b 0 h a b =∑（+8）w 648bk0)w 64a(k0+8k1)w 6464bk1 =7a 0 =∑(7b 0h a b =∑（+8）w 648bk0)w 64a(k0+8k1) =7a 0 =∑(7b 0h a b =∑（+8）w 8bk0)w 64a(k0+8k1), [w 6464bk1=1] =7a 0 =∑(7b 0f b =∑（）w 8bk0)w 64a(k0+8k1) =7a 0=∑(7b 0 f b =∑（）w 8bk0w 64ak0 )w 648ak1) =7a 0=∑(7 b 0f b =∑（）w 8bk0w 64ak0 )w 8ak1) w 8bk0旋转因子是1、-1、 等，可以用加减等运算实现，只有w 64ak0要乘旋转因子。 7b 0f b =∑（） w 8bk0用一个蝶8运算。 F(k0)= 7 b 0f b =∑（）w 8bk0,f(b)=h(a+2b), 红色旋转因子还未找到处理方法，使第二级也变为pe8运算 18W =134689222222 ------＝＋＋＋＋＋ 汇编下的FFT算法实现 ;----------------------------------------------------------- ; 快速付里叶变换子程序 ; ; 入口 : 一维四字节浮点数数组的首地址 ; ; ; 以 2 为基β的值,信号抽样的点数 = 2^β; ; ; 出口 : 在原数组的位置保存结果的值; ; ; 补充说明 : ; ; 1. 该子程序所需 RAM 的容量为 2^β*12 字节,另外需要 ; ; 少量的堆栈作为临时变量的存储空间; ; ; 2. 所需 RAM 空间以输入的首地址为基址,向增加的方向 ; ; 扩展; ; ; 应用举例 : ; ; PUSH #0A000H ; 数组首地址压栈 ; PUSH #6 ; β值压栈 ; CALL FFT ; ;----------------------------------------------------------- PROC FFT FFT: LD FFT_CX,2[SP] LD FFT_AX,#0001H ; SHL FFT_AX,FFT_CL ; 计算采样点数 PUSH FFT_AX ; 将采样点数压栈 SHL FFT_AX,#1 ; LD FFT_BX,6[SP] ; 根据采样点数, ADD FFT_BX,FFT_AX ; 计算出其余 PUSH FFT_BX ; 三个被占用 ADD FFT_BX,FFT_AX ; 空间的首地址; PUSH FFT_BX ; 并依次将首地址 SHL FFT_AX,#1 ; ADD FFT_BX,FFT_AX ; 压栈; PUSH FFT_BX ; LD FFT_CX,6[SP] LD FFT_FX,#2 LOOP1: CLR FFT_AX SUB FFT_EX,FFT_CX,#1 LD FFT_BX,FFT_EX LD FFT_DX,10[SP] LOOP2: FFT算法实现实验报告 辛旸 PB10210006 实验目的 1、加深对快速傅里叶变换的理解。 2、掌握FFT 算法及其程序的编写。 3、掌握算法性能评测的方法。 实验内容 1.编写自己的FFT算法： 代码如下： function [ X ] = Sampling( x,N ) %myFFT实现FFT时域取样算法 % 输出：生成FFT序列X(k),输入：欲变换序列x（n），FFT变换长度N（可缺省）（1） if ~exist('N','var'); %检查是否有变换长度N输入 （2） N=length(x); %若无，则令N等于序列长度 （3） end （4） if N （18） w=zeros(1,N); （19） for m=0:1:N-1; %确定反序序列k1和正序序列k的关系 （20） k=m; （21） for n=i-1:-1:0; %从高位开始依次将各位移至反序位 （22） k1(m+1)=k1(m+1)+fix(k/(2^n))*(2^(i-1-n)); （23） k=rem(k,2^n); （24） end; （25） end （26） for l=1:1:N; （27）X(k1(l)+1)=x(l); %生成反序序列 （28） w(l)=exp(-1i*2*pi/N*(l-1)); %生成旋转因子 （29） end （30） for l=0:1:i-1; %控制FFT运算级数 （31） for m=1:1:N; %每一级中有N/2个蝶形运算 （32） if rem((m-1),2^(l+1))<2^l; %找到蝶形运算的上半部分（33） b=X(m)+X(m+2^l)*w(2^(i-1-l)*rem((m-1),2^l)+1); %将结果暂存至b （34）X(m+2^l)=a(m)-X(m+2^l)*w(2^(i-1-l)*rem((m-1),2^l)+1); （35）X(m)=b; %实现原位运算 （36） end （37） end （38） end 2.选择实验1中的典型信号序列验证算法的有效性： 为方便比较两个算法，编写了myCompare函数计算两种算法的运行时间，并绘制频谱曲线 代码如下： function [ t1,t2,e ] = myCompare( x,N ) %myCompare函数：比较自己编写的算法与系统自带算法的差异 % 输入：与变换信号序列x（n）和欲变换长度N % 输出：自己编写的函数的执行时间t1，系统自带函数的执行时间t2，两者计算序列的差异平方和e tic; X1=myFFT(x,N); t1=toc; tic; X2=fft(x,N); 基2的DIT蝶形算法源代码及注释如下： /************FFT***********/ //整个程序输入和输出利用同一个空间x[N]，节约空间 #include DSP课程设计报告 设计题目：256点FFT 院系：计算机科学学院 专业：自动化 年级：2008级 姓名： 学号： 指导教师： 页脚内容0 2011年11 月28日 256点FFT的实现 一、设计目的 1、加深对DFT算法原理和基本性质的理解； 2、熟悉FFT的算法原理和FFT子程序的算法流程和应用； 3、学习用FFT对连续信号和时域信号进行频谱分析的方法； 4、学习DSP中FFT的设计和编程思想； 5、学习使用CCS的波形观察器观察波形和频谱情况； 二、设计内容 给定256 采样点,求频谱,统计运行时间并在PC 上显示。 三、设计原理 快速傅里叶变换（FFT）是一种高效实现离散傅里叶变换（DFT）的快速算法，是数字信号处理中最为重要的工具之一，它在声学，语音，电信和信号处理等领域有着广泛的应用。 页脚内容1 快速傅里叶变换FFT 旋转因子WN 有如下的特性。 对称性：WNk+N/2=-WNk 周期性：WNn(N-k)=WNk(N-n)=WN-nk 利用这些特性，既可以使DFT中有些项合并，减少了乘法积项，又可以将长序列的DFT分解成几个短序列的DFT。FFT就是利用了旋转因子的对称性和周期性来减少运算量的。 FFT的算法是将长序列的DFT分解成短序列的DFT。例如：N为偶数时，先将N点的DFT分解为两个N/2点的DFT，使复数乘法减少一半：再将每个N/2点的DFT分解成N/4点的DFT，使复数乘又减少一半，继续进行分解可以大大减少计算量。最小变换的点数称为基数，对于基数为2的FFT算法，它的最小变换是2点DFT。 一般而言，FFT算法分为按时间抽取的FFT（DITＦＦＴ）和按频率抽取的ＦＦＴ（DIF FFT）两大类。ＤIF FFT算法是在时域内将每一级输入序列依次按奇／偶分成２个短序列进行计算。而DIF FFT算法是在频域内将每一级输入序列依次奇／偶分成２个短序列进行计算。两者的区别是旋转因子出现的位置不同，得算法是一样的。在DIF FFT算法中，旋转因子出现在输入端，而在DIF FFT算法中它出现在输入端。假定序列x(n)的点数N是2的幂，按照DIF FFT算法可将其分为偶序列和奇序列。 偶序列：x(2r)=x1(r) 奇序列：x(2r+1)=x2(r) 其中：r=0,1,2,…,N/2-1 则x(n)的DFT表示为 页脚内容2 毕业设计 [论文] 题目：FFT算法研究及基2-FFT算法的C语言实现学院：电气与信息工程学院 专业：电气工程及其自动化 姓名：XXX 学号：XXXXXX 指导老师：XXX 完成时间：2015年06月01日 摘要 离散傅立叶变换（DFT）常常用于计算信号处理。DFT算法可以得到信号的频域特性，因为该算法在计算上是密集的，长时间的使用时，计算机不能实时进行信号处理。所以DFT被发现之后的相当长时间内是没被应用到实际的项目。到了二十世纪六十年代中期一种新的计算方法被研究者发现出来，它就是FFT。FFT 并不是一种新的获取频域特征的方式，而是离散傅里叶变换的一种快速实现算法。 数字信号处理在当今科技发展中发展很迅速，不但是在传统的通信领域，其他方面也经常用到。利用快速傅里叶变换，实现了信号频域的变换处理。对于信号的处理，往往会和数学中的算法联系到一起。如果处理得当，将会对气象，地理信息等的发展，起到举足轻重的作用，同时对世界其他领域的发展有很大的促进作用。 关键词: FFT算法，C语言，编译实现 Abstract Discrete Fourier Transform (DFT) is often used to calculate the signal processing to obtain frequency domain signals. DFT algorithm can get the frequency domain characteristics of the signal, because the algorithm is computationally intensive, long-time use, the computer is not conducive to real-time signal processing. So DFT since it was discovered in a relatively long period of time is not to be applied to the actual projects until a new fast discrete Fourier calculation method --FFT is found in discrete Fourier transform was able to actually project has been widely used. FFT is not a new way to get the frequency domain, but the discrete Fourier transform of a fast algorithm. Fast Fourier Transform (FFT) is a digital signal processing important tool that the time domain signal into a frequency-domain signal processing. matched filtering has important applications. FFT is a discrete Fourier transform (DFT) is a fast algorithm, it can be a signal from the time domain to the frequency domain. Some signals in the time domain is not easy to observe the characteristics of what is, but then if you change the signal from the time domain to the frequency domain, it is easy to see out of. This design is required to be familiar with the basic principles of FFT algorithm, based on the preparation of C language program to achieve FFT algorithm real number sequence. Keywords: FFT algorithm, C language compiler to achieve DSP课程设计报告 设计题目256点FFT 院系：计算机科学学院 专业：__________ 自动化 年级：2008级 姓名：__________________________ 学号：__________________________ 指导教师： _______________________ 2011年11月28日 256 点FFT 的实现 一、设计目的 1、加深对DFT算法原理和基本性质的理解； 2、熟悉FFT的算法原理和FFT子程序的算法流程和应用； 3、学习用FFT对连续信号和时域信号进行频谱分析的方法； 4、学习DSP中FFT的设计和编程思想； 5、学习使用CCS的波形观察器观察波形和频谱情况； 二、设计内容 给定256 采样点,求频谱,统计运行时间并在PC 上显示。 三、设计原理 快速傅里叶变换(FFT是一种高效实现离散傅里叶变换( DFT的快速算法，是数字信 号处理中最为重要的工具之一，它在声学，语音，电信和信号处理等领域有着广泛的应用。 快速傅里叶变换FFT 旋转因子WN有如下的特性。 对称性：WNk+N/2=-WNk 周期性：WNn(N-k)=WNk(N-n)=WN-nk 利用这些特性，既可以使DFT中有些项合并，减少了乘法积项，又可以将长序列的DFT 分解成几个短序列的DFT。FFT就是利用了旋转因子的对称性和周期性来减少运算量的。 FFT的算法是将长序列的DFT分解成短序列的DFT＞例如：N为偶数时，先将N点的DFT 分解为两个N/2点的DFT,使复数乘法减少一半：再将每个N/2点的DFT分解成N/4点的DFT，使复数乘又减少一半，继续进行分解可以大大减少计算量。最小变换的点数称为基数，对于基数为2的FFT 算法，它的最小变换是2点DFT。 一般而言，FFT算法分为按时间抽取的FFT( DIT FFT)和按频率抽取的FFT ( DIF FFT) 两大类。D IF FFT算法是在时域内将每一级输入序列依次按奇/偶分成2个短序列进行计算。而DIF FFT算法是在 频域内将每一级输入序列依次奇/偶分成2个短序列进行计算。两者的区别是旋转因子出现的位置不同，得算法是一样的。在DIF FFT算法中，旋转因子出现 在输入端，而在DIF FFT 算法中它出现在输入端。 假定序列x(n)的点数N是2的幕，按照DIF FFT算法可将其分为偶序列和奇序列。 偶序列：x(2r)=x1(r) 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 FFT的C语言算法实现 程序如下： /************FFT***********/ #include void divi(complex ,complex ,complex *);/*复数除法*/ void output(); /*输出结果*/ complex x[N], *W;/*输出序列的值*/ int size_x=0;/*输入序列的长度，只限2的N次方*/ double PI; int main() { int i,method; system("cls"); PI=atan(1)*4; printf("Please input the size of x:\n"); /*输入序列的长度*/ scanf("%d",&size_x); printf("Please input the data in x[N]:(such as:5 6)\n"); /*输入序列对应的值*/ for(i=0;i实验二 FFT算法的MATLAB实现
按时间抽取的基2FFT算法分析及MATLAB实现
FFT的定点DSP实现
FFT-C快速傅里叶变换超级详细的原代码
用matlab实现fft算法
实验2FFT算法实现
(完整word版)stm32F103进行FFT算法教程
128点FFT算法设计方法
汇编语言实现的FFT算法
FFT算法实现实验报告
fft算法代码注释及流程框图
DSP课程设计报告(256点FFT的实现)
FFT算法研究及基2-FFT算法的C语言实现
DSP课程设计报告(256点FFT的实现)
FFT算法C语言程序代码
FFT的C语言算法实现
相关主题
文本预览