1、功率谱估计功率谱-信号功率在频域的分布规律内容 5.1 确定信号的谱估计 5.2 平稳随机信号的功率谱估计 5.3 参数功率谱密度估计 5.4 基于子空间特征值分析的功率谱估计（高分辨谱估计）5.1 确定信号的谱分析连续时间信号首先通过一个低通（抗混叠）滤波器，然后采样连续时间信号首先通过一个低通（抗混叠）滤波器，然后采样得到离散时间信号，选择桢长为得到离散时间信号，选择桢长为N交叠为交叠为N0的采样数据，然后的采样数据，然后加窗，最后用加窗数据的一个合适长度的加窗，最后用加窗数据的一个合适长度的DFT作为谱估计作为谱估计(1)时域采样和抗混叠滤波2( )( )()( )jftjj mmX f
2、x t edtX ex meFTDTFT：(2) 频域采样和时域添零操作频域采样:DTFT-DFT1(2 / )0()( )( )( )jj mmNjN mkmX ex meX kx me补零操作10)/2()()(, 00),()(MmmkMjeemxkXMmNNmmxmx 通过对采样信号后面补零来得到更高密更高密度的频谱度的频谱。低密度离散频谱，在使用线性内插重建连续频谱时会产生误差。提高频谱密度可以减少这个误差。 (3) 能量泄漏，分辨率损失和加窗操作NR( )( )( )x nx n w n可用数据=完整数据矩形窗()R1()()()()()d2jjjjjNXeX eW eX e W
3、e 傅利叶变换：-(1)/2()( )jjNRW eAe矩形窗幅频特性矩形窗幅频特性分辨率损失 1(2 / )0( ) cos0.35cos0.40.25cos0.8( )( ) ( )NjN mkmx mmmmX kx mwme可以看出频率分辨率决定于数据窗的持续时间（采样点数-1）窗口函数对分辨率和泄漏的影响v频率分辨率频率分辨率：窗口函数频谱的主瓣宽度。v泄漏：旁瓣谱峰的水平（导致虚假谱峰出现） 因此要根据实际情况选择不同的窗口。矩形窗具有最窄矩形窗具有最窄的主瓣宽度。旁瓣振幅的减的主瓣宽度。旁瓣振幅的减小只能以降低分辨率为代价小只能以降低分辨率为代价。利用矩形窗口分辨两个频率，应当使两
4、个频率的绝对差值大于矩形窗频谱的主瓣宽度：|w1-w2|mainWin各种不同窗口的属性比较和实域频域特性窗口旁瓣峰值（db）近似主瓣宽度精确主瓣宽度矩形窗-134pi/(N-1)1.81pi/(N-1)汉宁窗-328pi/(N-1)5.01pi/(N-1)汉明窗-438pi/(N-1)6.27pi/(N-1)凯撒窗-A(A-8)pi/(2.285N-1)多尔夫-切比雪夫窗-A5.2 平稳随机信号的谱估计 功率谱为自相关函数的Fourier变换 Wiener-Khintchine 定理 自相关的估计值 估计值的均值与方差mfmjxxxxemrfP2)()(1)Wiener-Khintchine
5、 定理1936年222)()(*)()()(mfmjxxemxfXfXfXfPnxxmnxnxmnxnxEmr)()()()()(2222()*( )( ) ()( )()( )( )( )jfmjfmxxmmnjfnjf n mxxnnrm ex n x nm ex n ex nm eXf X fPf (2) 自相关估计1, 2 , 1 , 0, )()(|1)(1|0NmkmxkxmNmrmNkxxN | | 101 ( )( ) ()mxxkrmx k x kmNunbiasedBiasedNonnegative definiteSmaller variance- 平稳随机信号的谱估计估
6、计估计1：估计估计2：估计量的均值和方差)(|1)()(1)(1|0mrNmmkxkxENmrExxmNkxxkxxxxxxxxmkrmkrkrNmrVar)()()(1)( 2Mean:Variance:是一个渐进无偏5.2.2 非参数功率谱密度估计方法 周期图法 Bartlett法（平均多个周期图, 采用不同数据块） Welch 法 (平均多个周期图, 采用重叠的数据块) Blackman-Tukey 法 (周期图平滑)(1) 周期图法 Periodogram Sir Arthur Schuster in 1899 Mean Variance22102| )(|1)(1)(fXNemxNf
7、PNmfmjxxefmjNNmxxxxmrNmfXENfPE21)1(2)(|1)(1)(efmjmxxxxNmrfPE2)()(lim222sin2sin1)()(fNfNfPfPVarxxxx一个无偏但不一致的估计An unbiased but not a consistent estimate12(1)( )( )NjfmxxxxmNPfr me(2) Bartlett 法平均多个不同数据块的周期图估计结果KiiXXBxxfKfPP1)()(1)(Mean:dvvfNvfvNfEfEPPPXXiXXBXX2/12/12)()(sin)(sin)(1)()(An asymptoticall
8、y unbiasedVariance:222sin2sin1)(1)(fNfNfPKfPVarXXBXX1/k of the variance of the periodogram(3) Welch 法 采用有重叠的数据块Mean:dvvfNvfvNfEfEPPPXXiXXBXX2/12/12)()(sin)(sin)(1)()(An asymptotically unbiasedVariance:222sin2sin1)(1)(fNfNfKfVarPPXXBXX1/k of the variance of the periodogramWelch-Bartlett 方法(4) Blackma
9、n-Tukey 法 加窗的谱估计Mean:Variance:Uw/k of the variance of the periodogramUw相关窗口的能量fmjNNmxxBTXXemrmwfP21)1()( )()(1)1(2)()()(NNmfmjcxxBTXXemwmrfPdvvfWvPfPcXXBTXX2121)()()(2w( )( )BTXXXXUVar PfPfN相关窗谱为单位面积时，为渐进无偏估计Blackman-Tukey 法方法理论和实践的比较 对功率谱加窗平滑等价于对估计的自相关序列进行加窗。对窗口有一定要求：三角窗(Bartlett)和Parzen窗可以，但不适用于汉明
10、，汉宁或凯瑟窗。5.3 参数功率谱密度估计 Autoregressive (AR) Model（自回归） Moving-Averatge (MA) Model（移动平均） Autoregressive Moving-Average (ARMA) Model（自回归移动平均）估计步骤：估计步骤： (1) 信号建模信号建模 (2) 估计模型参数估计模型参数 (3) 计算功率谱计算功率谱本节内容 5.3.1 信号模型及其功率谱 5.3.2 AR,MA,ARMA 模型与 他们的功率谱 5.3.3 AR 模型的功率谱估计 Yule-Walker 方程 AR 模型特性 参数估计方法 5.3.4 MA模型的
11、功率谱估计 5.3.5 ARMA模型的功率谱估计平稳随机信号平稳随机信号白噪声白噪声线性移不变系统线性移不变系统有理函数模型有理函数模型:5.3.1 信号模型及其功率谱2222jjjxxjB eSeH eA e(1) Autoregressive-Moving Average (ARMA) 模型 000, 1qkkkpkkkb zB zH zaA za z 10pqkkkkx na x nkb u nk 22111xxB z BzSzH z HzA z Az输入输出关系:系统方程:功率谱密度:零零-极点模型极点模型 111111 11 real is 1 : 1 then , real is
12、If zHerHerHzHerznhrezerzrezzHzHnhjjjjjj 2122xxjjxxSzB z B zSeB e(2) Moving Average 模型 0qkMAkkHzB zb z 0qkkx nb u nk输入输出关系:系统方程:功率谱密度:全零点模型全零点模型全极点模型全极点模型 212211xxjxxjSzA z A zSeA e(3) Autoregressive 模型 011ARpkkkHzA za z 1pkkx na x n ku n输入输出关系:系统方程:功率谱密度:模型模型:模型模型:模型模型: pAR qMAqpARMA,21,paa21,qbb,1p
13、aa21,qbb模型参数模型参数 k 1100 xxpkpkxxnllRmEx n x nmEx na x nmku nma RmkEx n u nmEx n u nmEh l u nlu nmh l E u nl u nm 220 lh lmlhm 5.3.2 AR 模型参数估计模型参数估计 222212000 ,00,0 causal 0 ( 0 ): 0lim1,00,0,00,0,0 (1)0,0zpxxkxxkpkxxkhmE x n u nmhmmhmmhH zmE x n u nmmmRma Rmkmma Rmkam 初值定理 1)( , 1, 0 0,1)( 0, 00,0 ,
14、0,020020121apmmkmRaammkmRamkmRamkmRamRpkxxkpkxxkpkxxkpkxxkxx（1）Yule-Walker 方程 001011011021pxxxxxxxxxxxxxxxxxxaaRpRpRpRRRpRRR 00011011pxxxxxxxxxxxxaaRpRpRpRRR20 0pkxxka Rmkm求解上面方程式，即得到求解上面方程式，即得到AR模型参数模型参数5.3.3 MA和ARMA模型的功率谱估计ARMA 模型的模型的 Yule-Walker 方程方程 2101,0, 1, pq mkxxk mkkxxpkxxkaR m kb h kmqR m
15、aR m km q MA 模型的模型的 Yule-Walker 方程方程 , 0, 1, 0,02qmqmbbmRmqkkmkxx（*1） , qmxxxxmqSzRm zmqMA 模型的功率谱估计 不需要估计模型参数不需要估计模型参数bk, 只需根据已知数据估计出只需根据已知数据估计出|m|q的自相关函数，代入上式计算即可。的自相关函数，代入上式计算即可。 （*2）( )x n()/ ()Bz Az( )e n( )A z( )v n( )B z(1) 首先计算首先计算AR模型参数。模型参数。（*1）式）式(2) 利用利用AR模型对模型对x(n)滤波得到滤波得到v(n)(3) 利用利用v(n
16、)计算计算MA模型功率谱。模型功率谱。（*2）式）式211vvxxpjkkkSSa eARMA模型的功率谱估计5.4 基于子空间特征值分析的功率谱估计目标信号目标信号： 已知在白噪声中的已知在白噪声中的M个负指数序列和的个负指数序列和的N 个个采样值采样值 11, 0,1iMMjniciiix ns nv nA ev nnN, 1,ijciiAA eiM 和 需要估计需要估计. where are uncorrelated random variables that uniformly distributed over the intervali2 , 0 nvi例例1: 白噪声中单个复指数序
17、列白噪声中单个复指数序列Signal autocorrelation matrixNoise autocorrelation matrix 111111111 12c11121 111 112, 0, 1,1;011,1011, EE jnjcccTTj NjTcHHHHvHvx nA ev nnNAAenxxx Nneenvvv NnAnnnnAnnnnPnnPnnXSVSVxevxevRxxeevveeIRRReeRI TNjjeen11111e nnPH111eeRS11111111121211NjNjNjjNjjeeeeeePSR信号矢量：信号矢量：信号的自相关矩阵：信号的自相关矩阵：
18、因为矩阵因为矩阵 的秩为的秩为1，所以仅有一个非零特征值，所以仅有一个非零特征值SR 的非零特征值的非零特征值:SR111111111111111: e0eIReeeeeeeeRSSNPNPNPPPHHNvvv,32-对应的特征向量对应的特征向量- 的非零特征值的非零特征值SR 是厄米共轭矩阵，所以其他的特征向量是厄米共轭矩阵，所以其他的特征向量与与e1 正交。正交。 SRNiiH, 3, 2 , 01veIRV2v噪声自相关矩阵噪声自相关矩阵222222 000000vvvvvvdiagVR是个满秩矩阵是个满秩矩阵VR设设 为信号自相关矩阵为信号自相关矩阵 的特征值，则的特征值，则iSR22
19、1max22222 vvviiviiviiiviiviNPvvvvvRvIRvRSSX 的特征值的特征值 :XR- 的特征值的特征值XR- 的最大特征值的最大特征值XR- 的其他特征值的其他特征值 XR的特征向量与的特征向量与 的相同，为的相同，为XRSRiv从从 的特征值和特征向量中提取信号的特征值和特征向量中提取信号 参参数的计算步骤数的计算步骤:1. 对自相关矩阵对自相关矩阵 进行特征值分解。其最大进行特征值分解。其最大特征值等于特征值等于 ，其他特征值等于，其他特征值等于2. 使用这些特征值计算功率使用这些特征值计算功率 和噪声方差和噪声方差XR nsXR21maxvNP2v1P2v计
20、算步骤minmax1min21NPv3. 从最大特征值所对应特征矢量从最大特征值所对应特征矢量 确定信号确定信号频率频率 例如，例如，1maxv 1argmax1v 111111max121maxjTNjjjeveeeev令令 为自相关矩阵为自相关矩阵 的噪声特征矢量，即具有的噪声特征矢量，即具有特征值特征值 的一个特征矢量；并且令的一个特征矢量；并且令 为特征矢为特征矢量量 的第的第 i个成份。个成份。频率估计方程频率估计方程 :ivXR2v nviiv 2210v11iHNnnjijienvePe正交条件正交条件: NiiH, 3, 2 , 01ve频率估计方程求取不同频率点处的上述方程值
21、。求取不同频率点处的上述方程值。分母在分母在 处将趋于处将趋于0。因此频率方程在。因此频率方程在 处将趋于无穷大。这样，理论上讲频率方程的处将趋于无穷大。这样，理论上讲频率方程的峰值位置可以用来估计指数序列的频率。然而，峰值位置可以用来估计指数序列的频率。然而，由于这个方法仅使用了一个特征向量，因此可由于这个方法仅使用了一个特征向量，因此可能对于矩阵能对于矩阵 的估计误差比较敏感。我们可以的估计误差比较敏感。我们可以使用对所有噪声特征矢量的加权平均来代替单使用对所有噪声特征矢量的加权平均来代替单个特征矢量。个特征矢量。11XR 221011vjiNHjniinPevn ee1221vjNiHi
22、iiPee例例2: 白噪声中两个复指数序列白噪声中两个复指数序列 22121 ; ; 2, 1 ,1, 1, 0 21vjicinjcnjcieAANnnveAeAnxi 2211122221 1122221 11222E E , HHHHHHvHHvnnAnnAnnnnPnnPnnPnnPnnXSVRxxeeeevveeeeIReeeeRI2, 1212 , PPdiagvHPeeEIEPERX为更精确描述上面分解，可以使用为更精确描述上面分解，可以使用矩阵形式：矩阵形式： 1 12 2ccnAnAnnxeevP1,P2分别为第一个和第二个复正弦波的功率。分别为第一个和第二个复正弦波的功率。
23、令令 和和 为矩阵为矩阵 的特征向量和特征值，的特征向量和特征值，并且把特征值按照降序排列：并且把特征值按照降序排列：ivXR因因 , 所以所以IRRSX2v2vixi 22220vIRvvvvvRvIRvRSSSXiiiviiviiiviivixixNxx21为特征值为特征值 ofiSR 由于信号自相关矩阵由于信号自相关矩阵 秩为秩为2，所以，所以只有两个非零特征值，并且他们都大于只有两个非零特征值，并且他们都大于零（因零（因 为非负定）。这样为非负定）。这样 矩阵矩阵 的的特征向量和特征值可以分为两个部分：特征向量和特征值可以分为两个部分：SRSRxR第一部分包含大于第一部分包含大于 的两
24、个特征值和对应的特征向量（称的两个特征值和对应的特征向量（称为信号特征向量）。两个向量张成一个子空间为信号特征向量）。两个向量张成一个子空间信号子信号子空间空间。第二部分包含那些等于第二部分包含那些等于 的两个特征值和对应的特征向量的两个特征值和对应的特征向量（称为噪声特征向量）。噪声特征向量张成一个（称为噪声特征向量）。噪声特征向量张成一个N-2维子空维子空间间噪声子空间噪声子空间。上面的定义有一点误导：因为噪声成份同时影响信号子空间上面的定义有一点误导：因为噪声成份同时影响信号子空间和噪声子空间和噪声子空间2v2v 因因 厄米共轭，特征向量厄米共轭，特征向量 相互正相互正交。因此，信号空间
25、和噪声空间是正交交。因此，信号空间和噪声空间是正交子空间。也就是说，对信号子空间中的子空间。也就是说，对信号子空间中的任一向量任一向量 和噪声子空间中的任一向量和噪声子空间中的任一向量 有下面成立：有下面成立：XRuv0vuHiv不像单个复指数序列的例不像单个复指数序列的例子，这里信号特征向量通子，这里信号特征向量通常不再等于常不再等于 和和 . 1e2e1v2viv1e2e2e1e然而，然而， 和和 同样位于由同样位于由 和和 张成的信号张成的信号子空间内。由于信号子空间和噪声子空间正子空间内。由于信号子空间和噪声子空间正交，那么交，那么 和和 同样与噪声特征矢量同样与噪声特征矢量 正正交。
26、（交。（i2）NiNiiHiH, 4, 3 , 0, 4, 3 , 021veveNiiHijieP32v1e我们仍然可以使用上面的频率方程得到对我们仍然可以使用上面的频率方程得到对两个频率值两个频率值的估计。的估计。通用情况通用情况: 一个广义平稳过程，在白噪一个广义平稳过程，在白噪声中包含声中包含M个不同的复指数序列个不同的复指数序列MTMvHTNjjjivHiiMiiPPdiagMieeePiii1121221, 2, 1 1PeeEIEPEReIeeRRRXVSX M个线性独个线性独立的向量立的向量 信号向量组成的信号向量组成的NM矩阵矩阵关于各个信号能量的对角阵关于各个信号能量的对角
27、阵这里这里 为矩阵为矩阵 的特征值。的特征值。 由于矩阵由于矩阵 的秩为的秩为 M ，所以，所以 前前M个特征值将大于个特征值将大于 ，后，后 个特征值将等于个特征值将等于 。因矩阵因矩阵 的特征值为的特征值为XR2vixiiSRXR2vMN2vSR矩阵矩阵 的特征值和特征向量可以分为两的特征值和特征向量可以分为两个部分：个部分：XR2. 噪声特征向量噪声特征向量MivixiM, 1 , ;,21vvNMivxiNM, 1 , ;,21vv1. 信号特征向量信号特征向量假设特征向量已经被模归一化假设特征向量已经被模归一化, 我们可以我们可以以下面形式对矩阵以下面形式对矩阵 进行分解：进行分解：
28、221112121 , ,vviNMMHHHiiNMivHiiMiviHiiNixidiagdiagMNNMNVSVSVVVSSSXXvvVvvVVVVVRvvvvvvRXR 所有信号向量所有信号向量 都位于信号子都位于信号子空间内。由信号子空间和噪声子空间的正空间内。由信号子空间和噪声子空间的正交性可以推出，信号向量交性可以推出，信号向量 正交于任何一正交于任何一个噪声特征向量：个噪声特征向量：Mee,1ie , 1 , 1 , 0NMkMikHive 信号的频率值可以使用频率估计方程信号的频率值可以使用频率估计方程进行估计：进行估计：NMiiHijeP12v1e 有两种不同类型的频率估计方
29、法基于以有两种不同类型的频率估计方法基于以上的频率估计方程上的频率估计方程:(1) Pisarenko 谐波分解谐波分解(2) MUltiple SIgnal Classification（MUSIC）(2) 基于频率估计方程的方法思想：信号频率值可以从自相关矩阵思想：信号频率值可以从自相关矩阵的对应于最小特征值的特征向量处估的对应于最小特征值的特征向量处估计得到。计得到。Pisarenko Harmonic Decomposition ( PHD 方法方法 ) PHD 方法缺点在于对于噪声敏感（由方法缺点在于对于噪声敏感（由于仅使用了一个特征向量），这限制了它于仅使用了一个特征向量），这限制
30、了它的广泛使用的广泛使用. 假设假设: 1. 信号中复指数序列的数目信号中复指数序列的数目M为已知为已知2. 个自相关序列的采样已知或者可个自相关序列的采样已知或者可以被估计出来以被估计出来1M当不知道复指数序列的确切数目时，使用这个方法需要格外小当不知道复指数序列的确切数目时，使用这个方法需要格外小心。心。 对于一个对于一个M+1M+1 维的自相关矩维的自相关矩阵阵 ，噪声子空间的维数显然为，噪声子空间的维数显然为1，噪声，噪声子空间是被对应于最小特征值的特征向量子空间是被对应于最小特征值的特征向量 所张成。所张成。 将与每一个信号向量正将与每一个信号向量正交：交：XRMievMkjkHii
31、, 1 , 00minminveminvminv这样对这个特征向量系数的傅利叶变换这样对这个特征向量系数的傅利叶变换在每一个复指数序列的频率点在每一个复指数序列的频率点 处取值为处取值为0.i1 , Mi MkjkjekveV0minmin MiekvMkjkHii, 1 , 00minminve相应的，噪声矢量的相应的，噪声矢量的z变换具有变换具有M个零点个零点在单位圆上在单位圆上 MkjkMkzezkvzVk110minmin1 与求取与求取 的零点相似，也可使用的零点相似，也可使用这是频率估计方程的一个特殊形式，这是频率估计方程的一个特殊形式， with and . zVmin2min1
32、veHjPHDeP1 MN11MjPHDeP 峰值点所对应的频率被作峰值点所对应的频率被作为复指数序列的频率估计为复指数序列的频率估计.jPHDeP尽管写为功率谱的形式，尽管写为功率谱的形式， 被叫做伪谱被叫做伪谱（或特征谱）。因为它不包含任何关于复指（或特征谱）。因为它不包含任何关于复指数序列或者噪声功率的信息。数序列或者噪声功率的信息。如何估计噪声和复指数序列的功率呢？如何估计噪声和复指数序列的功率呢？假设假设: 信号子空间的特征向量信号子空间的特征向量 已经被规范化即已经被规范化即,1Mvv1iHivvMiiii, 1 , vvRXMiiiHiiiHi, 1 , vvvRvX功率估计功率
33、估计对下式两边都左乘一个矢量对下式两边都左乘一个矢量得到得到MieVPMiPPMiviMkjikviMkiHkkiiMkvHkkkHiiiHiiiHik, 1 , , 1 , , 1 , 21221212vevIeevvvvRvX2222121222222222212121212 , 1 , 22121vMvvMjMjMjMjjjjjjviMkjikPPPeVeVeVeVeVeVeVeVeVMieVPMkMMkEquation *注意，除注意，除P1,P2, ,PM之外，其他参数如各个频率，以及噪声之外，其他参数如各个频率，以及噪声方差都已经求得。求解这个方程即可得到各个复指数序列的功方差都已
34、经求得。求解这个方程即可得到各个复指数序列的功率。率。Step 1: 对于给定的一个白噪声中对于给定的一个白噪声中M个复指个复指数序列的随机过程，找到其自相关矩阵数序列的随机过程，找到其自相关矩阵 的最小特征值的最小特征值 和对应的特征向量和对应的特征向量 。Step 2: 令白噪声功率为最小特征值令白噪声功率为最小特征值 。令复指数序列频率等于特征向量令复指数序列频率等于特征向量 的的z变变换换xminminvXRxmin计算步骤 kMkzkvzV0minmin最靠近单位圆的最靠近单位圆的M个零点的角度个零点的角度或者下面频率估计方程的或者下面频率估计方程的M个峰点所对应的个峰点所对应的频率频率2min1veHjPHDePStep3: 计算复指数序列的功率。求解现行方计算复指数序列的功率。求解现行方程组程组 equation*
