当前位置：文档库 › 用极大似然法进行参数估计

用极大似然法进行参数估计

北京工商大学

《系统辨识》课程

上机实验报告

（2014年秋季学期）

专业名称：控制工程

上机题目：极大似然法进行参数估计

专业班级：

2015年1 月

一实验目的

通过实验掌握极大似然法在系统参数辨识中的原理和应用。

二实验原理

1 极大似然原理

设有离散随机过程}{k V 与未知参数θ有关，假定已知概率分布密度)(θk V f 。如果我们得到n 个独立的观测值,21,V V …n V ,，则可得分布密度)(1θV f ，)(2θV f ，…,)(θn V f 。要求根据这些观测值来估计未知参数θ，估计的准则是观测值{}{k V }的出现概率为最大。为此，定义一个似然函数

)

()()(),,,(2121θθθθn n V f V f V f V V V L = （1.1）

上式的右边是n 个概率密度函数的连乘，似然函数L 是θ的函数。如果L 达到极大值，}{k V 的出现概率为最大。因此，极大似然法的实质就是求出使L 达到极大值的θ的估值∧

θ。为了便于求∧

θ，对式（1.1）等号两边取对数，则把连乘变成连加，即 ∑==n

i i

V f L 1

)(ln ln θ （1.2）

由于对数函数是单调递增函数，当L 取极大值时，lnL 也同时取极大值。求式（1.2）

对θ的偏导数，令偏导数为0，可得

0ln =??θL

（1.3）

解上式可得θ的极大似然估计ML ∧

θ。

2 系统参数的极大似然估计

Newton-Raphson 法实际上就是一种递推算法，可以用于在线辨识。不过它是一种依每L

次观测数据递推一次的算法，现在我们讨论的是每观测一次数据就递推计算一次参数估计值得算法。本质上说，它只是一种近似的极大似然法。

设系统的差分方程为 )()()()()(1

k k u z b k y z a ξ+=-- （2.1）式中

111()1...n n a z a z a z ---=+++

1101()...n n b z b b z b z ---=+++

因为)(k ξ是相关随机向量，故（2.1）可写成

)()()()()()(1

11k z c k u z b k y z a ε---+= （2.2）

式中

)()()(1

k k z c ξε=- （2.3） n

n z

c z c z c ---+++= 1

1)( （2.4）

)(k ε是均值为0的高斯分布白噪声序列。多项式)(1-z a ，)(1-z b 和)(1-z c 中的系数

n n c c b b a a ,,,,,10,1和序列)}({k ε的均方差σ都是未知参数。

设待估参数

n a a 1[=θ n b b 0 ]T

n c c 1 （2.5）并设)(k y 的预测值为

+-+++-----=∧

∧∧∧∧)()()()1()(01n k u b k u b n k y a k y a k y n n

)()1(1n k e c k e c n -++-∧

∧

（2.6）式中)(i k e -为预测误差；i a ∧，i b ∧，i c ∧

为i a ，i b ，i c 的估值。预测误差可表示为

+-+-???--=-=∑∑=∧

=∧

∧

)()()()()()(01

i k u b i k y a k y k y k y k e n i i n i i

-+++-+++=??

?--∧-∧∧-∧-∧=∧∑)()()()1()(1

10111k u z b z b b k y z a z a i k e c n n n n n i i )()(2

1k e z c z

c z c n n -∧

-∧

+++ （2.7）

或者

)()1(1

1k e z c z c n

n -∧

-∧

+++ =-+++-∧

-∧

)()1(1

1k y z a z a n

)()(1

10k u z b z b b n

n -∧

-∧

∧

+++ （2.8）因此预测误差{})(k e 满足关系式

)()()()()()(1

k u z b k y z a k e z c -∧

-∧

-= （2.9）式中

n n z a z a z a -∧

-∧

+++= 1

1)( n n z b z b b z b -∧

-∧

∧-∧

+++= 1

101

)( n n z c z c z c -∧

-∧

+++= 1

1)(

假定预测误差)(k e 服从均值为0的高斯分布，并设序列{})(k e 具有相同的方差2

σ。因

为{})(k e 与)(1

-∧

z c ，)(1

-∧

z a 和)(1

-∧

z b 有关，所以2

σ是被估参数θ的函数。为了书写方便，把式（2.9）写成

)()()()()()(1

11k u z b k y z a k e z c ----= （2.10）

-------++-+= )1()1()()1()()(101k u b k u b n k y a k y a k y k e n

,2,1),()1()(1++=------n n k n k c k e c n k u b n n （2.11）

或写成

)()()()()(1

i k e c i k u b i k y a k y k e n

i i

n i i

-----+

=∑∑∑=== （2.12）

令k=n+1,n+2,…,n+N,可得)(k e 的N 个方程式，把这N 个方程式写成向量-矩阵形式

θN N N Y e Φ-= （2.13）式中

????????????+++=)()2()1(N n y n y n y Y N ,????????????+++=)()2()1(N n e n e n e e N ,??

??????

???????????

?=n n b b a a 01θ

??????-+-+--=Φ)1()1()

(N n y n y n y N )()2()1(N y y y --- )()2()1(N n u n u n u +++ )()2()1(N u u u

)1()1()(-++N n e n e n e ?

???

?)()2()1(N e e e 因为已假定{})(k e 是均值为0的高斯噪声序列，高斯噪声序列的概率密度函数为

]

)(21exp[)

2(122

12m y f --

πσ （2.14）

式中y 为观测值，2

σ和m 为y 的方差和均值，那么

)](21exp[)

2(122

12k e f σ

πσ-

（2.15）对于)(k e 符合高斯噪声序列的极大似然函数为

)21exp()

2(1)]}()2()1([21exp{)

2(1])([])2([])1([])(,),2(),1([),(22

222

T N N

N N e e N n e n e n e N n e f n e f n e f N n e n e n e L Y L σ

πσσ

πσθθθθσθ-

=++++++-

+++=+++=

（2.16）

或

]2)()(exp[)

2(1),(2

2σθθπσσθΦ-Φ--=

N T N N

N Y Y Y L （2.17）

对上式（2.17）等号两边取对数得

N T N N

T N N

N e e N N e e Y L 2222

221ln 22ln 2)21exp(ln )

2(1ln

),(ln σ

σπσ

πσσθ---=-

+= （2.18）

或写为

∑++=---=N n n k N k e N N Y L 1

222

)(21ln 22ln 2),(ln σσπσθ （2.19）求),(ln σθN Y L 对2

σ的偏导数，令其等于0，可得

0)(212),(ln 1

422=+-=??∑++=N n n k N k e N Y L σσσσθ （2.20）

则

J N

k e N k e N N n n k N n n k 2

)(212)(112122

===∑∑++=++=∧

σ （2.21）式中

∑++==N n n k k e J 1

)(21 （2.22）

2σ越小越好，因为当方差2σ最小时，)(2k e 最小，即残差最小。因此希望2σ的估值取最

小

J N

min 2

∧σ （2.23）因为式（2.10）可理解为预测模型，而e(k)可看做预测误差。因此使式（2.22）最小就是使误差的平方之和最小，即使对概率密度不作任何假设，这样的准则也是有意义的。因此可按J 最小来求n n c c b b a a ,,,,,10,1的估计值。

由于e(k)式参数n n c c b b a a ,,,,,10,1的线性函数，因此J 是这些参数的二次型函数。求使),(ln σθN Y L 最大的∧

θ，等价于在式（2.10）的约束条件下求∧

θ使J 为最小。由于J 对

i c 是非线性的，因而求J 的极小值问题并不好解，只能用迭代方法求解。求J 极小值的常用

迭代算法有拉格朗日乘子法和牛顿-拉卜森法。下面介绍牛顿-拉卜森法。整个迭代计算步骤如下：

(1)确定初始的0∧θ值。对于0∧

θ中的n b b a a ,,,0,1可按模型

)()()()()(1

k u z b k y z a k e -∧

-∧

-= （2.24）用最小二乘法来求，而对于0∧

θ中的

n c c ,1可先假定一些值。

(2)计算预测误差

)()()(k y k y k e ∧

-= （2.25）给出

∑++==N n n k k e J 1

)

(21

并计算

∑++=∧=

n n k k e N 1

)

σ （2.26）

(3)计算J 的梯度θ??J

和海赛矩阵 2

2θ

??J

,有 θ

θ??=??∑++=)()(1k e k e J N n n k （2.27）式中

???????=??n a k e a k e k e )()()(1

θ n b k e b k e ????)()(0 T

n c k e c k e ???

????)()(1

--------++-+??

=??)()1()()()1()([)(101n k u b k u b k u b n k y a k y a k y a a k e n n i

i )]()1(1n k e c k e c n ----

i i a n k e c a k e c a k e c i k y ?-?--?-?-?-?--=)

()2()1()(21 （2.28）即

i n

j j

i a j k e c i k y a k e ?-?--=??∑=)

()()(1

（2.29）同理可得

i n

j j

i b j k e c i k u b k e ?-?---=??∑=)

()()(1 （2.30） i n j j

i c j k e c i k e c k e ?-?---=??∑=)

()()(1

（2.31）将式（2.29）移项化简，有

n j j

i n j j i a j k e c a j k e c a k e i k y ?-?=?-?+??=-∑∑==)

()()()(01 （2.32）因为

j z k e j k e -=-)()( （2.33）

由)(j k e -求偏导，故

i a z k e a j k e ??=

?-?-)()( （2.34）将（2.34）代入（2.32），所以

j n

j j i i j n j j i n

j j z c a k e a z k e c a j k e c i k y -=-==∑∑∑??=??=?-?=-00

0)()()()( （2.35） n n z c z c z c ---+++= 1111)(

所以得

)()

()

i k y a k e z c i

-=??- （2.36）

同理可得（2.30）和（2.31）为 )()

