编程语言
首页 > 编程语言> > 现代法谱估计(1)Yule Walker 方程法MATLAB及Python实现

现代法谱估计(1)Yule Walker 方程法MATLAB及Python实现

作者:互联网

原理

AR模型的系统函数可以表示为:
在这里插入图片描述
如果在白噪声 激励下模型的输出为x(n),则模型输入、输出关系的时域表达式为:
在这里插入图片描述
此式为AR模型的差分方程。将白噪声 激励AR模型产生的输出x(n)叫做AR过程。
根据相关卷积定理,若y(n)=x(n)*h(n),则有
在这里插入图片描述
即卷积的相关等于相关的卷积。如果对上式两边求傅里叶变换,根据维纳辛钦定理和相关定理,有
在这里插入图片描述
即输出自功率谱等于输入自功率谱与系统能量谱的乘积。
根据谱分解定理,任何平稳随机信号x(n)都可以看成是由高斯白噪声激励一个因果稳定的可逆系统H(z)产生的输出,如下图所示
在这里插入图片描述
因此上述公式可以变换为
在这里插入图片描述
上式说明,平稳随机信号x(n)的功率谱可以用模型H(z)的参数来表示。
对于p阶AR模型的系统函数
在这里插入图片描述
可以看出其由p+1个待定参数:a1-ap和G。
系统输出x(n)可表示为:
在这里插入图片描述
上式可以解释为:用n时刻之前的p个值的线性组合来预测n时刻的值x(n),预测误差为Gw(n) 。在均方误差最小准则下,组合系数的选择应该使预测误差的均方误差最小。令均方误差为e(n),有
在这里插入图片描述
则有
在这里插入图片描述
由于
在这里插入图片描述
可以计算出
在这里插入图片描述
在最小均方误差(MMSE)准则下,要使预测值最佳地逼近x(n),参数的选择应使
在这里插入图片描述
也即
在这里插入图片描述
可得到
在这里插入图片描述

在这里插入图片描述
由上式可得p个方程,写成矩阵式为
在这里插入图片描述
上式用到自相关函数的偶对称性质。由这p个方程,可以求出p个参数ai。有了这些参数就可以根据自相关函数和参数ai求解系统增益G。联立可得
在这里插入图片描述
上式即为AR模型的正则方程,也叫Yule-Walker方程。
利用以上求出的p+1个参数,我们可以可以估计出该随机信号的功率谱。

程序及结果

(出于维护版权原因,此次只放截图)
MATLAB
程序:
在这里插入图片描述
在这里插入图片描述
结果:
在这里插入图片描述
Python
程序:
在这里插入图片描述
在这里插入图片描述
在这里插入图片描述
在这里插入图片描述
在这里插入图片描述
在这里插入图片描述
结果:
在这里插入图片描述

分析

由上图可见,我给程序输入的N为256,取的阶数p为128,信号中f1=0.1,f2=0.13,在图中我们可以看到对应的位置出现了峰值,由于给定幅度大小不同,故峰值大小有很大差别。与经典法谱估计不同的是,图中只出现了一半,而不是像经典法中的出现四个两两对称的峰,这是由于程序中使用了freqz函数(Python中为scipy.signal.freqz)求其频率响应,MATLAB中我们可以打这个函数的帮助手册:

h :Complex, n-element frequency response vector.
h:复n元频率响应矢量。
w: Frequency vector of length n, in radians/sample. w consists of n points equally spaced around the upper half of the unit circle (from 0 to π radians/sample).
w:长度n的频率矢量,单位为弧度/样本。w由n个点组成,这些点平均分布在单位圆的上半部分(从0到π弧度/样本)。

可见w的范围是0到π弧度,所以用freqz函数求出的频率响应只有0到π部分。由于h是复频率响应矢量,因此要求其功率谱要对频率响应求模平方。然后我利用了寻峰函数来得出峰值对应的横坐标(Python中是自己写了一个寻峰估计频率的方法),程序中将所有的峰值及其索引分别放入到两个数组中,且根据峰值大小降序排列,确定最大的两个峰值和其索引位置,然后由索引位置除以2N(freqz函数得出的是单位圆上半部分,频率范围对应为0-0.5,而对应的序列长度为N,因此根据比例关系可知,索引值/N=频率值/0.5也即,频率值=索引值/2N)。

(因原博客是以word形式写的,CSDN不支持Mathtype公式,故部分地方直接采用了截图形式)

标签:误差,Yule,函数,Python,模型,谱估计,峰值,AR,频率响应
来源: https://blog.csdn.net/qq_44628230/article/details/107004723