滤波与功率谱估计

1、清 华 大 学数字信号处理期末作业2013 年 1 月第一题 掌握去噪的方法1.1 题目描述MATLAB 中的数据文件 noisdopp 含有噪声，该数据的抽样频率未知。调出该数据，用你学过的滤波方法和奇异值分解的方法对其去噪。要求： 1尽可能多地去除噪声，而又不损害原信号； 2给出你去噪的原理与方法；给出说明去噪效果的方法或指标； 3形成报告时应包含上述内容及必要的图形，并附上原程序。1.2 信号特性分析MATLAB所给noisdopp信号极其频域特征如图1.1、图1.2。图1.1 含有噪声的noisdopp信号图1.2 noisdopp信号频域特性其中横坐标f采用归一化频率，即未知抽样频率

2、Fs对应2（与滤波器设计时参数一致）。信号基本特性是一个幅值和频率逐渐增加的正弦信号叠加噪声，噪声为均匀的近似白噪声，没有周期等特点。因为噪声信号能量在全频带均匀分布，滤波器截止频率过低则信号损失大，过高则噪声抑制小，认为频谱中含有毛刺较多的部分即为信噪比较小的部分，滤除这部分可以达到较好的滤波效果。先给定去噪效果的评定指标。信号开始阶段频率较高（如图1.3，红圈为信号值），一周期内采样点45个，即信号归一化频率达到0.40.5（Fs=2），难以从频域将这部分信号同噪声分离，滤波后信号损失较大，故对前128点用信噪比考察其滤波效果，定义：其中，为原nosidopp信号，为滤波后信号。SNR越大

3、表示滤除部分能力越小，可以反映滤波后信号对原信号的跟踪能力，对前128点主要考察SNR，越大滤波器性能越好。对128点以后的点，信号能量集中在低频部分，滤波后效果不仅仅需考察SNR，同时考察滤波后信号的平滑程度。定义平滑指标1r越小表示滤波后信号相比原信号更加平滑，滤波效果更好。图1.3 信号前100点本文采用整体滤波、分时段滤波和奇异值分解的方法分别对信号进行去噪，并采用信噪比和平滑指数两个指标评定去噪效果。1.3 滤波器设计根据信号的频谱，设计IIR滤波器和FIR滤波器对信号去噪，用重叠相加法对信号进行卷积运算（每128点截为16段）。根据信号的频谱，均采用低通滤波方法去噪，将滤波器截止频

4、率设为0.2。1.3.1 IIR滤波器滤波设计巴特沃斯低通滤波器，给定参数，由buttord得到阶数。用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 滤波后信

5、号图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点长度

6、的滤波器。得到滤波器幅频响应和冲击响应如图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滤波器，得

7、到的滤波效果为：前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 低通滤波器幅频特性具有很好的线性相位，

8、同样，输出信号存在N/2=17点的延时，得到的滤波后信号如图1.16。与原信号比较如图1.17。图1.16 滤波后信号图1.17 滤波后信号与原信号前128点信噪比SNR为6.167，后896点信号SNR为13.987，全信号平滑指标为0.0655，与窗函数法得到的滤波效果类似。对全信号采用同样的滤波器参数滤波，因为没有考虑信号频率随时间的变化，不能得到更好的滤波效果。1.4 分时段滤波用短时傅里叶变化（）估计信号的时频特性大致如图1.18。图1.18 Nosidopp信号的时频特性同样将信号分为16段，每段64点，对每段信号做离散傅里叶变化求其频谱，得到幅值最大点频率为，每两点间，因为信号的

9、包络为低频分量，需要保留，故设计为截止频率的低通滤波器。采用切比雪夫一致逼近法，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个分量反应了信号低频成分，而较小的分量反

