现代信号处理功率谱估计课件

七、最大熵谱估计1、利用最大熵的原则外推自相关函数2、最大熵谱估计与AR模型谱估计的等价性八、最大似然谱估计1、最小方差谱估计2、最大熵谱与最大似然谱估计的关系九、特征分解法谱估计1、正弦波用退化AR模型表示2、白噪声中正弦波组合用一特殊的ARMA模型表示3、特征分解法谱估计功率谱估计七、最大熵谱估计功率谱估计1十、Prony谱分析法1、利用最大熵的原则外推自相关函数2、最大熵谱估计与AR模型谱估计的等价性十一、多重信号分类MUSIC1、最小方差谱估计2、最大熵谱与最大似然谱估计的关系十二、特征分解法谱估计1、波束形成器2、特征子空间分析3、MUSIC算法及其改进功率谱估计十、Prony谱分析法功率谱估计2一、最大熵谱估计

1.利用最大熵的原则外推自相关函数

按照Shannon对熵的定义，

当随机变量X取离散值时，熵的定义为

(4.6.1)式中pi是出现状态i的概率。当X取连续值时，熵的定义为

(4.6.2)一、最大熵谱估计1.利用最大熵的原则外推自相关3式中,p(x)是X的概率密度函数，对于离散随机序列,概率密度函数用联合概率密度函数代替。显然,熵代表一种不确定性，最大熵代表最大的不确定性，或者说最大的随机性。下面我们研究对于有限的自相关函数值不作任何改变，对于未知自相关函数用最大熵原则外推，即不作任何附加条件的外推方法。假设x(n)是零均值正态分布的平稳随机序列，它的N维高斯概率密度函数为式中

式中,p(x)是X的概率密度函数，对于离散随机序列,概率4按照（4.6.2）式，x(n)信号的熵为

（4.6.3）

