第12章_1_参数模型功率谱估计
参数模型功率谱估计
LSI系统的输入、输出关系：
差分方程 卷积关系 以上两式是LSI系统的时域表示，无论对确 定性信号还是随机信号都成立。现假定输入、 输出是平稳随机信号（输入是白噪声）。
转移函数的两 种表示形式， 独立于信号。
谱分解 的Z域 表示
待辨识 的参数。
Px (z)
u2H (z)H (z1)
u2
B( z ) B( z 1 ) A(z) A(z1)
求解方法：由下面的差分方程入手：
两边同乘 x(n m) ，求均值
p
ak Ex(n m k)x(n) k 1
Eu(n m)x(n)
x(n) 和
u(n) 的
互相关
卷积 关系
因果 系统
结果1： 结果2：
结合 起来
正则方程 （Normal Eq.）
rx (0) rx (1) rx (1) rx (0)
k
可以得到使 最小的 1,L , p 及 min 。
不求导，使用正交原理：
E{x(n m)[x(n) xˆ(n)]} 0, m 1, 2, L , p
e(n)
p
rx (m) k rx (m k), m 1, 2, L , p k 1
Wiener-Hopf Eq.
min E{x(n)[x(n) xˆ(n)]}
u(n)
x(n)
AR模型
1 A(z)
白化滤 x(n)
波器
A(z)
e(n)
x(n)
线性预 测器
1 A(z)
xˆ(n)
e(n)
Yule-Walker 方程的快速计算
－Levinson-Durbin快速算法： 要求解的参数：
ap (1), ap (2),L , ap ( p), 2(min @p )
功率谱估计教材
1 ˆxx (m) r N
ˆ ( w) P BT
N |m|1
n 0
x ( n ) x ( n m)
m m
ˆxx (m) e jwm r
自相关法
由于在估计x的自相关函数时，数据的长度为N， 因此估计的自相关函数r̂xx(m)的长度为2N-1点：
ˆxx (m) r ˆxx (m) r 0
功率谱估计
--非参数估计方法
功率谱估计
经典功率谱估计（非参数法）
自相关法 周期图法
参数谱估计（参数法）
AR、MA、ARMA模型
经典谱估计法-自相关法
自相关法-BT(Blackman-Tukey提出)
随机信号的一个样本数据为[x(0),x(1),…,x(N-1)]，长 度为N。 先根据样本数据估计自相关函数r̂xx(m)，再利用FFT 变换，得到功率谱的估计PBT(w)。
m
jwm ˆ E[rxx (m)]e
窗函数法
则自相关函数的变化：
1 ˆxx (m)] E[r N
n
E[ x(n)v(n) x(n m)v(n m)]
1 E[ x(n) x(n m)] v(n)v(n m) N n
1 rxx (m) N
这样，功率谱估计为：
m N 1 m N 1
| m | N 1 else
ˆ ( w) P BT
jwm ˆ rxx (m) e
周期图法
相关法是利用样本数据对自相关函数进行估计， 进而估计功率谱密度，而周期图法则根据功率 谱密度的另一定义：
1 N 1 Pxx (w) lim E[ | x(n)e jwn |2 ] N N n 0
功率谱计算[解说]
![功率谱计算[解说]](https://img.taocdn.com/s3/m/3cec2a18f02d2af90242a8956bec0975f465a478.png)
功率谱计算功率谱估计在现代信号处理中是一个很重要的课题，涉及的问题很多。
在这里，结合matlab，我做一个粗略介绍。
功率谱估计可以分为经典谱估计方法与现代谱估计方法。
经典谱估计中最简单的就是周期图法，又分为直接法与间接法。
直接法先取N点数据的傅里叶变换（即频谱），然后取频谱与其共轭的乘积，就得到功率谱的估计；间接法先计算N点样本数据的自相关函数，然后取自相关函数的傅里叶变换，即得到功率谱的估计.都可以编程实现，很简单。
在matlab中，周期图法可以用函数periodogram实现。
但是周期图法估计出的功率谱不够精细，分辨率比较低。
因此需要对周期图法进行修正，可以将信号序列x(n)分为n个不相重叠的小段，分别用周期图法进行谱估计，然后将这n段数据估计的结果的平均值作为整段数据功率谱估计的结果。
还可以将信号序列x(n)重叠分段，分别计算功率谱，再计算平均值作为整段数据的功率谱估计。
这2种称为分段平均周期图法，一般后者比前者效果好。
加窗平均周期图法是对分段平均周期图法的改进，即在数据分段后，对每段数据加一个非矩形窗进行预处理，然后在按分段平均周期图法估计功率谱。
相对于分段平均周期图法，加窗平均周期图法可以减小频率泄漏，增加频峰的宽度。
welch法就是利用改进的平均周期图法估计估计随机信号的功率谱，它采用信号分段重叠，加窗，FFT等技术来计算功率谱。
与周期图法比较，welch法可以改善估计谱曲线的光滑性，大大提高谱估计的分辨率。
matlab中，welch法用函数psd实现。
调用格式如下：[Pxx,F] = PSD(X,NFFT,Fs,WINDOW,NOVERLAP)X：输入样本数据NFFT：FFT点数Fs:采样率WINDOW：窗类型NOVERLAP，重叠长度现代谱估计主要针对经典谱估计分辨率低和方差性不好提出的，可以极大的提高估计的分辨率和平滑性。
可以分为参数模型谱估计和非参数模型谱估计。
参数模型谱估计有AR模型，MA模型，ARMA模型等；非参数模型谱估计有最小方差法和MUSIC法等。
功率谱估计
功率谱估计引言：对信号和系统进行的分析研究、处理有两类方法：一类是在时域内进行，维纳滤波、卡尔曼滤波以及自适应滤波等都属于时域处理方法；另一类方法是频域研究方法。
对于确定性信号，傅里叶变换是在频率分析研究的理论基础，但是在实际生活中大多数信号是随机信号，而随机信号的傅里叶变换是不存在的，在实际应用中，通常通过采集和观测平稳随机过程的一个抽样序列的一段（有限个）数据，根据这有限个已知的数据来估计随机过程的功率谱问题来对随机信号进行分析，这即是频率谱估计。
功率谱估计是数字信号处理的主要内容之一，主要研究信号在频域中的各种特征，目的是根据有限数据在频域内通过用某种有效的方法来估计出其功率谱密度，从而得出信号、噪声及干扰的一些性质来，提取被淹没在噪声中的有用信号。
功率谱估计就是通过信号的相关性估计出接受到信号的功率随频率的变化关系，实际用途有滤波，信号识别（分析出信号的频率），信号分离，系统辨识等。
谱估计技术是现代信号处理的一个重要部分，还包括空间谱估计，高阶谱估计等。
按照Weiner —Khintchine 定理，随机信号的功率谱和其自相关函数服从傅里叶变换关系，可以得出功率谱的一个定义，如公式（1）所示：()jwm m xx jw xx e m re P -∞-∞=∑=)( 公式（1）对于平稳随机信号，服从各态历经性，集合平均可以用时间平均来代替，可以推出功率谱的另一定义。
如公式（2）所示：()])(121[2lim ∑-=-∞→+=N N n jwn N jw xx e n x N E e P 公式（2）频率谱估计主要分为经典谱估计和现代谱估计，经典谱估计是将数据工作区外的未知数据假设为零,相当于数据加窗,主要方法有相关法和周期图法;现代谱估计是通过观测数据估计参数模型再按照求参数模型输出功率的方法估计信号功率谱,主要是针对经典谱估计的分辨率低和方差性能不好等问题提出的,应用最广的是AR 参数模型。
参数法功率谱估计
参数法功率谱估计一、信号的产生（一）信号组成在本实验中，需要事先产生待估计的信号，为了使实验结果较为明显，我产生了由两个不同频率的正弦信号（频率差相对较大）和加性高斯白噪声组成的信号。
（二）程序N=1024;n=0:N-1;xn=2*cos(2*pi*0.2*n)+ cos(2*pi*0.213*n)+randn(1,1024);这样就产生了加有白噪声的两个正弦信号其波形如下0100200300400500600-8-6-4-2246810(a) 两个正弦信号与白噪声叠加的时域波形二、参数模型法功率谱估计（一）算法原理简介1．参数模型法是现代谱估计的主要内容，思路如下：① 假定所研究的过程)(n x 是由一个白噪声序列)(n 激励一个因果稳定的可逆线性系统)(z H 的输出；② 由已知的)(n x ，或其自相关函数)(m r x 估计)(z H 的参数；③ 由)(z H 的参数来估计)(n x 的功率谱。
2．自回归模型，简称AR 模型，它是一个全极点的模型。
“自回归”的含义是：该模型现在的输出是现在的输入和过去p 个输出的加权和。
此模型可以表现为以下三式：① ∑=+--=p k k n u k n x a n x 1)()()(；② ∑=-+==p k kk z a z A z H 111)(1)(；③ 2121)(∑=-+=p k jwkk jw x e a e P σ。
3．AR 模型的正则方程建立了参数k a 和)(n x 的自相关函数的关系，公式如下：=)(m r x ∑=--p k x k k m r a 1)( 1≥m 时，=)(m r x 21)(σ+-∑=k r a pk x k 0=m 时。
（二）两种AR 模型阶次的算法1．Yule-Walker 算法(自相关法)（1）算法主要思想Yule-Walker 算法通过解Yule-Walker 方程获得AR 模型参数。
从低阶开始递推，直到阶次p ，给出了在每一个阶次时的所有参数。
功率谱估计浅谈讲解
功率谱估计浅谈摘要：介绍了几种常用的经典功率谱估计与现代功率谱估计的方法原理，并利用Matlab对随机信号进行功率谱估计，对两种方法做出比较，分别给出其优缺点。
关键词：功率谱；功率谱估计；经典功率谱估计；现代功率谱估计前言功率谱估计是从频率分析随机信号的一种方法,一般分成两大类:一类是经典谱估计;另一类是现代谱估计。
由于经典谱估计中将数据工作区以外的未知数据假设为零，这相当于数据加窗，导致分辨率降低和谱估计不稳定。
现代谱估计则不再简单地将观察区外的未知数据假设为零，而是先将信号的观测数据估计模型参数，按照求模型输出功率的方法估计信号功率谱，回避了数据观测区以外的数据假设问题。
周期图、自相关法及其改进方法(Welch)为经典(非参数)谱估计方法, 其以相关和傅里叶变换为基础,对于长数据记录较适用,但无法根本解决频率分辨率低和谱估计稳定性的问题,特别是在数据记录很短的情况下,这一问题尤其突出。
以随机过程的参数模型为基础的现代参数法功率谱估计具有更高的频率分辨率和更好的适应性,可实现信号检测或信噪分离,对语音、声纳雷达、电磁波及地震波等信号处理具有重要意义,并广泛应用于通信、自动控制、地球物理等领域。
在现代参数法功率谱估计方法中,比较有效且实用的是AR模型法，Burg谱估计法,现代谱估计避免了计算相关,对短数据具有更强的适应性,从而弥补了经典谱估计法的不足,但其也有一些自身的缺陷。
下面就给出这两类谱估计的简单原理介绍与方法实现。
经典谱估计法经典法是基于传统的傅里叶变换。
本文主要介绍一种方法：周期图法。
周期图法由于对信号做功率谱估计，需要用计算机实现，如果是连续信号，则需要变换为离散信号。
下面讨论离散随机信号序列的功率谱问题。
连续时间随机信号的功率谱密度与自相关函数是一对傅里叶变换对，即：()()j x x S R e d +∞-Ω-∞Ω=⎰τττ若()x R m 是()x R Ω的抽样序列，由序列的傅里叶变化的关系，可得()()j j n x x m S e R m e ωω∞-=-∞=∑即()j x S e ω与()x R m 也是一对傅里叶变换对。
第14章_参数模型功率谱估计_2
ˆ (n) e(n) x(n) x
对同样一组数据，我们可以实现双向预测：
x(n p) x(n p 1)
x(n 1)
x ( n)
Forward Prediction
前向预测 误差序列
误差功率
Backward Prediction
对同一组数据 的后向预测
后向预测 误差序列
前、后向预测误差序列有如下的关系：
m 1, 2,
,p
初始条件 反射 系数
上述关系引出了线性预测中的Lattice结构。这 一结构在现代谱估计、语音信号处理中有着 重要的应用。
上述的关系还是集总平均。对实际的信号： 单个样本有限长，求均值要简化，对
取代
的范围
N点数据，前向预测误差序列范围
Burg算法 Marple算法
14.7 MA模型
b(0) 1
再推导一步，有：
非线性方程组
MA模型的正则方程
从谱估计的角度，MA模型等效于经典法中 的间接法，所以分辨率低。因此，MA模型 用于谱估计无优势。但，MA模型： 1. 常用于系统辨识；
2. ARMA模型中包含了MA部分。
求解算法：由于MA模型的正则方程是非线性方 程，所以人们提出了很多的求解算法，如谱分解、 基于迭代的方法、基于高阶AR模型近似的方法。 后者最好用，基础是Wold分解定理。 对 x(n) 建立 一个无穷阶 的AR模型 令其等效为 模型 于是有：
i 1
H i
Wp w I
Rp S p Wp
秩为p 1
相关矩阵的 分解：信号 部分和噪声 部分
Sp
VV
i 1 i i
p 1
H
i
功率谱估计
W(n)为零均值方差为1的AWGN，n=1,2,3……,128
1.1周期图法：
我们知道随机信号的功率谱和自相关函数是一对傅式变换对：
而自相关函数定义为：
对于平稳随机过程，并由功率谱的偶函数特性得：
实际得到的随机信号只能是它的一个样本的片断，因此只能用有限长的样本序列来估计功率谱，这相当于用一个有限宽度（N）的窗函数 去乘样本序列，于是有（用离散频率K代替ω）：
title('周期图法');
xlabel('Hz');
ylabel('dB/Hz');
window1=hamming(128);
noverlap=20; %数据20%的重叠
[Pxx1,f]=pwelch(xn,window1,noverlap,nfft,Fs,'onesided');
plot_Pxx1=10*log10(Pxx1);
仿真结果：
2.现代功率谱估计
现代功率谱估计即参数谱估计方法是通过观测数据估计参数模型再按照求参数模型输出功率的方法估计信号功率谱。主要是针对经典谱估计的分辨率低和方差性能不好等问题提出的。主要方法有最大熵谱分析法（AR模型法）、Pisarenko谐波分解法、Prony提取计点法、Prony谱分解法以及Carpon最大似然法。其中AR模型应用较多，具有代表性。常用的模型有ARMA模型、AR模型、MA模型。
这就是用样本序列片断的DFT来估计功率谱的式子。由于加了矩形窗，使得这种直接的周期图估计平滑性、一致性和分辨率不能满足实际要求，因此有必要对上式作一些修改，这些修改主要有两种方法：
1.分段平均：即将长度为N的数据分成L段（允许有重叠），分别求出每一段的功率谱，然后即以平均。这样L个平均的方插笔每个随机变量的单独方差小L倍。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
零阶预 测器的 误差等 于信号 的功率
0 rx (0) rx (0) Ex(n)x(n 0)
1
rx (0) 2
Px
(e
j
)d
递
km
m1 k 1
am1(k
)rx
(m
k)
rx
(m)
m1
推
公
am (k) am1(k) kmam1(m k)
式
k 1, 2, , m 1
m m1[1 km2 ]
现在希望用它们预测 x(n)
x(n p) x(n p 1)
x(n 1) x(n)
p
xˆ(n) k x(n k) k 1
线性预测
e(n) x(n) xˆ(n)
误差序列
E e2(n)
E
x(n)
p
k x(n
k
2 )
k 1
均方误差
令：
0, k 1, 2, , p
k
可以得到使 最小的 1, , p 及 min 。
线性预测器的误差序列等效于激励AR模型 的白噪声序列；
由LP的含意，因此AR模型也可以看作是在 最小平方意义上对数的拟合；
由于LP包含了对数据的外推，因此，对应的 谱估计所用数据的范围比实际的应有扩展， 因此可以提高分辨率。
u(n)
x(n)
AR模型
1 A(z)
白化滤 x(n)
波器
A(z)
e(n)
机变量
连续型随机变量
Burg最大熵谱估计的思路是： 已知某随机信号自相关函数 rx (m) 的 p 1 个值
rx (0), rx (1), , rx ( p) ，现希望以这 p 1 个值对
m p 的自相关函数予以外推。外推的方法很多， Burg的准则是：外推后的自相关函数对应的时 间序列具有最大的熵，即是最随机的。
总效果：
紧随 的峰值
紧跟 谱的峰值
4. AR谱的统计性质 AR谱估计的方差反比于 的长度N和SNR
上式等效于Yule-Walker 方程，对同样的模型
系数 a1, , ap ，因此必有
当
时，可以用下式外推：
外推后的
对应AR谱，因此AR谱有较
高的分辨率。而经典谱估计中无外推，即：
分辨 率低
注意到AR模型自相关函数的匹配：
设想：如果阶次
, 则AR谱对应的自相关函
数完全等于信号的自相关函数，AR谱等于真谱。
保证：r(0) rx (0)
r( p) rx ( p)
的递推方法很多。
所以
x(n)
很多
最大熵功率谱
r(m) 1
P
()e jmd
2 MEMS
原则：x(n)是所有各种可能外推所对应的时间序
列中最随机的，即含有最大信息量（熵）。再假
定 x(n) 是高斯的。在这三个条件制约下，有：
PMEMS() PAR ()
第12章 参数模型功率谱估计
12.1 平稳随机信号的参数模型 12.2 AR模型的正则方程与参数计算 12.3 AR模型谱估计的性质与阶次选择 12.4 AR模型的稳定性与信号建模 12.5 关于线性预测 12.6 AR模型系数的求解算法 12.7 MA模型 12.8 ARMA模型 12.9 Pisarenko谐波分解与MUSIC 算法
（b) p=10; (c) p=20; (d) p=30
最大熵谱估计: Burg 于 1975年博士论文。 Maximum Entropy Spectral Estimation，MESE）
关于熵：
设信源由 X x1, x2, , xM 这 M 个 事件组成：
x 产生 i 的概率是 P(xi )
min ：最小
预测误差功率
p
rx (m) k rx(m k), m 1, 2, , p k 1
线性预测的Wiener-Hopf Eq.
注意到：对同一信号 x(n) ，都使用其 rx (m)
得到了两组方程：
来自AR模型： Yule-Walk 方程
来自LP： Wiener-Hopf
方程
结论：对同一信号，二者是相同的，即
当
的真 实功率谱
AR谱
有：
AR模型自 相关函数 匹配性质
增加 ，等效地扩大了
相等的部分
所以，理论上：我们可用一个全极点模型来近
似已知谱
，达到任意精度。
由：
（1）全局跟随性质（global)
在
内紧随
（2）局部跟随性质（local)
因为均值为1，所以 在 上下波动
情况多 情况少
从对整个 积分的贡 献来考虑
k ak k 1, 2, , p
min
2
一个 p 阶AR模型的系数可用来构成一个 p 阶的线性预测器，反之亦然。并且：
p
由于 e(n) x (n) xˆ(n) x(n) k x(n k)
p
k 1
u(n) x(n) ak x(n k)
k 1
所以
u(n) e(n) 等效的概念
上面等效的含意是： min 应等于AR模型激励白噪声的功率 2 。
2)
a2
0
rx (0)
a
p
0
p 1:
rx (0) rx (1)
rx (1) rx (0)
1
a1(1)
1
0
a1(1) rx (1) / rx (0) k1
1 rx (0) rx2 (1) / rx (0) rx (0)[1 a12(1)]
1 0[1 k12 ]
LSI系统的输入、输出关系：
差分方程 卷积关系 以上两式是LSI系统的时域表示，无论对确 定性信号还是随机信号都成立。现假定输入、 输出是平稳随机信号（输入是白噪声）。
转移函数的两 种表示形式， 独立于信号。
谱分解 的Z域 表示
待辨识 的参数。
Px (z)
u2H (z)H (z1)
u2
B( z ) B( z 1 ) A(z) A(z1)
AR谱 可以证明：
AR谱对应的自相 关函数
AR模型自 相关函数 匹配性质
证明： 由
Note : H (z) 1/ A(z)
两边做DTFT反变换：
左边 ph(k) ph(0) p (k)
p
右边 a(m) ra (m) a(k)ra (m k) k 0 p
有： ra (m) a(k)ra (m k) p (k) m 0 k 1
x(n)
线性预 测器
1 A(z)
xˆ(n)
e(n)
Yule-Walker 方程的快速计算
－Levinson-Durbin快速算法： 要求解的参数：
ap (1), ap (2), , ap ( p), 2(min p )
思路：
利用Toeplitz 矩阵特点，由低阶 高阶
am (k) : m 1, 2, , p
利用Yule-Walker 方程，可求解出AR模型参数：
a1, a2, , ap , 2
于是模型可以构造，可以实现功率谱估计。
为了深入了解AR模型的特点，现探 讨另外一个问题，即线性预测问题：
提法：设 x(n) 在 n 时刻之前的 p 个数据
x(n p), x(n p 1), , x(n 1) 已知
AR（Auto—Regressive，自回归）模型
若：
并假定：b0 1
则：
全动平均）模型
若：
则：
全 零 点 模 型
ARMA（Auto-Regressive MovingAverage,自回归－移动平均） 模型
如果：
ai : i 1 ~ p 不全 bi : i 1 ~ q 为零
rx
(1)
rx (0)
rx (1)
rx
(2)
rx (1)
rx (0)
rx ( p) rx ( p 1) rx ( p 2)
rx ( p) 1 2
rx ( p 1)
a1
0
rx ( p
2)
a2
0
rx (0)
a
p
0
Toeplitz 自相关阵
2
R
a
Op
又称 YuleWalker 方程
12.1 平稳随机信号的参数模型
经典谱估计： 分辨率低（受窗函数长度的限制）; 方差性能不好； 方差和分辨率之间的矛盾。
对平稳信号建模: 用于功率谱估计:提高分辨率，减小方差； 也可用于信号的特征提取，预测，编码及 数据压缩 等。
从功率谱估计的角度，对平稳信号建模的步骤：
步骤1
假定所研究的平稳过程 x(n) 是由一白噪声 序列 u(n) 激励一线性系统所产生的输出；
2 Makhoul J. Linear Prediction: a tutorial review. Proc. IEEE, 62(April):561-580,1975
3 Kay S M. Modern Spectrum Estimation: Theory and Application. 1988
递推过程中，要始终保持：
P 阶AR模型（LP）有三组参数： 都是 p+1 个
可互相导出，请给出它们互相导出的公式。 AR模型
基于AR模型谱估计的实现：
步骤1
由
估计
步骤2
解Yule-walker方程，得估计的模型参数
步骤3
尚需离 散化
实际计算：
aˆp1 aˆN1 0
离散谱，用FFT计算
12.3 AR模型谱估计的性质
不求导，使用正交原理：
E{x(n m)[x(n) xˆ(n)]} 0, m 1, 2, , p
e(n)