10、应高频噪声，将较小的特征值置零再恢复信号，得到去噪的信号。奇异值分解方法滤波需要考虑两个问题一维信号构成用于奇异值分解的矩阵的方法如何去除特征值。本文采用教材式（9.5.18）给定的方式构造矩阵，对长度N的信号x(n)，定义矩阵X1为：得到的矩阵为方阵可最大程度利用其特征值，对nosidopp信号，N=1024，故选择L=513，M=512构造矩阵。利用svd分解得到特征值，如图1.21。同时得到矩阵U，V。图1.22 矩阵特征值可以看出特征值存在明显的转折点，将转折点（第30点）以后的特征值值为0。恢复信号矩阵：此时得到的矩阵将不再对称，按上式中对应位置得到滤波后信号为：得到的滤波后信号如图

11、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。对噪声去除减少，对原信号特征保持较好。构造矩阵时，

12、减小矩阵的行数、增加矩阵列数，或采用文献2中提出的构造矩阵的方法：不能使得到的特征值拐点给清晰。1.6 小结在以信噪比SNR和输出信号平滑指标两者作为滤波效果评定标准下，采用切比雪夫一致逼近法分时段设计低通滤波器得到了最好的滤波效果，最优指标为：前128点信噪比SNR为8.144，后896点信号SNR为13.734，全信号平滑指标为0.0773。直观观察奇异值分解去噪后滤波效果最佳。本文去噪方式存在的不足有给定的评定滤波效果的标准不够合理，缺乏标准信号的情况下，信噪比作为标准难以区分估计误差和噪声，给定的标准不能完全反应直观观察出的滤波效果缺乏对于信号所分段数L的研究奇异值分解方法中，对矩阵构

13、造方式对滤波效果的研究不足。第一题滤波器去噪部分MATLAB程序见附录I，奇异值分解去噪部分MATLAB程序见附录II。第二题 掌握功率谱估计的方法第2章2.1 试验数据的产生2.1.1 产生x(n)第一步，产生均值为0，功率为0.12的白噪声序列。利用教材所给出的式(1.10.2)，可以得到给定功率的白噪声序列，即rand函数产生的500点均匀分布的白噪声序列，再乘以常数得到。第二步，通过MATLAB函数sos2tf得到5个FIR子系统级联组成的FIR系统传递函数的系数firb和fira，firb即为FIR系统的冲击响应h(n)。将第一步得到的白噪声信号分别与h(n)卷积，即得到通过FIR系

14、统后的输出，并截取其中256点数据得到信号。第三步，产生四个复正弦信号，相加，再与v(n)相加，得到x(n)。四个复正弦信号叠加得到的信号幅值和相位图如图2.1，x(n) 幅值和相位图如图2.2。四个复正弦信号叠加信号是周期信号，周期为276。图2.1 复正弦信号波形图2.2 x(n)波形2.1.2 求信号真是功率谱第一步，计算白噪声信号功率谱。第二步，计算FIR系统频率响应。第三步，计算v(n)功率谱。由教材式(10.3.2)给出的随机信号通过线性系统理论，得到功率谱为求v(n)功率谱，先计算自相关函数在计算功率谱，得到第四步，求复正弦信号功率谱。设S(n)为4个复正弦信号的和。S(n)是确

15、定信号，且是功率信号，相求自相关函数再求功率谱按教材式(10.4.4)，用时域傅里叶变换模平方再除以时间长度得到的功率谱为S(n)是周期信号，故上式化为一周期N点内平均，取02等分取4096点转换为DFT若信号S(n)存在截断，信号长度N，傅里叶变换长度M，此时S(n)功率谱为得到信号x(n)的功率谱为取02等分取4096点，离散化的功率谱取上式中N=256，M=4096，得到真实功率谱如图2.3。图2.3 信号x真是功率谱2.1.3 计算给定频率点信噪比求出复正弦信号在四个频率点处的功率。根据教材式(3.2.32)，幅度和频率分别为的复正弦信号的功率为用离散的功率谱计算得到的结果相同，即复正

16、弦信号在功率集中在处，幅度为。噪声信号在处功率为。故在频率处的信噪比为信噪比都比较低。2.2 对x(n)做功率谱估计2.2.1 周期图法估计x(n)的功率谱周期图法估计随机信号功率谱的基本原理是根据的N点观察值做傅里叶变换，得到，然后取其平方，再除以N，作为对其真是功率谱的估计，周期图谱法估计的功率谱为令在02等分取4096点，得到离散化的功率谱对256点做4096点DFT，得到估计的曲线如图2.4。图2.4 周期图法估计x功率谱得到比较好的功率谱估计值，256点实际分辨率为，故四个频率点的谱线均可分开，能量较小处也有所反应。2.2.2 Welch方法改进周期图W点汉明窗表达式为每长度M=64

