1、ch64现代功率谱估计1.方差性能不好，不是方差性能不好，不是Px(W W)的一致估计的一致估计2.平滑周期图和平均周期图改善了周期图的方差平滑周期图和平均周期图改善了周期图的方差性能，但却降低了谱分辨率和增大了偏差。性能，但却降低了谱分辨率和增大了偏差。3.可能使短序列的功率谱估计出现错误的结果可能使短序列的功率谱估计出现错误的结果将观测数据以外的数据一律视为零，与实际不符。将观测数据以外的数据一律视为零，与实际不符。假设信号是白噪声通过假设信号是白噪声通过LTI系统产生的。由观测数系统产生的。由观测数据估计据估计LTI系统模型的参数。最后由系统模型的参数。最后由LTI系统模型的参数得系统模
2、型的参数得出功率谱。出功率谱。hh hkyk22j)e()(W WW WHPy 输入白噪声的自相关函数输入白噪声的自相关函数 2nnR h h (AutoRegressive model)(111)(0zAzazHpnnn (Moving Average model)qlllzbzH11)(AutoRegressive-Moving Average model)()(1)(00zAzBzazbzHpnnnqlll 若输入白噪声的功率谱若输入白噪声的功率谱2)(W Wh h P则输出序列的功率谱为则输出序列的功率谱为2j22j)e()()e()(W Wh hW W W WW WHPHPx 若能确
3、定模型中各参数若能确定模型中各参数和和就可以求得功率谱就可以求得功率谱Px(W W)描述p阶AR模型的差分方程为)(1knkynakyppnh将上式两端同乘以ykm再求数学期望 )(1kmkyEnkymkyEnakymkyEppnh h )(1mRnmRnamRyyppnyhmRynmRy mRyh h)(1mRnmRnamRyyppnyh h 由于输入是白噪声信号 Rh hn=2 n mRmhmRyh hh h对于因果系统,p阶AR模型的自相关函数与模型参数的关系 pmmnmRnamRyppny,1,0)(21 Yule-Walker(Y-W)方程2mh Y-W 方程的矩阵表示pmmnmRn
4、amRyppny,1,0)(21 例:p=3 时的Y-W 方程 000)3()2()1(101231012210132102333 aaaRRRRRRRRRRRRRRRRyyyyyyyyyyyyyyyy 若已知若已知Ryn，由，由Y-W方程解出各参数方程解出各参数a3(1),a3(2),a3(3),2。则可由则可由AR模型参数获得功率谱模型参数获得功率谱Py(W W)的估计值。的估计值。一阶一阶Y-W方程的解方程的解0)1(10 1 1 0211aRRRRyyyy解此方程得 0 1)1(1yyRRa)1(1 0121aRRyy)1(1(021aRy 二阶二阶Y-W方程的解方程的解00)2()1
5、(10 1 2 1 0 1 2 1 02222aaRRRRRRRRRyyyyyyyyy 1 0 1 20)2(2222yyyyyRRRRRa211 1)1(2yyRaR 1 02 1 1 0)1(222yyyyyyRRRRRRa)1()2()1(121aaa)2(2)1(1 02222aRaRRyyy2122)2(1 a p阶阶Y-W方程的递推解方程的递推解21111)()(pyppnypnpRnapRpa)1,2,1()()()()(11pnnpapananapppp2122)(1 ppppa(1)计算自相关函数的估计值计算自相关函数的估计值(2)求解一阶模型参数关函数的估计值求解一阶模型参
6、数关函数的估计值0 1)1(1yyRRa)1(1(02121aRy L-D算法估计功率谱的步骤算法估计功率谱的步骤(4)求出功率谱估计求出功率谱估计21j2ARe)(1)(pnnpnaPWW(3)由递推算法求解由递推算法求解p阶模型参数阶模型参数 L-D算法估计功率谱的步骤算法估计功率谱的步骤21111)()(pyppnypnpRnapRpa)1,2,1()()()()(11pnnpapananapppp2122)(1 ppppax：进行功率谱估计的输入有限长序列；p:模型的阶数NFFT：DFT的点数；Fs：绘制功率谱曲线的抽样频率，默认值为1；Pxx：功率谱估计值；f：Pxx值所对应的频率点
7、Pxx,f=pyulear(x,p,NFFT,Fs)例：例：利用利用L-D算法进行谱估计算法进行谱估计一序列含有白噪声和两个频率的余弦信号，一序列含有白噪声和两个频率的余弦信号，)4.0cos()3.0cos(kkkkxh h 利用利用L-D算法估计该序列的功率谱。算法估计该序列的功率谱。N=128;p=40;NFFT=2048;Fs=2;n=0:N-1;randn(state,0);x=cos(0.3*pi*n)+cos(0.4*pi*n)+randn(size(n);P,f=periodogram(x,NFFT,2);Py,fy=pyulear(x,p,NFFT,2);subplot(21
8、1);plot(f,10*log(P);grid;title(Periodogram);axis(0 1-60 60);subplot(212);plot(fy,10*log(Py);grid;title(Yule);axis(0 1-60 60);00.10.20.30.40.50.60.70.80.91-60-40-200204060 Periodogram 00.10.20.30.40.50.60.70.80.91-60-40-200204060 Yule2212)2()1(1zaza例例:前向线性预测前向线性预测(p=2)2)2(1)1(22kyakyaky f2kykyke2)2(1
9、)1(22kyakyaky2)2(1)1(22f2kyakyakykemin2f2f2keEE 若若 2,1,0f2nnkykeE则估计的均方误差达到最小。则估计的均方误差达到最小。且且minf22f2f2kykeEkeEE2)2(1)1(22f2kyakyakyke 1f2kykeE0 1)2(0)1(1 22yyyRaRaR 2f2kykeE00)2(1)1(222yyyRaRaRf2f2kykeEE 2)2(1)1(022yyyRaRaR00)2()1(10 1 2 1 0 1 2 1 0f222EaaRRRRRRRRRyyyyyyyyy例例:后向线性预测后向线性预测(p=2)yk12
10、kyyk)2(1)1(2 22kyakyaky2 2b2kykyke)2(1)1(222kyakyaky)()()(b2bzYzEzA2122)1()2(zzaa)(12zAz由由1,0,0b2nnkykeE2minb22b2b2kykeEkeEE可得可得00)2()1(10 1 2 1 0 1 2 1 0b222EaaRRRRRRRRRyyyyyyyyy比较前向比较前向预测系统和后预测系统和后向向预测系统的预测系统的方程后可得方程后可得 系数系数a2(1),a2(2)满足的方程相同满足的方程相同 2b2f2 EE Ab(z)=z 2A(z)在计算相关函数估计时，对在计算相关函数估计时，对N个
11、观测数据以个观测数据以外的数据作零的假设，故谱估计误差较大。外的数据作零的假设，故谱估计误差较大。直接从观测的数据利用线性预测器的前向和直接从观测的数据利用线性预测器的前向和后向预测的总均方误差之和为最小的准则来估计后向预测的总均方误差之和为最小的准则来估计反射系数，进而通过反射系数，进而通过L-D算法的递推公式求出算法的递推公式求出AR模型的优化参数。模型的优化参数。2阶前向预测误差 2)2(1)1(22f2 kyakyakyke1阶后向预测误差)1(11b1kyakyke L-D算法的递推公式)1()2()1()1(1212aaaa 2)2(1)1()2()1(2121f2 kyakyaa
12、akyke 1)1(2)2(1)1(121 kyakyakyaky 1b12f1f2 keKkekeK2=a2(2)2阶预测器的反射系数阶预测器的反射系数 一般地 1b1f1f keKkekepppp同理可得后向预测误差的递推公式 1f1b1bkeKkekepppp Kp=ap(p)为p阶预测器的反射系数。1b1f1f keKkekepppp 1f1b1bkeKkekepppp 前向和后向均方预测误差的总和为前向和后向均方预测误差的总和为 2b2f1kekeEppNpkp 由由0221pbpbppfpfpNpkppKkekeKkekeKE可得可得 1 1221211111kekekekeKbp
13、fpNpkbpfpNpkp(1)确定初始条件确定初始条件b0f0kykeke 121020kyNNk (2)从从p=1开始迭代计算开始迭代计算:计算计算AR模型参数模型参数1122b12f11b1f11 kekekekeKppNpkppNpkp递推递推p阶均方误差阶均方误差2122)1(pppK(3)递推高一阶前、后向预测误差递推高一阶前、后向预测误差 111keKkekebppfpfp 111keKkekefppbpbp(4)若阶数小于若阶数小于p，则阶数加，则阶数加1，回到步骤，回到步骤(2)进行下进行下一次迭代，直到达到预定阶数一次迭代，直到达到预定阶数p。(5)估计功率谱估计功率谱21
14、j2ARe)(1)(pnnpnaPWWx：进行功率谱估计的输入有限长序列；p:模型的阶数NFFT：DFT的点数；Fs：绘制功率谱曲线的抽样频率，默认值为1；Pxx：功率谱估计值；f：Pxx值所对应的频率点N=512;NFFT=1024;Fs=2;p=40;n=0:N-1;randn(state,0);x=cos(0.3*pi*n)+cos(0.32*pi*n)+randn(size(n);P,f=pyulear(x,p,NFFT,2);Pw,f2=pburg(x,p,NFFT,2);subplot(211);plot(f,10*log(P);grid;title(L-D);axis(0 1-30 60);subplot(212);plot(f2,10*log(Pw);grid;title(Burg);axis(0 1-30 60);00.10.20.30.40.50.60.70.80.91-200204060L-D00.10.20.30.40.50.60.70.80.91-200204060Burg00.10.20.30.40.50.60.70.80.91-200204060L-D00.10.20.30.40.50.60.70.80.91-200204060Burg
