Y-W方程的L-D递推算法和伯格递推算法的比较
合集下载
相关主题
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
伯格(Burg)递推算法步骤
(1) 确定初始条件
e [ k ] e [ k ] y[ k ]
f 0 b 0
2 0
1 N
N 1 k 0
y 2 [k ]
(2) 从p=1开始迭代计算: 计算AR模型参数
Kp 2 e fp 1 [ k ]e b p 1 [ k 1]
谱估计结果——p=40,N=128
Periodogram 60 40 20 0 -20 -40 -60 0 0.1 0.2 0.3 0.4 0.5 Yule 60 40 20 0 -20 -40 -60 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0.6 0.7 0.8 0.9 1
[Pxx,f] = pburg (x,p,NFFT,Fs) x：进行功率谱估计的输入有限长序列； p: 模型的阶数 NFFT：DFT的点数； Fs ：绘制功率谱曲线的抽样频率，默认值为1； Pxx：功率谱估计值； f：Pxx值所对应的频率点
利用Burg法进行谱估计程序
N=512;NFFT=1024;Fs=2;p=40; n=0:N-1;randn('state',0); x=cos(0.3*pi*n)+cos(0.32*pi*n)+randn(size(n)); [P,f]=pyulear(x,p,NFFT,2); [Pw, f2]=pburg(x,p,NFFT,2); subplot(211);plot(f,10*log(P));grid;title('L-D'); axis([0 1 -30 60]); subplot(212);plot(f2,10*log(Pw));grid;title('Burg '); axis([0 1 -30 60]);
k p
N 1
N 1
e pf [ k ] K p
2e b p [k ]
eb p [k ] K p
} 0
可得
Kp
2 e pf 1 [ k ]e b p 1 [ k 1]
k p N 1 k p
{e
f p 1
[k ] e
2
b p 1
[ k 1] }
2
伯格(Burg)递推算法
K p eb p 1[k 1]
f K p e p1[k ]
1]
(4) 若阶数小于p，则阶数加1，回到步骤(2)进行下 一次迭代，直到达到预定阶数p。 (5) 估计功率谱
ˆ ( ) P AR
2
1 a p ( n ) e jn
n 1 p 2
Burg算法估计频谱的MATLAB函数
*
* AR模型参数与自相关函数的关系
Y-W 方程的矩阵表示
R y [ m ] a p ( n ) R y [ m n ] 2 [ m ] m 0 ,1,, p
n 1
p
例: p=3 时的Y-W 方程
R y [ 0] R y [1 ] R y [ 2] R y [ 3]
AR模型阶数p=50 的谱估计结果
L-D 60 40 20 0 -20 0 0.1 0.2 0.3 0.4 0.5 Burg 60 40 20 0 0.6 0.7 0.8 0.9 1
-20
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
AR模型阶数p=80 的谱估计结果
R y [1 ] R y [0] R y [1 ] R y [ 2]
R y [ 2] R y [1 ] R y [ 0] R y [1 ]
R y [ 3 ] 1 2 R y [ 2 ] a 3 ( 1 ) 0 R y [1] a 3 ( 2 ) 0 R y [ 0 ] 0 a 3 ( 3)
(1) 计算自相关函数的估计值 (2) 求解一阶模型参数关函数的估计值
a1 (1) R y [1] R y [ 0]
12 R y [ 0](1 a1 (1) )
2
Y-W方程的L-D递推算法
L-D算法估计功率谱的步骤 (3)由递推算法求解p阶模型参数
2 p 1 a p (n) a p1 (n) a p ( p)a p1 ( p n)
K2=a2(2) 2阶预测器的反射系数
伯格(Burg)递推算法
预测误差的递推公式
一般地
efp[k] efp1[k] K peb ] p 1[k 1
同理可得后向预测误差的递推公式
b f eb [ k ] e [ k 1 ] K e p p 1 p p 1[k]
Kp=ap(p)为p阶预测器的反射系数。
解此方程得
a1 (1)
2 1
R y [1] R y [ 0]
2
Ry [0] Ry [1]a1 (1) R y [ 0](1 a1 (1) )
Y-W方程的L-D递推算法
二阶Y-W方程的解
2 R y [ 0] R y [1] R y [ 2] 1 2 R y [1] R y [ 0] R y [1] a 2 (1) 0 R y [ 2] R y [1] R y [ 0] 0 a ( 2 ) 2
2 [1 a p ( p ) ] p 1 2 p 2
a p ( p)
R y [ p ] a p 1 ( n ) R y [ p n ]
n 1
p 1
(n 1,2,, p 1)
(4) 求出功率谱估计
ˆ ( ) P AR
p
2
2 n 1
1 a p ( n ) e jn
若已知Ry[n] ，由Y-W方程解出各参数a3(1), a3 (2),
a3 (3), 2。则可由AR模型参数获得功率谱Py()的估计值。
Y-W方程的L-D递推算法
一阶Y-W方程的解
R y [ 0] R y [1] 1 12 R [1] R [ 0] y y a1 (1) 0
k p N 1 k p N 1
2 { e fp 1 [ k ]2 e b [ k 1 ] } p 1
递推p阶均方误差
2 2 2 ( 1 K ) p p p 1
伯格(Burg)递推算法
伯格(Burg)递推算法步骤 (3) 递推高一阶前、后向预测误差
f ep [k ] b e p [k ] f e p1[k ] b e p1[k
利用L-D算法估计频谱的MATLAB函数
[Pxx,f] = pyulear(x,p,NFFT,Fs)
x：进行功率谱估计的输入有限长序列； p: 模型的阶数 NFFT：DFT的点数； Fs ：绘制功率谱曲线的抽样频率，默认值为1； Pxx：功率谱估计值； f：Pxx值所对应的频率点
例：利用L-D算法进行谱估计
L-D 60 40 20
0
-20 0 0.1 0.2 0.3 0.4 0.5 Burg 60 40 20 0 -20 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0.6 0.7 0.8 0.9 1
伯格(Burg)递推算法
L-D算法缺点： 在计算相关函数估计时，对N个观测数据以 外的数据作零的假设，故谱估计误差较大。 伯格(Burg)递推算法基本思想：
直接从观测的数据利用线性预测器的前向和 后向预测的总均方误差之和为最小的准则来估计 反射系数，进而通过L-D算法的递推公式求出AR 模型的优化参数。
一序列含有白噪声和两个频率的余弦信号，
x[ k ] cos( 0.3 πk ) cos( 0.4 πk ) [ k ]
利用L-D算法估计该序列的功率谱。
N=128;p=40; NFFT=2048;Fs=2; n=0:N-1;randn('state',0); x=cos(0.3*pi*n)+cos(0.4*pi*n)+randn(size(n)); [P,f]=periodogram(x,[],NFFT,2) ; [Py,fy]=pyulear(x,p,NFFT,2); subplot(211);plot(f,10*log(P));grid;title(' Periodogram '); axis([0 1 -60 60]); subplot(212);plot(fy,10*log(Py));grid;title(' Yule'); axis([0 1 -60 60]);
f e2 [k] y[k] a1(1) a2(2)a1(1)y[k 1] a2(2) y[k 2]
y[k] a1(1) y[k 1] a2(2) y[k 2] a1(1) y[k 1]
f f b e2 [k] e1 [k] K2e1 [k 1]
伯格(Burg)递推算法
前向预测误差的递推公式 f 2阶前向预测误差 e2 [k] y[k] a2(1) y[k 1] a2(2) y[k 2]
1阶后向预测误差
L-D算法的递推公式
b e1 [k] y[k 1] a1(1) y[k]
a2(1) a1(1) a2(2)a1(1)
n 1 2 p 1 p 1
a p ( p)
a p (n) a p1 (n) a p ( p)a p1 ( p n)
2 [1 a p ( p ) ] p 1 2 p 2
(n 1,2,, p 1)
Y-W方程的L-D递推算法
L-D算法估计功率谱的步骤
伯格(Burg)递推算法
预测误差滤波器的格形结构
efp[k] efp1[k] K peb ] p 1[k 1
b f eb [ k ] e [ k 1 ] K e p p 1 p p 1[k]
f y[k ] e0 [k ]
e1f [k ]
K1
e2f [k ]
K2
e pf [k ]
2 2 Ry [ 0] R y [1]
2 [ 1 a ( 2 ) ] Ry [0] Ry [1]a2 (1) Ry [ 2]a2 ( 2) 2 1
Y-W方程的L-D递推算法
p阶Y-W方程的递推解
R y [ p ] a p 1 ( n ) R y [ p n ]
Kp
z 1
b e0 [k ]
K1
z 1
K2
b e2 [k ]
z 1
Kp
eb p [k ]
e1b [k ]
伯格(Burg)递推算法
反射系数Kp 的确定
前向和后向均方预测误差的总和为
Ep
由
N 1 k p
2 {efp[k]2 eb [ k ] } p
E p K p
{2e pf [ k ]
a2 ( 2)
2 R y [0]R y [ 2] R y [1]
R [ 0 ] R [1]
2 y 2 y
R y [ 2] a1 (1) R y [1]
12
a1 (1) a 2 ( 2) a1 (1)
2
a2 (1)
2 2
R y [ 0] R y [1] R y [1] R y [ 2]
(1) 确定初始条件
e [ k ] e [ k ] y[ k ]
f 0 b 0
2 0
1 N
N 1 k 0
y 2 [k ]
(2) 从p=1开始迭代计算: 计算AR模型参数
Kp 2 e fp 1 [ k ]e b p 1 [ k 1]
谱估计结果——p=40,N=128
Periodogram 60 40 20 0 -20 -40 -60 0 0.1 0.2 0.3 0.4 0.5 Yule 60 40 20 0 -20 -40 -60 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0.6 0.7 0.8 0.9 1
[Pxx,f] = pburg (x,p,NFFT,Fs) x：进行功率谱估计的输入有限长序列； p: 模型的阶数 NFFT：DFT的点数； Fs ：绘制功率谱曲线的抽样频率，默认值为1； Pxx：功率谱估计值； f：Pxx值所对应的频率点
利用Burg法进行谱估计程序
N=512;NFFT=1024;Fs=2;p=40; n=0:N-1;randn('state',0); x=cos(0.3*pi*n)+cos(0.32*pi*n)+randn(size(n)); [P,f]=pyulear(x,p,NFFT,2); [Pw, f2]=pburg(x,p,NFFT,2); subplot(211);plot(f,10*log(P));grid;title('L-D'); axis([0 1 -30 60]); subplot(212);plot(f2,10*log(Pw));grid;title('Burg '); axis([0 1 -30 60]);
k p
N 1
N 1
e pf [ k ] K p
2e b p [k ]
eb p [k ] K p
} 0
可得
Kp
2 e pf 1 [ k ]e b p 1 [ k 1]
k p N 1 k p
{e
f p 1
[k ] e
2
b p 1
[ k 1] }
2
伯格(Burg)递推算法
K p eb p 1[k 1]
f K p e p1[k ]
1]
(4) 若阶数小于p，则阶数加1，回到步骤(2)进行下 一次迭代，直到达到预定阶数p。 (5) 估计功率谱
ˆ ( ) P AR
2
1 a p ( n ) e jn
n 1 p 2
Burg算法估计频谱的MATLAB函数
*
* AR模型参数与自相关函数的关系
Y-W 方程的矩阵表示
R y [ m ] a p ( n ) R y [ m n ] 2 [ m ] m 0 ,1,, p
n 1
p
例: p=3 时的Y-W 方程
R y [ 0] R y [1 ] R y [ 2] R y [ 3]
AR模型阶数p=50 的谱估计结果
L-D 60 40 20 0 -20 0 0.1 0.2 0.3 0.4 0.5 Burg 60 40 20 0 0.6 0.7 0.8 0.9 1
-20
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
AR模型阶数p=80 的谱估计结果
R y [1 ] R y [0] R y [1 ] R y [ 2]
R y [ 2] R y [1 ] R y [ 0] R y [1 ]
R y [ 3 ] 1 2 R y [ 2 ] a 3 ( 1 ) 0 R y [1] a 3 ( 2 ) 0 R y [ 0 ] 0 a 3 ( 3)
(1) 计算自相关函数的估计值 (2) 求解一阶模型参数关函数的估计值
a1 (1) R y [1] R y [ 0]
12 R y [ 0](1 a1 (1) )
2
Y-W方程的L-D递推算法
L-D算法估计功率谱的步骤 (3)由递推算法求解p阶模型参数
2 p 1 a p (n) a p1 (n) a p ( p)a p1 ( p n)
K2=a2(2) 2阶预测器的反射系数
伯格(Burg)递推算法
预测误差的递推公式
一般地
efp[k] efp1[k] K peb ] p 1[k 1
同理可得后向预测误差的递推公式
b f eb [ k ] e [ k 1 ] K e p p 1 p p 1[k]
Kp=ap(p)为p阶预测器的反射系数。
解此方程得
a1 (1)
2 1
R y [1] R y [ 0]
2
Ry [0] Ry [1]a1 (1) R y [ 0](1 a1 (1) )
Y-W方程的L-D递推算法
二阶Y-W方程的解
2 R y [ 0] R y [1] R y [ 2] 1 2 R y [1] R y [ 0] R y [1] a 2 (1) 0 R y [ 2] R y [1] R y [ 0] 0 a ( 2 ) 2
2 [1 a p ( p ) ] p 1 2 p 2
a p ( p)
R y [ p ] a p 1 ( n ) R y [ p n ]
n 1
p 1
(n 1,2,, p 1)
(4) 求出功率谱估计
ˆ ( ) P AR
p
2
2 n 1
1 a p ( n ) e jn
若已知Ry[n] ，由Y-W方程解出各参数a3(1), a3 (2),
a3 (3), 2。则可由AR模型参数获得功率谱Py()的估计值。
Y-W方程的L-D递推算法
一阶Y-W方程的解
R y [ 0] R y [1] 1 12 R [1] R [ 0] y y a1 (1) 0
k p N 1 k p N 1
2 { e fp 1 [ k ]2 e b [ k 1 ] } p 1
递推p阶均方误差
2 2 2 ( 1 K ) p p p 1
伯格(Burg)递推算法
伯格(Burg)递推算法步骤 (3) 递推高一阶前、后向预测误差
f ep [k ] b e p [k ] f e p1[k ] b e p1[k
利用L-D算法估计频谱的MATLAB函数
[Pxx,f] = pyulear(x,p,NFFT,Fs)
x：进行功率谱估计的输入有限长序列； p: 模型的阶数 NFFT：DFT的点数； Fs ：绘制功率谱曲线的抽样频率，默认值为1； Pxx：功率谱估计值； f：Pxx值所对应的频率点
例：利用L-D算法进行谱估计
L-D 60 40 20
0
-20 0 0.1 0.2 0.3 0.4 0.5 Burg 60 40 20 0 -20 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0.6 0.7 0.8 0.9 1
伯格(Burg)递推算法
L-D算法缺点： 在计算相关函数估计时，对N个观测数据以 外的数据作零的假设，故谱估计误差较大。 伯格(Burg)递推算法基本思想：
直接从观测的数据利用线性预测器的前向和 后向预测的总均方误差之和为最小的准则来估计 反射系数，进而通过L-D算法的递推公式求出AR 模型的优化参数。
一序列含有白噪声和两个频率的余弦信号，
x[ k ] cos( 0.3 πk ) cos( 0.4 πk ) [ k ]
利用L-D算法估计该序列的功率谱。
N=128;p=40; NFFT=2048;Fs=2; n=0:N-1;randn('state',0); x=cos(0.3*pi*n)+cos(0.4*pi*n)+randn(size(n)); [P,f]=periodogram(x,[],NFFT,2) ; [Py,fy]=pyulear(x,p,NFFT,2); subplot(211);plot(f,10*log(P));grid;title(' Periodogram '); axis([0 1 -60 60]); subplot(212);plot(fy,10*log(Py));grid;title(' Yule'); axis([0 1 -60 60]);
f e2 [k] y[k] a1(1) a2(2)a1(1)y[k 1] a2(2) y[k 2]
y[k] a1(1) y[k 1] a2(2) y[k 2] a1(1) y[k 1]
f f b e2 [k] e1 [k] K2e1 [k 1]
伯格(Burg)递推算法
前向预测误差的递推公式 f 2阶前向预测误差 e2 [k] y[k] a2(1) y[k 1] a2(2) y[k 2]
1阶后向预测误差
L-D算法的递推公式
b e1 [k] y[k 1] a1(1) y[k]
a2(1) a1(1) a2(2)a1(1)
n 1 2 p 1 p 1
a p ( p)
a p (n) a p1 (n) a p ( p)a p1 ( p n)
2 [1 a p ( p ) ] p 1 2 p 2
(n 1,2,, p 1)
Y-W方程的L-D递推算法
L-D算法估计功率谱的步骤
伯格(Burg)递推算法
预测误差滤波器的格形结构
efp[k] efp1[k] K peb ] p 1[k 1
b f eb [ k ] e [ k 1 ] K e p p 1 p p 1[k]
f y[k ] e0 [k ]
e1f [k ]
K1
e2f [k ]
K2
e pf [k ]
2 2 Ry [ 0] R y [1]
2 [ 1 a ( 2 ) ] Ry [0] Ry [1]a2 (1) Ry [ 2]a2 ( 2) 2 1
Y-W方程的L-D递推算法
p阶Y-W方程的递推解
R y [ p ] a p 1 ( n ) R y [ p n ]
Kp
z 1
b e0 [k ]
K1
z 1
K2
b e2 [k ]
z 1
Kp
eb p [k ]
e1b [k ]
伯格(Burg)递推算法
反射系数Kp 的确定
前向和后向均方预测误差的总和为
Ep
由
N 1 k p
2 {efp[k]2 eb [ k ] } p
E p K p
{2e pf [ k ]
a2 ( 2)
2 R y [0]R y [ 2] R y [1]
R [ 0 ] R [1]
2 y 2 y
R y [ 2] a1 (1) R y [1]
12
a1 (1) a 2 ( 2) a1 (1)
2
a2 (1)
2 2
R y [ 0] R y [1] R y [1] R y [ 2]