式中det(Rxx(N))表示矩阵Rxx(N)的行列式，由上式表明为使熵最大，要求det(Rxx(N)最大。按照（4.6.2）式，x(n)信号的熵为（4.6.3）式5若已知N+1个自相关函数值rxx(0),rxx(1),…,rxx(N)，下面用最大熵方法外推rxx(N+1)。设rxx(N+1)确实是信号自相关函数的第N+2个值，根据自相关函数的性质，由N+2个自相关函数组成的矩阵为（4.6.4）

若已知N+1个自相关函数值rxx(0),rx6它必须是非负定的矩阵，

即

（4.6.5）

将行列式展开，det(Rxx(N+1))是rxx(N+1)的二次函数，该二次函数系数的符号是：(-1)1+N+2(-1)1+N+1=-1，且det(Rxx(N+1))对rxx(N+1)的二次导数是-2det［Rxx(N-1)］，它是负值，负值表示det(Rxx(N+1))对rxx(N+1)的一次导数是减函数，det(Rxx·(N+1))作为rxx(N+1)的函数，凹口向下，那么只有一个最大值。为选择rxx(N+1)使det(Rxx(N+1)最大，解下列方程：(4.6.6)它必须是非负定的矩阵，即（4.6.5）将行列式展开，d7用数学归纳法，得到

（4.6.7）

上式是rxx(N+1)的一次函数，可以解出rxx(N+1)。继续再将rxx(N+1)代入Rxx(N+2)和det(Rxx(N+2))中，求det(Rxx(N+2))对rxx(N+2)的最大值，得到rxx(N+2);以此类推，可推出任意多个其它自相关函数值，而不必假设它们为零，这就是最大熵谱估计的基本思想。用数学归纳法，得到（4.6.7）上式是rxx(N+1)的82.最大熵谱估计与AR模型谱估计的等价性

我们已经知道AR模型信号自相关函数与模型参数服从Yule-Walker方程，即

将m≥1的情况写成矩阵形式：

m>0m=02.最大熵谱估计与AR模型谱估计的等价性我们已经9式中ai是AR模型系数,i=1,2,3,…,N,。在AR模型中,列写齐次方程式,可得（4.6.8）

及式中ai是AR模型系数,i=1,2,3,…,N,。10利用N个参数,由齐次方程组即可解得a1,a2,…,aN值,再将得到的参数值代入(4.6.8)式,并将它整理成行列式:可以看出AR模型得到的结果与按最大熵外推rxx(N+1)得到的结果一致，这就证明了当x(n)为高斯分布时的最大熵谱估计与AR模型法是等价的。上式（4.6.8)是rxx(N+1)的一次函数，由此可解得rxx(N+1)。再用类似的方法求得rxx(N+2)，rxx(N+3)，┄，然后确定功率谱估计。利用N个参数,由齐次方程组即可解得a1,a2,…,aN值,再11最大熵谱估计用下式计算信号功率谱：

（4.6.9）

最大熵谱估计用下式计算信号功率谱：（4.6.9）12二、最大似然谱估计１、最小方差谱估计最大似然谱估计是用一个FIR滤波器实现，该滤波器对所关心频率的正弦信号，可以无失真地通过，而对于其它频率的信号，让其频响尽可能地小，亦即将它们尽可能地滤除。此时,滤波器输出的均方值，就作为信号的功率谱估计。设实信号用x(n)表示，FIR滤波器系统函数用A(z)表示：

输出y(n)为

(4.6.10)二、最大似然谱估计输出y(n)为(4.6.10)13式中

输出信号的均方值为

(4.6.11)上式中T表示转置,H表示共轭转置,Rp=E［XXT］是Toeplith自相关矩阵,为求,必须先求FIR滤波器的系数。求这些系数的原则是：在所关心频率ωi处，信号x(n)无失真地通过，即在ωi处的传输函数为1：

式中输出信号的均方值为(4.6.11)上式中T表示转置14式中

(4.6.12)另外一个原则是在ωi附近的频率分量尽量衰减掉，即ω≠ωi处,滤波器输出y(n)的均方差最小,即(4.6.11)式最小,此时作为信号x(n)的功率谱估计 。因此,最大似然谱估计称为最小方差谱估计更为合适，但由于习惯也可以仍称为最大似然谱估计。在以上原则下，使方差最小的滤波器系数和分别为［30］、［31］

式中(4.6.12)另外一个原则是在ωi附近的频率分量尽15应该指出,此时并不是真正意义上的信号功率谱,只是描述了信号功率谱的相对强度。

应该指出,此时并不是真正意义上的信号功率谱,只是描16

2、最大熵谱与最大似然谱估计的关系伯格证明了最大熵谱PMEM与最大似然谱PMLM估计的关系从上式可知最大似然谱估计相当于从最大熵谱估计的最低分辨率到最高分辨率的平均，所以最大熵谱估计的分辨率要比最大似然谱估计的分辨率高。但最大似然谱具有更大的统计稳定性，对模型阶数的依赖性要小于最大熵谱估计。另外在最大熵谱估计中提到，它的最大缺点就是求得最佳频率成分后，其相应的振幅值并不代表原来的振幅值，尚须用其他办法来近似确定。通过其他两位同学的介绍，我们知道，频谱估计中，振幅谱常用傅立叶变换（传统法）求得，功率谱可通过振幅谱的平方求得，另外也可通过自相关函数的傅立叶变换求得。随机信号一般只作功率谱估计，所以功率谱估计在谱估计中占有重要地位。它的主要缺点是失去了相位信息，因此光靠功率谱是无法恢复信号的。2、最大熵谱与最大似然谱估计的关系伯格证明了最大熵谱PME17三、特征分解法谱估计

4.7.1正弦波用退化AR模型表示无论是实正弦波还是复正弦波，都可以用一个退化AR模型表示，设P个实正弦波组成的信号用下式表示：

（4.7.1）

式中，初相位θi是在区间（-π，π）均匀分布的随机变量，首先分析下面的三角恒等式：

-π＜ω＜π

三、特征分解法谱估计4.7.1正弦波用退化AR模型表示18令x(n)=sin(ωn+θ),则上式变为

（4.7.2）

将上式进行Z变换，得到

（4.7.3）

这样(4.7.2)式的特征多项式为

（4.7.4）

令x(n)=sin(ωn+θ),则上式变为（4.7.2）19上式的两个根分别是:z1=ejω,z2=e-jω,它们共轭成对，且模为1。由这两个根可以确定正弦波的频率。对比AR模型的系统函数，可以把正弦波信号用一个特殊的AR（2）模型表示，括弧中的2表示模型是二阶的。该AR模型的激励白噪声方差趋于0，极点趋于单位圆。通常称为退化的AR模型。这一模型系数有两个，即2cosω和1，（4.7.2）式是模型的差分方程。

上式的两个根分别是:z1=ejω,z2=e-jω,它们共轭成20对于P个实正弦波,特征多项式是

（4.7.5）

上式是z-1的2P阶多项式，可以表示为

（4.7.6）

注意上式中的系数ak(k=1,2,3,…,2P),必须保证它的根共轭成对。考虑到根共轭成对，也可表示为

（4.7.7）

对于P个实正弦波,特征多项式是（4.7.5）上式是z-21这样由(4.7.6)式，P个正弦波组合的模型用下面2P阶差分方程描述

（4.7.8）

对于复正弦波情况，P个复正弦波组成的信号是

（4.7.9）

用一个退化的AR（p）模型表示的差分方程为

（4.7.10）

这样由(4.7.6)式，P个正弦波组合的模型用下面2P阶差分22其特征多项式为

(4.7.11)其根为

1≤i≤P

注意这里的根不是共轭成对出现的。总结以上P个正弦波组合是一个退化的AR(2P)过程，独立参量个数为P个；P个复正弦波的组合是退化的AR（P）过程，独立参量个数仍为P个。实正弦过程相应的退化AR过程的阶数比复正弦情况的阶数高1倍。

其特征多项式为(4.7.11)其根为1≤i≤P注意这234.7.2白噪声中正弦波组合用一特殊的ARMA模型表示白噪声中正弦波组合的信号为

（4.7.12）

式中，w(n)为白噪声,且

将(4.7.12)式中x(n)的用AR(2P)表示，即将(4.7.8)式带入(4.7.12)式中，

得到

（4.7.13）

4.7.2白噪声中正弦波组合用一特殊的ARMA模型表示（424将（4.7.12）式中的n用n-i代替,x(n-i)=y(n-i)-w(n-i)再将上式带入(4.7.13)式,得到

(4.7.14)上式可以看成一个特殊的ARMA(2P,2P)模型，括弧中的两个2P分别表示ARMA模型系统函数分子和分母的阶次。它与一般的ARMA模型比较，有三方面不同：

将（4.7.12）式中的n用n-i代替,x(n-i)=y(25(1)它的AR部分和MA部分具有相同的参数，它们存在共同的因子；(2)由于特征多项式(4.7.6)式的根的模为1，故AR部分特征多项式不满足平稳性条件，MA部分特征多项式也不满足可逆性条件；（3）AR部分的y(n)=x(n)+w(n)，y(n)是含白噪声的观测值，而通常为信号的x(n)不含白噪声。(1)它的AR部分和MA部分具有相同的参数，264.7.3特征分解法谱估计这种特殊的ARMA模型结构，不能用一般的ARMA模型结构求解。下面介绍特征分解技术。将（4.7.14）式写成矩阵形式：

YTA=WTA

(4.7.15)式中

4.7.3特征分解法谱估计YTA=WTA(4.7.15)27用向量Y左乘(4.7.15)式,并取数学期望,得到

E［YYT］A=E［YW

T］A

(4.7.16)式中

用向量Y左乘(4.7.15)式,并取数学期望,得到E［28将上面关系式带入(4.7.16)式,得到(4.7.17)式中,σw2是自相关函数Ryy的特征值;A是对应σw2的特征矢量。

由于x(n)与w(n)互不相关，由（4.7.12）式求y(n)的自相关函数,得到

式中Rxx是x(n)的自相关函数,可以推导出在特征方程(4.7.17)式中,σw2是Ryy的最小特征值,且Ryy的阶数超过(2P+1)(2P+1)时，σw2就是其多重最小特征值。这一结论为我们寻求向量A提供了重要的依据。

(4.7.18)将上面关系式带入(4.7.16)式,得到(4.7.17)29当特征向量A求出后，就可以通过解特征方程，解出各个根，求出各正弦波的频率值，特征方程为

（4.7.19）

该方程有2P个根，这些根在z平面的单位圆上，它们是

i=1,2,3,…,2P

（4.7.20）

式中的ωi即是正弦波的频率。当各正弦频率由(4.7.19)式求出后,各正弦波的功率也可以求出,对于白噪声中P个实正弦组合信号,它的自相关函数为

（4.7.21）

当特征向量A求出后，就可以通过解特征方程，解出各个根30m≠0(4.7.22)式中

qi是频率为ωi的正弦波的幅度，Pi是其功率。将(4.7.22)式写成矩阵形式:FP=r(4.7.23)式中

(4.7.24)m≠0(4.7.22)式中qi是频率为ωi的正弦波的幅31PT=［P1,P2,…,PP］

(4.7.25)(4.7.26)由于各正弦波的频率已求出,矩阵F,r已知,由(4.7.23)式解出各正弦波的功率或幅度。

最后噪声功率由下式求出:(4.7.27)以上就是皮萨论科谱分解法的全过程。PT=［P1,P2,…,PP］(4.7.25)(432Prony谱分析法

Prony曾在1975年提出用指数函数的一个线性组合去描述等间距采样数据的数学模型，其方法并不是通常意义下的谱估计技术。这里介绍有Kay和Marple（1981年）提出的可用于估计非有理式谱密度的扩充的Prony方法。基本的Prony法是一种内插法。它是采用等间隔复指数值的线性组合来拟合观测数据。每个复指数包含幅度和指数因子两部分，因此p个复制数值就含有2p个代求因子，若用p个指数项去逼近2p个数据，这种逼近可以精确实现。Prony谱分析法Prony曾在1975年提出用指数33

当用p个指数去逼N(N>2p)个数据时，可用最小二乘法准则来实现Prony方法的谱估计但这是一种非线性最小二乘法问题，直接求解相当困难，常用迭代法求解。经过适当的扩展后，Prony方法可用来计算有理式功率谱密度，本节介绍这种扩展Prony方法。扩展Prony方法采用的数学模型为一组p个具有任意幅值、相位、频率与衰减因子的复指数线型组合来拟合观测数据，其离散时间的函数形式为Prony谱分析法当用p个指数去逼N(N>2p)个数据时，可用最小34现代信号处理功率谱估计课件35现代信号处理功率谱估计课件36现代信号处理功率谱估计课件37现代信号处理功率谱估计课件38现代信号处理功率谱估计课件39现代信号处理功率谱估计课件40现代信号处理功率谱估计课件41现代信号处理功率谱估计课件42现代信号处理功率谱估计课件43现代信号处理功率谱估计课件44多重信号分类MUSIC

--现代谱估计的应用多重信号分类MUSIC

--现代谱估45

MUSIC为：MultipleSignalClassification即多重信号分类。MUSIC方法是一种估计空间参数的现代谱计算方法，它将功率谱推广为空间谱，是最早问世的子空间方法。这种算法是通过对数据协方差矩阵进行本征分解获得信号空间谱估计的。MUSIC为：46

当空间存在多个信号源时，常常用要对些空间信号进行分离，以便跟踪或检测我们感兴趣的空间信号，抑制那些被认为是干扰的空间信号。为此，需要使用天线阵列对多个空间信号进行接收。对天线阵列接收的空间信号所进行的分析与处理统称为阵列信号处理。

当空间存在多个信号源时，常常用要对些空间信号进行分离，以47

空间谱是阵列信号处理中的一个重要概念。时域频谱表示信号在各个频率上的能量分布，而空间谱则是表示信号在空间各个方向上的能量分布，它能够对多个空间信号进行识别（即分类）。空间谱是阵列信号处理中的一个重要概念。时域频谱表示信号在48波束形成器特征子空间分析MUSIC算法及其改进波束形成器491.波束形成器

最佳波束形成器的基本原理是使来自非期望波达方向（DOA－directionofarrival）的任何干扰所贡献的功率最小，但是又能够保持“在观测方向上的信号功率不变”，也就是说通过将个阵元输出进行加权求和，在一时间内将天线阵列波束“导向”到一个方向上，通过对期望信号得到最大输出功率即给出波达方向的估计，因此这个也可以看成为一个尖锐的空间带通滤波器。MUSIC算法作为一种经典的高分辨的波达方向估计算法，与传统测向方法相比，它具有很高的分辨率1.波束形成器最佳波束形成器的基本原理是使来自非期望波50

在空间传播的多个信号源（窄带信号），现利用一天线阵列对这些信号进行接收。每根天线为一个阵元，它们都是全向天线，即无任何方向性。这里假定各个阵元等间距地直线排列。如下图：在空间传播的多个信号源（窄带信号），现利用一天线阵列对这51窄带包络变化缓慢，所以等距线阵各阵元接收到的同一信号的包络相同。假设为远场信号，则到达各阵元的方向角相同，用表示，称为波达方向（角），定义为信号到达阵元的直射线与阵列法线方向之间夹角。对一般的远场信号而言，同一信号到达不同的阵元存在一个波程差，这个波程差导致各接收阵元间的相位差。以上图的阵元1作为基准点，令信号传播延迟引起的相位差为，则由上图易知，波达方向与相位差之间的关系为因此，信号在第k个阵元上的接收信号为窄带包络变化缓慢，所以等距线阵各阵元接收到的同52

假设阵列有m个阵元组成，则信号到达各阵元的相位差所组成的向量为称为信号的方向向量。若共有p个信号位于远场（p<m），则在第k个阵元上的接收信号为其向量形式可以写成假设阵列有m个阵元组成，则信号到达各阵元的相位差53其中分别为m*p的方向矩阵和p*1的信号向量。具有3.9.6结构的是Vandermonde矩阵，当各不相同时，则列线性独立。在阵列信号处理中，一次采样称为一次快拍。假定每个阵元上共观测到N次快拍的接收信号,…,。波束形成问题就是依靠这些观测值，求得某个期望信号的波达方向。其中54解决问题：先假定信号与各阵元上的观测噪声独立统计，并且各观测噪声具有相同的方差。此时考察第k个信号的接收，那其它信号全部视为干扰信号，应抑制。我们知道一个阵列经过加权求和后，可以将阵列的方向增益聚焦在一个方向上，相当于形成一个波束。为此设计权向量

解决问题：55得到输出信号则N次快拍阵列输出的平均功率为其中

得到输出信号56

我们的目的是保证来自某个确定方向的信号能正确被接收，而其它方向的信号或干扰被完全抑制。用数学表示就是，在保证所需方向的信号输出为为常数的条件下，使阵列的输出功率最小化，即上面的第二个式子表明只有第k个信号的波达方向可以做为期望抽取。我们的目的是保证来自某个确定方向的信号能正确被接收，而57

下面用Largange常数法求解上述优化问题。为此构造目标函数对求偏导得，从中得到使输出能量最小化的最佳波束形成器此式代入，可知所以，这一波束形成器是Capon于1969年提出的，称为最小方差无畸变响应波束形成器。下面用Largange常数法求解上述优化问题。为此构造目58

将的结果代入平均功率的表达式，则但在实际应用中由于不能确切知道信号的方向，所以只能通过扫描的方法得到曲线这就是Capon定义的“空间谱”，它将峰值对应的定为p个信号的波达方向。它是著名的极大似然谱估计，也称为最小方差法。将的结果代入平均功率的表达式，则592.特征子空间分析窄带远场信号DOA数学模型为

对此模型，我们做以下假设：假设1：对于不同的值，向量相互线性独立假设2：加性噪声向量的每个元素都是零均值的白噪声，它们不相关，并且有相同的方差假设3：矩阵非奇异，即rank（P）＝pMUSIC是以特征子空间分析为基础的2.特征子空间分析窄带远场信号DOA数学模型为MUSIC是以60在假设1～3下，由模型的表达式很容易得到

由于，则为Herimitian矩阵。令其特征分解为,式中，U为特征向量构成的酉矩阵。

在假设1～3下，由模型的表达式很容易得到61

将R的特征值按照从小到大的顺序排列，则前p个与信号有关，其数值大于，从p+1开始的特征值完全决定于噪声，其数值等于。很自然地，前p个主特征值称为信号特征值，其余m－p个次特征值称为噪声特征值。根据信号特征值和噪声特征值，可以将特征矩阵U的列向量分为两个部分，即U=[S|G]将R的特征值按照从小到大的顺序排列，则前p个与信号有关623.MUSIC算法及其改进

由于是R的噪声特征对，故有用G右乘，有所以，进而有,也就是说。等价有,显然，当时，将上述等价式改写成标量形式，可以定义一种类似功率谱的函数：3.MUSIC算法及其改进由于是R的噪声特征对63

上式去峰值的p个w值就可以得到p个信号的波达方向。上式称为噪声子空间法，而下式则称为信号子空间法。在实际应用中，通常将w划分为数百个等间距的单位，得到，然后将每个值代入功率谱式以求得最大的峰值所对应的w值。另外，执行MUSIC算法式选择噪声子空间还是信号子空间方式，决定于G和S中哪个有更小的唯数。除了计算量有所不同外，两种方式没有本质区别。上式去峰值的p个w值就可以得到p个信号的64MUSIC算法步骤：步骤1：计算样本自相关矩阵的特征值分解以得到主特征值和次特征值，看p和m－p哪个更小，若p小则选择信号子空间法，否则选择噪声子空间法，并存储其相应的特征向量。步骤2：计算MUSIC谱P(wi)，其中，网格可取作2*pi*0.001等步骤3:找出P（w）的p个峰值，它们就是待求的MUSIC估计值。MUSIC算法步骤：65优点：实践证明，MUSIC算法具有较好的性能和较高的效率，能提供高分辨率及渐近无偏的到达角估计.用特征值的大小不同，将协方差矩阵的特征空间划分成正交的信号子空间和噪声子空间，概念直观、清晰，可以高分辨地处理多个同时信号，可以提取信号数量，可以高精度和高一致性地估计信号参量等等。优点：66缺点：MUSIC算法建立在不相干信号模型的基础上的，对于有相干信号存在的情况，该算法将会失效。而对广泛使用相干技术的现代雷达而言。如果两个目标的多普勒频移相等，从目标返回的回波信号将是完全相关的。下面就将探讨解相干的MUSIC算法。缺点：67

由以上的讨论可知，MUSIC算法在理想条件下具有良好的性能，但在信号相干时算法的性能变得很坏。若信号相干，则信号子空间的“扩散”到噪声子空间，这会导致某些相干源的矢量与噪声子空间不完全正交，从而无法正确估计信号源方向。目前关于解相干的处理基本有两类：降唯和非降唯。这里我们只讨论降唯处理中的基于空间平滑算法。空间平滑算法在一般情况下只适用于均匀线阵，空间平滑MUSIC是利用子阵平滑恢复数据协方差矩阵。由以上的讨论可知，MUSIC算法在理想条件下具有良好的68

空间平滑的基本思想是将等距线阵分为若干个相重叠的子阵列，若各子阵列的流形相同，则各子阵列的协方差矩阵可以进行平均运算。空间平滑技术的原理图如下：假设m为阵元数，N 为信号源数将均匀线阵m个阵元分成如图所示的M个子阵列，每个子阵列的阵元数为p，即有m＝M+p-1。空间平滑的基本思想是将等距线阵分为若干个相重叠的子阵列，69

如上图所示，取第一个子阵为参考子阵，则对于第k个子阵有数据模型其中，，A为m*N维的方向向量，也就是参考子阵的方向向量。表示矩阵D的k-1次幂。于是该子阵数据协方差矩阵为P为前面所定义的信号的自相关矩阵。如上图所示，取第一个子阵为参考子阵，则对于第k个子阵有70

前向空间平滑MUSIC方法对满秩方差的恢复是通过求各子阵协方差矩阵的均值来实现的，即取前向平滑修正的协方差矩阵为其中，定理：如果子阵阵元数目大于等于信号源数，则当子阵列数目大于等于信号源数时前向空间平滑数据协方差矩阵是满秩的。前向空间平滑MUSIC方法对满秩方差的恢复是通过求各子71

空间平滑的缺点是：阵列的有效孔径减小了，因为子阵列比原阵列小。然而，尽管存在这一空径损失，空间平滑变换减轻了所有子空间估计技术的局限性，并将保留一维谱搜索的计算有效性。空间平滑的缺点是：阵列的有效孔径减小了，因为子阵列比原72基于空间平滑的MUSIC算法：步骤1:有阵列的接收数据得到数据协方差矩阵步骤2:利用本部分介绍的方法对R进行修正步骤3:利用修正后的协方差矩阵进行MUSIC谱估计，找出极大值对应的信号方向。基于空间平滑的MUSIC算法：73附：为A的共轭转置，容易证明一实的方阵,若，称称其为正交矩阵.一复值方阵，若，称其为酉矩阵。则正交矩阵实际上就是实的酉矩阵。附：74附：若矩阵,则称其为Hermitian矩阵。若A非奇异、为Hermitian矩阵并可以成则A有唯一的分解为得到上式，必须，，也就是说一个Hermitian矩阵有实的特征根，其特征向量组成正交集，当然Q就是酉矩阵。因此Hemitian矩阵可用酉矩阵来对角化。附：75七、最大熵谱估计1、利用最大熵的原则外推自相关函数2、最大熵谱估计与AR模型谱估计的等价性八、最大似然谱估计1、最小方差谱估计2、最大熵谱与最大似然谱估计的关系九、特征分解法谱估计1、正弦波用退化AR模型表示2、白噪声中正弦波组合用一特殊的ARMA模型表示3、特征分解法谱估计功率谱估计七、最大熵谱估计功率谱估计76十、Prony谱分析法1、利用最大熵的原则外推自相关函数2、最大熵谱估计与AR模型谱估计的等价性十一、多重信号分类MUSIC1、最小方差谱估计2、最大熵谱与最大似然谱估计的关系十二、特征分解法谱估计1、波束形成器2、特征子空间分析3、MUSIC算法及其改进功率谱估计十、Prony谱分析法功率谱估计77一、最大熵谱估计

1.利用最大熵的原则外推自相关函数

按照Shannon对熵的定义，

当随机变量X取离散值时，熵的定义为

(4.6.1)式中pi是出现状态i的概率。当X取连续值时，熵的定义为

(4.6.2)一、最大熵谱估计1.利用最大熵的原则外推自相关78式中,p(x)是X的概率密度函数，对于离散随机序列,概率密度函数用联合概率密度函数代替。显然,熵代表一种不确定性，最大熵代表最大的不确定性，或者说最大的随机性。下面我们研究对于有限的自相关函数值不作任何改变，对于未知自相关函数用最大熵原则外推，即不作任何附加条件的外推方法。假设x(n)是零均值正态分布的平稳随机序列，它的N维高斯概率密度函数为式中

式中,p(x)是X的概率密度函数，对于离散随机序列,概率79按照（4.6.2）式，x(n)信号的熵为

（4.6.3）

式中det(Rxx(N))表示矩阵Rxx(N)的行列式，由上式表明为使熵最大，要求det(Rxx(N)最大。按照（4.6.2）式，x(n)信号的熵为（4.6.3）式80若已知N+1个自相关函数值rxx(0),rxx(1),…,rxx(N)，下面用最大熵方法外推rxx(N+1)。设rxx(N+1)确实是信号自相关函数的第N+2个值，根据自相关函数的性质，由N+2个自相关函数组成的矩阵为（4.6.4）

若已知N+1个自相关函数值rxx(0),rx81它必须是非负定的矩阵，

即

（4.6.5）

将行列式展开，det(Rxx(N+1))是rxx(N+1)的二次函数，该二次函数系数的符号是：(-1)1+N+2(-1)1+N+1=-1，且det(Rxx(N+1))对rxx(N+1)的二次导数是-2det［Rxx(N-1)］，它是负值，负值表示det(Rxx(N+1))对rxx(N+1)的一次导数是减函数，det(Rxx·(N+1))作为rxx(N+1)的函数，凹口向下，那么只有一个最大值。为选择rxx(N+1)使det(Rxx(N+1)最大，解下列方程：(4.6.6)它必须是非负定的矩阵，即（4.6.5）将行列式展开，d82用数学归纳法，得到

（4.6.7）

上式是rxx(N+1)的一次函数，可以解出rxx(N+1)。继续再将rxx(N+1)代入Rxx(N+2)和det(Rxx(N+2))中，求det(Rxx(N+2))对rxx(N+2)的最大值，得到rxx(N+2);以此类推，可推出任意多个其它自相关函数值，而不必假设它们为零，这就是最大熵谱估计的基本思想。用数学归纳法，得到（4.6.7）上式是rxx(N+1)的832.最大熵谱估计与AR模型谱估计的等价性

我们已经知道AR模型信号自相关函数与模型参数服从Yule-Walker方程，即

将m≥1的情况写成矩阵形式：

m>0m=02.最大熵谱估计与AR模型谱估计的等价性我们已经84式中ai是AR模型系数,i=1,2,3,…,N,。在AR模型中,列写齐次方程式,可得（4.6.8）

及式中ai是AR模型系数,i=1,2,3,…,N,。85利用N个参数,由齐次方程组即可解得a1,a2,…,aN值,再将得到的参数值代入(4.6.8)式,并将它整理成行列式:可以看出AR模型得到的结果与按最大熵外推rxx(N+1)得到的结果一致，这就证明了当x(n)为高斯分布时的最大熵谱估计与AR模型法是等价的。上式（4.6.8)是rxx(N+1)的一次函数，由此可解得rxx(N+1)。再用类似的方法求得rxx(N+2)，rxx(N+3)，┄，然后确定功率谱估计。利用N个参数,由齐次方程组即可解得a1,a2,…,aN值,再86最大熵谱估计用下式计算信号功率谱：

（4.6.9）

最大熵谱估计用下式计算信号功率谱：（4.6.9）87二、最大似然谱估计１、最小方差谱估计最大似然谱估计是用一个FIR滤波器实现，该滤波器对所关心频率的正弦信号，可以无失真地通过，而对于其它频率的信号，让其频响尽可能地小，亦即将它们尽可能地滤除。此时,滤波器输出的均方值，就作为信号的功率谱估计。设实信号用x(n)表示，FIR滤波器系统函数用A(z)表示：

输出y(n)为

(4.6.10)二、最大似然谱估计输出y(n)为(4.6.10)88式中

输出信号的均方值为

(4.6.11)上式中T表示转置,H表示共轭转置,Rp=E［XXT］是Toeplith自相关矩阵,为求,必须先求FIR滤波器的系数。求这些系数的原则是：在所关心频率ωi处，信号x(n)无失真地通过，即在ωi处的传输函数为1：

式中输出信号的均方值为(4.6.11)上式中T表示转置89式中

(4.6.12)另外一个原则是在ωi附近的频率分量尽量衰减掉，即ω≠ωi处,滤波器输出y(n)的均方差最小,即(4.6.11)式最小,此时作为信号x(n)的功率谱估计 。因此,最大似然谱估计称为最小方差谱估计更为合适，但由于习惯也可以仍称为最大似然谱估计。在以上原则下，使方差最小的滤波器系数和分别为［30］、［31］

式中(4.6.12)另外一个原则是在ωi附近的频率分量尽90应该指出,此时并不是真正意义上的信号功率谱,只是描述了信号功率谱的相对强度。

应该指出,此时并不是真正意义上的信号功率谱,只是描91

2、最大熵谱与最大似然谱估计的关系伯格证明了最大熵谱PMEM与最大似然谱PMLM估计的关系从上式可知最大似然谱估计相当于从最大熵谱估计的最低分辨率到最高分辨率的平均，所以最大熵谱估计的分辨率要比最大似然谱估计的分辨率高。但最大似然谱具有更大的统计稳定性，对模型阶数的依赖性要小于最大熵谱估计。另外在最大熵谱估计中提到，它的最大缺点就是求得最佳频率成分后，其相应的振幅值并不代表原来的振幅值，尚须用其他办法来近似确定。通过其他两位同学的介绍，我们知道，频谱估计中，振幅谱常用傅立叶变换（传统法）求得，功率谱可通过振幅谱的平方求得，另外也可通过自相关函数的傅立叶变换求得。随机信号一般只作功率谱估计，所以功率谱估计在谱估计中占有重要地位。它的主要缺点是失去了相位信息，因此光靠功率谱是无法恢复信号的。2、最大熵谱与最大似然谱估计的关系伯格证明了最大熵谱PME92三、特征分解法谱估计

4.7.1正弦波用退化AR模型表示无论是实正弦波还是复正弦波，都可以用一个退化AR模型表示，设P个实正弦波组成的信号用下式表示：

（4.7.1）

式中，初相位θi是在区间（-π，π）均匀分布的随机变量，首先分析下面的三角恒等式：

-π＜ω＜π

三、特征分解法谱估计4.7.1正弦波用退化AR模型表示93令x(n)=sin(ωn+θ),则上式变为

（4.7.2）

将上式进行Z变换，得到

（4.7.3）

这样(4.7.2)式的特征多项式为

（4.7.4）

令x(n)=sin(ωn+θ),则上式变为（4.7.2）94上式的两个根分别是:z1=ejω,z2=e-jω,它们共轭成对，且模为1。由这两个根可以确定正弦波的频率。对比AR模型的系统函数，可以把正弦波信号用一个特殊的AR（2）模型表示，括弧中的2表示模型是二阶的。该AR模型的激励白噪声方差趋于0，极点趋于单位圆。通常称为退化的AR模型。这一模型系数有两个，即2cosω和1，（4.7.2）式是模型的差分方程。

上式的两个根分别是:z1=ejω,z2=e-jω,它们共轭成95对于P个实正弦波,特征多项式是

（4.7.5）

上式是z-1的2P阶多项式，可以表示为

（4.7.6）

注意上式中的系数ak(k=1,2,3,…,2P),必须保证它的根共轭成对。考虑到根共轭成对，也可表示为

（4.7.7）

对于P个实正弦波,特征多项式是（4.7.5）上式是z-96这样由(4.7.6)式，P个正弦波组合的模型用下面2P阶差分方程描述

（4.7.8）

对于复正弦波情况，P个复正弦波组成的信号是

（4.7.9）

用一个退化的AR（p）模型表示的差分方程为

（4.7.10）

这样由(4.7.6)式，P个正弦波组合的模型用下面2P阶差分97其特征多项式为

(4.7.11)其根为

1≤i≤P

注意这里的根不是共轭成对出现的。总结以上P个正弦波组合是一个退化的AR(2P)过程，独立参量个数为P个；P个复正弦波的组合是退化的AR（P）过程，独立参量个数仍为P个。实正弦过程相应的退化AR过程的阶数比复正弦情况的阶数高1倍。

其特征多项式为(4.7.11)其根为1≤i≤P注意这984.7.2白噪声中正弦波组合用一特殊的ARMA模型表示白噪声中正弦波组合的信号为

（4.7.12）

式中，w(n)为白噪声,且

将(4.7.12)式中x(n)的用AR(2P)表示，即将(4.7.8)式带入(4.7.12)式中，

得到

（4.7.13）

4.7.2白噪声中正弦波组合用一特殊的ARMA模型表示（499将（4.7.12）式中的n用n-i代替,x(n-i)=y(n-i)-w(n-i)再将上式带入(4.7.13)式,得到

(4.7.14)上式可以看成一个特殊的ARMA(2P,2P)模型，括弧中的两个2P分别表示ARMA模型系统函数分子和分母的阶次。它与一般的ARMA模型比较，有三方面不同：

将（4.7.12）式中的n用n-i代替,x(n-i)=y(100(1)它的AR部分和MA部分具有相同的参数，它们存在共同的因子；(2)由于特征多项式(4.7.6)式的根的模为1，故AR部分特征多项式不满足平稳性条件，MA部分特征多项式也不满足可逆性条件；（3）AR部分的y(n)=x(n)+w(n)，y(n)是含白噪声的观测值，而通常为信号的x(n)不含白噪声。(1)它的AR部分和MA部分具有相同的参数，1014.7.3特征分解法谱估计这种特殊的ARMA模型结构，不能用一般的ARMA模型结构求解。下面介绍特征分解技术。将（4.7.14）式写成矩阵形式：

YTA=WTA

(4.7.15)式中

4.7.3特征分解法谱估计YTA=WTA(4.7.15)102用向量Y左乘(4.7.15)式,并取数学期望,得到

E［YYT］A=E［YW

T］A

(4.7.16)式中

用向量Y左乘(4.7.15)式,并取数学期望,得到E［103将上面关系式带入(4.7.16)式,得到(4.7.17)式中,σw2是自相关函数Ryy的特征值;A是对应σw2的特征矢量。

由于x(n)与w(n)互不相关，由（4.7.12）式求y(n)的自相关函数,得到

式中Rxx是x(n)的自相关函数,可以推导出在特征方程(4.7.17)式中,σw2是Ryy的最小特征值,且Ryy的阶数超过(2P+1)(2P+1)时，σw2就是其多重最小特征值。这一结论为我们寻求向量A提供了重要的依据。

(4.7.18)将上面关系式带入(4.7.16)式,得到(4.7.17)104当特征向量A求出后，就可以通过解特征方程，解出各个根，求出各正弦波的频率值，特征方程为

（4.7.19）

该方程有2P个根，这些根在z平面的单位圆上，它们是

i=1,2,3,…,2P

（4.7.20）

式中的ωi即是正弦波的频率。当各正弦频率由(4.7.19)式求出后,各正弦波的功率也可以求出,对于白噪声中P个实正弦组合信号,它的自相关函数为

（4.7.21）

当特征向量A求出后，就可以通过解特征方程，解出各个根105m≠0(4.7.22)式中

qi是频率为ωi的正弦波的幅度，Pi是其功率。将(4.7.22)式写成矩阵形式:FP=r(4.7.23)式中

(4.7.24)m≠0(4.7.22)式中qi是频率为ωi的正弦波的幅106PT=［P1,P2,…,PP］

(4.7.25)(4.7.26)由于各正弦波的频率已求出,矩阵F,r已知,由(4.7.23)式解出各正弦波的功率或幅度。

最后噪声功率由下式求出:(4.7.27)以上就是皮萨论科谱分解法的全过程。PT=［P1,P2,…,PP］(4.7.25)(4107Prony谱分析法

Prony曾在1975年提出用指数函数的一个线性组合去描述等间距采样数据的数学模型，其方法并不是通常意义下的谱估计技术。这里介绍有Kay和Marple（1981年）提出的可用于估计非有理式谱密度的扩充的Prony方法。基本的Prony法是一种内插法。它是采用等间隔复指数值的线性组合来拟合观测数据。每个复指数包含幅度和指数因子两部分，因此p个复制数值就含有2p个代求因子，若用p个指数项去逼近2p个数据，这种逼近可以精确实现。Prony谱分析法Prony曾在1975年提出用指数108

当用p个指数去逼N(N>2p)个数据时，可用最小二乘法准则来实现Prony方法的谱估计但这是一种非线性最小二乘法问题，直接求解相当困难，常用迭代法求解。经过适当的扩展后，Prony方法可用来计算有理式功率谱密度，本节介绍这种扩展Prony方法。扩展Prony方法采用的数学模型为一组p个具有任意幅值、相位、频率与衰减因子的复指数线型组合来拟合观测数据，其离散时间的函数形式为Prony谱分析法当用p个指数去逼N(N>2p)个数据时，可用最小109现代信号处理功率谱估计课件110现代信号处理功率谱估计课件111现代信号处理功率谱估计课件112现代信号处理功率谱估计课件113现代信号处理功率谱估计课件114现代信号处理功率谱估计课件115现代信号处理功率谱估计课件116现代信号处理功率谱估计课件117现代信号处理功率谱估计课件118现代信号处理功率谱估计课件119多重信号分类MUSIC

--现代谱估计的应用多重信号分类MUSIC

--现代谱估120

MUSIC为：MultipleSignalClassification即多重信号分类。MUSIC方法是一种估计空间参数的现代谱计算方法，它将功率谱推广为空间谱，是最早问世的子空间方法。这种算法是通过对数据协方差矩阵进行本征分解获得信号空间谱估计的。MUSIC为：121

当空间存在多个信号源时，常常用要对些空间信号进行分离，以便跟踪或检测我们感兴趣的空间信号，抑制那些被认为是干扰的空间信号。为此，需要使用天线阵列对多个空间信号进行接收。对天线阵列接收的空间信号所进行的分析与处理统称为阵列信号处理。

当空间存在多个信号源时，常常用要对些空间信号进行分离，以122

空间谱是阵列信号处理中的一个重要概念。时域频谱表示信号在各个频率上的能量分布，而空间谱则是表示信号在空间各个方向上的能量分布，它能够对多个空间信号进行识别（即分类）。空间谱是阵列信号处理中的一个重要概念。时域频谱表示信号在123波束形成器特征子空间分析MUSIC算法及其改进波束形成器1241.波束形成器

最佳波束形成器的基本原理是使来自非期望波达方向（DOA－directionofarrival）的任何干扰所贡献的功率最小，但是又能够保持“在观测方向上的信号功率不变”，也就是说通过将个阵元输出进行加权求和，在一时间内将天线阵列波束“导向”到一个方向上，通过对期望信号得到最大输出功率即给出波达方向的估计，因此这个也可以看成为一个尖锐的空间带通滤波器。MUSIC算法作为一种经典的高分辨的波达方向估计算法，与传统测向方法相比，它具有很高的分辨率1.波束形成器最佳波束形成器的基本原理是使来自非期望波125

在空间传播的多个信号源（窄带信号），现利用一天线阵列对这些信号进行接收。每根天线为一个阵元，它们都是全向天线，即无任何方向性。这里假定各个阵元等间距地直线排列。如下图：在空间传播的多个信号源（窄带信号），现利用一天线阵列对这126窄带包络变化缓慢，所以等距线阵各阵元接收到的同一信号的包络相同。假设为远场信号，则到达各阵元的方向角相同，用表示，称为波达方向（角），定义为信号到达阵元的直射线与阵列法线方向之间夹角。对一般的远场信号而言，同一信号到达不同的阵元存在一个波程差，这个波程差导致各接收阵元间的相位差。以上图的阵元1作为基准点，令信号传播延迟引起的相位差为，则由上图易知，波达方向与相位差之间的关系为因此，信号在第k个阵元上的接收信号为窄带包络变化缓慢，所以等距线阵各阵元接收到的同127

假设阵列有m个阵元组成，则信号到达各阵元的相位差所组成的向量为称为信号的方向向量。若共有p个信号位于远场（p<m），则在第k个阵元上的接收信号为其向量形式可以写成假设阵列有m个阵元组成，则信号到达各阵元的相位差128其中分别为m*p的方向矩阵和p*1的信号向量。具有3.9.6结构的是Vandermonde矩阵，当各不相同时，则列线性独立。在阵列信号处理中，一次采样称为一次快拍。假定每个阵元上共观测到N次快拍的接收信号,…,。波束形成问题就是依靠这些观测值，求得某个期望信号的波达方向。其中129解决问题：先假定信号与各阵元上的观测噪声独立统计，并且各观测噪声具有相同的方差。此时考察第k个信号的接收，那其它信号全部视为干扰信号，应抑制。我们知道一个阵列经过加权求和后，可以将阵列的方向增益聚焦在一个方向上，相当于形成一个波束。为此设计权向量

解决问题：130得到输出信号则N次快拍阵列输出的平均功率为其中

得到输出信号131

我们的目的是保证来自某个确定方向的信号能正确被接收，而其它方向的信号或干扰被完全抑制。用数学表示就是，在保证所需方向的信号输出为为常数的条件下，使阵列的输出功率最小化，即上面的第二个式子表明只有第k个信号的波达方向可以做为期望抽取。我们的目的是保证来自某个确定方向的信号能正确被接收，而132

下面用Largange常数法求解上述优化问题。为此构造目标函数对求偏导得，从中得到使输出能量最小化的最佳波束形成器此式代入，可知所以，这一波束形成器是Capon于1969年提出的，称为最小方差无畸变响应波束形成器。下面用Largange常数法求解上述优化问题。为此构造目133

将的结果代入平均功率的表达式，则但在实际应用中由于不能确切知道信号的方向，所以只能通过扫描的方法得到曲线这就是Capon定义的“空间谱”，它将峰值对应的定为p个信号的波达方向。它是著名的极大似然谱估计，也称为最小方差法。将的结果代入平均功率的表达式，则1342.特征子空间分析窄带远场信号DOA数学模型为

对此模型，我们做以下假设：假设1：对于不同的值，向量相互线性独立假设2：加性噪声向量的每个元素都是零均值的白噪声，它们不相关，并且有相同的方差假设3：矩阵非奇异，即rank（P）＝pMUSIC是以特征子空间分析为基础的2.特征子空间分析窄带远场信号DOA数学模型为MUSIC是以135在假设1～3下，由模型的表达式很容易得到

由于，则为Herimitian矩阵。令其特征分解为,式中，U为特征向量构成的酉矩阵。

在假设1～3下，由模型的表达