数字信号处理课程论文
Xxxxxxxxx
从经典谱估计到参数模型功率谱估计
经典的两种功率谱估计方法是由Fourier变换作为理论基础的，一种是19世纪末由舒适特直接利用Fourier级数去拟合某类信号时提出的周期图法，又称直接法，即对给定的有限长且实的观测序列x（n）（0 n N 1），作Fourier变换
X(w) x(n)e
n 0
N 1
jwn
（1）
可得到周期图的计算公式为
IN(w)
12
X(w) （2） N
另一种是由布莱克曼—图基运用平稳过程的维纳—辛欽定理，由采样数据序列间接得到的功率谱估计法，又称间接法。即对平稳随机信号作N次观测，得到N个观测值x（0），x（1），。。。x（N-1）；然后作自相关估计：
(m) 1 RN
N
N m 1
n 0
x(n)x(n m)
m M 1(M N) （3）
(m)作离散Fourier变换得到功率谱估计 对RN
(w) S
可记：
M 1m (M 1)
(m) jwmRNe
M N （4）
IN(w)
1
X(w)X*(w) N
m (N 1)
N 1
(m) jwm （5） RNe
(m)（易证它是有偏估计）式（5）表示周期图就是自相关估计R在M=N时的Fourier变换，N
同时表明当M=N时，直接法与间接法所得结果相同。
由于不管观测数据中N多大，直接、间接方法得出的估计均非理想的估计，他们的主要缺陷为：（1）弱信号被强信号的旁瓣淹没；（2分辨率约为数据长度的倒数，导致分辨率不高；（3）频谱存有旁瓣，导致出现泄漏现象而使主瓣失真。
现代谱估计方法是以随机过程的参数模型为基础的，早在1967年，美国学者伯格（Burg）提出了最大熵谱估计方法，其基本思想是：已知有限（N）个相关序列函数的功率谱，估计值是某一最不可预测或最随机序列的功率谱，则该序列的前N个相关函数值与待估序列的已知相关函数值相等。用信息论观点即最不确定或最随机事件的熵最大，故起名为最大熵估计。Burg方法利用了LD算法的优点，使算法稳定且计算量不大，但发现用于正弦波信号
的功率谱估计时，有谱线分裂和频谱偏移现象，后来者提出了前、后向AR参数估计法、Marple法、迭代法、数据自适应加权法等种种改进算法，以使谱线分裂现象消除并提高频率分辨率。1971年Van der Bos证明了一维最大熵估计与自回归谱估计的等效性，1973年Pisarenko提出的谐波分析法为谱估计问题的讨论开辟了新方向。1979年Schmidt公开发表了被后人称为MUSIC算法的多信号分类算法，可以认为是Pisarenko方法的推广。此外还有诸如Capon的最大使然谱估计方法、Prony的谱密度与谱线估计方法、Shore的最小交叉熵谱估计方法以及Lp范数最小化的谱估计方法等等。目前，国内外学者仍在不断努力，致力于构造出更好、更有效、更实用的方法，不断推动现代谱估计方法的向前发展。
参数模型化的思路：
a． 假定所研究的过程x（n）是由一个输入序列u（n）激励一个线性系统H（z）的输出，
如图1-1所示。 b． 由已知的x（n），或其自相关函数rx(m)来估计H（z）的函数。 c． 由H（z）的参数来估计x（n）的功率谱。
u(n)x(n)
H(z)
图1-1 参数模型
在上图所示，H(z)是一个因果的线性移不变离散时间系统，当然，它应该是稳定的，其单位抽样响应h（n）是确定的。输出序列x（n）可以是平稳的随机序列，也可以是稳定性的时间序列。若x（n）是确定性的，那么u（n）是一个冲激序列，若x（n）是随机的，那么u（n）应是一个白噪声序列。
不论x（n）是确定信号还是随机信号，对上图的线性系统，u（n）和x（n）之间总有如下的输入、输出关系：
x(n) akx(n k) bku(n k) （6）
k 1
k 0
pq
及
x(n) h(k)u(n k) （7）
k 0
对式子（6）和（7）两边分别取Z变换，并假定b0 1，可得
H(z)
式中：
B(z)
（8） A(z)
A(z) 1 akz k （9-a）
k 1
p
B(z) 1 bkz k （9-b）
k 1
q
H(z) h(k)z k （9-c）
k 0
为了保证H（z）是一个稳定的且是最小相位系统，A（z），B（z）的零点都应在单位圆内。
假定u（n）是一个方差为 的白噪声序列，由随机信号通过线性系统的理论可知，输出序列x（n）的功率谱：
2
Px(ejw)
B(e)B(e)
A(e)A(e)
*
jw2
2jw*jw
jw
2B(ejw)
A(e)
jw
2
2
（10）
这样，如果激励白噪声的方差 及模型参数a1，a2.。。。，ap；b1，b2，。。。bq已知，那么由上式可求出序列x（n）的功率谱。
现对（6）式分三种情况来讨论：
①如果b1,b2, ,bq全为零，那么（6）（8）（10）式分别变为：
x(n) akx(n k) u(n) （11）
k 1
p
H(z)
1
A(z)
11 akz k
k 1p
（12）
Px(e)
jw
2
ake jwk
k 1p
2
（13）
此三式给出的模型称为自回归模型，简称AR模型，它是一个全极点的模型。自回归的定义是：该模型现在的输出时现在的输入和过去p个输出的加权和。 ②如果a1,a2, ,ap全为零，那么（6）（8）（10）式子分别变为：
x(n) bku(n k) u(n) bku(n k),b0 1 （14）
k 0
k 1
H(z) B(z) 1 bkz k （15）
k 1
q
Px(ejw) 2 bke jwk （16）
k 1
q
2
此三式给出的模型称为移动平均模型，简称MA模型，它是一个全零点的模型。
③若a1,a2, ,ap；b1,b2, ,bq不全为零，则（6）式给出的模型是自回归移动平均模型，
简称ARMA模型。显然，ARMA模型是一个既有零点、又有极点的模型
AR，MA，ARMA是功率谱估计中最主要的参数模型，由后面的讨论知，AR模型的正则方程是一组线性方程，而MA和ARMA模型是非线性方程。由于AR模型具有一系列好的性能，因此，是被研究最多并获得管饭应用的一种模型。本文将较为详细地讨论AR模型参数的计算、谱的性能及与其他算法（线性预测和最大熵谱估计等）的关系。
AR模型的正则方程与参数计算：假定w（n）、x（n）都是实平稳的随机信号，w
（n）为白噪声，方差为 ，现在我们希望建立AR模型的参数a1,a2, ,ap和x（n）的自相关函数的关系，也即AR模型的正则方程。 按定义
2
xx(m) E x(n)x(n m) （17）
将式（11）的关系代入上式，得
p
xx(m) E x(n) akx(n k m) (n m)
k 1
ak xx(m k) E x(n) (n m)
k 1
p
（18）
按式（11），x（n）只与 (n)相关而与 (n m)(m 1)无关，故式（18）中的第二项为
m 0 0,
E x(n) (n m) E[( akx(n k) (n)) (n m)] E (n) (n m) 2
m 0
代入式（18）
（19）
p
ak xx(m k) , m 0 1
xx(m) kp
a ( k) 2,m 0 kxx k 1
或
p
xx(m) ak xx(m k) 0 , m 0 k 1 p
(m) a ( k) 2, m 0
xxkxx
k 1
（20）
将m=1， ，p分别代入式（19）并写成矩阵形式，得
xx( 1) xx( 2) xx( (p-1)) a1 xx(1) xx(0)
(1) a (2) (0) ( 1) ( (p-2))xxxxxx xx 2 xx
a (p 1) (p 2) (p 3) (0) (p) p xx xxxxxx xx
如果利用自相关函数的偶对称性，r（m）=r（-m），则式子可写为：
2
xx( 1) xx( 2) xx( p) xx(0) 1
(1) a (0) ( 1) ( (p-1))xxxxxxxx 1 0
xx(p 1) xx(p 2) xx(0) xx(p) 0 ap
（21）
2
（22）
式（21）及式（22）就是Yule-Walker方程。由式（21）或(22)可解得各ak及 。
设x（n）在n时刻之前的p个数据{x(n-p),x(n-p+1),…,x(n-1)}已知，我们希望利用这p个数据来预测n时刻的值x（n），现在我们用线性预测的方法实现：
(n) akx(n k) （23） x
k 1
p
(n)和真实值x（n）之间的误差e（n）记预测值x，则
(n) （24） e(n) x(n) x
因此总的预测误差总功率为：
E{e(n)} E{[x(n) akx(n k)]2 （25）
2
k 1
p
根据正交原理，为求使得p最小的ak值，则：
rx(m) akrx(m k) （26）
k 1
p
将式（22）与式（26）比较可见，如果二式中的
xx( )是代表同一个随机信号 x(n) 的自相
关函数，则式（20）与式（26）的区别仅仅在于式（26）中的l=1,2, , p(式（26）中的l相当于式（20）中的m)，而式（20）中的m可以取任何大于零的整数，包括1,2, , p，也包括p的值。如果我们取式（20）的m=1,2, , p的p个方程(即式(22)，解得ak(k=1,2, , p)的p个参量，则这p个ak必将等于利用式(26)的p个方程解得的p个apk，即
ak apk k 1, 2, , p
（27）
我们再将式(26)与式(20)作比较，由于二式的又有
xx( )是代表同一个 x(n) 的自相关函数，
ak apk
，故有
2 E e2(n) min
（28）
于是，式(26)全同于式20)。
式(26)与(20)说明，如果一个随机序列x(n)的形成系统是一个p阶的AR模型(即x(n)可由一白噪激励一p阶AR模型产生)，形成系统的传递函数为
H(z)
11 akz k
k 1p
（29）
其功率谱密度为
Pxx( ) H( )
2
2
2
akz k
k 1
z ej
p
2
（30）
2
Ee(n)min，
则其中的ak与 分别等于x(n)的p阶线性预测器的系数apk与最小均方误差
2
因此AR模型法又称为线性预测AR模型法。这种估计功率谱密度的方法可归结为求AR模型的各ak参数的问题，也可以归结为求在最小均方误差准则下的线性预测器的系数apk的问题。
最大熵谱估计法的基本思想及其与线性预测AR模型法的等价性
在前面中已经讨论到经典法(BT PSD法)用已知的有限个(N个)自相关函数序列的估值求功率谱估计时，是将此有限个估值以外的自相关序列的数据认为是零，因而得不到好的分辨率。J.P.Burg于1967年提出的MESE法与此不同，它是基于将已知的有限长度的自相关序列以外的数据用外推法求得，而不是把它们当作是零。如果假设已知
xx(0), xx(1), , xx(N), 问题在于按什么原则外推 xx(N 1), xx(N 2), 。在保证自相关
函数的Toeplitz矩阵是正定的情况下有无穷多种外推法，Burg认为外推的自相关函数应使时间序列表现出最大熵，因此把Burg提出的这种方法称之为最大熵谱估计法。
所谓熵是代表一种不定度，最大熵为最大不定度，即它的时间序列最随机，而它的PSD应是最平伏(最白色)。按Shannon对熵的定义，当x的取值为离散的时，熵H定义为
H pilnpi
i
（31）
这里pi为出现状态i的概率。当x的取值为连续的时，熵被定义为
H p(x)lnp(x)dx E ln(p(x))
（32）
这里p(x)为概率密度函数。当我们处理时间序列由联合概率密度函数维高斯分布为
x0, x1, , xN的实现问题时，概率密度应
p(x0, x1, , xN)代替。设 xn 为零均值，高斯分布的随机过程。一
p(x)
221
e x/2
2
（33）
其中
2 x2p(x)dx
（34）
N维高斯分布为
p(x1, x2, , xN) (2 )
其中
-N/2
1 1
(det (N))exp X (N) X
2
（35）
x1
x X 2
xN
2 11 12 1N 2 21222N (N)
2
N2 NN N1
xx( 1) xx( N) xx(0) xx(0) xx( (N-1)) xx(1)
xx(N) xx(N 1) xx(0)
2 ii Var[xi] E xi2 xx(0)
(36)
这里
ij Cov[xixj] Exixj xx(j i)
(37)
下面我们先求一维高斯分布的信号的熵，然后推广到N维。得
x2 2H p(x)ln p(x) dx p(x) ln2 dx2 2
ln2
2
x p(x)dx
2
p(x)dx
2 2
因为
22
p(x)dx 1及xp(x)dx
代入上式得到一维高斯分布的熵为
H ln2 2
1
ln2 2 ln ln2 2e2
(38)
同理可求得N维高斯分布信号的熵为
H ln(2 e)N/2(det )1/2
(39)
式中det 代表矩阵 的行列式，要使熵H最大，就要求det 最大。
如果已知 xx(0), xx(1), , xx(N), 现欲求得 xx(N 1)。由于自相关函数的矩阵必是正定的，故矩阵 (N 1)的行列必大于零，即
det (N 1)
xx(0) xx(1) xx(N 1) xx(1) xx(0) xx(N)
xx(N 1) xx(N) xx(0)
0
(40)
d (N 1)
为了得到最大熵，要求det (N) 最大，为此用 xx(N 1)对上式微分，使d xx(N 1)求得使det (N 1) 最大的 xx(N 1)，满足下列方程：
0
，
xx(1) xx(0) xx(N 1)xx(2) xx(1) xx(N 2)
xx(N 1) xx(N) xx(1)
0
(41)
上式是 xx(N 1)的一次函数，由此式可解出 xx(N 1)。于是又可以此 xx(N 1)为已知，再用类似方法求得 xx(N 2)，以此类推。这样每步都按最大熵的原则外推后一个自相关序列的值，可以外推到任意多个而不必认为它们是零。这就是最大熵谱估计法的基本思想。 可以证明这种按最大熵外推自相关函数的结果与AR模型是等价的。为了证明这一点，将式(20)，即
xx(m) ak xx(m k) 0 m 0
k 1
N
(20)
中的m分别用1, 2, , N+1代入，写成下列方程组：
xx(1) a1 xx(0) aN xx(N 1) 0 xx(2) a1 xx(1) aN xx(N 2) 0
xx(N) a1 xx(N 1) aN xx(0) 0
(42) (43)
xx(N 1) a1 xx(N) aN xx(1) 0
式（42)即式(20)只是把阶数p写成N。如果我们从式(42)的N个线性方程中解得的N个AR参数a1, a2, , aN值，代入式(5.88)并将其整理成行列式的形式，即可得
xx(1) xx(0) xx(N 1)xx(2) xx(1) xx(N 2)
xx(N 1) xx(N) xx(1)
0
(44)
注意，从AR模型得到的式(44)与按最大熵外推 xx(N 1)得到的式（43)完全相同，这就证明了当x(n)为高斯分布时最大熵谱估计法与AR模型法是等价的。
AR模型谱估计的性质： 1. AR谱的平滑特性
由于AR模型是一个有理分式，因而估计出的谱要比经典法的谱平滑。
2. AR谱的分辨率
由信号的时宽-宽带积可知，长度为N的信号，若抽样间隔为Ts，那么由DFT作谱分析时，其分辨率粗略地为fs/N。总之，经典谱估计的分辨率反比于使用的信号的长度。现代谱估计的分辨率可以不受此限制，这是因为，对给定的数据x（n），n=0,1，。。。，N-1，虽然其估计出的自相关函数也是有限长，即m=-（N-1）~（N-1），但现代谱估计的一些方法隐含着数据和自相关函数的外推，使其可能的长度超过给定的长度。例如，AR模型是在最小均方意义上对给定的数据的拟合，即：
(n) akx(n k) （23） x
k 1
p
(n)可能达到的长度是从0至（N-1+p） (n)代替x（n）这样x。在此之外，若用x，还可
以继续外推。
3. AR谱的匹配性质
jw
记Pe(e)为误差序列e（n）的功率谱，则
Px(e) PAR(ejw)
jw
Pe(ejw)A(e)
jw
2
（45）
min
A(e)
jw
2
（46）
比较上面两式可以看出，如果信号x（n）的功率谱Px(ejw)被AR模型的谱PAR(ejw)来
jw模拟，则误差序列e（n）的功率谱P，e(e)将由一个常数 min来模拟。具有这样谱的e（n）
对确定性信号来说，是一个冲激函数，对随机信号来说，将是一个白噪声。从这里也可以看出A(z)的白化作用。
记e（n）的自相关函数为re(m)，则
1
E{e(n)} re(0)
2
21
P(e)dw x
2
jw
Px(ejw) （47）
将（45）式带入（47）式得到有：
min
2
Px(ejw)
dw （48）
PAR(ejw)
式中p是误差序列e（n）的均方值，若令上式中的p最小，同样可得到正则方程。这样，我们对AR模型或线性预测，又可以在频域给出新的解释：给定一个随机信号x（n）的功率谱，我们希望用PAR(e)来模拟它，使两者比值的积分为最小。由前面可知p的最小值
jw
为pmin，所以（48）式变为：
12
Px(ejw)
dw 1 （49） jw
PAR(e)
由（49）式，我们可以得出PAR(ejw)和Px(ejw)相匹配的一些性质。从整体上看PAR(ejw)将均匀地和Px(ejw)相跟随。因为均值为1，所以PAR(ejw)将在Px(ejw)的上下波动，即在有的区域PAR(ejw)>Px(ejw)，而在另外的区域PAR(ejw)<Px(ejw)。但是，从对（49）式积分的贡献来说，PAR(ejw)/Px(ejw)<1的贡献要比PAR(ejw)/Px(ejw)>1的贡献大。总的一句话，PAR(ejw)是Px(ejw)的包络的一个好近似。从AR模型自身的特点也可以理解这一点。因为H（z）是一个全极点的模型，因而它易于表示谱峰，而不易于表示谱谷。总之，在整个频率范围内，PAR(ejw)和Px(ejw)相跟随，但在每一个局部处，它跟随Px(ejw)的峰点要比跟随谷点的程度好。
4. AR谱的统计特性
严格地分析AR谱的方差比较困难，目前尚未有一个解析表达式。粗略地讲，AR谱的方差反比于数据x（n）的长度N和信噪比SNR。
5. AR模型谱估计方法的不足
其一：AR谱的分辨率和求AR模型时所使用信号的信噪比SNR有着密切的关系。由于
2
零点的存在使谱的动态范围减小，从而降低了分辨率。 m越大，即y（n）的信噪比SNR
越小，谱的分辨率降低得越明显。
其二：如果x（n）是含有噪声的正弦信号，在应用时发现，谱峰的位置易受x（n）的初相位的影响，且在有的算法中还可能出现“谱线分裂”的现象，即在本来应只有一个谱线的位置分裂出两个谱线。
其三：谱估计的质量受到阶次p的影响。P选的过低，谱太平滑，反应不出谱峰；p选的过大，可能会产生虚假的峰值。
例题12.11: 一．
求出信号x（n）和真实功率谱Px(ejw)
1. 利用rand产生一组白噪声数据，题目中要求的白噪声长度是N=256，为了保证产生
的数据较好地近似白噪声，数据的长度需要取得长一些，此处取4096。程序取
2=0.12。
2. 将噪声u（n）通过系统H（z），得到输出v（n）；则v（n）的真实功率谱可以求出，
即Pv
2
2
H(ejw)
3. 求出三个正弦信号。已知三个正弦信号的归一化频率分别是f1=0.1，f2=0.25，
f3=0.26，其幅度可由下式求出 SNRi 10log10(Ai/2 ),i 1,2,3，由此式分别求出A1，A2，A3
2
2
4. 求出所要的信号x（n） x(n) v(n)
Asin(2 fn)
i
i
i 1
jw
2
3
3
5．求出真实功率谱，正弦信号的功率谱是两条垂直线，所以可以直接叠加
Px(e) H(e)
jw
2
i 1
2
Ai2[ (f f') (f f')]
上图为信号x（n）的真实功率谱
二．1.自相关法求解AR模型估计其功率谱，模型阶次分别为
p=8,11,30
p=8
p=11
p=30
2.用Burg方法重复上述步骤，p=5,7,9
p=5
p=7
p=9
比较自相关法和Burg法，明显看出，无论是在谱的分辨率方面还是在模型的阶次选择方面，Burg算法都有明显的优点。在阶次为8和11时，自相关不能将f=0.25和f=0.26分开，直到30后，才稍微将两个谱峰分开。相比之下，Burg算法在阶次等于5的时候就将两个谱峰明显分开，具有明显的优点。 附程序： clear all; close all;
%题目给出的参数------------------------------------------------------- N=256;
P=0.12;
h=[1 -0.1 0.09 0.648]; %FIR滤波器系数
SNR=[10 50 50]; F=[0.1 0.25 0.26]; %三个信噪比和归一化频率
%步骤一：产生白噪声----------------------------------------------- a=sqrt(P);
u=rand(1,20*N); u=u-mean(u);u=u*a;
%步骤二：滤波处理----------------------------------------------------- v1=filter(h,1,u);
v=v1(1000:1000+N-1);
%步骤三：产生正弦序列---------------------------------------------------- A=sqrt(2*P*10.^(SNR/10)) %根据信噪比确定三个正弦函数的幅度 x1=A(1)*sin(2*pi*F(1)*[0:N-1]); x2=A(2)*sin(2*pi*F(2)*[0:N-1]); x3=A(3)*sin(2*pi*F(3)*[0:N-1]);
%步骤四：叠加信号，得到信号x--------------------------------------------- x=v+x1+x2+x3; save testdat x;
%步骤五：求出信号x的真实功率谱--------------------------------------------- M=4096;H=freqz(h,1,M,1);H=abs(H);H=H.*H;H=P*H; f=0:1/M:0.5-1/M;
power_y1=pi*A(1)*A(1)/2; power_y2=pi*A(2)*A(2)/2; power_y3=pi*A(3)*A(3)/2; f1=0.1;f2=0.25;f3=0.26; f1_n=round(f1*M); f2_n=round(f2*M); f3_n=round(f3*M);
H(f1_n)=H(f1_n)+power_y1; H(f2_n)=H(f2_n)+power_y2; H(f3_n)=H(f3_n)+power_y3; PSD=10*log10(H/max(H)+eps);
figure('color','w');plot(f,PSD(1:2048));grid on;
load testdat; N=2048;
% 使用自相关法得到功率谱估计； [xpsd,F]=pyulear(x,8,N,1); pmax=max(xpsd); xpsd=xpsd/pmax;
xpsd=10*log10(xpsd+0.000001); figure('color','w');subplot(321);
plot(F,xpsd);grid on; %title('pyulear 8');
[xpsd,F]=pyulear(x,11,N,1); pmax=max(xpsd); xpsd=xpsd/pmax;
xpsd=10*log10(xpsd+0.000001); subplot(323);
plot(F,xpsd);grid on; %title('pyulear 11');
[xpsd,F]=pyulear(x,30,N,1); pmax=max(xpsd); xpsd=xpsd/pmax;
xpsd=10*log10(xpsd+0.000001); subplot(325);
plot(F,xpsd);grid on; % title('pyulear 30');
[xpsd,F]=pburg(x,5,N,1); pmax=max(xpsd); xpsd=xpsd/pmax;
xpsd=10*log10(xpsd+0.000001); subplot(322);
plot(F,xpsd);grid on; % title('pburg 5');
[xpsd,F]=pburg(x,7,N,1); pmax=max(xpsd); xpsd=xpsd/pmax;
xpsd=10*log10(xpsd+0.000001); subplot(324);
plot(F,xpsd);grid on; % title('pburg 7');
[xpsd,F]=pburg(x,9,N,1); pmax=max(xpsd); xpsd=xpsd/pmax;
xpsd=10*log10(xpsd+0.000001); subplot(326);
plot(F,xpsd);grid on; %title('pburg 9'); ...
