这是我第二次写课程设计报告，虽说写的不怎么样，但我还较满意。因为我享受了过程，几天忙着找资料，有很努力的去做，想想如果把这种认真放到做别的事情，就算最后没成功也会很高兴，比起不努力去做总是会有成效的。
《数字信号处理》
课程设计报告
信号的频谱分析及MATLAB实现
专 业： 通信工程 班 级： 通信09-2BF 组 次： 第5组 姓 名： 杨慧银 学 号： 14062300513
信号的频谱分析及MATLAB实现
张登奇，杨慧银
(湖南理工学院 信息与通信工程学院，湖南 岳阳 414006)
摘 要：DFT是在时域和频域上都已离散的傅里叶变换，适于数值计算且有快速算法，是利用计算机实现信号频谱分析的常用数学工具。文章介绍了利用DFT分析信号频谱的基本流程，重点阐述了频谱分析过程中误差形成的原因及减小分析误差的主要措施，实例列举了MATLAB环境下频谱分析的实现程序。通过与理论分析的对比，解释了利用DFT分析信号频谱时存在的频谱混叠、频谱泄漏及栅栏效应，并提出了相应的改进方法。
关键词：MATLAB；频谱分析；离散傅里叶变换；频谱混叠；频谱泄漏；栅栏效应
引 言
信号的频谱分析就是利用傅里叶分析的方法，求出与时域描述相对应的频域描述，从中找出信号频谱的变化规律，以达到特征提取的目的
［1］
。不同信号的傅里叶分析理论与方法，在有关专业书中都
有介绍，但实际的待分析信号一般没有解析式，直接利用公式进行傅里叶分析非常困难。DFT是一种时域和频域均离散化的傅里叶变换，适合数值计算且有快速算法，是分析信号的有力工具。下面以连续时间信号为例，介绍利用DFT分析信号频谱的基本流程，重点阐述频谱分析过程中可能存在的误差，实例列出MATLAB环境下频谱分析的实现程序。
1 分析流程
实际信号一般没有解析表达式，不能直接利用傅里叶分析公式计算频谱，虽然可以采用数值积分方法进行频谱分析，但因数据量大、速度慢而无应用价值。DFT在时域和频域均实现了离散化，适合数值计算且有快速算法，是利用计算机分析信号频谱的首选工具。由于DFT要求信号时域离散且数量有限，如果是时域连续信号则必须先进行时域采样，即使是离散信号，如果序列很长或采样点数太多，计算机存储和DFT计算都很困难，通常采用加窗方法截取部分数据进行DFT运算。对于有限长序列，因其频谱是连续的，DFT只能描述其有限个频点数据，故存在所谓栅栏效应。总之，用DFT分析实际信号的频谱，其结果必然是近似的。即使是对所有离散信号进行DFT变换，也只能用有限个频谱数据近似表示连续频谱；如果对离散信号进行了加窗处理，则会因截断效应产生吉伯斯现象；倘若是连续信号，则还会出现频谱混叠。但如果合理选择参数，分析误差完全可以控制在允许范围内，利用DFT分析信号的频谱在工程上是完全可行的
［2］
。分析信号频谱的基本流程如图1所示。
LPFx(t)x'(t)A/Dx(n)w(n)图1 信号频谱分析的基本流程 2 分析误差
利用DFT（实际是用FFT）对连续或离散信号进行频谱分析时，如果信号连续一般要进行采样和截
DFTxw(n)X(k) 断，即使信号离散也往往需要进行加窗截断。用有限的离散数据进行DFT变换，得到有限个DFT数据值，与原信号的频谱肯定不同，这种不同就是分析误差。下面按信号频谱分析的基本流程，分别介绍误差形成的原因及减小分析误差的主要措施，为实际分析过程中适当选择参数提供理论依据。 2.1 混叠现象
对连续信号进行频谱分析时，先要对信号进行采样，理论上要求采样频率fs必须大于两倍信号的最高频率
［3］
。在满足采样定理条件下，采样序列的数字频谱能准确反映连续信号的模拟频谱，否则会
发生频谱混叠现象。严格地讲，实际信号的持续时间有限、频谱无限，为了尽可能减少频谱混叠，信号在采样之前一般都要进行预滤波处理。预滤波也不可能是理想低通，所以频谱混叠不可避免。在实际工作中，为了减小频谱混叠的影响，可通过适当提高防混叠滤波器的指标和适当增大采样频率来实现，采样频率常取信号最高频率的2.5～3倍。各类连续信号采样频率的选取原则与方法可参考文［2］。 2.2 截断效应
利用计算机对离散序列或连续信号的采样序列进行DFT运算时，往往要进行截断，即将离散序列进行加窗处理。对离散序列的加窗实际上是将离散序列与窗函数相乘，加窗后信号的频谱是加窗前信号的频谱与窗函数频谱的卷积，造成截断后信号的频谱与截断前信号的频谱不同，这就是所谓截断效应。截断效应对频谱分析的影响主要表现在两个方面：
1．频谱泄漏 原序列经截断后，频谱会向两边展宽，通常称这种展宽为泄漏。频谱泄漏使频谱变模糊，分辨率变差，泄漏程度和窗函数幅度谱主瓣宽度有关。窗型一定，窗口越长，主瓣越窄，频谱泄漏越小。窗口长度一定，矩形窗主瓣最窄，频谱泄漏最小，但其旁瓣的幅度最大。
2．谱间干扰 对原序列截断，频谱不仅会向附近展宽，还会形成许多旁瓣，引起不同频率间的干扰，简称谱间干扰。特别是强信号谱的旁瓣可能湮没弱信号的主谱或误认为是另一假信号的主谱线。矩形窗的旁瓣幅度大，谱间干扰严重。相对而言，布莱克曼窗的旁瓣幅度比矩形窗小，谱间干扰小，但其主瓣过渡带宽，分辨率差。
采样频率或采样周期是在满足混叠误差前提下选取的，当采样频率或采样周期确定后，适当增加窗口长度有利于减小截断效应。工程上，可用试探法确定窗口长度M ，即将M加倍，分别进行DFT运算，直到相邻两个长度的计算结果接近，取长度较小的M ，这样既可满足截断效应要求，又可使存储单元最小且运算速度最快。如对频率分辨率F0有要求，则窗口长度M可取fs或大于且接近该值的2的整
F0数幂。在窗口长度一定情况下，如果希望引起频谱扩展的过渡带窄，可选矩形窗，但其旁瓣大，谱间干扰严重。若选用布莱克曼窗，旁瓣幅度小，谱间干扰相对较小，但主瓣过渡带更宽，分辨率会进一步下降
［4］
。
2.3 栅栏效应
对加窗后的序列进行DFT运算时，DFT长度必须大于或等于加窗序列的长度，否则会作自动截短处理。实际的DFT运算一般采用FFT算法，其长度取大于或等于加窗序列的2的整数幂，不足进行补零处理，得到的DFT值是对加窗序列的连续谱进行等间隔取样的结果。这就好比通过一个有很多缝隙的栅栏去观察一个连续频谱，很多地方会被栅栏挡住，故称栅栏效应。在加窗序列的尾部补零可使频谱的取样点更密，相当于加密了栅栏的缝隙，使原来看不到的谱分量可能看得到，减小了栅栏效应，但由于被观察的连续谱并没有发生变化，故频率分辨率并没有提高，最多只能说可视分辨率提高了信号的长度来实现。
［5］
。要提高信
号的频率分辨率，选择主瓣窄的截断窗有一定的改善，但谱间干扰会更严重，根本上只能通过增加原始
3 分析实例
对信号进行频谱分析时，由于信号不同，傅里叶分析的频率单位也可能不同，频率轴有不同的定标方式。为了便于对不同信号的傅里叶分析进行对比，这里统一采用无纲量的归一化频率单位，即模拟频率对采样频率归一化；模拟角频率对采样角频率归一化；数字频率对2π归一化；DFT的k值对总点数归一化。同时，为了便于与理论值进行对比，理解误差的形成和大小，这里以确定信号的幅度谱分析为例进行分析说明。假设信号为：x(t)?eu(t)，分析过程：首先利用CTFT公式计算其模拟频谱的理论
?t值；然后对其进行等间隔理想采样，得到x(n)序列，利用DTFT公式计算采样序列的数字连续频谱理论值，通过与模拟频谱的理论值对比，理解混叠误差形成的原因及减小误差的措施；接下来是对x(n)序列进行加窗处理，得到有限长加窗序列xw(n)，再次利用DTFT公式计算加窗后序列xw(n)的数字连续频谱，并与加窗前x(n)的数字连续频谱进行对比，理解截断误差形成的原因及减小误差的措施；最后是对加窗序列进行DFT运算，得到加窗后序列xw(n)的DFT值，它是对xw(n)数字连续频谱进行等间隔采样的采样值，通过对比，理解栅栏效应及DFT点数对栅栏效应的影响。利用MATLAB实现上述分析过程的程序如下：
clc;close all;clear;
%CTFT程序，以x(t)=exp(-t) t>=0 为例 %利用数值运算计算并绘制连续信号波形 L=4, %定义信号波形显示时间长度 fs=4,T=1/fs; %定义采样频率和采样周期
t_num=linspace(0,L,100);%取若干时点,点数决定作图精度 xt_num=exp(-1*t_num);%计算信号在各时点的数值 subplot(3,2,1);plot(t_num,xt_num),%绘信号波形 xlabel('时间(秒)'),ylabel('x(t)'),%加标签 grid,title('(a) 信号时域波形'),%加网格和标题 %利用符号运算和数值运算计算连续信号幅度谱的理论值 syms t W %定义时间和角频率符号对象 xt=exp(-1*t)*heaviside(t),%连续信号解析式 XW=fourier(xt,t,W),%用完整调用格式计算其傅氏变换 %在0两边取若干归一化频点，点数决定作图精度 w1=[linspace(-0.5,0,50),linspace(0,1.5,150)];
XW_num=subs(XW,W,w1*2*pi*fs);%利用置换函数求频谱数值解 mag1=abs(XW_num);%计算各频点频谱的幅值
subplot(3,2,2);plot(w1,mag1),%绘制归一化频率幅值谱线 xlabel('频率(*2*pi*fs)rad/s'),ylabel('幅度'),%加标签 grid,title('(b) 连续信号幅频理论值'),%加网格和标题 %DTFT程序，以x(n)=exp(-nT) n>=0 为例 %利用数值运算计算并绘制离散信号图形 N=L*fs+1;n_num=0:N-1;%生成信号波形采样点
xn_num=exp(-1*n_num*T);%计算信号理想采样后的序列值 subplot(3,2,3);stem(n_num,xn_num,'b.'),%绘序列图形 xlabel('n'),ylabel('x(n)'),%加标签
grid,title('(c) 理想采样图形'),%加网格和标题 %利用符号运算和数值运算计算离散信号幅度谱的理论值 syms n z w %定义符号对象 xn=exp(-n*T), %定义离散信号
Xz=ztrans(xn,n,z),%用完整调用格式计算其Z变换 %利用复合函数计算序列傅里叶变换的解析解 Z=exp(j*w);Hejw=compose(Xz,Z,z,w); Hejw_num=subs(Hejw,w,w1*2*pi);%求频谱数值解