功率谱估计报告

1功率谱估计一、实验内容信号为两个正弦信号加高斯白噪声，各正弦信号的信噪比均为10dB，长度为N，信号频率分别为1f和2f，初始相位120，取1/0.2sff，2/sff取不同数值：0.3，0.25。sf为采样频率。分别用Levinson递推法和Burg法进行功率谱估计，并分析改变数据长度、模型阶数对谱估计结果的影响。二、实验原理（一）Levinson递推法：自相关法——列文森（Lenvison）递推法是已知信号观测数据，估计功率谱。它的出发点是选择AR模型参数使预测误差功率最小。假设信号xn的数据区在01nN范围，有P个预测系数，N个数据经过冲激响应为0,1,piaip的滤波器，输出预测误差en的长度为Np，因此有预测误差功率为211200111NpNpppinnienxnaxniNNen的长度长于数据的长度，上式中数据xn在01nN以外补充零点，相当于对无穷长的信号加窗处理，会引入误差。上式对系数pia的实部和虚部求微分使预测误差功率最小，得12ˆˆˆˆ0111ˆˆˆˆ1022ˆˆˆˆ120pxxxxxxxxpxxxxxxxxppxxxxxxxxarrrprarrrprarprprrp（Yule-Walker方程）式中自相关函数采用有偏自相关估计，即2Levinson-Durbin算法：是一种按阶次递推的算法。它以0AR和1AR模型参数作为初始条件，计算2AR模型参数；再用2AR模型参数计算3AR模型参数，k阶模型参数由k-1阶模型参数计算得到。一直计算出ARp模型参数为止。一阶AR模型1p的Yule-Walker方程为2111011100xxxxxxxxrrrra由该方程解出1110xxxxrar2211,110xxar然后令2,3,4,p，以此类推，可以得到一般递推公式如下：pk称为反射系数，1pk。22212p，随着阶数增加，预测误差功率将减少或不变。由k=1开始递推，递推到k=p，依次得到各阶模型参数，2221112122212,,,,,,,,pppppaaaaaa3AR模型的各个系数及模型输入白噪声方差求出后，信号功率谱用下式计算222211/1pjwjwjwixxwwpiipeHeae这种方法递推效率高，当阶数变化时，无需从头计算。但需要预先估计出信号自相关函数，当观测数据长度较短时，估计误差较大，会出现谱峰频率偏移和谱线分裂；如数据很长，估计自相关函数较准确。（二）Burg递推法：Levinson-Durbin递推法需要由观测数据估计自相关函数，这是它的缺点。而伯格递推法则由信号观测数据直接计算AR模型参数。伯格递推法利用Levinson-Durbin递推公式，导出前向预测误差与后向预测误差，并按照使它们最小的原则求出pk，从而实现不用估计自相关函数，直接用观测数据得出结果。Burg递推法思想：借助格型预测误差滤波器，求前向、后向预测误差平均功率，选择pk使其最小，求出pk。之后，再利用Levinson-Durbin递推法求模型参数和输入噪声方差。设信号xn的观测数据区间：01nN,前向、后向预测误差功率分别用pf和pb表示，预测误差平均功率用p表示，公式分别为211NfpfpnpenNp2101NbpbpnenNp12ppfpb前向、后向观测误差公式分别为41pfppkknxnxnkae1pbppkknxnpxnpkae01pa上式中，信号项的自变量最大的是n,最小的是n-p，为了保证计算范围不超出给定的数据范围，在pf和pb计算公式中，选择求和范围为：1pnN。为求预测误差平均功率p最小时的反射系数pk，令0ppk,将前、后向预测误差的递推公式代入得11121211211NfbppnppNfbppnpnneeknneeBurg递推法求AR模型参数的递推公式总结：(1)12010,00NxxxxnrxnrN⌒⌒(2)00fbnxnenxne0,1,2,10,1,2,1nNnN…，…，(3)11121211211NfbppnppNfbppnpnneeknnee(4)211pppk(5),1,1,ppipippiaaka1,2,1ip…，(6),pppka(7)111111ffbppppfbbppppnnkneeennkneee1,2,1,1,2,2nppNnpppN…，…，5三、实验结果及分析(一）成原始信号1222sin()sin()ssfnfnsff，观测信号xsw，这里取1sf，10.2f，20.3f，30N，20M。(二)Levenson递推法1、取1sf，10.2f，20.3f或20.25f，阶数20M不变，实验不同数据长度对功率谱估计的影响1.1)信号长度30N,1sf，10.2f，20.3f,阶数20M功率谱估计1.2)信号长度30N,1sf，10.2f,20.25f,阶数20M功率谱估计62、2.1)信号长度150N,1sf，10.2f，20.3f,阶数20M功率谱估计2.2)信号长度150N,1sf，10.2f，20.25f,阶数20M功率谱估计３、3.1)信号长度2000N,1sf，10.2f，20.3f,阶数20M功率谱估7计3.2)信号长度2000N,1sf，10.2f，20.25f,阶数20M功率谱估计分析：由以上三个实验对比，可以看出当观测数据长度较短时，估计误差较大，会出现谱峰频率偏移与谱线分裂；当数据很长时，估计自相关函数较准确，但计算量较大。2、取1sf，10.2f，20.3f，信号长度100N不变，实验不同模型阶数对功率谱估计的影响1）阶数M=282）阶数M=43）阶数M=894)阶数M=16分析：由以上几个实验对比，可以看出当阶次较低，会使谱估计产生偏移，降低分辨率；当阶次越高，分辨率越高；当阶次太高，会使估计误差加大，谱峰分裂。(三)Burg递推法1、取1sf，10.2f，20.3f或20.25f，阶数20M不变，实验不同数据长度对功率谱估计的影响1.1)信号长度30N,1sf，10.2f，20.3f,阶数20M功率谱估计10Elapsedtimeis0.513413seconds.1.2)信号长度30N,1sf，10.2f,20.25f,阶数20M功率谱估计Elapsedtimeis0.493875seconds.2、2.1)信号长度150N,1sf，10.2f，20.3f,阶数20M功率谱估计11Elapsedtimeis0.505008seconds.2.2)信号长度150N,1sf，10.2f，20.25f,阶数20M功率谱估计Elapsedtimeis0.465947seconds.3、3.1)信号长度2000N,1sf，10.2f，20.3f,阶数20M功率谱估计Elapsedtimeis0.583838seconds.3.2)信号长度2000N,1sf，10.2f，20.25f,阶数20M功率谱估计12Elapsedtimeis0.586020seconds.分析：由以上三个实验对比，可以看出当观测数据长度较短时，估计误差较大，会出现谱峰频率偏移与谱线分裂；当数据很长时，估计自相关函数较准确，但计算量较大。频率越靠近的谱估计，需要的阶数越高。３、取1sf，10.2f，20.3f，信号长度900N不变，实验不同模型阶数对功率谱估计的影响Ｍ=4M=8Elapsedtimeis0.493989seconds.Elapsedtimeis0.347844seconds.M=32M=6413Elapsedtimeis0.365664seconds.Elapsedtimeis0.403261seconds.M=256M=512Elapsedtimeis6.183147seconds.Elapsedtimeis6.275853seconds分析：由以上六个实验对比，可以看出当阶次较低，会使谱估计产生偏移，降低分辨率；当阶次越高，分辨率越高；当阶次太高，会使估计误差加大，谱峰分裂。当Ｍ值取得过大，运行时间过久，以致无法达到效果，即功率谱Ｐ无法得出。四、实验总结本次试验采用分别用Levinson递推法和Burg递推法进行功率谱估计，并分析改变数据长度、模型阶数对谱估计结果的影响。通过实验，学习了Levinson递推法和Burg递推法的基本原理和一般流程，和如何选择14AR模型的阶次，并使用Matlab语言，编写源代码，完成实验过程。在实验过程中，分别设计了不同信号长度、不同的AR模型的阶次和不同频率组合而成的4组实验，并在实验后，分析对比了实验结果。五、源代码（1）Levenson递推法clearall;cleartic;%产生信号fs=1;%设采样频率为1N=100;%数据长度改变数据长度会导致分辨率的变化；f1=0.2*fs;%第一个sin信号的频率,f1/fs=0.2f2=0.3*fs;%第二个sin信号的频率，f1/fs=0.2或者0.3M=16;%滤波器阶数的最大取值,超过则认为代价太大而放弃L=2*N;%有限长序列进行离散傅里叶变换前，序列补零的长度n=1:N;s=sin(2*pi*f1*n/fs)+sin(2*pi*f2*n/fs);%s为原始信号x=awgn(s,10);%x为观测信号，即对原始信号加入白噪声，信噪比10dB%画出原始信号和观测信号figure(1);subplot(2,1,1);15plot(s,'b'),xlabel('时间'),ylabel('幅度'),title('原始信号s');%axis([0150-33]);grid;subplot(2,1,2);plot(x,'r'),xlabel('时间'),ylabel('幅度'),title('观测信号x');%axis([0150-33]);grid;%计算自相关函数rxx=xcorr(x,x,M,'biased');%计算有偏估计自相关函数，长度为-M到M，共2M+1r0=rxx(M+1);%r0为零点上的自相关函数，相对于-M，第M+1个点为零点R=rxx(M+2:2*M+1);%R为从1到第M个点的自相关函数矩阵%Levinson递推算法%确定矩阵大小a=zeros(M,M);FPE=zeros(1,M);%FPE：最终预测误差，用来估计模型的阶次var=zeros(1,M);%求初值a(1,1)=-R(1)/r0;%一阶模型参数16var(1)=(1-(abs(a(1,1)))^2)*r0;%一阶方差FPE(1)=var(1)*(M+2)/(M);%递推forp=2:Msum=0;fork=1:p-1%求a(p,p)sum=sum+a(p-1,k)*R(p-k);enda(p,p)=-(R(p)+sum)/var(p-1);fork=1:p-1%求a(p,k)a(p,k)=a(p-1,k)+a(p,p)*a(p-1,p-k);endvar(p)=(1-a(p