清华大学《数字信号处理》期末作业2013 年 1 月第一题掌握去噪的方法1.1 题目描述MATLAB 中的数据文件noisdopp 含有噪声，该数据的抽样频率未知。

调出该数据，用你学过的滤波方法和奇异值分解的方法对其去噪。

要求：1．尽可能多地去除噪声，而又不损害原信号；2．给出你去噪的原理与方法；给出说明去噪效果的方法或指标；3．形成报告时应包含上述内容及必要的图形，并附上原程序。

1.2 信号特性分析MATLAB所给noisdopp信号极其频域特征如图1.1、图1.2。

图1.1含有噪声的noisdopp信号图1.2 noisdopp 信号频域特性其中横坐标f 采用归一化频率，即未知抽样频率Fs 对应2（与滤波器设计时参数一致）。

信号基本特性是一个幅值和频率逐渐增加的正弦信号叠加噪声，噪声为均匀的近似白噪声，没有周期等特点。

因为噪声信号能量在全频带均匀分布，滤波器截止频率过低则信号损失大，过高则噪声抑制小，认为频谱中含有毛刺较多的部分即为信噪比较小的部分，滤除这部分可以达到较好的滤波效果。

先给定去噪效果的评定指标。

信号开始阶段频率较高（如图1.3，红圈为信号值），一周期内采样点4~5个，即信号归一化频率达到0.4~0.5（Fs=2），难以从频域将这部分信号同噪声分离，滤波后信号损失较大，故对前128点用信噪比考察其滤波效果，定义：22()10lg (()())k k x k SNR y k x k =-∑∑其中，()x k 为原nosidopp 信号，()y k 为滤波后信号。

SNR 越大表示滤除部分能力越小，可以反映滤波后信号对原信号的跟踪能力，对前128点主要考察SNR ，越大滤波器性能越好。

对128点以后的点，信号能量集中在低频部分，滤波后效果不仅仅需考察SNR，同时考察滤波后信号的平滑程度。

定义平滑指标[1]22((1)())((1)())kky k y krx k x k+-=+-∑∑r越小表示滤波后信号相比原信号更加平滑，滤波效果更好。

图1.3信号前100点本文采用整体滤波、分时段滤波和奇异值分解的方法分别对信号进行去噪，并采用信噪比和平滑指数两个指标评定去噪效果。

1.3 滤波器设计根据信号的频谱，设计IIR滤波器和FIR滤波器对信号去噪，用重叠相加法对信号进行卷积运算（每128点截为16段）。

根据信号的频谱，均采用低通滤波方法去噪，将滤波器截止频率设为0.2。

1.3.1 IIR 滤波器滤波设计巴特沃斯低通滤波器，给定参数0.17,0.23,0.1,60p s p s R dB R dB ωω====，由buttord 得到阶数7N =。

用MATLAB 函数butter 设计滤波器，并得到70点频率响应如图1.4，滤波器冲击响应如图1.5。

图1.4 巴特沃斯低通滤波器幅频特性图1.5 巴特沃斯低通滤波器冲击响应采用重叠相加法计算nosidopp信号于巴特沃斯低通滤波器卷积结果。

nosidopp信号每128点分为一段，共16段，得到滤波后信号有1093点。

由图1.5可以得到，设计得到的巴特沃斯滤波器存在30点的延时，取卷积后信号31点开始的1024点为滤波后信号y。

得到滤波后信号如图1.6。

与原信号比较如图1.7。

图1.6滤波后信号图1.7滤波后信号与原信号直观看滤波后信号去除了大部分噪声，但开始高频部分信号损失较多（如图1.7），后尾部信号频率较低，滤波后噪声含量仍比较多。

图1.8信号起始部分前128点信噪比SNR为3.429，后896点信号SNR为12.715，全信号平滑指标为0.0642。

本题中设计的IIR滤波器的问题在于①不是线性相位，只能近似得到滤波后信号的延时，滤波后信号难以同原信号“对齐”②冲击响应无限长，有限长序列通过IIR系统得到的输出也是无限长，存在截断误差③分时段设计滤波器时计算较慢。

1.3.2 窗函数法设计FIR滤波器滤波选择性能较优的汉宁窗设计低通滤波器，截止频率0.2，N取34，得到35点长度的滤波器。

得到滤波器幅频响应和冲击响应如图1.9，图1.10图1.9低通滤波器幅频特性图1.10低通滤波器冲击响应该低通滤波器具有很好的线性相位，输出信号存在N/2=17点的延时，得到的滤波后信号如图1.11。

与原信号比较如图1.12。

图1.11滤波后信号图1.12滤波后信号与原信号直观感觉比巴特沃斯滤波器的滤波结果在信号尾部纹波减小，信号开始部分损失也比较大。

得到前128点信噪比SNR为6.448，后896点信号SNR为14.418，全信号平滑指标为0.0643。

由计算得到的SNR和平滑指标看，窗函数法得到的滤波器滤波效果比巴特沃斯滤波器要好。

采用矩形窗设计同样阶数的FIR滤波器，得到的滤波效果为：前128点信噪比SNR为10.303，后896点信号SNR为15.738，全信号平滑指标为0.184，如图1.13。

图1.13滤波后信号与原信号显然，全信号采用同样的滤波器参数滤波，信号高频分量保持较好，则会保留较多的噪声。

若在频域对信号直接加窗，相当于N=∞的窗函数滤波器，恢复信号有明显的吉布斯现象，如图1.14。

图1.14滤波后信号的吉布斯现象1.3.3 一致逼近法设计FIR滤波器滤波选择性切比雪夫一致逼近法设计低通滤波器，截止频率0.2，N取34，通带、阻带加权均取1，得到35点长度的滤波器。