()

i k u b k e z c i

--=??- （2.37） )()

()(1

i k e c k e z c i

--=??- （2.38）根据（2.36）构造公式

)(])([)]

([)

(1i k y j j i k y a j i k e z c j

-=---=?--?- （2.39）

将其代入（2.36），可得

j a k e z c a j i k e z c ??=?--?--)

()

()]([)

(11 (2.40）消除)(1

-z c 可得

)

1()()(a i k e a j i k e a k e j i ?+-?=

?+-?=?? （2.41）同理可得（2.37）和（2.38）式

)

()()(b i k e b j i k e b k e j i ?-?=

?+-?=?? （2.42）

)

1()()(c i k e c j i k e c k e j i ?+-?=

?+-?=?? （2.43）式（2.29）、式（2.30）和式（2.31）均为差分方程，这些差分方程的初始条件为0，可通过求解这些差分方程，分别求出e(k)关于n n c c b b a a ,,,,,10,1的全部偏导数，而这些偏导数分别为)}({k y ，)}({k u 和)}({k e 的线性函数。下面求关于θ的二阶偏导数，即

∑∑++=++=??+??????????=??N n n k N

n n k k e k e k e k e J T

22122

)

()()()(θθθθ （2.44）当∧

θ接近于真值θ时，e(k)接近于0。在这种情况下，式（2.44）等号右边第2项接近

于0，2

2θ

??J

可近似表示为 T

N n n k k e k e J ??

???????=??∑++=θθθ)()(122 （2.45）

则利用式（2.45）计算2

2θ??J

比较简单。

(4)按牛顿-拉卜森计算θ的新估值1∧

θ，有

021201)(∧

???

?????-??-=-∧

∧

θθθ

θθJ J （2.46）重复(2)至(4)的计算步骤，经过r 次迭代计算之后可得r ∧

θ，近一步迭代计算可得

r J J r r ∧?????

???-??-=-∧

+∧

θθθθθ121)(2

（2.47）如果

42212

10-∧

∧

+∧

<-r

r σσσ

（2.48）

则可停止计算，否则继续迭代计算。

式（2.48）表明，当残差方差的计算误差小于%01.0时就停止计算。这一方法即使在噪声比较大的情况也能得到较好的估计值∧

θ。

三实验内容

设SISO 系统差分方程为

)()2()1()2()1()(2121k k u b k u b k y a k y a k y ξ+-+-=-+-+

其中极大似然估计法默认()e k 为：

112()()()()(1)(2)e k C z k k c k c k ξξξξ-==+-+-

辨识参数向量为 θ

＝1[a 2a 1b 2b c 1 c 2]T

式中，)(k ξ为噪声方差各异的白噪声或有色噪声；u(k)和y(k)表示系统的输入输出变量。

四实验流程图

五代码实现

程序如下：

U=[1.147 0.201 -0.787 -1.584 -1.052 0.866 1.152 1.573 0.626...

0.433 -0.958 0.810 -0.044 0.947 -1.474 -0.719 -0.086 1.099...

1.450 1.151 0.485 1.633 0.043 1.326 1.706 -0.340 0.890...

0.433 -1.177 -0.390 -0.982 1.435 -0.119 -0.769 -0.899 0.882...

-1.008 -0.844 0.628 -0.679 1.541 1.375 -0.984 -0.582 1.609...

0.090 -0.813 -0.428 -0.848 -0.410 0.048 -1.099 -1.108 0.259...

-1.627 -0.528 0.203 1.204 1.691 -1.235 -1.228 -1.267 0.309...

0.043 0.043 1.461 1.585 0.552 -0.601 -0.319 0.744 0.829...

-1.626 -0.127 -1.578 -0.822 1.469 -0.379 -0.212 0.178 0.493...

-0.056 -0.1294 1.228 -1.606 -0.382 -0.229 0.313 -0.161 -0.810... -0.277 0.983 -0.288 0.846 1.325 0.723 0.713 0.643 0.463...

0.786 1.161 0.850 -1.349 -0.596 1.512 0.795 -0.713 0.453...

-1.604 0.889 -0.938 0.056 0.829 -0.981 -1.232 1.327 -0.681...

0.114 -1.135 1.284 -1.201 0.758 0.590 -1.007 0.390 0.836...

-1.52 -1.053 -0.083 0.619 0.840 -1.258 -0.354 0.629 -0.242...

1.680 -1.236 -0.803 0.537 -1.100 1.417 -1.024 0.671 0.688...

-0.123 -0.952 0.232 -0.793 -1.138 1.154 0.206 1.196 1.013...

1.518 -0.553 -0.987 0.167 -1.445 0.630 1.255 0.311 -1.726...

0.975 1.718 1.360 1.667 1.111 1.018 0.078 -1.665 -0.760...

1.184 -0.614 0.994 -0.089 0.947 1.706 -0.395 1.222 -1.351...

0.231 1.425 0.114 -0.689 -0.704 1.070 0.262 1.610 1.489...

-1.602 0.020 -0.601 -0.020 -0.601 -0.235 1.245 1.226 -0.204...

0.926 -1.297] %输入数据

Y=[0.086 2.210 0.486 -1.812 -3.705 -2.688 1.577 2.883 3.705...

1.642 0.805 -

2.088 0.946 -0.039 1.984 -2.545 -1.727 -0.231...

2.440

3.583 2.915 1.443 3.598 0.702 2.638 3.611 -0.168...

1.732 0.666

2.377 -0.554 -2.088 2.698 0.189 -1.633 -2.010...

1.716 -1.641 -1.885 1.061 -0.968

2.911

3.088 -1.629 -1.533...

3.030 0.614 -1.483 -1.029 -1.948 -1.066 -0.113 -2.144 -2.626...

0.134 -3.043 -1.341 0.338 2.702 3.813 -1.924 -2.813 -1.795...

3.002 1.027 1.027 2.755 3.584 1.737 -0.837 -0.617 1.703...

2.045 -2.886 -0.542 -2.991 -1.859

3.045 0.068 -0.375 0.451...

1.036 0.153 -0.474

2.512 -2.681 -0.954 -0.307 0.628 -0.270...

-0.277 0.983 -0.288 0.846 1.325 0.723 1.750 1.401 1.340...

0.916 1.396 2.446 2.103 2.432 -1.486 3.031 2.373 -0.763...

-0.752 -3.207 1.385 -1.642 -0.118 1.756 -1.613 -1.690 2.136... -1.136 -0.005 -2.210 2.331 -2.204 0.983 1.347 -1.691 0.595...

1.809 -

2.204 -2.330 -0.454 1.290 2.080 -1.990 -0.770 1.240...

-0.252 3.137 -2.379 1.206 1.221 -1.977 2.471 -1.680 1.148...

1.816 0.055 -1.856 0.269 -1.323 -

2.486 1.958 0.823 2.481...

2.209

3.167 -0.762 -2.225 -0.123 -2.786 1.026 2.843 1.071...

-3.317 1.514 3.807 3.388 3.683 -1.935 -1.935 0.309 -3.390...

-2.124 2.192 -0.855 -1.656 0.016 1.804 3.774 -0.059 2.371...

-2.322 -0.032 2.632 0.565 -1.460 -1.839 1.917 0.865 3.180...

3.261 -2.755 -0.536 -1.171 -0.905 -3.303 -0.834 2.490 3.039...

0.134 1.901]%输出数据

na=2;nb=1;nc=2;d=1;

nn=max(na,nc);

L=size(Y,2);

xiek=zeros(nc,1); %白噪声估计初值

yfk=zeros(nn,1); %yf(k-i)

ufk=zeros(nn,1); %uf(k-i)

xiefk=zeros(nc,1); %vf(k-i)

thetae_1=zeros(na+nb+1+nc,1); %参数估计初值

P=eye(na+nb+1+nc);

for k=3:L

%构造向量

phi=[-Y(k-1);-Y(k-2);U(k-1);U(k-2);xiek]; %组建h（k）

xie=Y(k)-phi'*thetae_1;

phif=[-yfk(1:na);ufk(d:d+nb);xiefk];

%递推极大似然参数估计算法

K=P*phif/(1+phif'*P*phif);

thetae(:,k)=thetae_1+K*xie;

P=(eye(na+nb+1+nc)-K*phif')*P;

yf=Y(k)-thetae(na+nb+2:na+nb+1+nc,k)'*yfk(1:nc); %yf(k)

uf=U(k)-thetae(na+nb+2:na+nb+1+nc,k)'*ufk(1:nc); %uf(k)

xief=xie-thetae(na+nb+2:na+nb+1+nc,k)'*xiefk(1:nc); %vf(k)

%更新数据

thetae_1=thetae(:,k);

for i=nc:-1:2

xiek(i)=xiek(i-1);

xiefk(i)=xiefk(i-1);

end

xiek(1)=xie;

xiefk(1)=xief;

for i=nn:-1:2

yfk(i)=yfk(i-1);

ufk(i)=ufk(i-1);

end

yfk(1)=yf;

ufk(1)=uf;

end

thetae_1

figure(1)

plot([1:L],thetae(1:na,:));

xlabel('k'); ylabel('参数估计a');

legend('a_1','a_2'); axis([0 L -2 2]);

figure(2)

plot([1:L],thetae(na+1:na+nb+1,:));

xlabel('k'); ylabel('参数估计b');

legend('b_1','b_2'); axis([0 L -1.5 2]);

figure(3)

plot([1:L],thetae(na+nb+2:na+nb+nc+1,:));

xlabel('k'); ylabel('参数估计c');

legend('c_1','c_2'); axis([0 L -2 2]);

六实验结果

实验结果如下所示。 thetae_1 = -0.1222 0.0365 1.7962 0.0004 0.0655 -0.0018

可知，估计值a 1=-0.1222，a 2=0.0365，b 1=1.7962，b 2= 0.0004，c 1= 0.0655，c 2= -0.0018

参数估计a

图1 A 参数波形变化情况

参数估计b

参数估计c

图3 C 参数波形变化情况

图10-11 图10-12有162个数据，结果如下： thetae_1 = -0.7095 0.0489 1.8688 -0.9850 -0.1768 -0.0005

可知，估计值a 1=-0.7095，a 2=0.0489，b 1= 1.8688，b 2= -0.9850，c 1= -0.1768，c 2= -0.0005

参数估计a

参数估计b

图2 B 参数波形变化情况

参数估计c

图3 C 参数波形变化情况

六实验总结

通过此次实验，对极大似然法的原理和方法进行了研究和对其matlab 的仿真，可发现其离线辨识效果很好，但是其在线辨识的计算量很大，实验数据不同，出现的效果也不一样，递推的极大似然法呈现病态，所得参数估计取值很不稳定，以至于不能采用。

功率谱估计方法的比较

功率谱估计方法的比较摘要：本文归纳了信号处理中关键的一种分析方法, 即谱估计方法。概述了频谱估计中的周期图法、修正的协方差法和伯格递推法的原理，并且对此三种方法通过仿真做出了对比。关键词：功率谱估计；AR 模型；参数引言：谱估计是指用已观测到的一定数量的样本数据估计一个平稳随机信号的谱。由于谱中包含了信号的很多频率信息，所以分析谱、对谱进行估计是信号处理的重要容。谱估计技术发展渊源很长，它的应用领域十分广泛，遍及雷达、声纳、通信、地质勘探、天文、生物医学工程等众多领域，其容、方法都在不断更新，是一个具有强大生命力的研究领域。谱估计的理论和方法是伴随着随机信号统计量及其谱的发展而发展起来的，最早的谱估计方法是建立在基于二阶统计量, 即自相关函数的功率谱估计的方法上。功率谱估计的方法经历了经典谱估计法和现代谱估计法两个研究历程，在过去及现在相当长一段时间里，功率谱估计一直占据着谱估计理论里的核心位置。经典谱估计也成为线性谱估计，包括BT 法、周期图法。现代谱估计法也称为非线性普估计，包括自相关法、修正的协方差法、伯格（Burg ）递推法、特征分解法等等。原理：经典谱估计方法计算简单，其主要特点是谱估计与任何模型参数无关，是一类非参数化的方法。它的主要问题是：由于假定信号的自相关函数在数据的观测区间以外等于零，因此估计出来的功率谱很难与信号的真实功率谱相匹配。在一般情况下，经典法的渐进性能无法给出实际功率谱的一个满意的近似，因而是一种低分辨率的谱估计方法。现代谱估计方法使用参数化的模型，他们统称为参数化功率谱估计，由于这类方法能够给出比经典法高得多的频率分辨率，故又称为高分辨率方法。下面分别介绍周期图法、修正的协方差法和伯格递推法。修正的协方差法和伯格递推法采用的模型均为AR 模型。（1）周期图法周期图法是先估计自相关函数, 然后进行傅里叶变换得到功率谱。假设随机信号x(n)只观测到一段样本数据，n=0, 1, 2, …, N-1。根据这一段样本数据估计自相关函数，如公式（1）对（1）式进行傅里叶变换得到（2）式。 ∑--=+=1||0 *) ()(1 )(?m N n xx m n x n x N m r

功率谱估计

功率谱估计及其MATLAB仿真詹红艳（201121070630控制理论与控制工程）摘要:从介绍功率谱的估计原理入手分析了经典谱估计和现代谱估计两类估计方法的原理、各自特点及在Matlab中的实现方法。关键词:功率谱估计;周期图法;AR参数法;Matlab Power Spectrum Density Estimation and the simulation in Matlab Zhan Hongyan （201121070630Control theory and control engineering） Abstract:Mainly introduces the principles of classical PSD estimation and modern PSD estimation,discusses the characteristics of the methods of realization in Matlab.Moreover,It gives an example of each part in realization using Matlab functions. Keywords:PSDPstimation,Periodogram method,AR Parameter method,Matlab 1引言现代信号分析中,对于常见的具有各态历经的平稳随机信号,不可能用清楚的数学关系式来描述,但可以利用给定的N个样本数据估计一个平稳随机信号的功率谱密度叫做功率谱估计(PSD)。它是数字信号处理的重要研究内容之一。功率谱估计可以分为经典功率谱估计(非参数估计)和现代功率谱估计(参数估计)。功率谱估计在实际工程中有重要应用价值,如在语音信号识别、雷达杂波分析、波达方向估计、地震勘探信号处理、水声信号处理、系统辨识中非线性系统识别、物理光学中透镜干涉、流体力学的内波分析、太阳黑子活动周期研究等许多领域,发挥了重要作用。 Matlab是MathWorks公司于1982年推出的一套高性能的数值计算和可视化软件,人称矩阵实验室,它集数值分析、矩阵运算、信号处理和图形显示于一体,构成了一个方便的、界面友好的用户环境,成为目前极为流行的工程数学分析软件。也为数字信号处理进行理论学习、工程设计分析提供了相当便捷的途径。本文的仿真实验中,全部在Matlab6.5环境下调试通过;随机序列由频率不同的正弦信号加高斯白噪声组成。 2经典功率谱估计经典功率谱估计是将数据工作区外的未知数据假设为零,相当于数据加窗。经典功率谱估计方法分为:相关函数法(BT法)、周期图法以及两种改进的周期图估计法即平均周期图法和平滑平均周期图法,其中周期图法应用较多,具有代表性。 1.1相关函数法(BT法) 该方法先由序列x(n)估计出自相关函数R(n),然后对R(n)进行傅立叶变换,便得到x(n)的功率谱估计。当延迟与数据长度相比很小时,可以有良好的估计精度。 Matlab代码示例1: Fs=500;%采样频率 n=0:1/Fs:1;

参数法功率谱估计

参数法功率谱估计一、信号的产生（一）信号组成在本实验中，需要事先产生待估计的信号，为了使实验结果较为明显，我产生了由两个不同频率的正弦信号（频率差相对较大）和加性高斯白噪声组成的信号。（二）程序 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 -2 2 4 6 8 10 (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 k k z a z A z H 111)(1)(； ③ 212 1)(∑=-+=p k jwk k 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 p k x k 0=m 时。（二）两种AR 模型阶次的算法 1．Yule-Walker 算法(自相关法) （1）算法主要思想 Yule-Walker 算法通过解Yule-Walker 方程获得AR 模型参数。从低阶开始递推，直到阶次p ，给出了在每一个阶次时的所有参数。公式如下： ① 11 11/])()()([--=-∑+--=m m k x x m m m r k m r k a k ρ； ② )()()(11k m a k k a k a m m m m -+=--；

功率谱密度估计方法的MATLAB实现

功率谱密度估计方法的MATLAB实现在应用数学和物理学中，谱密度、功率谱密度和能量谱密度是一个用于信号的通用概念，它表示每赫兹的功率、每赫兹的能量这样的物理量纲。在物理学中，信号通常是波的形式，例如电磁波、随机振动或者声波。当波的频谱密度乘以一个适当的系数后将得到每单位频率波携带的功率，这被称为信号的功率谱密度（power spectral density, PSD）或者谱功率分布（spectral power distribution, SPD）。功率谱密度的单位通常用每赫兹的瓦特数（W/Hz）表示，或者使用波长而不是频率，即每纳米的瓦特数（W/nm）来表示。信号的功率谱密度当且仅当信号是广义的平稳过程的时候才存在。如果信号不是平稳过程，那么自相关函数一定是两个变量的函数，这样就不存在功率谱密度，但是可以使用类似的技术估计时变谱密度。信号功率谱的概念和应用是电子工程的基础，尤其是在电子通信系统中，例如无线电和微波通信、雷达以及相关系统。因此学习如何进行功率谱密度估计十分重要，借助于Matlab工具可以实现各种谱估计方法的模拟仿真并输出结果。下面对周期图法、修正周期图法、最大熵法、Levinson递推法和Burg法的功率谱密度估计方法进行程序设计及仿真并给出仿真结果。以下程序运行平台：Matlab R2015a（8.5.0.197613）一、周期图法谱估计程序 1、源程序 Fs=100000; %采样频率100kHz N=1024; %数据长度N=1024 n=0:N-1; t=n/Fs; xn=sin(2000*2*pi*t); %正弦波,f=2000Hz Y=awgn(xn,10); %加入信噪比为10db的高斯白噪声 subplot(2,1,1); plot(n,Y) title('信号') xlabel('时间');ylabel('幅度');

参数法功率谱估计

二、参数模型法功率谱估计（一）算法原理简介 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 k k z a z A z H 111) (1 )(； ③ 2 12 1)(∑=-+= p k jwk k 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 p k x k 0=m 时。

现代信号处理经典的功率谱估计

《现代信号处理》姓名：李建强学号：201512172087 专业：电子科学与技术作业内容：在MATLAB平台上对一个特定的平稳随机信号进行经典功率谱估计和现代功率谱估计的比较一、前言功率谱估计是信息学科中的研究热点,在过去的30多年里取得了飞速的发展。在许多工程应用中，它能给出被分析对象的能量随频率的分布情况。平滑周期图是一种计算简单的经典方法，它的主要特点是与任何模型参数无关，但估计出来的功率谱很难与信号的真是功率谱相匹配。与周期图方法不同，现代谱估计主要是针对经典谱估计(周期图和自相关法)的分辨率低和方差性能不好的问题而提出的。其使用参数化的模型，能够给出比周期图方法高得多的频率分辨率。其内容极其丰富，涉及的学科和领域也相当广泛，按是否有参数大致可分为参数模型估计和非参数模型估计，前者有AR模型、MA模型、ARMA模型、PRONY指数模型等;后者有最小方差方法、多分量的MUSIC方法等。二、总体概述本次实验分别使用经典的功率谱估计（如周期图法）与AR模型法对某一特定的平稳随机信号进行其功率谱估计，由图像得到信号的频率。利用MATLAB平台，直观形象地观察并比较二者估计效果的区别，以便于加深对功率谱估计的理解和掌握。三、具体的实现步骤 1、经典法功率谱估计周期图法又称直接法，它是从随机信号x(n)中截取N长的一段，把它视为能量有限的

真实功率谱的估计的一个抽样。 1.1、实现步骤（1）、模拟系统输出参数x(n)=A*sin(2πf1*n)+B*sin(2πf2*n)，包括序列长度N（128或512或1024，加性高斯白噪声（AGWN）功率一定，设置A，B，f1，f2，n的值。（2）、应用周期图法（不加窗）对信号的功率谱密度进行估计，使用直接法在MATLAB 平台上进行编程实现。（3）、输出相应波形图，进行观察，记录。 1.2 MATLAB源代码实现 clear all; %清除工作空间所有之前的变量 close all; %关闭之前的所有的figure clc; %清除命令行之前所有的文字 n=1:1:128; %设定采样点n=1-128 f1=0.2; %设定f1频率的值0.2 f2=0.213; %设定f2频率的值0.213 A=1; %取定第一个正弦函数的振幅 B=1; %取定第一个正弦函数的振幅 a=0; %设定相位为0 x1=A*sin(2*pi*f1*n+a)+B*sin(2*pi*f2*n+a); %定义x1函数，不添加高斯白噪声x2=awgn(x1,3); %在x1基础上添加加性高斯白噪声，信噪比为3,定义x2函数temp=0; %定义临时值，并规定初始值为0 temp=fft(x2,128); %对x2做快速傅里叶变换 pw1=abs(temp).^2/128; %对temp做经典功率估计

功率谱估计浅谈汇总

功率谱估计浅谈摘要：介绍了几种常用的经典功率谱估计与现代功率谱估计的方法原理，并利用Matlab对随机信号进行功率谱估计，对两种方法做出比较，分别给出其优缺点。关键词：功率谱；功率谱估计；经典功率谱估计；现代功率谱估计前言功率谱估计是从频率分析随机信号的一种方法,一般分成两大类:一类是经典谱估计;另一类是现代谱估计。由于经典谱估计中将数据工作区以外的未知数据假设为零，这相当于数据加窗，导致分辨率降低和谱估计不稳定。现代谱估计则不再简单地将观察区外的未知数据假设为零，而是先将信号的观测数据估计模型参数，按照求模型输出功率的方法估计信号功率谱，回避了数据观测区以外的数据假设问题。周期图、自相关法及其改进方法(Welch)为经典(非参数)谱估计方法, 其以相关和傅里叶变换为基础,对于长数据记录较适用,但无法根本解决频率分辨率低和谱估计稳定性的问题,特别是在数据记录很短的情况下,这一问题尤其突出。以随机过程的参数模型为基础的现代参数法功率谱估计具有更高的频率分辨率和更好的适应性,可实现信号检测或信噪分离,对语音、声纳雷达、电磁波及地震波等信号处理具有重要意义,并广泛应用于通信、自动控制、地球物理等领域。在现代参数法功率谱估计方法中,比较有效且实用的是AR模型法，Burg谱估计法,现代谱估计避免了计算相关,对短数据具有更强的适应性,从而弥补了经典谱估计法的不足,但其也有一些自身的缺陷。下面就给出这两类谱估计的简单原理介绍与方法实现。经典谱估计法经典法是基于传统的傅里叶变换。本文主要介绍一种方法：周期图法。周期图法由于对信号做功率谱估计，需要用计算机实现，如果是连续信号，则需要变换为离散信号。下面讨论离散随机信号序列的功率谱问题。连续时间随机信号的功率谱密度与自相关函数是一对傅里叶变换对，即：

随机信号的功率谱估计方法

数字信号处理II ——随机信号的功率谱估计方法

一、实验目的 1．利用自相关函数法和周期图法实现对随机信号的功率谱估计。 2．观察数据长度、自相关序列长度、信噪比、窗函数、平均次数等对谱估计的分辨率、稳定性、主瓣宽度和旁瓣效应的影响。 3．学习使用FFT 提高谱估计的运算速度。 4．体会非参数化功率谱估计方法的优缺点。二、实验原理与方法假设信号()x n 为平稳随机过程，其自相关序列定义为: {}*()()()m E x n x n m φ+ (0.1) 其中{}E ?表示取数学期望，{}* ?表示取共轭。根据定义，()x n 的功率谱密度()P w 与自相关序列()m φ存在如下关系： ()()j m m P m e ωωφ+∞ -=-∞ = ∑ (0.2) 1 ()()2j m m P e d π ωπφωωπ - = ? (0.3) 然而，实际中我们很难得到准确的自相关序列()m φ，只能通过随机信号的一段样本序列来估计信号的自相关序列，进而得到信号的功率谱估计。目前常用的线性谱估计方法有两种：自相关函数法和周期图方法，本实验将对这两种方法分别予以讨论。 1．自相关函数法假设已知随机信号()x n 的N 个观测样本，则其自相关序列可以用下式进行估计: ||1 *0 1 ?()()()||1||N m n m x n x n m m N N m φ --==+≤--∑ (0.4) 当仅使用长度为2M-1的自相关序列时，对其进行傅立叶变换即可得到功率谱估计如下:

11 ??()()M j m m M P m e ωωφ --=-+=∑ (0.5) 其中M 为加窗长度，Re ()c M W m 为矩形窗函数，定义如下: Re 1,||()0,||c M m M W m m M

周期图法估计功率谱

随机信号谱估计方法的Matlab实现摘要：功率谱估计是随机信号分析中的一个重要内容。从介绍功率谱的估计原理入手分析经典谱估计和现代谱估计两类估计方法的原理、各自特点及在Matlab中的实现方法。经典功率谱估计的方差大、谱分辨率差，分辨率反比于有效信号的长度，但现代谱估计的分辨率不受此限制。给出了功率谱估计的应用。关键词：功率谱估计；周期图法；AR参数法； 1 引言在一般工程实际中，随机信号通常是无限长的，例如，传感器的温漂，不可能得到无限长时间的无限个观察结果来获得完全准确的温漂情况，即随机信号总体的情况，一般只能在有限的时间内得到有限个结果，即有限个样本，根据经验来近似地估计总体的分布。有时，甚至不需要知道随机信号总体地分布，而只需要知道其数字特征，如均值、方差、均方值、相关函数、功率谱的比较精确的情况即估计值。功率谱估计(PSD)是用有限长的数据估计信号的功率谱，它对于认识一个随机信号或其他应用方面都是重要的，是数字信号处理的重要研究内容之一。功率谱估计可以分为经典谱估计(非参数估计)和现代谱估计(参数估计)。 2 .平均周期图法和平滑平均周期图法对于周期图的功率谱估计, 当数据长度N 太大时, 谱曲线起伏加剧, 若N 太小, 谱的分辨率又不好,因此需要改进。两种改进的估

计法是平均周期图法和平滑平均周期图法。 (1)Bartlett 法: Bartlett 平均周期图的方法是将N 点的有限长序列x(n)分段求周期图再平均。 Matlab 代码示例1: fs=600; n=0:1/fs:1; xn=cos(2*pi*20*n)+3*cos(2*pi*90*n)+randn(size(n)); nfft=512; window=hamming(nfft); %矩形窗 noverlap=0;%数据无重叠 p=0.9;%置信概率 [Pxx,Pxxc]=psd(xn,nfft,fs,window,noverlap,p); index=0:round(nfft/2- 1); k=index*fs/nfft; plot_Pxx=10*log10(Pxx(index+1)); plot_Pxxc=10*log10(Pxxc(index+1)); figure(1) plot(k,plot_Pxx); figure(2) plot(k,[plot_Pxx plot_Pxx- plot_Pxxc

功率谱估计方法综述

功率谱估计方法综述：简介：随机信号的持续时间是无限长的，因此随机信号的总能量是无限的，因而随机过程的任意一个样本寒暑都不满足绝对可积条件，所以其傅里叶变换不存在。尽管随机信号的总能量是无限的，但其平均功率却是有限的，因此，要对随机信号的频域进行分析，应从功率谱出发进行研究才有意义。信号的功率谱密度描述随机信号的功率在频域随频率的分布。功率谱估计(PSD)是用有限长的数据来估计信号的功率谱，即利用给定的N个样本数据估计一个平稳随机信号的功率谱密度。背景：功率谱估计在实际工程中有重要应用价值，如在语音信号识别、雷达杂波分析、波达方向估计、地震勘探信号处理、水声信号处理、系统辨识中非线性系统识别、物理光学中透镜干涉、流体力学的内波分析、太阳黑子活动周期研究等许多领域，发挥了重要作用。功率谱估计方法主要分为2大类：非参数化方法（又称经典功率谱估计）和参数化方法（又称现代功率谱估计）。非参数化方法有相关函数法(BT法)、周期图法、平均周期图法、平滑平均周期图法等；而参数化谱估计有R模型法、移动平均模型法(简称MA模型法)、自回归移动平均模型法(简称ARMA模型法)、最大熵谱分析法(AR模型法)、Pisarenko谐波分解法、Prony 提取极点法、Prony谱线分解法以及capon最大似然法等，由于涉及许多复杂数学计算，在此未作详细数学推导，以下介绍几种常用的功率谱估计方法一、非参数化方法（经典法）经典功率谱估计是将数据工作区外的未知数据假设为零，相当于数据加窗。 1、自相关法又称相关函数法(BT法)，根据维纳—辛钦定理：平稳随机过程的自相关函数和功率谱函数是一傅里叶变换对，对于平稳随机信号来说，其相关函数是确定性函数，故其功率谱也是确定的．这样可由平稳随机离散信号的有限个离散值，求出自相关函数，然后作Fourier变换，得到功率谱。由于随机序列{X(n)}的自相关函数R(n)=E[X(n)X(n+m)]定义在离散点m上，设取样间隔为错误！未找到引用源。，可将随机序列的自相关函数用连续时间函数表示为错误！未找到引用源。等式两边取傅里叶变换，则随机序列的功率谱密度错误！未找到引用源。错误！未找到引用源。错误！未找到引用源。 BT法是先估计自相关函数Rx(m)(|m|=0,1,2…,N-1),然后再经过离散傅里叶变换求的功率谱密度的估值错误！未找到引用源。。即错误！未找到引用源。其中错误！未找到引用源。可有式得到。

burg法实现功率谱估计

用burg 法实现功率谱估计参数模型法是现代谱估计中的主要内容，AR 模型参数的求解有三种方法：自相关法、Burg 递推算法和改进协方差法。Burg 算法不是直接估计AR 模型的参数，而是先估计反射系数Km,再利用Levinson 关系式求得AR 模型的参数。Burg 算法采用的数据加窗方法是协方差法，不含有对已知数据段之外的数据做人为的假设。 1.其原理如下： Burg 算法是使前向预测误差和后向预测误差均方误差之和最小来求取Km 的，它不对已知数据段之外的数据做认为假设。计算m 阶预测误差的递推表示公式如下： x(n) (n)(n)(n)1)-(n (n)1) -(n (n)(n)0f 0f 1-m m 1-b m 1-m f 1-m m e e e e ==+=+=e k e e k e b b m b m f 求取反射系数的公式如下： } 1)]-(n [(n)]{[1)] -(n (n)[2-2b 1-m 2f 1-m b 1-m f 1-m m e e e e +=E E k 对于平稳随机过程，可以用时间平均代替集合平均，因此上式可写成： [][][][]{}p ,2,1,1)-(n (n)1)-(n (n)2-1-2 1-2 1-1-m n 1-1-，?=+=∑∑==m N m n b m f m N b m f m m e e e e k 这样便可求得AR 模型的反射系数。将m 阶AR 模型的反射系数和m-1阶AR 模型的系数代入到Levinson 关系式中，可以求得AR 模型其他的p-1个参数。 Levinson 关系式如下： 1-m 1,2,i i),-(m (i)(i)1-m 1-m m ，?=+=a k a a m m 阶AR 模型的第m+1个参数G ， ρm 2 G =, 其中ρm 是预测误差功率，可由递推公式

功率谱估计

功率谱估计浅谈

————————————————————————————————作者：————————————————————————————————日期：