山东大学《生物医学信号处理》实验3 用FFT作谱分析

实验3 用FFT作谱分析、实验目的（）进一步加深DFT算法原理和基本性质的理解（因为FFT只是DFT的一种快速算法，所以FFT的运算结果必然满足DFT的性质）（）熟悉FFT算法原理及子程序的应用。（）掌握用FFT对连续信号和时域离散信号进行频谱分析的基本方法。了解可能出现的分析误差和原因，以便在实际中正确应用FFT。、实验原理如果用FFT对模拟信号进行谱分析，首先要把模拟信号转换成数字信号，转换时要求知道模拟信号的最高截止频率，以便选择满足采样定理的采样频率。一般选择采样频率是模拟信号中最高频率的34倍。另外要选择对模拟信号的观测时间，如果采样频率和观测时间确定，则采样点数也确定了。这里观测时间和对模拟信号进行谱分析的分辨率有关，最小的观测时间和分辨率成倒数关系。要求选择的采样点数和观测时间大于它的最小值。用FFT作谱分析时，要求做FFT的点数服从的整数幂，这一点在上面选择采样点数时可以考虑满足，即使满足不了，可以通过在序列尾部加完成。如果要进行谱分析的模拟信号是周期信号，最好选择观测时间是信号周期的整数倍。如果不知道信号的周期，要尽量选择观测时间长一些，以减少截断效应的影响。用FFT对模拟信号作谱分析是一种近似的谱分析。首先一般模拟信号（除周期信号外）的频谱是连续频谱，而用FFT作谱分析得到的是数字谱，因此应该取FFT的点数多一些，用它的包络作为模拟信号的近似谱。另外，如果模拟信号不是严格的带限信号，会因为频谱混叠现象引起谱分析的误差，这种情况下可以预先将模拟信号进行预滤，或者尽量将采样频率取高一些。一般频率混叠发生在折叠频率附近，分析时要注意因频率混叠引起的误差。最后要注意一般模拟信号是无限长的，分析时要截断，截断的长度和分辨率有关，但也要尽量取长一些，取得太短因截断引起的误差会很大。举一个极端的例子，一个周期性正弦波，如果所取观察时间太短，例如取小于一个周期，它的波形和正弦波相差太大，肯定误差很大，但如果取得长一些，即使不是周期的整倍数，这种截断效应也会小一些。、实验步骤及内容（）复习DFT的定义、性质和用DFT作谱分析的有关内容。（）复习FFT算法原理与编程思想。（）编制信号产生程序，产生以下典型信号供谱分析用：式中频率自己选择；（）分别以变换区间，对进行FFT，画出相应的幅频特性曲线。（）分别以变换区间，对，进行FFT，画出相应的幅频特性曲线。（）分别以变换区间，对进行FFT，画出相应的幅频特性曲线。（）分别对模拟信号选择采样频率和采样点数。对，周期，频率自己选择，采样频率，观测时间，采样点数用计算。对，选择采样频率，采样点数为，。（）分别将模拟信号转换成序列，用，表示，再分别对它们进行FFT，并画出相应的幅频特性曲线。4、实验用MATLAB函数介绍fft(); figure(); plot(); stem(); abs();title(); xlabel(); ylabel(); text(); hold on; axis(); grid on; subplot(); sin(); cos(); 等。、思考题（）在N=8时，和的幅频特性会相同吗? 为什么? N=16呢?（）如果周期信号的周期预先不知道，如何用FFT进行谱分析?、实验报告要求（）简述实验目的及实验原理。（）编程实现各实验内容，列出实验清单及说明。（）将实验结果和理论分析结果进行比较，分析说明误差产生的原因以及用FFT作谱分析时有关参数的选择方法。并总结实验所得的主要结论。（）简要回答思考题。例程参考用DFT对连续信号作谱分析。已知xa(t)=cos(200*pi*t)+sin(100*pi*t)+cos(50*pi*t);选取不同的截取长度Tp，观察用DFT进行频谱分析时存在的截取效应（频谱泄漏和谱间干扰）。在计算机上用DFT对模拟信号进行谱分析时，只能以有限大的采样频率fs对模拟信号采样。对有限点样本序列（等价于截取模拟信号一段进行采样）作DFT变换得到模拟信号的近似频谱clear;close all;fs=400;T=1/fs;Tp=0.04;N=Tp*fs;N1=N,4*N,8*N;%三种长度0.04s 4*0.04s 8*0.04s%矩形窗截断for m=1:3 n=1:N1(m); xn=cos(200*pi*n*T)+sin(100*pi*n*T)+cos(50*pi*n*T); Xk=fft(xn,4096); fk=fs*0:4095/4096; subplot(3,2,2*m-1); plot(fk,abs(Xk)/max(abs(Xk); if m=1 title('矩形窗截断'); end end %加海明窗截断for m=1:3 n=1:N1(m); wn=hamming(N1(m); xn=(cos(200*pi*n*T)+sin(100*pi*n*T)+cos(50*pi*n*T).*wn' Xk=fft(xn,4096); fk=fs*0:4095/4096; subplot(3,2,2*m) plot(fk,abs(Xk)/max(abs(Xk); if m=1 title('hamming窗截断'); end 3