功率谱估计的参数模型方法---AR模型谱估计：自相关法&协方差法

信号的线性模型（零极点模型/Pole-Zero Modeling）
- （1）AR模型（自回归/全零点模型）
- （2）MA模型（移动平均/全极点模型）
- （3）ARMA模型（自回归/全零点模型）
自相关法（Yule-Walk 法）
协方差法
Reference

第一篇自己的博客，这篇文章将参考文献[1]中的AR模型方法与参考文献[2]中的方法联系起来，算是个笔记吧。看到这个代码有点懵acm.m，这样简单的问题都思考了好一会儿了，LaTex也不怎么会用。。。我的知识真是太匮乏了

信号的线性模型（零极点模型/Pole-Zero Modeling）

常系数线性差分方程：
Y[n]+a1Y[n−1]+⋯+apY[n−p]=b0X[n]+b1X[n−1]+⋯+bqX[n−q]Y[n] + a_1Y[n-1]+\cdots+a_pY[n-p]=b_0X[n]+b_1X[n-1]+\cdots+b_qX[n-q]Y[n]+a1Y[n−1]+⋯+apY[n−p]=b0X[n]+b1X[n−1]+⋯+bqX[n−q]
其中，X[n]X[n]X[n]是输入随机序列，Y[n]Y[n]Y[n]是输出序列。

（1）AR模型（自回归/全零点模型）

如果则如果q=0q=0q=0（右端仅有一项）
Y[n]+a1Y[n−1]+⋯+apY[n−p]=b0X[n]Y[n] + a_1Y[n-1]+\cdots+a_pY[n-p]=b_0X[n]Y[n]+a1Y[n−1]+⋯+apY[n−p]=b0X[n]
记为AR(p)AR(p)AR(p)，称Y[n]Y[n]Y[n]为ppp阶自回归模型

（2）MA模型（移动平均/全极点模型）

如果则如果p=0p=0p=0（左端仅有一项）
Y[n]=b0X[n]+b1X[n−1]+⋯+bqX[n−q]Y[n] =b_0X[n]+b_1X[n-1]+\cdots+b_qX[n-q]Y[n]=b0X[n]+b1X[n−1]+⋯+bqX[n−q]
记为MA(q)MA(q)MA(q)，称Y[n]Y[n]Y[n]为qqq阶移动平均模型

（3）ARMA模型（自回归/全零点模型）

如果则如果p,q0p,q0p,q0记为ARMA(p，q)ARMA(p，q)ARMA(p，q)，称Y[n]Y[n]Y[n]为(p,q)(p,q)(p,q)阶自回归移动平均模型

常系数线性差分方程对应的系统函数H(z)H(z)H(z)为
H(z)=b0+b1z−1+⋯+bqz−q1+a1z−1+⋯+apz−p)H(z)=\frac{b_0+b_1z^{-1}+\cdots+b_qz^{-q}}{1+ a_1z^{-1}+\cdots+a_pz^{-p}})H(z)=1+a1z−1+⋯+apz−pb0+b1z−1+⋯+bqz−q)
将随机信号建模为零均值单位方差白噪声激励H(z)H(z)H(z)产生：假设输入随机序列X[n]X[n]X[n]是零均值单位方差白噪声，功率谱为常数111，自相关函数为单位脉冲函数δ(k)\delta(k)δ(k),则输出功率谱为
P(ejω)=∣b0+b1e−jω+⋯+bqe−jqω∣2∣1+a1e−jω+⋯+ape−jpω∣2P(e^{j\omega})=\frac{|b_0+b_1e^{-j\omega}+\cdots+b_qe^{-jq\omega}|^2}{|1+ a_1e^{-j\omega}+\cdots+a_pe^{-jp\omega}|^2}P(ejω)=∣1+a1e−jω+⋯+ape−jpω∣2∣b0+b1e−jω+⋯+bqe−jqω∣2
一旦模型选定，下一步是由给定数据估计参数模型ai,bja_i,b_jai,bj,对于AR模型需要估计b0,a1,⋯,apb_0,a_1,\cdots,a_pb0,a1,⋯,ap。基于最小均方误差准则的四种参数估计方法：自相关法（Yule-Walk 法）、协方差法、修正协方差法、Burg法，后两种方法最小化前向加反向预测误差的平方和来求解参数。
argmaxakξMSEarg\;\underset{a_k}{max}\xi_{MSE}argakmaxξMSE
其中ξMSE=E{∣e(n)∣2}\xi_{MSE}=E\{|{e(n)|}^2\}ξMSE=E{∣e(n)∣2},e(n)=x(n)−x^(n)\quad e(n)=x(n)-\widehat x(n)e(n)=x(n)−x(n)