滤波器幅频响应如图1.15。

图1.15低通滤波器幅频特性具有很好的线性相位，同样，输出信号存在N/2=17点的延时，得到的滤波后信号如图1.16。

与原信号比较如图1.17。

图1.16滤波后信号图1.17滤波后信号与原信号前128点信噪比SNR 为6.167，后896点信号SNR 为13.987，全信号平滑指标为0.0655，与窗函数法得到的滤波效果类似。

对全信号采用同样的滤波器参数滤波，因为没有考虑信号频率随时间的变化，不能得到更好的滤波效果。

1.4 分时段滤波用短时傅里叶变化（128,(128),127,2nfft window hanning overlap Fs ====）估计信号的时频特性大致如图1.18。

图1.18 Nosidopp 信号的时频特性同样将信号分为16段，每段64点，对每段信号做离散傅里叶变化求其频谱，得到幅值最大点频率为max f ，每两点间2/640.0313f ∆==，因为信号的包络为低频分量，需要保留，故设计max 3f f +∆为截止频率的低通滤波器。

采用切比雪夫一致逼近法，N 取34，通带、阻带加权均取1，输出信号存在N/2=17点的延时。

得到的滤波后信号如图1.19。

与原信号比较如图1.20。

图1.19滤波后信号图1.20滤波后信号与原信号图1.21滤波后信号与原信号前100点前128点信噪比SNR为8.144，后896点信号SNR为13.734，全信号平滑指标为0.0773，得到了兼顾信号高频分量和去噪的效果，滤波效果比不分时段要好。

1.5 奇异值分解方法滤波奇异值分解滤波的基本原理是通过矩阵分解方法，得到信号的特征值，特征值中较大的k个分量反应了信号低频成分，而较小的分量反应高频噪声，将较小的特征值置零再恢复信号，得到去噪的信号。

奇异值分解方法滤波需要考虑两个问题①一维信号构成用于奇异值分解的矩阵的方法②如何去除特征值。

本文采用教材式（9.5.18）给定的方式构造矩阵，对长度N的信号x(n)，定义矩阵X1为：1(0)(1)(1)(1)(2)()(1)()(1)x x x Mx x x Mx L x L x N-⎡⎤⎢⎥⎢⎥=⎢⎥⎢⎥--⎣⎦XLLM M M ML得到的矩阵为方阵可最大程度利用其特征值，对nosidopp信号，N=1024，故选择L=513，M=512构造矩阵。

利用svd分解得到特征值，如图1.21。

同时得到矩阵U，V。

图1.22 矩阵特征值可以看出特征值存在明显的转折点，将转折点（第30点）以后的特征值值为0。

恢复信号矩阵：1112111111ˆˆˆ(0)(1)(1)ˆˆˆ(1)(2)()ˆˆˆˆ(1)()(1)T x x x M x x x M x L x L x N -⎡⎤⎢⎥⎢⎥==⎢⎥⎢⎥--⎣⎦X USV L L M M M M L 此时得到的矩阵1ˆX 将不再对称，按上式中对应位置得到滤波后信号ˆ()y x n =为：1ˆˆ()(),0,1,...,1k kxi x i i N k ==-∑ 得到的滤波后信号如图1.22。

与原信号比较如图1.23。

图1.23 滤波后信号图1.24滤波后信号与原信号前128点信噪比SNR为6.145，后896点信号SNR为13.827，全信号平滑指标为0.0708。

直观观察滤波后信号质量好与滤波器滤波后的信号。

特别对信号前100点，如图1.24，直观观察要好于同样信噪比下滤波器去噪的效果。

图1.25滤波后信号与原信号前100点若采用将特征值中小于均值的数置零的方式滤波，得到滤波后信号如图1.26。

图1.26滤波后信号前128点信噪比SNR为10.83，后896点信号SNR为15.357，全信号平滑指标为0.315。

对噪声去除减少，对原信号特征保持较好。

构造矩阵时，减小矩阵的行数、增加矩阵列数，或采用文献[2]中提出的构造矩阵的方法：1(0)(1)(1)(1)(2)()(1)(0)(2)x x x Mx x x Mx N x x N-⎡⎤⎢⎥⎢⎥=⎢⎥⎢⎥--⎣⎦XLLM M M ML不能使得到的特征值拐点给清晰。

1.6 小结在以信噪比SNR和输出信号平滑指标两者作为滤波效果评定标准下，采用切比雪夫一致逼近法分时段设计低通滤波器得到了最好的滤波效果，最优指标为：前128点信噪比SNR为8.144，后896点信号SNR为13.734，全信号平滑指标为0.0773。

直观观察奇异值分解去噪后滤波效果最佳。

本文去噪方式存在的不足有①给定的评定滤波效果的标准不够合理，缺乏标准信号的情况下，信噪比作为标准难以区分估计误差和噪声，给定的标准不能完全反应直观观察出的滤波效果②缺乏对于信号所分段数L 的研究③奇异值分解方法中，对矩阵构造方式对滤波效果的研究不足。

第一题滤波器去噪部分MATLAB 程序见附录I ，奇异值分解去噪部分MATLAB 程序见附录II 。

第二题 掌握功率谱估计的方法2.1 试验数据的产生2.1.1 产生x(n)第一步，产生均值为0，功率为0.12的白噪声序列12(),()u n u n 。

利用教材所给出的式(1.10.2)，可以得到给定功率的白噪声序列，即rand 函数产生的500点均匀分布的白噪声序列，再乘以常数12*12*0.12=1.2a P ==得到12(),()u n u n 。

第二步，通过MATLAB 函数sos2tf 得到5个FIR 子系统级联组成的FIR 系统传递函数的系数firb 和fira ，firb 即为FIR 系统的冲击响应h(n)。