《预测误差格型滤波器及伯格(精).docx》由会员分享，可在线阅读，更多相关《预测误差格型滤波器及伯格(精).docx（12页珍藏版）》请在金锄头文库上搜索。
1、输入：x(), = 0,1,2,3,1,阶数IP图3.3.5 Burg递推法流程图plot (PCC,k)；运行结果CDP)磐憩阱SN=1000时，周期图的结果N=100时，周期图的结果400204060801 001 20N=50时，周期图的结果作业：改变被分析数据的长度、AR模型的阶数、信噪比等参数，比较各种谱估计方法的特点!然后总结%采样频率设置（Hz）%信号的数据长度%产生高斯白噪声信号%信噪比%噪声信号的功率%由信噪比计算正弦信号的振幅%正弦序列%有用信号%正弦信号的功率%计算实际信躁比%输入信号=有用信号+高斯白噪声close all;clear;clc;fs= 10000：fl=
2、fs*0.2;f2=fs*0.25;N=100;t=O:l/fs:(N-l)/fs;randnCslate,0);wn=randn( 1 ,N);SNR=10;Pwn=sum(abs(wn). A2)/N;A1 =sqrt(2*Pwn* 10A(SNR/l 0);A2=AI;snl =A 1 *sin(2*pi*fl *t);sn2=A2*sin(2*pi*f2*t);sn=snl+sn2;Psn l=sum(abs(sn 1 ).A2)/N;Psn2=sum(abs(sn2).A2)/N;SNRreal I = 10*log 10(Psnl/Pwn);SNRreal2= 10*log 10(
3、Psn2/Pwn);xn=sn+wn;cf(I,:)=xn;eb(l,:)=xn;rou( 1 )=sum(abs(xn).A2)/N;%p=ladd 1=0;add2=0;for n=2:N,add 1 =add 1 +ef( l,n)*eb(l,n-1);add2=add2+abs(ef( 1 ,n)A2+abs(eb( I .n- 1)A2; end a( 1J )=-2*add l/add2;rou(2)=( 1 -abs(a( 1,1 )A2)*rou( 1);sigma( 1 )=rou(2);for n=3:N, ef(2,n)=ef( 1 ,n)+a( 1.1 )*eb( 1
4、,n-1); cb(2,n)=cb(l,n- l)+a(l,l)*ef(l,n);end %p=2% add 1=0;% add2=0;% for n=3:N,%add 1 =addl +ef(2,n)*eb(2,n-1);%add2=add2+abs(cf(2,n)A2+abs(cb(2,n-1 )A2;% end%a(2,2)=-2*addl/add2;% rou(3)=( 1 -abs(a(2,2)A2)*rou(2);% sigma(2)=rou(3);%a(2,l)=a(l,l)+a(2,2)*a(l,l);% for n=4:N,%ef(3,n)=ef(2,n)+a(2,2)*eb
5、(2,n-1);%eb(3,n)=eb(2,n-1 )+a(2,2)*ef(2,n);% end%迭代IP=20;for p=2:IP,add 1=0;add2=0;for n=p+l:N,add 1 =add 1 +ef(p,n)*eb(p,n-1); add2=add2+abs(ef(p,n)A2+abs(eb(p,n- 1)A2;enda(p.p)=-2*add l/add2;for n=p+2:N,ef(p+1 ,n)=ef(p,n)+a(p,p)*eb(p,n-1); eb(p+l,n)=eb(p,n-1 )+a(p,p)*ef(p,n); endrou(p+1 )=( 1 -abs
6、(a(p,p)A2)*rou(p);sigma(p)=rou(p+l);fbr k=l:p-l,a(p,k)=a(p-1 ,k)+a(p,p)*a(p- l,p-k);endendnfft=IO24;ifmod(nfft,2)=0Plength 二 nfft/2+1;elseif mod(nfft,2)= 1Plength=(nfift+l)/2;endF=linspace(0,fs/2,Plength);%频率坐标W=F*2*pi/fs;%数字角频率for omiga=l:Plengthadd3(omiga)=0;for 1=1 :p add3(omiga)=add3(oiniga)+a(p
7、,l)*cxp(-j*W(omiga)*l);endPxx(omiga)=sigma(p)/abs( 1 +add3(omiga)A2;endsubplot(211);plot(F,l 0*logl O(Pxx);grid on;xlabelC 频率(Hz);ylabel(助率谱幅度(dB);tille(编程实现Burg法功率谱估计，)；PxxB,FBJ=pburg(xn,p,nfft,fs);subplot(212);plot(FB,10*logl0(PxxB);grid on;xlabelC 频率(Hz);ylabel(功率谱幅度(dB);title(，命令实现Burg法功率谱估计)；勰实
8、现Burg法功率谱估计4020010004000500020003000癖(Hz)命令实现Burg法功率谱估计-200i丁 /-r、一I n)1/ ) 厂./ / / JI -Xrrrr-20-401000400020003000呻Hz)mpfill热&mpfi岫1藤棕5000600clear;clc;%Burg Data leng N=64clear;clc;N=1024;n=O:N-l;randn ( state*, 0);wn=5*randn(l, N);xn=sqrt(20)*sin(2*pi*0. 1*n)+sqrt(20)*sin(2*pi*0. 4*n)+wn; ef (1, :
9、)=xn; eb(l, :)=xn;tempsum=0;for i=l:l:Ntempsum=tempsum+xn(i). *xn(i);endthegma(1)=tempsum/N;%p=lsum1=0;sum2=0;for n=2:N,suml=suml+ef (1, n)*eb(l, nl);sum2=sum2+abs(ef(1, n)*2+abs(eb(l, n-1) 2;enda(l, l)=-2*suml/sum2;thegma= (l-abs(a(l, 1)2)*thegma(l);for n=3:N,ef (2, n)=ef (1, n)+a(l, l)*eb(l, n-1);
10、eb(2, n) =eb(l, n-l)+a(l, l)*ef (1, n);endP=l；IP=14;while pIPP=P+1；suml=0;sum2=0;for i=p+l:Nsuml=suml+ ef(p, i)*cb(p, il);sum2=sun)2+ef (p, i)*ef (p, i)+eb(p, i-l)*eb(p, i-1);a(p, p)=-2*suml/sum2;endthegma(p+l) = l-(abs(a(p, p) 2*thegma(p);Power Spectrum (dB)forTI-p+NNef (P+L iTef (pi)+a(pprebpill)
11、eb(p+L iTebpi)+a(ppTefpi) “enda(?k)na(PIL k)+a(pP)*1(PIL plk)- end end bnpa(p:.L zeros?NIP)- pxxuabs (Lhegma)2)PXXUPXX(1:N、2)一 pxxo*oglo(pxx)- 彳 ON nA(lz、2)“ 次 subplotp卜一)- por (n/M pxx)； xlabel( Frequency，)- ylabelC Power spectrum (dBy )-Lir+oc Burg DasQng NH64P) grid onmp) ErlMOQds MOCLmp) ErlMOQd
12、s MOCLwn=20*randn(l, N);例3. 5.1对一个含有白噪声和两个不同频率正弦信号的随机信号x()的功率谱进行分析。 解：设两个正弦信号玉()、易()的频率分别为。15和0. 2HZ,幅度分别为A, B；白噪声信号 noise(n)的方差为C,采样频率 = 1 HZ,则随机信号工()可描述为： fs=l;N=1000;t=0:l/fs:N;xl=A*sin(2*0. 15*pi*t); x2=B*sin(2*0. 2*pi*t);noise=C*randn(size(t); x=xl+x2+noisc;当参数C、A、B取不同值时，随机信号的波形如图3. 5. 1所示。51 p
13、 11 1 /I I f I Ji A |1.1 1A=B=1 , C=0, 演声信号)50100150200250300350400IPh hillI、h |L (L 1,(.IImH R fJ l| I Cl ill nfl ply Unruhs1 *.A=B= C= 1)50100150200250300350400L A i i A |.A=B=0.5, C=1)50100150200250300350400VAwn/仰r版.心L册卅M祈/W A=B=0.1, C=1)50100150200250300350400图3. 5.1 不同信噪比信号图3. 5. 2分别给出了基于burg递推
14、方法，用不同阶数的模型对不同信噪比信号的功率谱进 行估计的结果。图3. 5.2功率谱估计结果用两种不同的方法(YULEAR EQ和BURG和周期图法等)获得的结果比较clc;clear;A=l;B=l;C=l;fs二 1;N=100;%N=1000;%观测数据的点数不同，则分析结果有差异%从而说明了各种方法的不同t=O:l/fs:N;xl二A*sin(2*0. 15*pi*L);x2二B*sin(2*0. 2*pi*t);noise=C*randn (size (t);x=xl+x2+noise;P=5;P=15nfft=1024PxxB, FB=pburg(x, P, nfft, fs);figure(l);plot(FB, 10*logl0(PxxB);grid on;xlabel C 频率(Hz);ylabelC功率谱幅度(dB);titleC命令实现Burg法功率谱估计);hold on
