1、1.1 市场与市场营销市场与市场营销1.2 我国汽车市场的发展与现状我国汽车市场的发展与现状复习思考题复习思考题实验14 快速傅里叶变换(FFT)一、实验目的一、实验目的(1)加深对快速傅里叶变换(FFT)基本理论的理解。(2)了解使用快速傅里叶变换(FFT)计算有限长序列和无限长序列信号频谱的方法。(3)掌握用MATLAB语言进行快速傅里叶变换时常用的子函数。二、实验涉及的二、实验涉及的MATLAB子函数子函数1.fft功能：功能：一维快速傅里叶变换(FFT)。调用格式：调用格式：yfft(x)；利用FFT算法计算矢量x的离散傅里叶变换，当x为矩阵时，y为矩阵x每一列的FFT。当x的长度为2
2、的幂次方时，则fft函数采用基2的FFT算法，否则采用稍慢的混合基算法。yfft(x，n)；采用n点FFT。当x的长度小于n时，fft函数在x的尾部补零，以构成n点数据；当x的长度大于n时，fft函数会截断序列x。当x为矩阵时，fft函数按类似的方式处理列长度。2.ifft功能：功能：一维快速傅里叶逆变换(IFFT)。调用格式：调用格式：yifft(x)；用于计算矢量x的IFFT。当x为矩阵时，计算所得的y为矩阵x中每一列的IFFT。yifft(x，n)；采用n点IFFT。当length(x)n时，将x截断，使length(x)n。3.fftshift功能：功能：对fft的输出进行重新排列，将
3、零频分量移到频谱的中心。调用格式：调用格式：yfftshift(x)；对fft的输出进行重新排列，将零频分量移到频谱的中心。当x为向量时，fftshift(x)直接将x中的左右两半交换而产生y。当x为矩阵时，fftshift(x)同时将x的左右、上下进行交换而产生y。三、实验原理三、实验原理1.用用MATLAB提供的子函数进行快速傅里叶变提供的子函数进行快速傅里叶变换换从理论学习可知，DFT是唯一在时域和频域均为离散序列的变换方法，它适用于有限长序列。尽管这种变换方法是可以用于数值计算的，但如果只是简单的按照定义进行数据处理，当序列长度很大时，则将占用很大的内存空间，运算时间将很长。快速傅里叶
4、变换是用于DFT运算的高效运算方法的统称，FFT只是其中的一种。FFT主要有时域抽取算法和频域抽取算法，基本思想是将一个长度为N的序列分解成多个短序列，如基2算法、基4算法等，大大缩短了运算的时间。MATLAB中提供了进行快速傅里叶变换(FFT)的子函数，用fft计算DFT，用ifft计算IDFT。例例14-1 已知一个长度为8点的时域离散信号，n10，n27，在n04前为0，n0以后为1。对其进行FFT变换，作时域信号及DFT、IDFT的图形。解解 程序如下：n10；n27；n04；nn1：n2；Nlength(n)；xn(nn0)0；%建立时域信号subplot(2，2，1)；stem(n
5、，xn)；title(x(n)；k0：N1；Xkfft(xn，N)；%用FFT计算信号的DFTsubplot(2，1，2)；stem(k，abs(Xk)；title(XkDFT(x(n)；xn1ifft(Xk，N)；%用IFFT计算信号的IDFTsubplot(2，2，2)；stem(n，xn1)；title(x(n)IDFT(Xk)；运行结果如图14-1所示。图图14-1 例14-1用FFT求有限长序列的傅里叶变换例例14-2 将例13-5已知的两个时域周期序列分别取主值，得到x11，1，1，0，0，0，x20，1，2，3，0，0，求时域循环卷积y(n)并用图形表示。解解 本例将例13-5使
6、用DFT处理的计算，改为用FFT和IFFT进行循环卷积。程序如下(作图程序部分省略)：xn10，1，2，3，0，0；%建立x1(n)序列xn21，1，1，0，0，0；%建立x2(n)序列Nlength(xn1)；n0：N1；k0：N1；Xk1fft(xn1，N)；%由x1(n)的FFT求X1(k)Xk2fft(xn2，N)；%由x2(n)的FFT求YkXk1.*Xk2；%Y(k)X1(k)X2(k)ynifft(Yk，N)；%由Y(k)的IFFT求y(n)ynabs(yn)运行结果如图13-5所示，与例13-5用DFT计算的结果一致。2.用用FFT计算有限长序列的频谱计算有限长序列的频谱1)基
7、本概念一个序号从n1到n2的时域有限长序列x(n)，它的频谱X(ejw)定义为它的离散傅里叶变换，且在奈奎斯特(Nyquist)频率范围内有界并连续。序列的长度为N，则Nn2n11。计算x(n)的离散傅里叶变换(DFT)得到的是X(ejw)的N个样本点X(ejwk)。其中数字频率为kd)N2k(k式中：dw为数字频率的分辨率；k取对应(N1)/2到(N1)/2区间的整数。在实际使用中，往往要求计算出信号以模拟频率为横坐标的频谱，此时对应的模拟频率为式中：D为模拟频率的分辨率或频率间隔；Ts为采样信号的周期，Ts1/Fs；定义信号时域长度LNTs。kD)L2k()NT2k(/Tsskk在使用FF
8、T进行DFT的高效运算时，一般不直接用n从n1到n2的x(n)，而是取的主值区间(n0，1，N1)的数据，经FFT将产生N个数据，定位在k0，1，N1的数字频率点上，即对应0，2p。如果要显示p，p范围的频谱，则可以使用fftshift(X)进行位移。(n)x2)频谱的显示及分辨率问题例例14-3 已知有限长序列x(n)1，2，3，2，1，其采样频率Fs10 Hz。请使用FFT计算其频谱。解解 MATLAB程序如下：Fs10；xn1，2，3，2，1；Nlength(xn)；D2*pi*Fs/N；%计算模拟频率分辨率kfloor(N1)/2：(N1)/2)；%频率显示范围对应p，pXfftshi
9、ft(fft(xn，N)；%作FFT运算且移位psubplot(1，2，1)；plot(k*D，abs(X)，o：)；%横轴化成模拟频率作幅度谱title(幅度频谱)；xlabel(rad/s)；subplot(1，2，2)；plot(k*D，angle(X)，o：)；%横轴化成模拟频率作相位谱title(相位频谱)；xlabel(rad/s)；程序运行结果：absX 0.3820 2.6180 9.0000 2.6180 0.3820angleX 1.2566 2.5133 KG*4/50 2.5133 1.2566运行结果如图14-2所示。图14-2 例14-3有限长序列的频谱由图14-2
10、可知，当有限长序列的长度N5时，频谱的频率样本点数也为5，如图上用“”表示的点位。频率点之间的间距非常大，即分辨率很低。即使使用了plot命令的插值功能，显示出的曲线仍是断断续续的，与真实曲线有较大的误差。改变分辨率的基本方法是给输入序列补零，即增加频谱的密度。注意，这种方法只是改善了图形的视在分辨率，并不增加频谱的细节信息。将上述有限长序列x(n)1，2，3，2，1末尾补0到N1000点，将程序改为：Fs10；N1000；xn1，2，3，2，1；Nxlength(xn)；xn1，2，3，2，1，zeros(1，NNx1)；D2*pi*Fs/N；%计算模拟频率分辨率kfloor(N1)/2：(
11、N1)/2)；%频率显示范围对应p，pXfftshift(fft(xn，N)；%作FFT运算且移位psubplot(1，2，1)；plot(k*D，abs(X)；%横轴化成模拟频率作幅度谱title(幅度频谱)；xlabel(rad/s)；subplot(1，2，2)；plot(k*D，angle(X)；%横轴化成模拟频率作相位谱title(相位频谱)；xlabel(rad/s)；此时程序执行的结果如图14-3所示。由图可以看出，图形的分辨率提高，曲线几乎是连续的频谱了。图14-3 将例14-2有限长序列末尾补0到N1000时的频谱3)实偶序列如何补0例例14-4 已知一个矩形窗函数序列为采样
12、周期Ts0.5 s，要求用FFT求其频谱。解解 由于该序列是一个实的偶序列，因而补0时需要仔细分析。假定按N32补0，则主值区域在n031，FFT的输入应为5n05n1x(n)5)ones(1,11),Nzeros(1,，6)，ones(1Xn即原来n=5：1的前五个点移到n27：31中去了。下面考虑分别用N=32，64，512，观察不同N值代入对频谱的影响。程序如下，Ts=0.5;C=32,64,512;%输入不同的N值forr0：2；NC(r1)；xnones(1，6)，zeros(1，N11)，ones(1，5)；%建立x(n)D2*pi/(N*Ts)；kfloor(N1)/2：(N1)
13、/2)；Xfftshift(fft(xn，N)；subplot(3，2，2*r1)；plot(k*D，abs(X)；%幅度频谱subplot(3，2，2*r2)；stairs(k*D，angle(X)；%相位频谱end注意：注意：此处相位频谱使用了stairs，因为该相位频谱变化率比较陡峭。程序执行结果如图14-4所示。图14-4 将例14-4有限长序列补0到N32、64、512时的频谱如果将x(n)的输入写成xnones(1，11)，zeros(1，N11)；%建立x(n5)相当于起点不是取自n0而是n5，计算的是x(n5)的频谱。幅度频谱不受影响，相位频谱引入一个线性相位5w，如图14-5
14、所示。图14-5 将有限长位移序列x(n5)补0到N512时的频谱3.用用FFT计算无限长序列的频谱计算无限长序列的频谱用FFT进行无限长序列的频谱计算，首先要将无限长序列截断成一个有限长序列。序列长度的取值对频谱有较大的影响，带来的问题是引起频谱的泄漏和波动。例例14-5 已知一个无限长序列为采样频率Fs20 Hz，要求用FFT求其频谱。0n00nex(n)n0.5解解 MATLAB程序如下：Fs20；C8，16，128；%输入不同的N值forr0：2；NC(r1)；n0：N1；xnexp(0.5*n)；%建立x(n)D2*pi*Fs/N；kfloor(N1)/2：(N1)/2)；Xffts
15、hift(fft(xn，N)；subplot(3，2，2*r1)；plot(k*D，abs(X)；axis(80，80，0，3)；subplot(3，2，2*r2)；stairs(k*D，angle(X)；axis(80，80，1，1)；end运行结果如图14-6所示。图14-6 将无限长序列截断为N8，16，128时的频谱由图14-6可见，N值取得越大，即序列保留得越长，曲线精度越高。例例14-6 用FFT计算下列连续时间信号的频谱，并观察选择不同的Ts和N值对频谱特性的影响。xa(t)e0.01t(sin2tsin2.1tsin2.2t)t0解解 该题选择了三个非常接近的正弦信号，为了将各
16、频率成分区分出来，在满足奈奎斯特定理的条件下确定采样周期，选择三组数据，分别是Ts0.5 s、0.25 s和0.125 s；再确定N值，分别选择N256和2048。观察不同Ts和N的组合对频谱的影响。程序如下：T00.5，0.25，0.125，0.125；%输入不同的Ts值N0256，256，256，2048；%输入不同的N值forr1：4；TsT0(r)；NN0(r)；%赋Ts和N值 n0：N1；D2*pi/(Ts*N)；%计算模拟频率分辨率 xaZK(exp(0.01*n*Ts).*(sin(2*n*Ts)sin(2.1*n*Ts)sin(2.2*n*Ts)；kfloor(N1)/2：(N
17、1)/2)；XaTs*fftshift(fft(xa，N)；r，Xa(1)%输出Xa(1)的数值，供误差计算用 subplot(2，2，r)；plot(k*D，abs(Xa)，k)；axis(1，3，1.1*min(abs(Xa)，1.1*max(abs(Xa)；end运行结果如图14-7所示。图14-7 用FFT计算三个很靠近的谐波分量的频谱图由图14-7可以得出以下结论：(1)N同样取256(如前三个图形)，当Ts越大时，时域信号的长度LNTs保留得越长，则分辨率越高，频谱特性误差越小；反之，则分辨率越低，频谱特性误差越大，甚至丢失某些信号分量。(2)Ts相同(如后两个图形)，当N越大时，
18、在0，2p范围内等间隔抽样点数越多，且时域信号的长度LNTs保留得越长，则分辨率越高，频谱特性误差越小；反之，当N越小时，在0，2p范围内等间隔抽样点数越少，则有可能漏掉某些重要的信号分量，称为栅栏效应。四、实验任务四、实验任务(1)阅读并输入实验原理中介绍的例题程序，观察输出的数据和图形，结合基本原理理解每一条语句的含义。(2)已知有限长序列x(n)1，0.5，0，0.5，1，1，0.5，0，要求：用FFT算法求该时域序列的DFT、IDFT的图形；假定采样频率Fs20 Hz，序列长度N分别取8、32和64，使用FFT来计算其幅度频谱和相位频谱。(3)已知一个无限长序列x(n)0.5n(n0)，采样周期Ts0.2 s，要求序列长度N分别取8、32和64，用FFT求其频谱。五、实验预习五、实验预习(1)认真阅读实验原理，明确本次实验任务，读懂例题程序，了解实验方法。(2)根据实验任务预先输入和运行例题程序，编写实验程序。(3)预习思考题：快速傅里叶变换(FFT)与离散傅里叶变换(DFT)有何联系？简述使用快速傅里叶变换(FFT)的必要性。六、实验报告六、实验报告(1)列写调试通过的实验程序，打印或描绘实验程序产生的图形和数据。(2)思考题：回答实验预习思考题。使用MATLAB语言提供的快速傅里叶变换(FFT)有关子函数，进行有限长和无限长序列频谱分析时，需注意哪些问题？
