摘要：信号参数估计是现代信号处理的重要研究内容之一，在频域中进行傅里叶变换研究信号，是研究确定性信号最简单且有效的手段，但在现代信号分析中，对于常见的随机信号，不可能用清楚的数学关系式来描述，其傅里叶变换更不存在，转而可以利用给定的若干个样本数据估计来估计信号的参数。本学期在导师的指导下我学习了这门课程，了解到相关的知
识，深刻体会了信号参数估计的理论基础。本文主要介绍我对信号参数估计中的现代谱估计的理解和有关体会。
关键字：参数估计；随机信号；谱估计
引言：
功率谱估计是随机信号处理的重要内容，其技术渊源很长，而且在过去的40余年中
获得了飞速的发展。涉及到信号与系统、随机信号分析、概率统计、矩阵代数等一系列的基础学科，广泛应用于人民的日常生活及军事、工业、农业活动中，是一个具有强大生命力的研究领域。现代谱估计的方法又大致可分为参数模型谱估计和非参数模型谱估计，前者有AR 模型、MA 模型、ARMA 模型、PRONY 模型等，后者有最小方差方法、多分量的MUSIC 方法等。
一 现代谱估计方法的发展 1.1功率谱研究的发展过程
功率谱估计是数字信号处理的主要内容之一，主要研究信号在频域中的各种特征，目的是根据有限数据在频域内提取被淹没在噪声中的有用信号。
现代谱估计的提出主要是针对经典谱估计（周期图和自相关法）的分辨率和方差性能不好的问题。1967 年，Burg 提出的最大嫡谱估计，即是朝着高分辨率谱估计所作的最有意义的努力。虽然，Bartlett 在 1948年，Parzem 于 1957 年都曾经建议用自回归模型做谱估计，但在 Burg 的论文发表之前，都没有引起注意。
现代谱估计的内容极其丰富，涉及的学科及应用领域也相当广泛，至今，每年都有大量的论文出现。非参数模型谱估计的特点是其模型不是用有限参数来描述，而直接由相关函数序列得到，这种方法能提高低信噪比时的谱分辨率。参数模型谱估计是先根据过程的先验信息或者一些假定，建立一个数学模型来表示所给定采样数据的过程，或者选择一个较好的近似实际模型，而后利用采样数据序列或者自相关序列，估计该模型的参数，最后把参数代入到该模型对应的理论功率谱表达式，得到所需要的谱估计。
1.2 功率谱估计应用及用途
功率谱估计有着极其广泛的应用，不仅在认识一个随机信号时，需要估计它的功率谱。它还被广泛地应用于各种信号处理中。在信号处理的许多场所，要求预先知道信号的功率谱密度(或自相关函数) 。
功率谱估计就是通过信号的相关性估计出接受到信号的功率随频率的变化关系，实际用途有滤波，信号识别（分析出信号的频率），信号分离，系统辨识等。谱估计技术是现代信号处理的一个重要部分，还包括空间谱估计，高阶谱估计等。维纳滤波、卡尔曼滤波，可用于自适应滤波，信号波形预测等（火控系统中的飞机航迹预判）。
二．现代谱估计
现代功率谱估计即参数谱估计方法是通过观测数据估计参数模型再按照求参数模型输出功率的方法估计信号功率谱, 主要是针对经典谱估计的分辨率低和方差性能不好等问题提出的。常用模型有 ARMA 模型、 AR 模型、 MA 模型,AR 模型的正则方程是一组线性方程，而MA 和ARMA 模型是非线性方程。由于AR 模型具有一系列良好的性能，因此被研究最多也得到最广泛的应用。
2.1 AR模型
AR 模型又称为自回归模型, 是一个全极点的模型, 可用如下差分方程来表示:
p
x (n ) =-∑a k x (n -k ) +w (n ) （2-1）
k =1
其中，p 是 A R模型的阶数， {a k }=l ， 2 ， …， p 是p 阶 A R模型的参数．将该模型记为AR( p) ，它的系统转移函数为
H (z ) =
X (z )
=W (z )
11+∑a k z
k =1p
-k
(2-2)
在功率谱估计中，若观测的数据 x(n) 是平稳随机过程，则该系统的输入w( n ) 也可认为是平稳的，因 而根据线性系统对平稳随机信号的响应理论可得观测数据的功率谱为
2
(w ) =|H () |σw e =P x
2
jw
p
2w
-j k w 2
(2-3)
|1+∑a k e
k =1
|
由式可知，利用A R模型进行功率谱估计的实质是求解模型系数 {a k } 和σ
2w
的
问 题．
将式 ( 1 )两端乘以x ( n-m ) 求平均 ( 数学期望) ，可以求得观测数据的A R( p) 模型参数与自相关函数的关系式为
p
⎧
⎪-∑a k R xx (m -k )
m >0⎪p k =1
2⎪
(m ) =-(m -k ) +m =0 (2-4) ⎨∑a k R xx σR xx w
⎪k =1
m
⎪R xx (-m ) ⎪⎩
可见， 阶 AR 模型输出的相关函数具有递推的性质， 因而选用 AR 模型进行谱估计只需较少的观测数据将式 ( 4 )写成矩阵形式得
R (1) R (p ) ⎤⎡1⎤⎡σ2⎤⎡R (0)
⎢⎥⎢w ⎥⎢R (1) ⎥R (0) R (p -1) ⎥⎢a 1⎥⎢0⎥⎢ (2-5) =
⎢ ⎥⎢ ⎥⎢ ⎥⎢⎥⎢⎥⎢⎥R (p ) R (p -1) R (0) ⎣⎦⎢⎣0⎥⎦⎣a p ⎥⎦⎢ 上式就是著名的Yule-Walker( Y—W) 方程．它表明，只要已知观测数据的自相关数，就能求出AR 模型参数{a k } 和
σ
2w
，进而按式 ( 3 )求得信号功率谱
的估值。
另外，从 AR 模型的差分方程式可知，该模型的现在输出值是它本身过去值的回归，这与预测器存在着一定的相似性，它们之间有着非常密切的关系，即它们的系统函数互为倒数，也就是说预测误差滤波器
A (z ) 就是 AR 信号模
p
型H(z)的逆滤波器．因此通过预测误差滤波器优化设计使预测均方误差最小就可求得A R信号模型的最优参数，即P 阶线性预测器的预测系数{a p (k ) }等于p 阶A R模型的系数{a k }，其最小均方预测误差E p 等于自噪声方差
σ
2w
。
因此， 根据上述的Y-W 方程以及 AR 模型与预测误差滤波器之间的关系，就可提取 A R模型参数．目前主要有三种：Levinson-Durbin 算法、Burg 算法和Marple 算法。这三种方法都可以用时间平均代替集合平均的最小平方准则推导得到。在本文中，我们主要采用TBurg 法来估计信号的模型参数，因此主要介绍一下Burg 法。
Levinson-Durbin 算法
L -D 递推算法是在满足前向预测均方误差最小的前提下，先求得观测数据的自相关函数，然后利用Y - W方程的递推性质求得模型参数，进而根据式 ( 3 )求得功率谱的估值．它是模型阶次逐次加大的一种算法，即先计算阶次m=l 时的预测系数{a m (k ) }=
a 1(1) 和σw 1，再计算m= 2 时的a 2(1) ，a 2(2) 和σw 2 按
2
2
22
此依次计算到 阶次m=p时的a p (1) ，a p (2) ， …a p (p ) 及σwp 当σp 满足精度要求时即可停止递推．递推公式为：
E
a (k ) =a (k ) +a (m ) a
m
m -1
m
a m (m ) =-
R (m ) +∑a m -k (k ) R (m -k )
k =1
m -1
(2-6)
m -1
m -1
(m -k ) k=1,2,……m-1 (2-7)
其中
2
⎡1-⎤⎡1-⎤ （2-8） ==(m ) =R (0) (k ) ∏a a E m σwm ⎢E m -1m k ⎥⎢⎥⎣⎦⎣⎦k =1
2
2
m
Burg 算法
Burg 算法的基本思想是直接从观测的数据利用线性预测器的前向和后向预测的总均方误差之和为最小的准则来估计反射系数，进而通过 L -D算法的递推公式求出 AR 模型优化参数．设观察到的N 个数据为x( 0 ),x(1)， …， x ( N—1 ) ，则具体算法如下： (1)取m=1,初始化：e 0(n ) =
2
f
e (n ) =x(n), n=0,1……N-1
b
1N -12
=R (0) =x (n ) ∑σw 0
N n =0
(2)计算反射系数：ρ=
m
n =m
∑e
N -1
-2∑e m -1(n ) e m -1(n -1)
n =m f m -1
N -t
f b
(n ) +
e
2
b m -1
(n )
2
（3）计算滤波器系数及预测误差功率：a m (m ) =ρm k=1,2,……m-1
a
m
(k ) =a m -1(k ) +ρ
2
m
a
m -1
(m -k )
E m =σwm =(1-ρ) E m -1
m
f
f
m
2
(4)递推高一阶前、后向预测误差：e m (n ) =e m -1(n ) +ρ
2
e
b m -1
(n -1)
f m -1
e m (n ) =e m -1(n -1) +ρ
b
b
m
e
(n )
把m 更新为m+1 ，重复②～④直到σp 满足要求。
Marple 算法
Marple 算法又称为不受约束的最小二乘法，它的主要思路是为了摆脱因采用递推运算对确定预测系数的约束，让每一预测系数 ( 模型参数) 的确定直接与前、
后向预测的总的平方误差最小 ( 最小二乘法) 联系起来．即令总的平方误差
εp =∑
n =p
N -1
{[e (n ) ]+[e (n ) ]}
f
2
b
2
m -1
m -1
2
2
⎡p ⎤N -1⎡p ⎤
=∑⎢∑a p (k ) x (n -k ) ⎥+∑⎢∑a p (k ) x (n +k -p ) ⎥
n =p ⎣k =0⎦n =p ⎣k =0⎦
N -1
a
p
(0) =1
最小．由上式可见，总的平方误差εp 是系数a p (k ) 的函数．若把εp 对各预测系数a p (k ) 而非单一地对a p (p ) =ρ求导数，并令其为零，就可以得到一组线性方
p
程． 解此方程组所得的a p (k ) 就是在最小平方误差准则下的最优预测系数．但由于方程组系数矩阵不是Toeplitz 型，所以不能利用 L — D 算法求解．为了减
少运算量，Ma r p l e 提出一种格型结构的高效递推算法。
2.1.1 AR 模型的稳定性及阶的确定
AR(p)模型稳定的充分必要条件是H(z)的极点(即A(z)的根) 都在单位圆内。稳定的AR(p)模型将具有以下的性质：
(1)H(z)的全部极点或A(z)的所有根都在单位圆内。 (2)自相关矩阵是正定的。
(3)激励信号的方差(能量) 随阶次的增加而递。 (4)反射系数的模恒小于1。
但是在实际应用中，levinson 算法的已知数据(自相关值) 是由 来估计的，有限字长效应有可能造成大的误差，致使估计出来的AR(p)参数所构成的A(z)的根跑到单圆上或单位圆外，从而使模型失去稳定。在递推计算的过程中如果出现这种情况，将致，或|即停止递推计算。
通常事先并不知道AR 模型的阶。阶选得太低，功率谱受到的平滑太厉害，平滑后的谱已经分辨不出真实谱中的两个峰了。阶选的太高，固然会提高谱估计的分辨率，但同时会产生虚假的谱峰或谱的细节。
一种简单而直观的确定AR 模型的阶的方法，是不断增加模型的阶，同时观察预测误差功率，当其下降到最小时，对应的阶便可选定为模型的阶。但是预测误差功率(或AR 模型激励源的方差) 是随着阶次的增加而单调下降的，因此，很难确定降到什么程度才合适。另一方面，应该注意到，随着阶次增加，模型参数的数目亦增多了，谱估计的方差会变大(表现在虚假谱峰的出现) 。因此，不能简
单地依靠观察预测误差功率的下降来确定模型的阶。与此相应的另一种简单方法是观察各阶模型预测误差序列的周期图，当它很接近于平坦(白色谱) 时即对应最佳的阶。
2.2 MA模型及功率谱估计
MA模型是一个全零点模型，其功率谱不易体现信号中的峰值，即分辨率较低。从谱估计的角度来看，MA 模型谱估计等效于经典谱估计中的自相关法，若单纯为了对一段有限长数据做谱估计，就没必要求解MA 模型了。但在系统分析于辨识以及ARMA 谱估计中都要用到MA 模型，因此仍有必要讨论MA 系数的求解方法。目前提出的方法主要有：1，谱分解法；2用高阶的AR 模型来近似MA 模型；3，最大似然法。以第二种方法为例讨论MA 模型参数的提取。
ⅰ有N 点数据x （n ）建立一个p 阶的AR 模型，p>>q，可用AR 模型参数的计算方法求出p 阶的AR 系数a p (k ) ，k=1,2…p； ⅱ利用a
(k )
p
，k=1,2…p进行线性预测，等效于一个q 阶的AR 模型，
再一次利用AR 系数的求解方法得到b （k ），k=1,2，…，q 。 进而通过两次求AR 系数可得MA 模型的系数。
2.3ARMA 模型
理论上可以证明，在一定条件下，A R MA和 M A模型都可以用一个阶次足够大的A R模型来近似，所以在实际使用中我们用一个阶次比较高的A R模型代替A R MA模型来进行功率谱估计，避免了求解非线性方程组带来的困难。
3基于AR 模型的噪声源识别
3.1噪声源识别的方法： 3.1.1 消去法
所谓消去法，就是首先把附带全部装备的试验对象(如汽车、发动机等) 在一定条件下测定其总和的工作噪声，然后去除其中的一部分装备或控制这部分噪声传出，再按同样试验条件测定去除部分装备后的试验对象的工作噪声，则去除这部分装备前后试验对象噪声交化值，即被视作被拆卸或控制传出装备噪声的方法。这一方法试验难度大，耗时耗力，由于对部分设备要求停机而不能进行在线测量，因而不能获得实验对象噪声分布的整体概念.
3.1.2 声强测量法
声强测量法是利用双点声压梯度的积分来近似空气质点的振动速度，并利用FFT 来实现声强的实时测量。采用声强测量法分析噪声可获得丰富的噪声信息，不但可测出噪声的大小，而且还能测出声能的流向。该方法对环境要求低。现场识别能力强．但是很强的背景声会对声强测试产生不良的影响，产生较大的误差，并且其测试精度较难控制。
3.1.3 相干函数法
相干函数法是根据相关理论，用互相关函数描述被测部件的振动信号与外界噪声信号之间的关系．根据结果判断噪声是否由该部件振动引起，从而判断出各个声源．在实际噪声测量中，常常是复杂的多输入单输出系统。而且输入噪声源又可能是相互依赖的，因此这种方法也不适合判断多输入单输出的系统
3.1.4 频谱分析法
频谱分析法，由于声源噪声的形成机理不同。使每个声源的嗓声特点有差异，能量分布的频谱范围有的位于全频带，有的位于中、低频或高频。因而，在了解各组成声源噪声频谱特点的情况下，测量出总噪声频谱。再取得各部分声源的频谱，并将其与总声源频谱比较，就可以分析出各组成部分噪声贡献幅度。从而找出需控制的主要声源，以有效降低高能频率段的噪声。频谱分析法中，除应用幅值频谱外，还常应用功率谱。但是由于采用了传统的功率谱估计，所以存在频率分辨率低，旁瓣泄漏严重的缺点。
由于现代谱估计方法为信号建立了一个准确或近似的模型，从根本上摒弃了对数据序列加窗的隐含假设，因而能够更好地估计信号功率谱．利用现代谱方法估计功率谱时，不要求噪声源中只有一个信号为高斯分布，这与盲信号分离有很大区别(盲信号分离要求源信号中只能有一个信号为高斯分布) 。只要利用能量随频率的分布(即功率谱) 来求取噪声源的比重。现代谱估计的上述优点为噪声源识别提供了更为精确的方法。
3.2研究内容
由于缺少实际的噪声源，在此采用仿真的噪声信号来模拟噪声。我们知道AR 模型的参数与其阶次有关。阶选得太低，功率谱受到的平滑太厉害，平滑后的谱已经分辨不出真实谱中的两个峰了。阶选的太高，固然会提高谱估计的分辨率，但同时会产生虚假的谱峰或谱的细节。所以必须选择合适的AR 模型的阶次。本文采用选用Akaike(alC)准则作为确定模型阶的依据。将信号x1,x2按不同的比例混合在一起组成混合噪声以模仿生活中常见的噪声信号。
3.2.1仿真信号的选取 N=512;
fs=1; a1=5; a2=3; a3=2; w=2*pi/fs;
x1=a1*sin(0.5*w*(0:N-1))+a2*cos(w*(0:N-1)).^2+randn(1,N); x2=a2*sin(w*(0:N-1))+a2*cos(0.8*w*(0:N-1)).^2+randn(1,N); y=1:512; plot(y,x1); figure(2); plot(y,x2); save x1 x1; save x2 x2
;
首先判断AR 模型的阶数。这里应用的就是A1C 准则。我们用仿真信号x1为例进行了阶数测定，与模型阶的关系曲线如图所示。从图中可以发现模型的AIC 存在最小值，此时信号模型对应的阶数为12。所以选取12阶的模型。 clear;
load x1.mat ; x=x1;
N=length(x); M=20;
e1=zeros(M,N); e2=zeros(M,N); e1(1,1:N)=x; e2(1,1:N)=x; r=zeros(1,M); a=zeros(1,M); q1=zeros(1,M); q2=zeros(1,M); AIC=zeros(1,M);
a(1)=(1/N)*sum(x.^2);
q1(1)=sum(e1(1,2:N).*e2(1,1:N-1));
q2(1)=sum((e1(1,2:N).^2)+(e2(1,1:N-1).^2)); r(1)=2*q1(1)/q2(1); AIC(1)=log(a(1)); for k=2:M for n=k:N
e1(k,n)=e1(k-1,n)+((-1)*r(k-1)*e2(k-1,n-1)); e2(k,n)=(-1)*r(k-1)*e1(k-1,n)+e2(k-1,n-1); end
q1(k)=sum(e1(k,k+1:N).*e2(k,k:N-1));
q2(k)=sum((e1(k,k+1:N).^2)+(e2(k,k:N-1).^2)); r(k)=2*q1(k)/q2(k);
a(k)=(1-r(k-1)^2)*a(k-1); AIC(k)=log(a(k))+k*log(N)/N; end
[X,Y]=min(AIC); k=1:M;
plot(k,AIC); [X,Y]
3.2.3功率谱估计及噪声源识别
y1分别取10x1+x2，x1+10x2，x1+x2;我们尝试着去识别主噪声源
在进行信号混合时利用功率谱估计出哪个信号是主噪声源，求功率谱的方法如下：function argonglv(y1,y2,y3)
clc
clear
load x1.mat
load x2.mat
load x3.mat
y1=10*x1+x2;
k1=sqrt(sum(y1.^2));
y1=y1./k1
p1= arburg(y1,12);
la=length(p1);
ln=128-la;
ts=1/5000;
ab=zeros(1,ln);
bb=[p1 ab]
dw=1/128;
w=dw*[0:127]
f=fft(bb)
f1=abs(f);
f2=f1.^2;
f3=1./f2;
f4=10*log10(f3);
plot(w(1:128),f4(1:128),'-kX' )
hold on
figure (1);
y1=x1;
k1=sqrt(sum(y1.^2));
y1=y1./k1
p1= arburg(y1,12);
la=length(p1);
ln=128-la;
ts=1/5000;
ab=zeros(1,ln);
bb=[p1 ab]
dw=1/128;
w=dw*[0:127]
f=fft(bb)
f1=abs(f);
f2=f1.^2;
f3=1./f2
f4=10*log10(f3);
plot(w(1:128),f4(1:128),'-r*')
hold on
figure(1);
y1=x2;
k1=sqrt(sum(y1.^2));
y1=y1./k1
p1= arburg(y1,12);
la=length(p1);
ln=128-la;
ts=1/5000;
ab=zeros(1,ln);
bb=[p1 ab]
dw=1/128;
w=dw*[0:127]
f=fft(bb)
f1=abs(f);
f2=f1.^2;
f3=1./f2
f4=10*log10(f3);
plot(w(1:128),f4(1:128),':+')
legend('y1' , 'x1' , 'x2' );
(1) 当源信号的比重大于1时，我们可以明显地看出混合信号的功率谱与x1的功
率谱相接近，也就是混合信号以x1信号为主噪声源如下图示
(2)当源信号的比重小于1时，我们可以明显地看出混合信号的功率谱与x2相接近，也就是混合信号以源信号x2为主噪声源，功率谱如下图所示
(2) 当源信号x1与x2以l ：l 的比例混合时，我们难以确定哪个是主要源信号。
由上面的讨论知：我们可以利用现代谱估计的理论方法是可以对主要噪声源信号进行识别的，也促使我们利用真实的噪声信号进行进一步的研究。只要选用合适的参数模型阶数，求出信号的模型参数和功率谱图，然后利用曲线相似度和曲线关联度进一步求取噪声源信号在混合信号中的比重，即解决信号功率谱图中曲线的相似度问题。现代谱估计方法的噪声源识别作为信号处理方法的一种，克服了仪器测量和环境的影响，同时克服了传统功率谱估计的缺点，因而在噪声识别的仿真实验中取得了较好的结果。
摘要：信号参数估计是现代信号处理的重要研究内容之一，在频域中进行傅里叶变换研究信号，是研究确定性信号最简单且有效的手段，但在现代信号分析中，对于常见的随机信号，不可能用清楚的数学关系式来描述，其傅里叶变换更不存在，转而可以利用给定的若干个样本数据估计来估计信号的参数。本学期在导师的指导下我学习了这门课程，了解到相关的知
识，深刻体会了信号参数估计的理论基础。本文主要介绍我对信号参数估计中的现代谱估计的理解和有关体会。
关键字：参数估计；随机信号；谱估计
引言：
功率谱估计是随机信号处理的重要内容，其技术渊源很长，而且在过去的40余年中
获得了飞速的发展。涉及到信号与系统、随机信号分析、概率统计、矩阵代数等一系列的基础学科，广泛应用于人民的日常生活及军事、工业、农业活动中，是一个具有强大生命力的研究领域。现代谱估计的方法又大致可分为参数模型谱估计和非参数模型谱估计，前者有AR 模型、MA 模型、ARMA 模型、PRONY 模型等，后者有最小方差方法、多分量的MUSIC 方法等。
一 现代谱估计方法的发展 1.1功率谱研究的发展过程
功率谱估计是数字信号处理的主要内容之一，主要研究信号在频域中的各种特征，目的是根据有限数据在频域内提取被淹没在噪声中的有用信号。
现代谱估计的提出主要是针对经典谱估计（周期图和自相关法）的分辨率和方差性能不好的问题。1967 年，Burg 提出的最大嫡谱估计，即是朝着高分辨率谱估计所作的最有意义的努力。虽然，Bartlett 在 1948年，Parzem 于 1957 年都曾经建议用自回归模型做谱估计，但在 Burg 的论文发表之前，都没有引起注意。
现代谱估计的内容极其丰富，涉及的学科及应用领域也相当广泛，至今，每年都有大量的论文出现。非参数模型谱估计的特点是其模型不是用有限参数来描述，而直接由相关函数序列得到，这种方法能提高低信噪比时的谱分辨率。参数模型谱估计是先根据过程的先验信息或者一些假定，建立一个数学模型来表示所给定采样数据的过程，或者选择一个较好的近似实际模型，而后利用采样数据序列或者自相关序列，估计该模型的参数，最后把参数代入到该模型对应的理论功率谱表达式，得到所需要的谱估计。
1.2 功率谱估计应用及用途
功率谱估计有着极其广泛的应用，不仅在认识一个随机信号时，需要估计它的功率谱。它还被广泛地应用于各种信号处理中。在信号处理的许多场所，要求预先知道信号的功率谱密度(或自相关函数) 。
功率谱估计就是通过信号的相关性估计出接受到信号的功率随频率的变化关系，实际用途有滤波，信号识别（分析出信号的频率），信号分离，系统辨识等。谱估计技术是现代信号处理的一个重要部分，还包括空间谱估计，高阶谱估计等。维纳滤波、卡尔曼滤波，可用于自适应滤波，信号波形预测等（火控系统中的飞机航迹预判）。
二．现代谱估计
现代功率谱估计即参数谱估计方法是通过观测数据估计参数模型再按照求参数模型输出功率的方法估计信号功率谱, 主要是针对经典谱估计的分辨率低和方差性能不好等问题提出的。常用模型有 ARMA 模型、 AR 模型、 MA 模型,AR 模型的正则方程是一组线性方程，而MA 和ARMA 模型是非线性方程。由于AR 模型具有一系列良好的性能，因此被研究最多也得到最广泛的应用。
2.1 AR模型
AR 模型又称为自回归模型, 是一个全极点的模型, 可用如下差分方程来表示:
p
x (n ) =-∑a k x (n -k ) +w (n ) （2-1）
k =1
其中，p 是 A R模型的阶数， {a k }=l ， 2 ， …， p 是p 阶 A R模型的参数．将该模型记为AR( p) ，它的系统转移函数为
H (z ) =
X (z )
=W (z )
11+∑a k z
k =1p
-k
(2-2)
在功率谱估计中，若观测的数据 x(n) 是平稳随机过程，则该系统的输入w( n ) 也可认为是平稳的，因 而根据线性系统对平稳随机信号的响应理论可得观测数据的功率谱为
2
(w ) =|H () |σw e =P x
2
jw
p
2w
-j k w 2
(2-3)
|1+∑a k e
k =1
|
由式可知，利用A R模型进行功率谱估计的实质是求解模型系数 {a k } 和σ
2w
的
问 题．
将式 ( 1 )两端乘以x ( n-m ) 求平均 ( 数学期望) ，可以求得观测数据的A R( p) 模型参数与自相关函数的关系式为
p
⎧
⎪-∑a k R xx (m -k )
m >0⎪p k =1
2⎪
(m ) =-(m -k ) +m =0 (2-4) ⎨∑a k R xx σR xx w
⎪k =1
m
⎪R xx (-m ) ⎪⎩
可见， 阶 AR 模型输出的相关函数具有递推的性质， 因而选用 AR 模型进行谱估计只需较少的观测数据将式 ( 4 )写成矩阵形式得
R (1) R (p ) ⎤⎡1⎤⎡σ2⎤⎡R (0)
⎢⎥⎢w ⎥⎢R (1) ⎥R (0) R (p -1) ⎥⎢a 1⎥⎢0⎥⎢ (2-5) =
⎢ ⎥⎢ ⎥⎢ ⎥⎢⎥⎢⎥⎢⎥R (p ) R (p -1) R (0) ⎣⎦⎢⎣0⎥⎦⎣a p ⎥⎦⎢ 上式就是著名的Yule-Walker( Y—W) 方程．它表明，只要已知观测数据的自相关数，就能求出AR 模型参数{a k } 和
σ
2w
，进而按式 ( 3 )求得信号功率谱
的估值。
另外，从 AR 模型的差分方程式可知，该模型的现在输出值是它本身过去值的回归，这与预测器存在着一定的相似性，它们之间有着非常密切的关系，即它们的系统函数互为倒数，也就是说预测误差滤波器
A (z ) 就是 AR 信号模
p
型H(z)的逆滤波器．因此通过预测误差滤波器优化设计使预测均方误差最小就可求得A R信号模型的最优参数，即P 阶线性预测器的预测系数{a p (k ) }等于p 阶A R模型的系数{a k }，其最小均方预测误差E p 等于自噪声方差
σ
2w
。
因此， 根据上述的Y-W 方程以及 AR 模型与预测误差滤波器之间的关系，就可提取 A R模型参数．目前主要有三种：Levinson-Durbin 算法、Burg 算法和Marple 算法。这三种方法都可以用时间平均代替集合平均的最小平方准则推导得到。在本文中，我们主要采用TBurg 法来估计信号的模型参数，因此主要介绍一下Burg 法。
Levinson-Durbin 算法
L -D 递推算法是在满足前向预测均方误差最小的前提下，先求得观测数据的自相关函数，然后利用Y - W方程的递推性质求得模型参数，进而根据式 ( 3 )求得功率谱的估值．它是模型阶次逐次加大的一种算法，即先计算阶次m=l 时的预测系数{a m (k ) }=
a 1(1) 和σw 1，再计算m= 2 时的a 2(1) ，a 2(2) 和σw 2 按
2
2
22
此依次计算到 阶次m=p时的a p (1) ，a p (2) ， …a p (p ) 及σwp 当σp 满足精度要求时即可停止递推．递推公式为：
E
a (k ) =a (k ) +a (m ) a
m
m -1
m
a m (m ) =-
R (m ) +∑a m -k (k ) R (m -k )
k =1
m -1
(2-6)
m -1
m -1
(m -k ) k=1,2,……m-1 (2-7)
其中
2
⎡1-⎤⎡1-⎤ （2-8） ==(m ) =R (0) (k ) ∏a a E m σwm ⎢E m -1m k ⎥⎢⎥⎣⎦⎣⎦k =1
2
2
m
Burg 算法
Burg 算法的基本思想是直接从观测的数据利用线性预测器的前向和后向预测的总均方误差之和为最小的准则来估计反射系数，进而通过 L -D算法的递推公式求出 AR 模型优化参数．设观察到的N 个数据为x( 0 ),x(1)， …， x ( N—1 ) ，则具体算法如下： (1)取m=1,初始化：e 0(n ) =
2
f
e (n ) =x(n), n=0,1……N-1
b
1N -12
=R (0) =x (n ) ∑σw 0
N n =0
(2)计算反射系数：ρ=
m
n =m
∑e
N -1
-2∑e m -1(n ) e m -1(n -1)
n =m f m -1
N -t
f b
(n ) +
e
2
b m -1
(n )
2
（3）计算滤波器系数及预测误差功率：a m (m ) =ρm k=1,2,……m-1
a
m
(k ) =a m -1(k ) +ρ
2
m
a
m -1
(m -k )
E m =σwm =(1-ρ) E m -1
m
f
f
m
2
(4)递推高一阶前、后向预测误差：e m (n ) =e m -1(n ) +ρ
2
e
b m -1
(n -1)
f m -1
e m (n ) =e m -1(n -1) +ρ
b
b
m
e
(n )
把m 更新为m+1 ，重复②～④直到σp 满足要求。
Marple 算法
Marple 算法又称为不受约束的最小二乘法，它的主要思路是为了摆脱因采用递推运算对确定预测系数的约束，让每一预测系数 ( 模型参数) 的确定直接与前、
后向预测的总的平方误差最小 ( 最小二乘法) 联系起来．即令总的平方误差
εp =∑
n =p
N -1
{[e (n ) ]+[e (n ) ]}
f
2
b
2
m -1
m -1
2
2
⎡p ⎤N -1⎡p ⎤
=∑⎢∑a p (k ) x (n -k ) ⎥+∑⎢∑a p (k ) x (n +k -p ) ⎥
n =p ⎣k =0⎦n =p ⎣k =0⎦
N -1
a
p
(0) =1
最小．由上式可见，总的平方误差εp 是系数a p (k ) 的函数．若把εp 对各预测系数a p (k ) 而非单一地对a p (p ) =ρ求导数，并令其为零，就可以得到一组线性方
p
程． 解此方程组所得的a p (k ) 就是在最小平方误差准则下的最优预测系数．但由于方程组系数矩阵不是Toeplitz 型，所以不能利用 L — D 算法求解．为了减
少运算量，Ma r p l e 提出一种格型结构的高效递推算法。
2.1.1 AR 模型的稳定性及阶的确定
AR(p)模型稳定的充分必要条件是H(z)的极点(即A(z)的根) 都在单位圆内。稳定的AR(p)模型将具有以下的性质：
(1)H(z)的全部极点或A(z)的所有根都在单位圆内。 (2)自相关矩阵是正定的。
(3)激励信号的方差(能量) 随阶次的增加而递。 (4)反射系数的模恒小于1。
但是在实际应用中，levinson 算法的已知数据(自相关值) 是由 来估计的，有限字长效应有可能造成大的误差，致使估计出来的AR(p)参数所构成的A(z)的根跑到单圆上或单位圆外，从而使模型失去稳定。在递推计算的过程中如果出现这种情况，将致，或|即停止递推计算。
通常事先并不知道AR 模型的阶。阶选得太低，功率谱受到的平滑太厉害，平滑后的谱已经分辨不出真实谱中的两个峰了。阶选的太高，固然会提高谱估计的分辨率，但同时会产生虚假的谱峰或谱的细节。
一种简单而直观的确定AR 模型的阶的方法，是不断增加模型的阶，同时观察预测误差功率，当其下降到最小时，对应的阶便可选定为模型的阶。但是预测误差功率(或AR 模型激励源的方差) 是随着阶次的增加而单调下降的，因此，很难确定降到什么程度才合适。另一方面，应该注意到，随着阶次增加，模型参数的数目亦增多了，谱估计的方差会变大(表现在虚假谱峰的出现) 。因此，不能简
单地依靠观察预测误差功率的下降来确定模型的阶。与此相应的另一种简单方法是观察各阶模型预测误差序列的周期图，当它很接近于平坦(白色谱) 时即对应最佳的阶。
2.2 MA模型及功率谱估计
MA模型是一个全零点模型，其功率谱不易体现信号中的峰值，即分辨率较低。从谱估计的角度来看，MA 模型谱估计等效于经典谱估计中的自相关法，若单纯为了对一段有限长数据做谱估计，就没必要求解MA 模型了。但在系统分析于辨识以及ARMA 谱估计中都要用到MA 模型，因此仍有必要讨论MA 系数的求解方法。目前提出的方法主要有：1，谱分解法；2用高阶的AR 模型来近似MA 模型；3，最大似然法。以第二种方法为例讨论MA 模型参数的提取。
ⅰ有N 点数据x （n ）建立一个p 阶的AR 模型，p>>q，可用AR 模型参数的计算方法求出p 阶的AR 系数a p (k ) ，k=1,2…p； ⅱ利用a
(k )
p
，k=1,2…p进行线性预测，等效于一个q 阶的AR 模型，
再一次利用AR 系数的求解方法得到b （k ），k=1,2，…，q 。 进而通过两次求AR 系数可得MA 模型的系数。
2.3ARMA 模型
理论上可以证明，在一定条件下，A R MA和 M A模型都可以用一个阶次足够大的A R模型来近似，所以在实际使用中我们用一个阶次比较高的A R模型代替A R MA模型来进行功率谱估计，避免了求解非线性方程组带来的困难。
3基于AR 模型的噪声源识别
3.1噪声源识别的方法： 3.1.1 消去法
所谓消去法，就是首先把附带全部装备的试验对象(如汽车、发动机等) 在一定条件下测定其总和的工作噪声，然后去除其中的一部分装备或控制这部分噪声传出，再按同样试验条件测定去除部分装备后的试验对象的工作噪声，则去除这部分装备前后试验对象噪声交化值，即被视作被拆卸或控制传出装备噪声的方法。这一方法试验难度大，耗时耗力，由于对部分设备要求停机而不能进行在线测量，因而不能获得实验对象噪声分布的整体概念.
3.1.2 声强测量法
声强测量法是利用双点声压梯度的积分来近似空气质点的振动速度，并利用FFT 来实现声强的实时测量。采用声强测量法分析噪声可获得丰富的噪声信息，不但可测出噪声的大小，而且还能测出声能的流向。该方法对环境要求低。现场识别能力强．但是很强的背景声会对声强测试产生不良的影响，产生较大的误差，并且其测试精度较难控制。
3.1.3 相干函数法
相干函数法是根据相关理论，用互相关函数描述被测部件的振动信号与外界噪声信号之间的关系．根据结果判断噪声是否由该部件振动引起，从而判断出各个声源．在实际噪声测量中，常常是复杂的多输入单输出系统。而且输入噪声源又可能是相互依赖的，因此这种方法也不适合判断多输入单输出的系统
3.1.4 频谱分析法
频谱分析法，由于声源噪声的形成机理不同。使每个声源的嗓声特点有差异，能量分布的频谱范围有的位于全频带，有的位于中、低频或高频。因而，在了解各组成声源噪声频谱特点的情况下，测量出总噪声频谱。再取得各部分声源的频谱，并将其与总声源频谱比较，就可以分析出各组成部分噪声贡献幅度。从而找出需控制的主要声源，以有效降低高能频率段的噪声。频谱分析法中，除应用幅值频谱外，还常应用功率谱。但是由于采用了传统的功率谱估计，所以存在频率分辨率低，旁瓣泄漏严重的缺点。
由于现代谱估计方法为信号建立了一个准确或近似的模型，从根本上摒弃了对数据序列加窗的隐含假设，因而能够更好地估计信号功率谱．利用现代谱方法估计功率谱时，不要求噪声源中只有一个信号为高斯分布，这与盲信号分离有很大区别(盲信号分离要求源信号中只能有一个信号为高斯分布) 。只要利用能量随频率的分布(即功率谱) 来求取噪声源的比重。现代谱估计的上述优点为噪声源识别提供了更为精确的方法。
3.2研究内容
由于缺少实际的噪声源，在此采用仿真的噪声信号来模拟噪声。我们知道AR 模型的参数与其阶次有关。阶选得太低，功率谱受到的平滑太厉害，平滑后的谱已经分辨不出真实谱中的两个峰了。阶选的太高，固然会提高谱估计的分辨率，但同时会产生虚假的谱峰或谱的细节。所以必须选择合适的AR 模型的阶次。本文采用选用Akaike(alC)准则作为确定模型阶的依据。将信号x1,x2按不同的比例混合在一起组成混合噪声以模仿生活中常见的噪声信号。
3.2.1仿真信号的选取 N=512;
fs=1; a1=5; a2=3; a3=2; w=2*pi/fs;
x1=a1*sin(0.5*w*(0:N-1))+a2*cos(w*(0:N-1)).^2+randn(1,N); x2=a2*sin(w*(0:N-1))+a2*cos(0.8*w*(0:N-1)).^2+randn(1,N); y=1:512; plot(y,x1); figure(2); plot(y,x2); save x1 x1; save x2 x2
;
首先判断AR 模型的阶数。这里应用的就是A1C 准则。我们用仿真信号x1为例进行了阶数测定，与模型阶的关系曲线如图所示。从图中可以发现模型的AIC 存在最小值，此时信号模型对应的阶数为12。所以选取12阶的模型。 clear;
load x1.mat ; x=x1;
N=length(x); M=20;
e1=zeros(M,N); e2=zeros(M,N); e1(1,1:N)=x; e2(1,1:N)=x; r=zeros(1,M); a=zeros(1,M); q1=zeros(1,M); q2=zeros(1,M); AIC=zeros(1,M);
a(1)=(1/N)*sum(x.^2);
q1(1)=sum(e1(1,2:N).*e2(1,1:N-1));
q2(1)=sum((e1(1,2:N).^2)+(e2(1,1:N-1).^2)); r(1)=2*q1(1)/q2(1); AIC(1)=log(a(1)); for k=2:M for n=k:N
e1(k,n)=e1(k-1,n)+((-1)*r(k-1)*e2(k-1,n-1)); e2(k,n)=(-1)*r(k-1)*e1(k-1,n)+e2(k-1,n-1); end
q1(k)=sum(e1(k,k+1:N).*e2(k,k:N-1));
q2(k)=sum((e1(k,k+1:N).^2)+(e2(k,k:N-1).^2)); r(k)=2*q1(k)/q2(k);
a(k)=(1-r(k-1)^2)*a(k-1); AIC(k)=log(a(k))+k*log(N)/N; end
[X,Y]=min(AIC); k=1:M;
plot(k,AIC); [X,Y]
3.2.3功率谱估计及噪声源识别
y1分别取10x1+x2，x1+10x2，x1+x2;我们尝试着去识别主噪声源
在进行信号混合时利用功率谱估计出哪个信号是主噪声源，求功率谱的方法如下：function argonglv(y1,y2,y3)
clc
clear
load x1.mat
load x2.mat
load x3.mat
y1=10*x1+x2;
k1=sqrt(sum(y1.^2));
y1=y1./k1
p1= arburg(y1,12);
la=length(p1);
ln=128-la;
ts=1/5000;
ab=zeros(1,ln);
bb=[p1 ab]
dw=1/128;
w=dw*[0:127]
f=fft(bb)
f1=abs(f);
f2=f1.^2;
f3=1./f2;
f4=10*log10(f3);
plot(w(1:128),f4(1:128),'-kX' )
hold on
figure (1);
y1=x1;
k1=sqrt(sum(y1.^2));
y1=y1./k1
p1= arburg(y1,12);
la=length(p1);
ln=128-la;
ts=1/5000;
ab=zeros(1,ln);
bb=[p1 ab]
dw=1/128;
w=dw*[0:127]
f=fft(bb)
f1=abs(f);
f2=f1.^2;
f3=1./f2
f4=10*log10(f3);
plot(w(1:128),f4(1:128),'-r*')
hold on
figure(1);
y1=x2;
k1=sqrt(sum(y1.^2));
y1=y1./k1
p1= arburg(y1,12);
la=length(p1);
ln=128-la;
ts=1/5000;
ab=zeros(1,ln);
bb=[p1 ab]
dw=1/128;
w=dw*[0:127]
f=fft(bb)
f1=abs(f);
f2=f1.^2;
f3=1./f2
f4=10*log10(f3);
plot(w(1:128),f4(1:128),':+')
legend('y1' , 'x1' , 'x2' );
(1) 当源信号的比重大于1时，我们可以明显地看出混合信号的功率谱与x1的功
率谱相接近，也就是混合信号以x1信号为主噪声源如下图示
(2)当源信号的比重小于1时，我们可以明显地看出混合信号的功率谱与x2相接近，也就是混合信号以源信号x2为主噪声源，功率谱如下图所示
(2) 当源信号x1与x2以l ：l 的比例混合时，我们难以确定哪个是主要源信号。
由上面的讨论知：我们可以利用现代谱估计的理论方法是可以对主要噪声源信号进行识别的，也促使我们利用真实的噪声信号进行进一步的研究。只要选用合适的参数模型阶数，求出信号的模型参数和功率谱图，然后利用曲线相似度和曲线关联度进一步求取噪声源信号在混合信号中的比重，即解决信号功率谱图中曲线的相似度问题。现代谱估计方法的噪声源识别作为信号处理方法的一种，克服了仪器测量和环境的影响，同时克服了传统功率谱估计的缺点，因而在噪声识别的仿真实验中取得了较好的结果。