您好,欢迎访问三七文档
当前位置:首页 > 商业/管理/HR > 资本运营 > 有代码功率谱估计Levinson递推法Burg递推法随机信号处理
功率谱估计随机信号处理学号:姓名:随机信号处理实验报告实验三功率谱估计1实验内容信号为两个正弦信号加高斯白噪声,各正弦信号的信噪比均为10dB,长度为N,信号频率分别为1f和2f,初始相位120,取1/0.2sff,2/sff取不同数值:0.3,0.25。sf为采样频率。分别用Levinson递推法和Burg法进行功率谱估计,并分析改变数据长度、模型阶数对谱估计结果的影响。2实验原理2.1Levinson递推法:自相关法——列文森(Lenvison)递推法是已知信号观测数据,估计功率谱。它的出发点是选择AR模型参数使预测误差功率最小。假设信号xn的数据区在01nN范围,有P个预测系数,N个数据经过冲激响应为0,1,piaip的滤波器,输出预测误差en的长度为Np,因此有预测误差功率为211200111NpNpppinnienxnaxniNNen的长度长于数据的长度,上式中数据xn在01nN以外补充零点,相当于对无穷长的信号加窗处理,会引入误差。上式对系数pia的实部和虚部求微分使预测误差功率最小,得12ˆˆˆˆ0111ˆˆˆˆ1022ˆˆˆˆ120pxxxxxxxxpxxxxxxxxppxxxxxxxxarrrprarrrprarprprrp(Yule-Walker方程)式中自相关函数采用有偏自相关估计,即随机信号处理实验报告Levinson-Durbin算法:使一种按阶次递推的算法。它以0AR和1AR模型参数作为初始条件,计算2AR模型参数;再用2AR模型参数计算3AR模型参数,k阶模型参数由k-1阶模型参数计算得到。一直计算出ARp模型参数为止。一阶AR模型1p的Yule-Walker方程为2111011100xxxxxxxxrrrra由该方程解出1110xxxxrar2211,110xxar然后令2,3,4,p,以此类推,可以得到一般递推公式如下:pk称为反射系数,1pk。22212p,随着阶数增加,预测误差功率将减少或不变。由k=1开始递推,递推到k=p,依次得到各阶模型参数,2221112122212,,,,,,,,pppppaaaaaaAR模型的各个系数及模型输入白噪声方差求出后,信号功率谱用下式计算随机信号处理实验报告222211/1pjwjwjwixxwwpiipeHeae这种方法递推效率高,当阶数变化时,无需从头计算。但需要预先估计出信号自相关函数,当观测数据长度较短时,估计误差较大,会出现谱峰频率偏移和谱线分裂;如数据很长,估计自相关函数较准确。2.2Burg递推法:Levinson-Durbin递推法需要由观测数据估计自相关函数,这是它的缺点。而伯格递推法则由信号观测数据直接计算AR模型参数。伯格递推法利用Levinson-Durbin递推公式,导出前向预测误差与后向预测误差,并按照使它们最小的原则求出pk,从而实现不用估计自相关函数,直接用观测数据得出结果。Burg递推法思想:借助格型预测误差滤波器,求前向、后向预测误差平均功率,选择pk使其最小,求出pk。之后,再利用Levinson-Durbin递推法求模型参数和输入噪声方差。设信号xn的观测数据区间:01nN,前向、后向预测误差功率分别用pf和pb表示,预测误差平均功率用p表示,公式分别为211NfpfpnpenNp2101NbpbpnenNp12ppfpb前向、后向观测误差公式分别为1pfppkknxnxnkae随机信号处理实验报告1pbppkknxnpxnpkae01pa上式中,信号项的自变量最大的是n,最小的是n-p,为了保证计算范围不超出给定的数据范围,在pf和pb计算公式中,选择求和范围为:1pnN。为求预测误差平均功率p最小时的反射系数pk,令0ppk,将前、后向预测误差的递推公式代入得11121211211NfbppnppNfbppnpnneeknneeBurg递推法求AR模型参数的递推公式总结:(1)12010,00NxxxxnrxnrN⌒⌒(2)00fbnxnenxne0,1,2,10,1,2,1nNnN…,…,(3)11121211211NfbppnppNfbppnpnneeknnee(4)211pppk(5),1,1,ppipippiaaka1,2,1ip…,(6),pppka(7)111111ffbppppfbbppppnnkneeennkneee1,2,1,1,2,2nppNnpppN…,…,随机信号处理实验报告3实验结果及分析3.1原始信号1222sin()sin()ssfnfnsff,观测信号xsw这里取1sf,10.2f,20.3f,145N,20M。3.2Levenson递推法3.2.1取1sf,10.2f,20.3f或20.25f,阶数20M不变,实验不同数据长度对功率谱估计的影响1)信号长度35N随机信号处理实验报告信号长度N=35,阶数M=20的功率谱估计2)信号长度145N信号长度N=145,阶数M=20的功率谱估计3)信号长度2000N随机信号处理实验报告信号长度N=2000,阶数M=20的功率谱估计分析:由以上三个实验对比,可以看出当观测数据长度较短时,估计误差较大,会出现谱峰频率偏移与谱线分裂;当数据很长时,估计自相关函数较准确,但计算量较大。3.2.2取1sf,10.2f,20.3f,信号长度100N不变,实验不同模型阶数对功率谱估计的影响1)阶数M=22)阶数M=4随机信号处理实验报告3)阶数M=84)阶数M=16分析:由以上几个实验对比,可以看出当阶次较低,会使谱估计产生偏移,随机信号处理实验报告降低分辨率;当阶次越高,分辨率越高;当阶次太高,会使估计误差加大,谱峰分裂。3.3Burg递推法3.3.1取1sf,10.2f,20.3f或20.25f,阶数20M不变,实验不同数据长度对功率谱估计的影响1)信号长度35N信号长度N=35,阶数M=20的功率谱估计2)信号长度145N信号长度N=145,阶数M=20的功率谱估计3)信号长度2000N随机信号处理实验报告信号长度N=2000,阶数M=20的功率谱估计分析:由以上三个实验对比,可以看出当观测数据长度较短时,估计误差较大,会出现谱峰频率偏移与谱线分裂;当数据很长时,估计自相关函数较准确,但计算量较大。频率越靠近的谱估计,需要的阶数越高。3.3.2取1sf,10.2f,20.3f,信号长度900N不变,实验不同模型阶数对功率谱估计的影响1)阶数M=42)阶数M=8随机信号处理实验报告3)阶数M=164)阶数M=28随机信号处理实验报告分析:由以上几个实验对比,可以看出当阶次较低,会使谱估计产生偏移,降低分辨率;当阶次越高,分辨率越高;当阶次太高,会使估计误差加大,谱峰分裂。3.4实验总结本次试验采用分别用Levinson递推法和Burg递推法进行功率谱估计,并分析改变数据长度、模型阶数对谱估计结果的影响。通过实验,学习了Levinson递推法和Burg递推法的基本原理和一般流程,和如何选择AR模型的阶次,并使用Matlab语言,编写源代码,完成实验过程。在实验过程中,分别设计了不同信号长度、不同的AR模型的阶次和不同频率组合而成的4组实验,并在实验后,分析对比了实验结果。5源代码5.1Levenson递推法%clearall;Clear;tic;%产生信号fs=1;%设采样频率为1N=100;%数据长度改变数据长度会导致分辨率的变化f1=0.2*fs;%第一个sin信号的频率,f1/fs=0.2f2=0.3*fs;%第二个sin信号的频率,f1/fs=0.2或者0.3M=60;%滤波器阶数的最大取值,超过则认为代价太大而放弃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);随机信号处理实验报告plot(s,'b'),xlabel('时间'),ylabel('幅度'),title('原始信号s');grid;subplot(2,1,2);plot(x,'r'),xlabel('时间'),ylabel('幅度'),title('观测信号x');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;%一阶模型参数var(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,p)^2)*var(p-1);%求方差FPE(p)=var(p)*(M+1+p)/(M+1-p);%求最终预测误差end%确定AR模型的最佳阶数min=FPE(1);%求出FPE最小时对应的阶数p=1;随机信号处理实验报告fork=2:MifFPE(k)minmin=FPE(k);p=k;endend%功率谱估计W=0.01:0.01:pi;%功率谱以2*pi为周期,又信号为实信号,只需输出0到PI%即可;he=ones(1,length(W));%length()求向量的长度fork=1:phe=he+(a(p,k).*exp(-j*k*W));endPxx=var(p)./((abs(he)).^2);%功率谱函数%画出功率谱F=W*fs/(pi*2);%将角频率坐标换算成HZ坐标,便于观察figure;plot(F,abs(Pxx))xlabel('频率/Hz'),ylabel('功率谱P'),title(['AR模型的最佳阶数p='int2str(p)]);grid;toc;5.2Burg递推法clearall;cleartic;%产生信号fs=1;%设采样频率为1N=900;%数据长度改变数据长度会导致分辨率的变化;f1=0.2*fs;%第一个sin信号的频率,f1/fs=0.2f2=0.3*fs;%第二个sin信号的频率,f1/fs=0.2或者0.3M=300;%滤波器阶数的最大取值,超过则认为代
本文标题:有代码功率谱估计Levinson递推法Burg递推法随机信号处理
链接地址:https://www.777doc.com/doc-2319645 .html