随机信时号处理.doc

参数功率谱估计的实现及其与经典谱估计的比较 参数功率谱估计的实现及其与经典谱估计的比较1、 摘要频域分析是从频率角度对信号进行分析研究，对于确定性信号来说，通常使用傅立叶变换或者傅立叶技术展开得到频域表达，但对随机信号而言，由于其时域波形的随机性及能量无限，没有确定的时域表达式，无法用傅立叶变换直接将其变换到频域中去研究。表达随机信号通常使用概率密度函数，根据维纳-辛钦定理，广义随机过程的功率谱与自相关函数是一对傅立叶变换，所以自然想到可以用功率谱来研究随机过程的频域性质。对随机信号的功率谱估计方法通常分为两大类：经典谱估计和现代谱估计经典谱估计是基于维纳-辛钦定理，从自相关函数出发通过傅立叶变换得到功率谱，而现代谱估计则将随机信号看成白噪声通过一个滤波器的输出。现代谱估计就是通过记录的信号序列估计滤波器参数，从而得到其频率响应，最后通过得到其功率谱。从频率分辨率来看，经典谱估计效果一般不如参数谱估计好，而且参数谱估计的不同算法对不同的采样序列也有不同的效果，本文将对周期图法，基于L-D快速递推算法的Y-W法和Burg算法使用Matlab进行编程实现并做比较。2、 关键字功率谱估计周期图法L-D递推Y-W算法Burg算法3、 原理1. 经典谱估计广义平稳随机信号经典谱估计基于维纳-辛钦定理：广义平稳随机的自相关函数与其功率谱是一堆傅立叶变换。所以要求功率谱，只需由随机序列求出自相关函数然后进行傅立叶变换即可， m=0，1，2，N-1周期图法：虽然经典谱估计方法比较直观简单，但由于随机序列相当于对信号加床，所以求自相关函数后傅立叶变换的功率谱往往受到窗函数影响，不是信号真实谱，所以就产生了以下的现代谱估计。参数功率谱估计：根据谱分解定理，任何平稳随机信号都可以看成是由白噪声经过一个因果稳定可逆系统产生的输出，如上图。其系统函数全部极点都在z平面单位圆内部，将表示成，实际中，很多物理过程都可以近似表示成AR，ARMA，MA过程，其中AR过程为，也就是可以看成白噪声通过系统后的输出。如果根据有限的数据记录（看成系统输出）和已知输入（白噪声）确定出系统的参数，那么所要研究的随机过程就是白噪声通过已知系统之后的输出，于是可以用系统的参数来描述索要研究的物理过程，例如要分析实平稳随机信号的有理功率谱，根据平稳随机信号通过线性系统的性质可知：，其中是白噪声的功率谱（为常数），是的能量谱对线性时不变的稳定系统，单位脉冲响应是确定性能量信号，其能量谱是频谱幅度的平方。上式说明平稳随机信号的功率谱可以用模型的参数来表示。所以只要建立合适的信号模型，并估计出模型的参数，就可以分析随机信号的功率谱。上面仅提到了AR模型，是因为柯尔莫哥洛夫定理表明，任何ARMA或MA过程都可以用无限阶的AR过程表示，即使在建模式选择了不正确的模型，只要模型阶数足够高，均可以得到一个合理的近似表示。AR模型的正则方程和参数计算：正则方程，是在均方误差最小准则下建立的模型参数与相关函数的关系他是在二阶统计意义上的建模，由于高斯过程仅用二阶统计量就可以完全描述，所以这种基于相关函数的建模方法只适合于高斯过程。AR模型的参数可以借助解线性方程获得，而且MA和ARNA模型都可以用阶数足够高的AR模型表示，所以实际中通常选用AR模型。Yule-Walker方程：假定和都是实平稳随机信号，其中是方差为1的高斯白噪声序列，为AR过程，则可将系统输出表示为：即用n时刻之前p个值的线性组合来预测，预测误差为，在均方误差最小准则下，的选择应使均方值最小。据此准则有，经过推导得到AR模型的正则方程，即Y-W方程：直接解Y-W方程的计算量一般较大，这里使用L-D快速递推算法：其中，递推初始条件：，自相关函数求法： ，这样做是保证求取自相关函数只用到观测数据，而不对序列补零。从观测得到随机序列求得，在代入递推公式即可求得从而得到，调用freqz函数后得到，再由求得。Burg算法：Burg算法与L-W算法的区别在于它没有直接使用观测随机序列的自相关函数求取AR模型的参数，而实现估计反射系数。再利用Levinson关系式求AR模型参数。前向和后向预测误差递推公式为：反射系数估计准则是使前向预测误差功率。即得到 Burg算法同样不对未知数据做人为假设。该算法的优点是求得的AR模型保证稳定，另外，由于Burg算法不估算自相关函数，所以性能优于自相关法。尤其在短数据时具有较高分辨率。4、 实现程序clear all%信号产生，两个频率的正弦信号加白噪声t0=1; %采样间隔tt=128; %采样时长N=tt/t0; %采样点数-1M=27; %模型阶数t=0:t0:tt;f1=0.2; %信号频率f2=0.213;y=sin(2*pi*f1*t)+sin(2*pi*f2*t)+randn(1,N+1);plot(t,y)title('信号加高斯噪声')%经典谱估计F=fft(y,1024);F=(F.2)/1024;m=(0:1023)/(1024);figure;plot(m,10*log(abs(F);title('周期图法经典谱估计')xlabel('对采样频率归一化')ylabel('10lg(H)')w0=Freqval(abs(F),1/t0) %基于Y-W方程的参数谱估计%采用L-D快速递推算法实现R=zeros(1,N+1);for i=0:N %计算相关函数 for r=0:(N-i) R(i+1)=R(i+1)+y(r+1)*y(r+i+1); end R(i+1)=R(i+1)/(N+1);enda(1,1)=(-1)*R(2)/R(1);P(1)=R(1);P(2)=R(1)*(1-(a(1,1)2);p=1;while(abs(a(p,p)<1 && (p<=M)%此处p表示每次运算 p=p+1; %中最高阶数，反射系 X=0; %数k，而直接是a(p,p) for i=1:(p-1) X=X+a(p-1,i)*R(p-i+1); end a(p,p)=(-1)*(R(p+1)+X)/P(p); for i=1:p-1 a(p,i)=a(p-1,i)+a(p,p)*a(p-1,p-i); end P(p+1)=P(p)*(1-(a(p,p)2); h=p-1; %h表示模型最高阶，一般为Mendfor i=1:h A(i)=a(h,i);endG=P(p)(1/2);A=1 A;H=freqz(G,A,1024,'whole');figure;plot(m,20*log(abs(H);title('基于Y-W方程的参数谱估计')xlabel('对采样频率的归一化')ylabel('10lg(H)')w1=Freqval(abs(H),1/t0) %Burg算法实现信号参数谱估计p=1;fenzi=0;fenmu=0;ef=zeros(p,N+1);eb=zeros(p,N+1);ef(1,:)=y;eb(1,:)=y;for i=p:N fenzi=fenzi+ef(1,i+1)*eb(1,i); fenmu=fenmu+ef(1,i+1)2+eb(1,i)2;enda(1,1)=(-2)*fenzi/fenmu;while(a(p,p)<1 && p<=M) fenzi=0; fenmu=0; p=p+1; for i=1:N ef(p,i+1)=ef(p-1,i+1)+a(p-1,p-1)*eb(p-1,i); eb(p,i+1)=eb(p-1,i)+a(p-1,p-1)*ef(p-1,i+1); end for i=p:N fenzi=fenzi+ef(p,i+1)*eb(p,i); fenmu=fenmu+ef(p,i+1)2+eb(p,i)2; end a(p,p)=(-2)*fenzi/fenmu; endh=p-1;for i=1:h-1 a(h,i)=a(h-1,i)+a(h,h)*a(h-1,h-i); A(i)=a(h,i);endA(h)=a(h,h);A=1 A;G=R(1)(1/2);H=freqz(G,A,1024,'whole');figure;plot(m,20*log(abs(H);title('Burg算法');xlabel('对采样频率的归一化');ylabel('10lg(H)');w2=Freqval(abs(H),1/t0)%FreqVal函数用于测定算法检测到的信号频率，设计为监测两个function z=FreqVal(x,y)%x为函数序列，y为采样频率q=length(x)/2;w=zeros(1,2);z=zeros(1,2);for j=2:q-1 if x(j)>=x(j-1) && x(j)>=x(j+1) if x(j)>w(1) w(1)=x(j); z(1)=j; end endendfor j=2:q-1 if x(j)>=x(j-1) && x(j)>=x(j+1) if x(j)>w(2) if j=z(1) w(2)=x(j); z(2)=j; end end endendz=(z.*y)./(2*q);五、实验结果简析程序一次执行结果如下：W0=0.2002 0.2158w1 =0.1992 0.2139w2 =0.2139 0.1973从以上结果可以对三种方法做一比较，随即序列表达式为，从检测到两个正弦信号的频谱幅度来看周期图法结果约为10dB，而Y-W法两峰值在2530dB之间，明显高于前者，Burg法可以达到接近40dB。再从谱估计图两峰值图形比较，由于加窗效应，周期图法的宽度是最宽的，而Y-W法峰宽度比较窄，效果好于周期图法，相比之下，Burg法峰宽度最窄，效果在三者中最好。最后再看三种方法得到的精确估计值与真实值的比较，随机信号谐波频率真实值：，而w0 =0.2002 0.2158，w1 =0.1992 0.2139，w2 =0.2139 0.1973，Y-W法和Burg法对0.213的估计好于周期图法，对0.2的估计，周期图法好于Y-W法，而Burg法较差。当然，这只是一次实验的结果，要想对三种方法对随机信号频率精确估值效果进行研究，可以进行大量实验，再求取各自均值和方差等统计量进行研究。6 总结通过本次作业实验，对使用Y-W方法和和Burg算法对随机信号进行谱估计有了更深的，更清晰的认识，报告撰写过程中，需要不断阅读研究有关周期图、Y-W算法、L-D快速递推法、预测误差滤波器及Burg算法的基本原理和公式推导过程，加深了对这些理论原理和各种特点的认识。在实验程序编写过程中，通过对各种算法的具体实现，不仅掌握了随机信号现代谱估计Y-W法和Burg法的实现过程，也对matlab的使用更加熟悉，当然，实验程序还有很多可以优化的地方，但总的来说是一个进步。