信号频谱的误差分析及MATLAB实现
2014-07-20任淑萍王欣峰
任淑萍，王欣峰
（山西大学 工程学院 电子信息工程系，山西 太原 030013）
0 引言
在现代数字信号处理中，信号的频域分析比时域分析更加具有清晰的物理概念和深刻含义（如检测、滤波、压缩等方面），因而在信息技术领域得到广泛的应用。不同信号的傅里叶变换的理论分析与方法在相关专业教材中都有介绍，但实际需要的分析信号一般没有算术解析式，直接利用公式进行傅里叶分析成为难题。DFT（Discrete Fourier Transform）是一种时域和频域均离散化的傅里叶变换，适合数值计算，且有快速算法更适合实时处理信号，从而使信号的实时处理设备能够得以简化实现，因此DFT是现代信号谱分析的主要工具。在频谱分析中需要对无限长的时域信号进行加窗截断，而时域的加窗会造成频谱泄漏的现象，影响频率分辨率，从而影响频谱分析的精度。本文重点阐述频谱分析过程中可能存在的误差，以两个频率很近的正弦序列为例［1］，介绍了利用DFT分析信号频率分辨率的影响因素。
1 信号频谱分析流程
DFT要求信号时域离散且数据长度有限，而实际中只能在有限长的时间内观测和采集信号样本，且只能利用有限个数据样本来分析信号的频谱。如果时域连续，则须先进行时域采样，对于离散信号，如果序列太长或采样点数过多，计算机存储和DFT计算都很困难，这样就要采用加窗方法截取数据进行DFT运算。对于有限长序列，因其频谱是连续的，DFT只能描述其有限个频点数据，故存在栅栏效应。因此，用DFT分析实际信号的频谱，其结果必定是近似的。即使对所有离散信号进行DFT变换，也只能用有限个频谱数据近似表示连续频谱；倘若是连续信号，则还会出现频谱混叠；如果对离散信号进行了加窗处理，则会因截断效应产生吉伯斯现象。但如果合理选择参数，分析误差完全可以控制在允许范围内，可见利用DFT分析信号的频谱在工程上是完全可行的。
分析信号频谱的基本流程如图1 所示。连续信号x（t）首先经过低通滤波器（LPF）滤波后得到x′（t），再经过A／D变换成数字信号x（n），此时的序列x（n）依然是无限长的，要利用DFT进行分析，必须经过w（n）加窗截断变成有限长的序列xw（n），再经过DFT变换后成为序列的频谱X（k）。
图1 信号频谱分析的基本流程
2 用DFT对连续信号进行频谱分析出现的问题及对策
用有限的离散数据样本进行DFT分析，得到有限个频点的DFT数值，与原信号的频谱肯定会有不同，而这种不同就是分析误差。以下依据信号频谱分析的基本流程，介绍误差形成的原因及减小分析误差的主要措施，为实际分析过程中适当选择参数提供重要的理论依据［2］。
2.1 混叠现象
由于抽样后序列的频谱是原连续信号频谱以抽样频率为周期的周期延拓，对时域有限长信号，抽样频率不能满足抽样定理，因此会出现频谱混叠现象，使抽样后序列的频谱函数不能如实反映原连续信号的频谱，因此产生频谱分析误差。减弱频谱混叠有两种方法：一是提高抽样频率，即减小抽样时间间隔；二是对被抽样信号提前进行抗混叠滤波，将非带限信号变成带限信号。由于一般信号的高频部分是以大于频率一次方的倒数衰减，因而提高抽样频率对减小频谱混叠是有效的，这种方法的缺点是丢失了部分高频分量，优点就是有效地保护了低频分量，使低频分量不因抽样而得到干扰，同时有效减少了抽样点数。
2.2 截断效应
利用DFT处理的是有限时宽序列，若连续信号x（t）在时域无限长，则抽样后的序列x（n）也是无限长的，此时利用DFT是无能为力的。在实际应用中，必须将x（n）先截断处理，使之成为有限长序列，而截断就相当于在时域乘以一个矩形窗函数，数据突然截断，窗内数据并不改变。时域截断会引起频谱展宽和波动，同时会造成频谱混叠，使不同频率分量的频谱产生混叠失真。
在实际信号处理中，对信号的加窗截断是不可避免的，因此窗函数的使用也是不可避免的。为了减小截断效应，其一是增加截断的长度，使窗的宽度加宽；其二是不要使数据突然锐截止，即不要使用矩形窗，而要缓慢截断，使用其他缓变类型的窗函数，如三角窗、汉明窗、汉宁窗、布莱克曼窗等。窗函数的主瓣越窄、旁瓣越小且衰减的越快，泄漏就越小。
2.3 栅栏效应
对加窗截断后的序列进行DFT分析时，DFT长度必须大于或等于窗序列本身的长度，否则会作自动截断处理。实际上DFT运算一般采用FFT算法，其长度取大于或等于加窗序列的2的整数幂，不足长度时进行补零，得到的DFT函数值是对加窗序列的连续谱进行等间隔抽样的结果。这就形如通过一个有缝隙的栅栏去观察连续信号的频谱，很多地方会被栅栏挡住，故称栅栏效应。在加窗后序列的尾部补零可使频谱的取样点变密，减小栅栏缝隙，使原来看不到的谱分量能看得到，减小了栅栏效应，但由于被观察的连续谱并没有发生变化，因此频率分辨率并没有提高，只能说是可视分辨率或者说计算分辨率提高了。如要提高信号的频率分辨率，选择主瓣窄的截断窗会有一定的改善，但谱间干扰会更严重，而最有效的方法是增加原始信号的长度。
3 实例分析
加窗截断造成频谱泄漏的同时还会使频率分辨率降低，如何提高分析信号的频率分辨率，本节实例通过MATLAB仿真给出了很好的解释。“频率分辨率”是信号处理中的基本概念，是指能将信号中的两个频率靠得很近的谱峰区分开的能力。最通用的方法是用两个不同频率的正弦信号来研究分辨率的大小，能分辨的两个正弦信号的频率越靠近，表明其分辨率越高。
设序列x（n）＝cos（0.48πn）＋cos（0.52πn），序列有2个主要频率0.48π和0.52π。本文画出了其10点DFT、64点DFT及10点后序列补零至64点的DFT。
（1）10点序列DFT程序如下：
clear
n＝［0：1：9］；
xn＝cos（0.48＊pi＊n）＋cos（0.52＊pi＊n）；xn1＝xn（1：1：10）；
N＝10；
n＝［0：1：N－1］；
k＝［0：1：N－1］；
WN＝exp（－j＊2＊pi／N）；
nk＝n＇＊k；
WNnk＝WN.＾nk；
Xk＝xn1＊WNnk；
subplot（121）
stem（n，xn1）；
xlabel（＇n＇）；ylabel（＇x（n）＇）；
axis（［0 9 －2.5 2.5］）；
n1＝0：1：9；
subplot（122）
plot（n1，abs（Xk（1：1：10）））
xlabel（＇k＇）；ylabel（＇｜X（k）｜＇）；
10点DFT的正弦序列及频谱图形见图2 。
图2 10点的正弦序列及频谱图形
（2）64点DFT程序如下：
clear
n＝［0：1：9］；
xn＝cos（0.48＊pi＊n）＋cos（0.52＊pi＊n）；xn1＝［xn（1：1：10），zeros（1，54）］；
N＝64；
n＝［0：1：N－1］；
k＝［0：1：N－1］；
WN＝exp（－j＊2＊pi／N）；
nk＝n＇＊k；
WNnk＝WN.＾nk；
Xk＝xn1＊WNnk；
subplot（121）
stem（n，xn1）；
xlabel（＇n＇）；ylabel（＇x（n）＇）；
axis（［0 63－2.5 2.5］）；
n1＝0：1：63；
subplot（122）
plot（n1，abs（Xk（1：1：64）））
xlabel（＇k＇）；ylabel（＇｜X（k）｜＇）；
64点DFT的正弦序列及频谱图形见图3 。
（3）10点序列补零至64点DFT程序如下：
clear
n＝［0：1：63］；
xn＝cos（0.48＊pi＊n）＋cos（0.52＊pi＊n）；
xn1＝xn（1：1：64）；
N＝64；
n＝［0：1：N－1］；
k＝［0：1：N－1］；
WN＝exp（－j＊2＊pi／N）；
nk＝n＇＊k；
WNnk＝WN.＾nk；
Xk＝xn1＊WNnk；
subplot（121）
stem（n，xn1）；
xlabel（＇n＇）；ylabel（＇x（n）＇）；
axis（［0 63－2.5 2.5］）；
n1＝0：1：63；
subplot（122）
plot（n1，abs（Xk（1：1：64）））
xlabel（＇k＇）；ylabel（＇｜X（k）｜＇）；
10点序列补零后的正弦序列及频谱图形见图4 。
图3 64点的正弦序列及频谱图形
由上例可以看出，要想提高频谱的物理分辨率必须增加样本的长度，而补零的方法只是增加了计算分辨率，平滑了序列的连续谱，可以在一定程度上克服由于加窗导致的频谱泄漏，却没有提高DFT的频率分辨率。因为频率分辨率是由于加窗截断时窗谱的主瓣有一定宽度造成的，因此，分辨相邻很近的两个频率分量的能力（即频率分辨率）完全取决于窗函数窗谱主瓣的宽度，而不取决于进行DFT分析时所用的频点数。主瓣宽度是由输入数据的个数决定的，无论补多少个零都不能有效增加输入数据的个数，因此对分辨率不会产生影响。
4 结论
DFT解决了利用计算机来分析信号频谱的重点和难点。在实际频谱分析［3］过程中，会出现频谱混叠、频谱泄漏、谱间干扰和栅栏效应等。理解了这几种误差形成的原因和可能产生的不良后果，实际分析问题时，只要合理选择分析参数，就会使结果控制在工程所允许范围内。
图4 10点序列补零后的正弦序列及频谱图形
［1］姚天仁.数字信号处理［M］.北京：清华大学出版社，2011.
［2］桂志国，杨明.数字信号处理原理及应用［M］.北京：国防工业出版社，2011.
［3］张登奇，慧银.信号的频谱分析及 MATLAB实现［J］.湖南理工学院学报（自然科学版），2010（3）：29-33.