自相关法（Yule-Walk 法）

为求解aka_kak,取ξMSE\xi_{MSE}ξMSE对aka_kak求偏导等于000,化简得
∑l=1palrx(k−l)=−rx(k)k=1,2,⋯,p\sum_{l=1}^{p}a_lr_x(k-l)=-r_x(k)\quad k=1,2,\cdots,p∑l=1palrx(k−l)=−rx(k)k=1,2,⋯,p
写成矩阵形式为：
Rxa=−rx\boldsymbol{R_xa=-r_x}Rxa=−rx\quad参考文献[1]中的式(3.7)
此式被称为线性预测的正则方程（Normal equations），其中Rx\boldsymbol{R_x}Rx是数据x(n)x(n)x(n)的自相关矩阵，是N∗NN*NN∗N的Hermite、Toeplitz矩阵：
Rx=[rx(0)rx∗(1)⋯rx∗(p−1)rx(1)rx(0)⋯rx∗(p−2)⋯⋯⋯⋯rx(p−1)rx(p−2)⋯rx(0)]\boldsymbol{R_x}= \begin{bmatrix} r_x(0)&r_x^*(1)&\cdots&r_x^*(p-1)\\ r_x(1)&r_x(0)&\cdots&r_x^*(p-2)\\ \cdots&\cdots&\cdots&\cdots\\ r_x(p-1)&r_x(p-2)&\cdots&r_x(0)\\ \end{bmatrix}Rx=⎣⎢⎢⎡rx(0)rx(1)⋯rx(p−1)rx∗(1)rx(0)⋯rx(p−2)⋯⋯⋯⋯rx∗(p−1)rx∗(p−2)⋯rx(0)⎦⎥⎥⎤
其中：
rx(k)=1N∑n=kNx(n)x∗(n−k)r_x(k) = \frac{1}{N}\sum_{n=k}^{N}x(n)x^*(n-k)rx(k)=N1∑n=kNx(n)x∗(n−k)NNN为数据长度，whenn<0orn>N,x(n)=0when\quad n<0\quad or\quad n>N,x(n)=0whenn<0orn>N,x(n)=0
a=[a1,a2,⋯,ap]Ta={[a_1,a_2,\cdots,a_p]}^Ta=[a1,a2,⋯,ap]T
式中Rx\boldsymbol{R_x}Rx又可以写成
Rx=XHX\boldsymbol{R_x}=\boldsymbol X^H\boldsymbol XRx=XHX
X=[x(0)00⋯0x(1)x(1)0⋯0x(2)x(1)x(0)⋯0⋯⋯⋯⋯⋯x(N)x(N−1)x(N−2)⋯x(N−p+1)0x(N)x(N−1)⋯x(N−p+2)⋯⋯⋯⋯⋯000⋯x(N)]\boldsymbol X = \begin{bmatrix} x(0)&0&0&\cdots&0\\ x(1)&x(1)&0&\cdots&0\\ x(2)&x(1)&x(0)&\cdots&0\\ \cdots&\cdots&\cdots&\cdots&\cdots\\ x(N) &x(N-1) &x(N-2)&\cdots&x(N-p+1)\\ 0&x(N)&x(N-1)&\cdots&x(N-p+2)\\ \cdots&\cdots&\cdots&\cdots&\cdots\\ 0&0&0&\cdots&x(N)\\ \end{bmatrix}X=⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎡x(0)x(1)x(2)⋯x(N)0⋯00x(1)x(1)⋯x(N−1)x(N)⋯000x(0)⋯x(N−2)x(N−1)⋯0⋯⋯⋯⋯⋯⋯⋯⋯000⋯x(N−p+1)x(N−p+2)⋯x(N)⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎤
rx\boldsymbol{r_x}rx可以写成Rx=XHX2\boldsymbol{R_x}=\boldsymbol X^H\footnotesize{\boldsymbol X_2}Rx=XHX2，X2\footnotesize{\boldsymbol X_2}X2是X\boldsymbol{X}X的第二列。
因此正则方程可写为：
XHXa=−XHX2\boldsymbol X^H\boldsymbol X\boldsymbol a=-\boldsymbol X^H\footnotesize{\boldsymbol X_2}XHXa=−XHX2
即：
Xa=−X2\boldsymbol X\boldsymbol a=-\footnotesize{\boldsymbol X_2}Xa=−X2\quad参考文献[2]中的式(4.21)(4.81)(4.123)
利用自相关法估计参数的代码acm.m在这里