17、分一段，得到Bartlett法估计的每段功率谱表达式为其中平均后的功率谱作为最后估计的功率谱令在02等分取4096点，用DFT代替DTFT，可得到离散化的功率谱估计。取M=64，每段重叠32点，分7段，用上述Welch法得到的功率谱估计曲线如图2.5。图2.5 Welch法估计x功率谱功率谱较周期谱法更光滑，但是分辨率下降，不能区分0.23、0.24两处频率值的信号。64点汉明窗的主瓣宽度为，故不能分辨两个频率值，若到达0.01的分辨率，至少需400点的窗函数，超过信号长度，等同于周期谱法，故改进的周期谱法不能达到更好的功率谱估计效果。2.2.3 自相关法估计功率谱自相关法估计功率谱的基本原理

18、是，根据信号的N点观察值，计算自相关函数的估计值，再对做傅里叶变化得到估计的功率谱。复信号的自相关函数计算方法为再得到功率谱估计为根据上述公式，取M=63，令在02等分取4096点，用DFT代替DTFT，得到的离散的功率谱计算方式为得到估计的功率谱曲线如图2.6（已舍去负值）。图2.6 自相关法估计x功率谱谱线比周期谱法光滑，但不能分辨0.23、0.24频率处的两个谱线，减小M或使用窗函数可使谱线光滑，但将更加降低频率分辨率。2.2.4 AR模型估计功率谱的自相关法自相关法求解AR模型的基本原理是用信号的N点观察值计算得到的自相关函数估计值代替Yule-Walker方程中的，得到AR模型参数。

19、采用Levinson-Durbin递推算法计算。初始条件递推关系为用MATLAB函数levinson可以求解上述方程参数，得到p阶AR模型参数。自相关求解AR模型的最小预测误差随p阶数增加而减小，p为6时，最小预测误差77.15，p为15时最小预测误差为43.3，取p为15，得到AR模型系数，和最小预测误差，从而得到估计的功率谱N=4096离散得到得到估计的功率谱曲线如图2.7。图2.7 AR模型自相关算法估计x功率谱，p=15基本估计出四处频率点的复正弦信号，但不能分辨0.23和0.24两处频率的谱线，信噪比较低的-0.16频率点处的信号估计不佳，若降低阶数p，则根本估计不出该处的谱线，p取

20、10时得到的功率谱如图2.8。图2.8 AR模型自相关算法估计x功率谱，p=102.2.5 AR模型估计功率谱的Burg算法根据教所给的Burg算法计算。设定初始条件，由得到。由得到m=1时的参数。再由求得，由表达式求得。由求得m=2时参数。重复上述步骤，可得到m=p时参数。根据上述算法，在MATLAB中编写程序实现上述算法。得到p=18时的参数，并利用MATLAB函数arburg验证了计算结果的正确性。得到估计的功率谱如图2.9。图2.9 AR模型Burg算法估计x功率谱，p=18比自相关法有更好的频谱，可以分辨0.23和0.24频率处的谱线，-0.16处的信号也更加明显。但Burg算法在p

21、较小时（如10），得到的功率谱不仅不能分辨较近的频谱和幅度小的信号，且谱线的相对幅度与信号幅度对应发生变化，如图2.10。图2.10 AR模型Burg算法估计x功率谱，p=102.2.6 AR模型估计功率谱的Burg算法根据信息论准则利用自相关算法，k由1开始增加，通过levinson函数求解最小估计误差，计算，得到最小值，即为最合适的阶次p，此时，最小估计误差49.38。p为10时利用自相关算法的AR模型的功率谱估计结果如前2.8图。2.3 小结对本题中信号的功率谱估计不足在于产生的实验信号中的白噪声不是理想白噪声，不同次实验产生的功率不稳定，两者的独立性也不佳缺乏对功率谱估计效果的定量指标

