ARMA谱估计与系统辨识 清华大学《现代信号处理》讲义-张贤达
现代信号处理
根据处理对象和应用背景的不同而选择相应 的处理方法
8
课程主要内容
时域离散随机信号的分析 维纳滤波和卡尔曼滤波 自适应数字滤波器 功率谱估计 时频分析
9
第一章 时域离散随机信号的分析
1.1 1.2 1.3 1.4 1.5 1.6
均方值（二阶原点矩 ）：
方差（二阶中心矩 ）： 协方差：
D E[ X ]
2 2 2
x f ( x)dx
2
2 E[ X ]
x f ( x)dx
2
cov[ X , Y ] E[( X X )(Y Y )* ] E[ XY * ] E[ X ]E[Y ]*
Cxx ( X m , X n ) Rxx (m, n)
这种情况下， 自相关函数和自协方差函数没有什么区别。
23
互相关函数定义为
Rxy (m, n) E[ X Y ]
* m n
* xm yn f X m ,Yn ( xm , m, yn , n)dxmdyn
互协方差函数定义为
4
课程讨论的主要问题－1
对信号特性的分析
研究对象：确定性信号－>随机信号； 研究目的：提取信号中的有用信息； 主要内容：
随机信号的统计特性； 随机信号的参数建模； 功率谱估计（经典谱估计和现代谱估计）； 时频分析（短时傅立叶变换、维格纳变换、小波变换）
5
ARMA谱估计与系统辨识 清华大学《现代信号处理》讲义-张贤达
Ax b
b
A
1 x
0
-b A+-e Ez = 0 或 B + Dz = 0
扰动矩阵
总体最小二乘TLS: Total Least Squares
思想：寻求一个解z，使得
m
n1
1/ 2
2
dij min
i1 j1
定义代价函数
Σ
diag(
2 11
,
2 22
,
,
2 nn
)
主奇异值：p个大的奇异值（p个信号分量的能量） 次奇异值：其它小奇异值（扰动或误差的能量）
信号与噪声的分离：
准则一：归一化比值
v(k)
2 11
2 11
2 kk
1/
2
2 nn
1/
2
1
若阈值=0.995，v(k)>阈值的最小整数k定为矩阵A的“有效秩”。
其中：A(z) 1 a1z1 apz p B(z) 1 b1z1 bq zq
ARMA模型描述的线性时不变（LTI）系统
e(n) hi x(n)
传递函数：
H (z)
B(z) A( z )
hi z i
i
x(n) e(k )hnk e(n) hn k
bq1bq
c1
2b0bq
cq
非线性方程，MA参数辨识 （Newton－Raphson迭代）
协方差函数的Fourier变换
Px (z)
5.7 因果ARMA模型辨识 高斯有色噪声中的谐波恢复 清华大学《现代信号处理》讲义 -张贤达
∞
k 2
( j ) h ( j + n ) ∑ a (i ) h ( j + m i )
i =0
p
b( j + m )
∵ ∑ a (i ) x ( n i ) = ∑ b ( j ) e ( n j )
i=0 j =0
q
h
δ
修正Yule-Walker方程： 方程： 修正 方程
∑ a (i ) c
Φ 4 = Φ1 + Φ 2 + Φ 3
四阶累积量定义
c4 x (τ 1 ,τ 2 ,τ 3 ) = cum x * ( n ), x ( n + τ 1 ), x ( n + τ 2 ), x ( n + τ 3 ) = α (1)α (2)α (3)α * (4) e j (ω1τ1 +ω 2τ 2 +ω 3τ 3 ) +
τ1 = τ 2 = τ 3 = τ
c4 x (τ ,τ ,τ ) = ∑ α (i ) e jωiτ = c4 x (τ )
4
i =1 p
虚拟谐波： 虚拟谐波：
x ( n ) = ∑ α (i ) e j (ω i n +Φ i )
2 i =1
p
Rx (τ ) = ∑ α (i ) e jωiτ
f m +1 ( m, 0) ≠ 0
比检验 f ( m, 0) ≠ 0 更稳定
残差时间序列法： 残差时间序列法：
已辨识出， 假设 a (i ), i = 1, , p已辨识出，并滤波 y ( n ) = x ( n ) + v ( n )
y ( n ) = ∑ a (i ) y ( n i )
《现代信号处理》教学大纲
《现代信号处理》教学大纲适用专业：信息与通信工程、物联课程性质：学位课网工程、电子与通信学时数：32 学分数： 2课程号：M081001 开课学期：秋季第(1)学期大纲执笔人：何继爱大纲审核人：陈海燕一、课程的地位和教学目标现代信号处理作为信息类专业研究生的一门专业基础课，是在传统数字信号处理基础上，基于概率统计的思想，用数理统计、优化估计、线性代数和矩阵计算等工具，研究有限数据量的随机信号的分析与处理，且系统可能是时变、非线性的，它是近代才发展起来的前沿学科。
主要讨论基于信号模型分析和滤波的基本理论和基本方法；以现代谱估计和自适应滤波为核心内容，并介绍现代信号处理的新技术。
该课程为众多信号处理的应用领域打下基础，包括通信、声学、图像、雷达、声纳、生物医学等领域的信号处理。
本课程的知识目标是使学生牢固掌握现代信号处理一些最基本的理论、方法和应用，并能跟踪和学习新的理论、方法和技术；内容涉及随机信号统计分析、现代谱估计、自适应滤波器、时频分析与二次型时频分布、信号多速率变换、盲信分离和阵列信号处理方法等；建立现代信号处理的知识体系，对课程内容总体把握；具有一定的实验和模拟仿真的基本知识。
了解现代信号处理重要新技术的发展趋势，为从事信息与通信工程及相关电子系统的工程设计打下坚实的基础。
本课程的能力目标是通过课程的学习提高学生的分析计算方法、演绎推理方法和归纳法等基本数学处理方法；运用数学、物理及工程概念及方法发现问题、分析问题和解决问题的能力，以及理论与实际相结合的能力；能够触类旁通，提高学生的科学学习方法；掌握通信学科的信号分析与处理基本理论和技能，思路开阔，具有运用所学知识的能力、搜集和提炼信息的能力、团队合作能力、表达能力和创新能力等。
本课程的专业素质目标通过本课程的课堂学习、单元知识及章节总结、习题及专题研讨培养学生培养良好严谨的科学研究态度和正确的思维方法，使学生敢于提出问题、善于分析问题和解决问题的能力及具有团队合作精神。
清华大学《现代信号处理》课件
现代信号处理（离散随机信号处理）电子工程系本课程要讨论的主要问题：（1）对信号特性的了解随机信号（随机过程，时间序列––随机过程的一个实现）信号模型→参数估计→现代谱估计：参数化谱估计讨论信号模型及模型参数的估计问题，比较参数谱估计方法和周期图方法的优劣。
（2）对统计意义下最优滤波器设计的研究平稳条件下：Wiener滤波器理论非平稳条件下：Kalman滤波理论上的目标，实际算法可达到的最佳结果（3）对环境的自适应，具备“学习能力”的滤波算法自适应均衡、波束形成、线性自适应滤波器（4）更多信息的利用，挖掘（针对非高斯问题）线性系统、功率谱：二阶矩，高斯过程的完全刻划非线性、多谱：高阶量，循环平稳（5）对时间（空间）–––频率关系的适应性：全局特性与局域特性，小波变换，时频分析信号处理算法设计面向的几个主要因素n信噪比n先验知识n雷达n通信系统n电子对抗n对先验知识的利用：统计基础上的假设、学习过程n算法复杂性与性能要求的匹配性一些进展中的课题盲自适应信号处理序列贝叶斯估计、粒子滤波阵列信号处理等等与信号处理紧密关联的学科人工神经网络统计学习理论模式识别等等教材n张旭东，陆明泉：离散随机信号处理，2005年10月，清华大学出版社主要参考书①S. Haykin, Adaptive Filter theory, Third Edition, Prentice-Hall, 1996，//Fouth Edition 2001 (电子工业出版社均有影印本)①S.M. Kay, Modern Spectral Estimation: Theory & Application,Prentice-Hall, 1988①S.M. Kay, Fundamentals of Statistical Signal Processing: Estimation Theory, Prentice Hall PTR, 1993.①S. Mallat, A Wavelet Tour of Signal Processing, Academic press, 1998，Second Edition 1999①扬福生, 小波变换的工程分析与应用, 科学出版社, 2000.① D. G. Manolakis, et,al. Statistical and Adaptive Signal Processing, Mcgraw-Hall, 2000.①J. G. Proakis, et al. Algorithms for Statistical Signal Processing, Prentice hall, 2002①张贤达现代信号处理第2版清华大学出版社课程成绩n平时作业10％n2个Matlab作业40％（布置后2周内提交）n期末开卷考试50％1.1随机信号基础被噪声干扰的初相位是随机值的正弦波信号本质上均是随机的，但将信号作为随机信号处理，还是做为确定信号处理，与我们的应用目标和我们的先验知识有关，一般地，我们总是选择对应用有利的处理方式。
现代信号分析与处理技术_第1讲_参数估计方法
一、估计子的偏差和无偏估计
ˆ ˆ 1、θ 是θ 的无偏估计子：θ 满足
ˆ E (θ ) = θ ˆ ˆ ˆ 否则θ 是有偏估计子,估计的偏差为： b(θ ) = E (θ ) − θ
ˆ ˆ 2、θ 是θ 的渐近无偏估计子:若对所有θ , N → ∞ 时, b(θ ) → 0 .
1 N −1 ˆ 例 1、样本均值估计的无偏性： m x = ∑ xn N n =0 1 N −1 1 N −1 ˆ E [ m x ] == ∑ E[ xn ] = ∑ m x = m x 无偏估计 N n =0 N n =0
2
一般将式子右边的分母记着 I (θ ) ，称为 Fisher 信息量：
⎡ ∂ ⎤ I (θ ) = E ⎢ ln f ( x;θ ) ⎥ ⎣ ∂θ ⎦
2
Cramer-Rao不等式(对矢量参数的情况):(介绍)
若估计的参数是矢量 θ , 并将似然函数的对数表示为 L=lnf(x;θ), 则构造Fisher信息矩阵(p×p):
p列
⎡ r (0) r (1) ˆ = ⎢ r (1) r (0) Rx ⎢ ⎢ r (2) ( p) r (2) ( p − 1) x ⎣x
r ( p) ⎤ r ( p − 1) ⎥ = 1 XX T ⎥ N ⎥ (2) rx (0) ⎦
对r(1)(l)构造的自相关阵,没有上式的分解,所以不能保证半正定性.
例 2、样本方差估计的无偏性：
1 N −1 2 ˆ x = ∑ ( xn − m x ) 2 1) 均值 m x 已知时： σ N n=0 1 N −1 1 N −1 2 2 2 2 ˆ E [σ x ] = ∑ E [( xn − m x ) ] = ∑ σ x = σ x 无偏 N n=0 N n=0 1 N −1 2 2 ˆ ˆ ˆ 2) 均值取估计 m x 时： σ x = ∑ ( xn − m x ) N n =0 ˆ 记 m x = x 。由于各样本 xi 是独立同分布的,故有：
MUSIC方法_清华大学《现代信号处理》讲义_-张贤达
改进方法1: (求根MUSIC方法)
基本思想：Pisarenko谐波分解 (不需一维搜索)
a H ( )G 0
j
或
j ( m 1)
G H a( ) 0
T
a( ) 1, e , , e
z e j
p( z ) 1, z, , z
m 1 T
波束形成器：
w opt
1 H R xx a (d ) 1 H a(d )R xx a (d )
5. 改进的MUSIC方法
改进方法1：
ˆ ( ) a H ( )Ua P( ) H a ( )GG H a( )
p
ˆ 2 U
i 1
2 i
i
H s s 2 k k
观测空间 = 信号子空间 + 噪声子空间
特征值分解后，与大特征值对 应 与小特征值对 应
子空间的几何意义：
U S, G
H H H S S S S G H U U H S, G H I H G S G G G
S S I p , GH G Im p , G H S 0 S H G 0
Vandermonde矩 阵
j p e j ( m 1) p e 1
方向矩阵
满列秩 1 2 p
1 j1 e j ( m 1)1 e
1 e j2 e j ( m 1)2
2
加性噪声
2
1 lim N N
2
n 1
N
z (n) w H E x(n)x H (n) w
2
现代信号处理课件
P( H 0 ) H1 Lnl ( z ) Ln Ln ........( 1 28 ) H0 P( H1 )
则有 η＝1，Lnη＝0
21:20 24
§1-3最大后验概率准则 Maximum Posteriori Probability
称为最大后验概率准则，常简称为MAP准则。
即 p(z |H0) < p(z |H1)----(1-30) 时 判决为H1，否则判决为H0。 P(z | Hi), i=0, 1 为在给定观测值为z的条件下，Hi为真的概率， 此值为后验概率。
最大后验概率准则与最小总错误概率准则是等价的
21:20
26
例1: 设一个二元通信系统发送1V,0V的信号，受到2 为1/12w加性高斯噪声的干扰。系统发送1V 0V信号的 概率分别是0.6和0.4，代价分别为C00= -2, C01=8, C10= 6,
假设――所要检验的对象的可能情况或状态
检验――检测系统所做的判决过程
21:20 13
检测分类
二元检测：只有两种可能的假设
多元检测：有多个可能的假设 复合假设：信号是一随机过程的实现，其均 值或方差可处于某个数值范围内
序列检测：按取样观测值出现的次序进行处 理和判决
21:20 14
二元假设检验可能的情况
H0假设为真，判决H0（正确）；代价－C00 H1假设为真，判决H0（漏警）；代价－C01
H0假设为真，判决H1（虚警）；代价－C10 H1假设为真，判决H1（正确）；代价－C11
21:20 15
贝叶斯准则(Bayes)
代价、风险最小
源有两个输出，两个输出发生的概率已知，即先验概率已知P(H0), P(H1)分 别为假设H0和H1发生的概率。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
极点部分
极点的作用：决定系统的稳定性和因果性
A( z ) 0 即极点不在单位圆上
因果性：称x(n)是e(n)的因果函数，若
⑴ hi
i
⑵x(n) hi e(n i )
i 0
即因果系统要求极点在单位圆以内，A(z)的根|z|<1
零点的作用：决定系统的可逆性，即
1 z 1 的无数多项 中含有 A( z )
无限冲激响应(IIR)系统
白噪声中的AR过程： x(n) s(n) v(n)
v(n) ~ WN(0, v2)
B( z )
2
Px ( ) Ps ( ) Pv ( )
2 v2 w 2 A( z) A( z 1 ) z e j A( z)
n 0
N 1
2
谱窗
功率谱曲线平滑, 但分辨率下降
Px ( ) Rx (k ) w(k )e jkT
要提高分辨率，使用参数化的谱估计！ 经典谱估计：使用FFT的谱估计 现代谱估计：参数化谱估计
3.1 ARMA谱估计与系统辨识
平稳ARMA过程
离散随机过程 {x(n)} 服从线性差分方程：
ARMA功率谱估计的两种线性方法
Cadzow谱估计子
B( z ) B( z 1 ) N ( z ) N ( z 1 ) Px ( z ) 1 A( z ) A( z ) A( z ) A( z 1 )
2
N ( z) A( z 1 ) N ( z 1 ) A( z) 2 B( z) B( z 1 )
A( z ) x(n) B( z )e(n)
则功率谱
e(n)～N (0, 2 )
Px ( ) 2
其中
B( z ) A( z )
2
z e jw
B( z ) B( z 1 ) 2 A( z ) A( z 1 )
z e jw
A( z 1 ) 1 a1 z a p z p A* ( z ) B( z 1 ) 1 b1 z bq z q B* ( z )
若A(z)和B(z)无可对消公共因子，且 a p 0 ，则AR参数 a1 ,, a p 可由p个修正Yule-Walker方程唯一确定或辨识。
a R (l i) R (l )，
i 1 i x x
p
l q 1,, q p
Rx (q 1) R (q 2) x Rx (q p )
x ( n)
i
h e( n i )
i
Rx (k ) E{x( n) x* ( n k )}
* E hi e( n j ) hi e ( n k i ) i 0 j 0
hi h j E e(n j )e* (n k i )
i 0 i j 0 j
p
q
ห้องสมุดไป่ตู้
n
a R (l i) a h h
2 i 0 i x i 0 i j 0 j
j l i
2
h b
j 0 j
j l
a R (l i) 0，
i 0 i x
p
l q
修正Yule-Walker方程(MYW方程)
定理(AR参数的可辨识性)：
A = UΣV H 其中U为m m酉矩阵，V为n n酉矩阵。
酉矩阵： -1 = UH U
2 2 2 Σ diag(11, 22 ,, nn )
主奇异值：p个大的奇异值（p个信号分量的能量）
次奇异值：其它小奇异值（扰动或误差的能量）
信号与噪声的分离： 准则一：归一化比值
1 v(k ) 2 2 1/ 2 11 nn
x(n) ai x(n i) e(n) b j e(n j )
i 1 j 1
AR参数
MA参数
a x(n i) b e(n j)
i 0 i j 0 j
p
q
e(n) ~ N (0, 2 ) 后向移位算子：z j x(n) x(n j )
第三章
现代谱估计
清华大学自动化系 张贤达 zxd-dau@ 电话：62794875
经典谱估计
样本 直接法
假设已零均值化， 2k N
x(0), x(1),, x( N -1)
周期函数
N 1 n 0 jnT
X N ( ) x(n)e
Px ( )
i 0 j 0
hi h j 2 (k i j )
i 0 j 0
BBR公式： Rx (k )
2
h h
i 0
i ik
BBR公式 Rx (k )
p p
2
h h
i 0
i ik
a h( n i ) b ( n j ) b
2 11 2 kk 1/ 2
若阈值=0.995，v(k)>阈值的最小整数k定为矩阵A的“有效秩”。
准则二：使用归一化奇异值 kk kk ，且 11 1 11 kk <某个很小的阈值(0.05)的最小整数k定为有效秩。
AR阶数确定的信息量准则法
N p q 1 ˆ2 FPE ( p, q) wp N p q 1
p
p
k 0,1,, q
Kaveh谱估计子： Px ( )
ck z k k q
q
1 i 1 ai z
p
i
2
z e jw
ARMA功率谱密度的特例
A( z ) x(n) B( z )e(n)
特例一：MA过程
A( z ) 1
i
x(n) B( z )e(n)
H 1 ( z )
是否存在。
1 A( z ) H ( z ) B( z )
可逆性：称e(n)是x(n)的可逆函数，若
(1)存在序列 i ，并满足
(2)
i
i 0
i
——可逆系统的稳定性 ——可逆性条件
e( n ) i x ( n i )
ARMA过程的功率谱密度
Rx ( q 1 p) 1 0 Rx ( q 1) Rx ( q 2 p) a1 0 Rx (q p 1) Rx (q) a p 0 Rx (q)
又 其中
Px ( z )
k
Cx ( k ) z
k
(k ) z
k 0
k
(k ) z k
k 0
1 Cx (k ), (k ) 2 Cx (k ),
k 0 其他
则
i N ( z ) i 0 ni z (k ) z k A( z ) p ai z i k 0 p i 0
若构造：
Re
Rx (qe 1) R (q 2) x e Rx (qe M ) Rx (qe ) Rx (qe 1 pe ) Rx (qe 1) Rx (qe 2 pe ) Rx (qe M 1) Rx (qe M pe )
1 2 X N ( ) N
间接法
1 N 1 Rx (k ) x(n) x(n k ) N n 0
Px ( ) Rx (k )e jkT
n 0
N 1
周期图法
数据窗
有偏估计,平滑性差
加窗函数
1 Px ( ) N
N 1 k 0
x(n)c(n)e jnT
修正Yule-Walker方程
a x(n i) b e(n j)
i 0 p i j 0 q j
p
q
令e(n) (n) 令x(n) h(n)
n
a h( n i ) b ( n j ) b
i 0 i j 0 j
等价
高斯白噪 N (0, 2 )
x ( n)
k
e(k )h
nk
e(n) hn
满足ARMA模型的条件： (1)冲激响应系数必须绝对可求和： hk (系统稳定) (2)A(z)和B(z)无公共因子(p,q唯一)
k
(3)系统是物理可实现的(因果系统)
B( z ) H ( z) A( z )
p
ai z i ，比较系数得 两边同乘 i 0
nk ai (k i )
i 0
p
k 0,1,, p
ai 所以，Cadzow谱估计子的关键：估计AR阶数p和AR参数
Kaveh谱估计子
B( z ) B( z ) Px ( z ) 1 A( z ) A( z ) A( z ) A( z 1 )
i 1 p
线 谱 加性白噪声中的可预测过程： x(n) s(n) v(n)
x(n) ai x(n i) v(n) a j v(n j )
i 1 j 1
p
p
特殊的ARMA
所以：
白噪声中的AR过程 = ARMA过程 白噪声中的可预测过程 = 特殊的ARMA过程
B( z ) H ( z ) hi z 1 b1 z 1 bq z q A( z ) i