协方差法

Reference

1.杨绿溪. 现代数字信号处理[M]. 2007.
2.Hayes M H. Statistical Digital Signal Processing and Modeling[J]. 1996. Chapter 4, 式(4.40)\quad书中MATLAB代码在这里
3.陈明. 信息与通信工程中的随机过程(第2版)[M]. 2005.

功率谱估计的参数模型方法---AR模型谱估计：自相关法协方差法相关推荐

现代信号处理——AR模型谱估计
AR谱估计方法可归结为求解AR模型系数或线性预测器系数的问题. AR模型参数估计方法:信号预测误差最小原则(或预测误差功率最小) 自相关法(Levison递推法) Burg法协方差法修正协方差法( ...
java中arSigal_基于AR模型谱估计算法（Yule-Walker方法与Burg方法）的C++实现
头文件: /* * Copyright (c) 2008-2011 Zhang Ming (M. Zhang), zmjerry@163.com * * This program is free so ...
【时间序列模型】AR模型（原理剖析+MATLAB代码）
目录一.AR模型的原理剖析 1.1 AR模型原理 1.2 模型的参数估计 1.3 模型的定阶方法 (1)AIC (2)BIC 二.MATLAB代码三.AR相关论文参考文献: 前言时间序列分析方 ...
burg算法matlab求功率谱,基于AR模型的Burg算法功率谱估计.doc
基于AR模型的Burg算法功率谱估计三种功率谱估计方法性能研究 1.前言: 我们已经知道一个随机信号本身的傅里叶变换并不存在,因此无法像确定性信号一样用数字表达式来精确表达它,而只能用各种统计平均量 ...
matlab ar谱分析,AR模型功率谱估计的典型算法比较及MATLAB实现
CHINANEW TELECOMMUNICATIONS 功率谱估计是信息学科中的研究热点, 在过去的 30 多年里取得了飞速的发展.现代谱估计主要是针对经典谱估计( 周期图和自相关法) 的分辨率低和方 ...
ar模型功率谱估计 matlab,AR模型功率谱估计及Matlab实现
南昌大学实验报告学生姓名:学号:专业班级: 实验类型:□验证□综合□设计□创新实验日期:实验成绩: 一.实验名称基于AR模型的功率谱估计及Matlab实现二.实验目的 1.了解现代谱估计方法,深 ...
基于AR模型的功率谱matlab,基于AR模型法功率谱估计的Matlab实现
2005iF 12月笫 21卷第 6期武警工程学院学报 jOURNAL OF ENGG COLLEGE OF ARMED POLICE FORCE DeC．2005 Vo1．21 No．6 [计 ...
QuantitativeFinance：量化金融之金融时间序列分析之ES/ETS/GARCH模型的简介、Box-Jenkins方法-AR/MA/ARMA/ARIMA模型的简介及其建模四大步骤之详细攻略
QuantitativeFinance:量化金融之金融时间序列分析之ES/ETS/GARCH模型的简介.Box-Jenkins方法-AR/MA/ARMA/ARIMA模型的简介及其建模四大步骤之详细攻略 ...
AR模型在信号处理中的应用
本文目标:分析AR模型并求解AR模型的输出x(n)的功率谱. 1. AR模型概念观数字信号处理功率谱估计方法分经典功率谱估计和现代功率谱估计,现代功率谱估计以参数模型功率谱估计为代表,参数功率谱模型 ...

功率谱估计的参数模型方法---AR模型谱估计：自相关法协方差法