22、。2.12.2.4，2.2.6节相关MATLAB程序见附录III，2.2.5节相关MATLAB程序见附录VI。29参考文献 1陈强, 黄声享, 王韦. 小波去噪效果评价的另一指标J. 测绘信息与工程, 2008(5):13-14. 2周江, 杨浩. 分段滤波在心率变异性功率谱分析中的应用J. 重庆大学学报(自然科学版), 2001(4):95-97. 3海兰萍, 姜占才, 李振起. 一种改进的奇异值分解语音降噪算法J. 科技信息, 2012(6):156-158. 4迟华山, 王红星, 赵培洪, 等. 基于S变换的线性调频信号时频滤波J. 无线电通信技术, 2012(1):21-24. 5DE

23、LL. 常用图像去噪比较与分析Z. 20101. 6奇异值分解降噪的改进方法-张磊彭伟才原春晖刘彦Z. 7胡广书. 数字信号处理M. 第二版. 北京: 清华大学出版社, 2008. 8自适应线性调频信号时变滤波器在分数傅里叶域Z. 9姚文俊. 自相关法和Burg法在AR模型功率谱估计中的仿真研究J. 计算机与数字工程, 2007(10):32-34.10白噪声背景下正弦信号频率和功率估计的有效算法Z.11李军, 王凯. 三种现代功率谱估计方法性能研究J. 科技信息, 2010(30):554-555.附录I% 期末大作业第一题clear all;clc;close allload noisdo

