《随机信号处理》word版

《《随机信号处理》word版》由会员上传分享，免费在线阅读，更多相关内容在应用文档-天天文库。

1、参数功率谱估计的实现及其与经典谱估计的比较参数功率谱估计的实现及其与经典谱估计的比较一、摘要频域分析是从频率角度对信号进行分析研究，对于确定性信号来说，通常使用傅立叶变换或者傅立叶技术展开得到频域表达，但对随机信号而言，由于其时域波形的随机性及能量无限，没有确定的时域表达式，无法用傅立叶变换直接将其变换到频域中去研究。表达随机信号通常使用概率密度函数，根据维纳-辛钦定理，广义随机过程的功率谱与自相关函数是一对傅立叶变换，所以自然想到可以用功率谱来研究随机过程的频域性质。对随机信号的功率谱估计方法通常分为两大类：经典谱估计和现代谱估计经典谱估计是基于维纳-辛钦定理，

2、从自相关函数出发通过傅立叶变换得到功率谱，而现代谱估计则将随机信号看成白噪声通过一个滤波器的输出。现代谱估计就是通过记录的信号序列估计滤波器参数，从而得到其频率响应，最后通过得到其功率谱。从频率分辨率来看，经典谱估计效果一般不如参数谱估计好，而且参数谱估计的不同算法对不同的采样序列也有不同的效果，本文将对周期图法，基于L-D快速递推算法的Y-W法和Burg算法使用Matlab进行编程实现并做比较。二、关键字功率谱估计周期图法L-D递推Y-W算法Burg算法三、原理1.经典谱估计广义平稳随机信号经典谱估计基于维纳-辛钦定理：广义平稳随机的自相关函数与其功率谱是一堆傅

3、立叶变换。所以要求功率谱，只需由随机序列求出自相关函数然后进行傅立叶变换即可，m=0，1，2，……N-1周期图法：虽然经典谱估计方法比较直观简单，但由于随机序列相当于对信号加床，所以求自相关函数后傅立叶变换的功率谱往往受到窗函数影响，不是信号真实谱，所以就产生了以下的现代谱估计。参数功率谱估计：根据谱分解定理，任何平稳随机信号都可以看成是由白噪声经过一个因果稳定可逆系统产生的输出，如上图。其系统函数全部极点都在z平面单位圆内部，将表示成，实际中，很多物理过程都可以近似表示成AR，ARMA，MA过程，其中AR过程为，也就是可以看成白噪声通过系统后的输出。如果根据有限

4、的数据记录（看成系统输出）和已知输入（白噪声）确定出系统的参数，那么所要研究的随机过程就是白噪声通过已知系统之后的输出，于是可以用系统的参数来描述索要研究的物理过程，例如要分析实平稳随机信号的有理功率谱，根据平稳随机信号通过线性系统的性质可知：，其中是白噪声的功率谱（为常数），是的能量谱对线性时不变的稳定系统，单位脉冲响应是确定性能量信号，其能量谱是频谱幅度的平方。上式说明平稳随机信号的功率谱可以用模型的参数来表示。所以只要建立合适的信号模型，并估计出模型的参数，就可以分析随机信号的功率谱。上面仅提到了AR模型，是因为柯尔莫哥洛夫定理表明，任何ARMA或MA过程都

5、可以用无限阶的AR过程表示，即使在建模式选择了不正确的模型，只要模型阶数足够高，均可以得到一个合理的近似表示。AR模型的正则方程和参数计算：正则方程，是在均方误差最小准则下建立的模型参数与相关函数的关系他是在二阶统计意义上的建模，由于高斯过程仅用二阶统计量就可以完全描述，所以这种基于相关函数的建模方法只适合于高斯过程。AR模型的参数可以借助解线性方程获得，而且MA和ARNA模型都可以用阶数足够高的AR模型表示，所以实际中通常选用AR模型。Yule-Walker方程：假定和都是实平稳随机信号，其中是方差为1的高斯白噪声序列，为AR过程，则可将系统输出表示为：即用n时

6、刻之前p个值的线性组合来预测，预测误差为，在均方误差最小准则下，的选择应使均方值最小。据此准则有，经过推导得到AR模型的正则方程，即Y-W方程：直接解Y-W方程的计算量一般较大，这里使用L-D快速递推算法：其中，递推初始条件：，，自相关函数求法：，，这样做是保证求取自相关函数只用到观测数据，而不对序列补零。从观测得到随机序列求得，在代入递推公式即可求得从而得到，调用freqz函数后得到，再由求得。Burg算法：Burg算法与L-W算法的区别在于它没有直接使用观测随机序列的自相关函数求取AR模型的参数，而实现估计反射系数。再利用Levinson关系式求AR模型参数。

7、前向和后向预测误差递推公式为：反射系数估计准则是使前向预测误差功率。即得到Burg算法同样不对未知数据做人为假设。该算法的优点是求得的AR模型保证稳定，另外，由于Burg算法不估算自相关函数，所以性能优于自相关法。尤其在短数据时具有较高分辨率。四、实现程序clearall%信号产生，两个频率的正弦信号加白噪声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,

8、y)tit