24、pp %读入信号%信号特性x=noisdopp;nx=length(x);plot(x);axis tight;xlabel('n');ylabel('noisdopp');grid onplot(1:100,x(1:100),1:100,x(1:100),'ro');axis tight;xlabel('n');ylabel('noisdopp');grid onx_fft=fft(x);figuresubplot(2,1,1);plot(1:nx/2*2./nx,abs(x_fft(1:nx/2);axis t

25、ight;xlabel('f');ylabel('|X(k)|');subplot(2,1,2);plot(1:nx/2*2./nx,unwrap(angle(x_fft(1:nx/2);axis tight;xlabel('f');ylabel('phase X(k)');% 分段处理信号；每M点分一段M=128/2;% 设计滤波器 %FIR，窗函数法阶数为L，有L+1点，有L/2点的延时 L=34;fc=0.2; %Fs=2 b1=fir1(L,fc,hanning(L+1);figure;freqz(b1);figure;s

26、tem(b1);xlabel('n');ylabel('h(n)');grid on %FIR，切比学夫一致逼近法 F=0,fc-0.03,fc+0.03,1; A=1,1,0,0; W=1,1;b2=remez(L,F,A,W);freqz(b2) %将H(1)置为一% b2_mo=max(abs(fft(b2);% b2=remez(L,F,A,W)./b2_mo; %IIR，巴特沃斯滤波器 Niir,wniir=buttord(fc-0.03,fc+0.03,0.1,60);bb,ab=butter(Niir,fc); %bb,ab=bilinear(bb

27、s,abs,2); L_b3=2*L+1;b3=impz(bb,ab,L_b3);freqz(b3);figure;stem(b3);xlabel('n');ylabel('h(n)');grid on %FIR，窗函数法,吉布斯现象 nb51=1:L/2-1;nb52=L/2+1:L+1;wc_b5=fc*2*pi; b5=sin(nb51-L/2)*wc_b5)./(pi*(nb51-L/2),wc_b5./pi,sin(nb52-L/2)*wc_b5)./(pi*(nb52-L/2);% 用于重建信号y1=zeros(1,nx+length(b1)-1);

28、y2=zeros(1,nx+length(b2)-1);y3=zeros(1,nx+length(b3)-1);y4=zeros(1,nx+length(b1)-1);y5=zeros(1,nx+length(b5)-1);ylw=zeros(1,nx);for n=1:(nx-M)/(M-n_over)+1 xw=x(n-1)*(M-n_over)+1:(n-1)*(M-n_over)+M); %直接频域加窗 xw_fft=fft(xw); xl=xw_fft(1:5),zeros(1,length(xw)-10),xw_fft(length(xw)-4:length(xw); yl=iff

29、t(xl,'symmetric'); ylw(n-1)*(M-n_over)+1:(n-1)*(M-n_over)+M)=ylw(n-1)*(M-n_over)+1:(n-1)*(M-n_over)+M)+yl; %分段调整滤波频率,取频率最大值左右一定带宽的信号 max_fft,n_max_fft=max(xw_fft(1:M/2); %通带中心频率 f0=n_max_fft./M*2; %设置带宽 fw=3;df=1./M*2; if n_max_fft-fw<=0 %即为低通滤波器，FIR，切比学夫一致逼近法 F=0,f0,f0+fw*df,1; b4=remez(

30、L,F,A,W); else %F=0,f0-fw*df,f0-fw*df+df,f0+fw*df-df,f0+fw*df,1; %即为带通滤波器 F=0,f0,f0+fw*df,1;b4=remez(L,F,A,W); end %信号还原 yw5=conv(xw,b5); y5(n-1)*(M-n_over)+1:(n-1)*(M-n_over)+M+length(b5)-1)=y5(n-1)*(M-n_over)+1:(n-1)*(M-n_over)+M+length(b5)-1)+yw5; yw4=conv(xw,b4); y4(n-1)*(M-n_over)+1:(n-1)*(M-n_

31、over)+M+length(b4)-1)=y4(n-1)*(M-n_over)+1:(n-1)*(M-n_over)+M+length(b4)-1)+yw4; yw3=conv(xw,b3); y3(n-1)*(M-n_over)+1:(n-1)*(M-n_over)+M+length(b3)-1)=y3(n-1)*(M-n_over)+1:(n-1)*(M-n_over)+M+length(b3)-1)+yw3.' yw2=conv(xw,b2); y2(n-1)*(M-n_over)+1:(n-1)*(M-n_over)+M+length(b2)-1)=y2(n-1)*(M-n_

32、over)+1:(n-1)*(M-n_over)+M+length(b2)-1)+yw2; yw1=conv(xw,b1); %信号幅值不对 y1(n-1)*(M-n_over)+1:(n-1)*(M-n_over)+M+length(b1)-1)=y1(n-1)*(M-n_over)+1:(n-1)*(M-n_over)+M+length(b1)-1)+yw1; endfigure my1,ny_y1=max(b1);y1_yx=y1(ny_y1:length(y1);plot(y1_yx);axis tight;xlabel('n');ylabel('y')

33、;grid onfigure;plot(1:nx,x,1:nx,y1_yx(1:nx),'r-');axis tight;xlabel('n');ylabel('x,y');grid onmy2,ny_y2=max(b2);y2_yx=y2(ny_y2:length(y2);plot(y2_yx);axis tight;xlabel('n');ylabel('y');grid onfigure;plot(1:nx,x,1:nx,y2_yx(1:nx),'r-');axis tight;xlabel(

34、'n');ylabel('x,y');grid on%figure;%plot(ylw) % 直接频域加窗，等于时域用sinc函数截断，效果不好，吉布斯现象my3,ny_y3=max(b3);y3_yx=y3(ny_y3:length(y3);plot(y3_yx);axis tight;xlabel('n');ylabel('y');grid onplot(1:nx,x,1:nx,y3_yx(1:nx),'r-');axis tight;xlabel('n');ylabel('x,y

35、9;);grid onmy4,ny_y4=max(b4);y4_yx=y4(ny_y4:length(y4);plot(y4_yx);axis tight;xlabel('n');ylabel('y');grid onfigure;plot(1:nx,x,1:nx,y4_yx(1:nx),'r-');axis tight;xlabel('n');ylabel('x,y');grid onmy5,ny_y5=max(b5);y5_yx=y5(ny_y5:length(y5);% plot(y5_yx);axis ti

36、ght;xlabel('n');ylabel('y');grid on% plot(1:nx,x,1:nx,y5_yx(1:nx),'r-');axis tight;xlabel('n');ylabel('x,y');grid on%计算峰值信噪比% y_MSN=y3_yx(129:nx);% y_MSN128=y3_yx(1:128);% y_MSN=y1_yx(129:nx);% y_MSN128=y1_yx(1:128);% y_MSN=y5_yx(129:nx);% y_MSN128=y5_yx(1:128

37、);% y_MSN=y2_yx(129:nx);% y_MSN128=y2_yx(1:128);y_MSN=y4_yx(129:nx);y_MSN128=y4_yx(1:128);x_MSN128=x(1:128);x_MSN=x(129:nx);MSN=(y_MSN-x_MSN)*(y_MSN-x_MSN)'PSNR=10*log10(x_MSN*x_MSN'/MSN);MSN128=(y_MSN128-x_MSN128)*(y_MSN128-x_MSN128)'PSNR128=10*log10(x_MSN128*x_MSN128'/MSN128);%平滑尺度

38、% y_r1=y3_yx(2:nx);y_r2=y3_yx(1:nx-1);% y_r1=y1_yx(2:nx);y_r2=y1_yx(1:nx-1);% y_r1=y5_yx(2:nx);y_r2=y5_yx(1:nx-1);% y_r1=y2_yx(2:nx);y_r2=y2_yx(1:nx-1);y_r1=y4_yx(2:nx);y_r2=y4_yx(1:nx-1);x_r1=x(2:nx);x_r2=x(1:nx-1);r_phx=(x_r1-x_r2)*(x_r1-x_r2)'r_phy=(y_r1-y_r2)*(y_r1-y_r2)'r_ph=r_phy/r_phx

39、;% 用短时傅里叶变化估计信号时频特性% spectrogram(x,128,127,128,2); %计算自相关函数% rx=xcorr(x,x);subplot(3,1,1);plot(rx);xlabel('n');ylabel('r')% ry=xcorr(y,y);subplot(3,1,2);plot(ry);% ry1=xcorr(y1,y1);subplot(3,1,3);plot(ry1);%对不足一周期正弦信号FFT分析效果% nt=0:1/300:127/300;% yt=sin(2*pi*nt);% Eb=(b4*b4');% y

40、tl=conv(yt,b4);% plot(yt',ytl(1:length(yt)')%plot(abs(fft(yt)附录II% 奇异值分解方法降噪clear all;clc;close all;load noisdopp %读入信号%信号特性x=noisdopp;%plot(x);nx=length(x);% LxL=nx方法M=nx/2;X1=zeros(M+1,M);for i=1:M+1 for j=1:M X1(i,j)=x(i+j-1); endendU,S,V=svd(X1);plot(abs(diag(S);axis tight;xlabel('k&

41、#39;);grid on%滤波S1=zeros(i,j);% S1(1:30,1:30)=S(1:30,1:30); %衰减很明显s_u=sum(diag(S)./length(diag(S);s_n=find(diag(S)-s_u>0);s_nq=max(s_n);S1(1:s_nq,1:s_nq)=S(1:s_nq,1:s_nq); %均值下去除%重建,奇异值分解降噪的改进方法A=U*S1*V.'y=0*x;for i=1:nx ytemp=0; if i<=M for j=1:i ytemp=ytemp+A(j,i-j+1); end y(i)=ytemp/i;

42、else for j=1:nx-i+1 ytemp=ytemp+A(M+1-j+1,i-M+j-1); end y(i)=ytemp/(nx-i+1); endendplot(y);axis tight;xlabel('n');ylabel('y');grid on;figureplot(1:nx,x,1:nx,y,'r-');axis tight;xlabel('n');ylabel('x,y');grid on%计算峰值信噪比y_MSN=y(129:nx);y_MSN128=y(1:128);x_MSN128=

43、x(1:128);x_MSN=x(129:nx);MSN=(y_MSN-x_MSN)*(y_MSN-x_MSN)'PSNR=10*log10(x_MSN*x_MSN'/MSN);MSN128=(y_MSN128-x_MSN128)*(y_MSN128-x_MSN128)'PSNR128=10*log10(x_MSN128*x_MSN128'/MSN128);%平滑尺度y_r1=y(2:nx);y_r2=y(1:nx-1);x_r1=x(2:nx);x_r2=x(1:nx-1);r_phx=(x_r1-x_r2)*(x_r1-x_r2)'r_phy=(y_

44、r1-y_r2)*(y_r1-y_r2)'r_ph=r_phy/r_phx;附录III% 第二题clear all;clc;close allN=500; %信号长度%功率为0.12的白噪声产生方式(1.10.2)P=0.12;a=sqrt(12*P);%产生噪声序列u1=rand(1,N)*a;u2=rand(1,N)*a;%5个FIR子系统B1=1,1.98,0.98;B2=1,-1.98,0.98;B3=1,-1.8418,0.98;B4=1,-1.5,0.98;B5=1,-1.2727,0.98;B=B1;B2;B3;B4;B5;sos=B,ones(5,1),zeros(5,

45、2);%得到级联系统传递函数firb,fira=sos2tf(sos);v1=conv(u1,firb);v2=conv(u2,firb);nx=256;v=v1(100:100+nx-1)+1i*v2(100:100+nx-1);%产生复正弦信号exp_a=6,12,12,2;exp_f=0.12,0.23,0.24,-0.16;for i=1:length(exp_a) exp_n=0:nx-1; exp_x(i,1:nx)=exp_a(i)*exp(1i*2*pi*exp_f(i)*exp_n);endexp_s=sum(exp_x);x=exp_s+v;save('x_sv.m

46、at','x') %保存用于burg算法subplot(2,1,1);plot(abs(exp_s);axis tight;xlabel('n');ylabel('|exp|');grid onsubplot(2,1,2);plot(angle(exp_s);axis tight;xlabel('n');ylabel('phase exp/rad');grid on%画x波形图figuresubplot(2,1,1);plot(abs(x);axis tight;xlabel('n');yl

47、abel('|x|');grid onsubplot(2,1,2);plot(angle(x);axis tight;xlabel('n');ylabel('phase x/rad');grid onM=4096;nw=-M/2:M/2-1;fw=nw./M;w=nw*(2*pi./M);bt=zeros(1,length(w);for k=0:10 bt=firb(k+1).*exp(-1i*w*k)+bt;endPv=abs(bt).2*2*P;%Ps=1./length(w)*(abs(fft(exp_s,length(w).2);Ps=z

48、eros(1,length(w);for i=1:4 expndt=find(2*pi*exp_f(i)-w<0); expd=min(expndt); Ps(expd)=length(x).2./length(w)*exp_a(i).2;end% plot(Ps)Px=Pv+Ps;figure;plot(fw,Px);grid on;axis tight;xlabel('f');ylabel('Px');%计算信噪比SNR=10*log10(exp_a.2./2/P);%周期图法功率谱估计%求信号FFTXper=fftshift(fft(x,M);Ppe

49、r=abs(Xper).2./M;figureplot(fw,Pper);grid on;axis tight;xlabel('f');ylabel('Pper');%用Welch法估计mw=64;nxw=length(x);lw=(nxw-mw/2)/(mw/2);%得到汉明窗wnw=0:mw-1;wn_h=0.54-0.46*cos(2*pi*wnw/mw);Uw=1./mw*(wn_h*wn_h');%将x分段Pper_w=zeros(1,M);for i=1:lw xni=(x(i-1)*mw/2+1:(i+1)*mw/2).*wn_h; xni

50、_fft=fft(xni,M); Pwt=1./(M*Uw)*(abs(xni_fft).2); Pper_w=Pper_w+Pwt;endPper_w=1./lw*fftshift(Pper_w)*16;figure;plot(fw,Pper_w);grid on;axis tight;xlabel('f');ylabel('Pper Welch');%自相关法估计m_bt=63;n_bt=length(x);%计算自相关函数r_bt=zeros(1,m_bt);for i=0:m_bt-1 r_bt(i+1)=1./n_bt*sum(x(1:n_bt-i)'.'.*x(1+i:n_bt);endPbt=(2*r