预测误差格型滤波器及伯格精

1、结束图3.3.5 Burg递推法流程图close all;clear;clc;fs=10000;f1=fs*0.2;f2=fs*0.25;N=100;t=0:1/fs:(N-1)/fs;randn(state,0);wn=randn(1,N);SNR=10;Pwn=sum(abs(wn)，2)/N;A1=sqrt(2*Pwn*10A(SNR/10);A2=A1;sn1=A1*sin(2*pi*f1*t);sn2=A2*sin(2*pi*f2*t);sn=sn1+sn2;Psn1=sum(abs(sn1).A2)/N;Psn2=sum(abs(sn2)，2)/N; SNRreal1=10*log

2、10(Psn1/Pwn); SNRreal2=10*log10(Psn2/Pwn); xn=sn+wn;%采样频率设置（Hz）%信号的数据长度%产生高斯白噪声信号%信噪比%噪声信号的功率%由信噪比计算正弦信号的振幅%正弦序列%有用信号%正弦信号的功率%计算实际信躁比%输入信号=有用信号+高斯白噪声ef(1,:)=xn;eb(1,:)=xn;rou(1)=sum(abs(xn).A2)/N;%p=1add1=0;add2=0;for n=2:N,add1=add1+ef(1,n)*eb(1,n-1);add2=add2+abs(ef(1,n)A2+abs(eb(1,n-1)A2; enda(1,

3、1)=-2*add1/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,n-1);eb(2,n)=eb(1,n-1)+a(1,1)*ef(1,n);end % %p=2% add1=0;% add2=0;% for n=3:N,%add1=add1+ef(2,n)*eb(2,n-1);%add2=add2+abs(ef(2,n)A2+abs(eb(2,n-1)A2;% end% a(2,2)=-2*add1/add2;% rou(3)=(1-abs(a(2,2)A2

4、)*rou(2);% sigma(2)=rou(3);% a(2,1)=a(1,1)+a(2,2)*a(1,1);% for n=4:N,%ef(3,n)=ef(2,n)+a(2,2)*eb(2,n-1);%eb(3,n)=eb(2,n-1)+a(2,2)*ef(2,n);% end%迭代IP=20;for p=2:IP,add1=0;add2=0;for n=p+1:N,add1=add1+ef(p,n)*eb(p,n-1);add2=add2+abs(ef(p,n)A2+abs(eb(p,n-1)A2;enda(p,p)=-2*add1/add2;for n=p+2:N,ef(p+1,n)

5、=ef(p,n)+a(p,p)*eb(p,n-1);eb(p+1,n)=eb(p,n-1)+a(p,p)*ef(p,n);endrou(p+1)=(1-abs(a(p,p)A2)*rou(p);sigma(p)=rou(p+1);for k=1:p-1,a(p,k)=a(p-1,k)+a(p,p)*a(p-1,p-k);endendnfft=1024;if mod(nfft,2)=0Plength=nfft/2+1;elseif mod(nfft,2)=1Plength=(nfft+1)/2;endF=linspace(0,fs/2,Plength);%频率坐标W=F*2*pi/fs;%数字角

6、频率for omiga=1:Plengthadd3(omiga)=0;for l=1:padd3(omiga)=add3(omiga)+a(p,l)*exp(-j*W(omiga)*l);endPxx(omiga)=sigma(p)/abs(1+add3(omiga)A2;endsubplot(211);plot(F,10*log10(Pxx);grid on;xlabel(频率(Hz);ylabel(功率谱幅度(dB);title(编程实现Burg法功率谱估计);PxxB,FB=pburg(xn,p,nfft,fs);subplot(212);plot(FB,10*log10(PxxB);g

7、rid on;xlabel(频率(Hz);ylabel(功率谱幅度(dB);title(命令实现Burg法功率谱估计);编程实现Burg法功率谱估计40200-20500001000200030004000频率(Hz)命令实现Burg法功率谱估计0-20-40-60500001000200030004000频率(Hz)clear;clc;%Burg Data leng N=64 clear;clc;N=1024;n=0:N-1;randn(state,0);wn=5*randn(1,N);xn=sqrt(20)*sin(2*pi*0.1*n)+sqrt(20)*sin(2*pi*0.4*n)+

8、wn;ef(1,:)=xn;eb(1,:)=xn;tempsum=0;for i=1:1:Ntempsum=tempsum+xn(i).*xn(i);endthegma(1)=tempsum/N;%p=1sum1=0;sum2=0;for n=2:N,sum1=sum1+ef(1,n)*eb(1,n-1);sum2=sum2+abs(ef(1,n)A2+abs(eb(1,n-1)A2;enda(1,1)=-2*sum1/sum2;thegma(2)=(1-abs(a(1,1)A2)*thegma(1);for n=3:N,ef(2,n)=ef(1,n)+a(1,1)*eb(1,n-1);eb(

9、2,n)=eb(1,n-1)+a(1,1)*ef(1,n);endp=1;IP=14;while pIPp=p+1;sum1=0;sum2=0;for i=p+1:Nsum1=sum1+ ef(p,i)*eb(p,i-1);sum2=sum2+ef(p,i)*ef(p,i)+eb(p,i-1)*eb(p,i-1);a(p,p)=-2*sum1/sum2;endthegma(p+1)=1-(abs(a(p,p)A2*thegma(p);for i=p+2:Nef(p+1,i)=ef(p,i)+a(p,p)*eb(p,i-1);eb(p+1,i)=eb(p,i-1)+a(p,p)*ef(p,i);

10、endfor k=1:1:p-1a(p,k)=a(p-1,k)+a(p,p)*a(p-1,p-k);endendb=1,a(p,:),zeros(1,N-p-1);Pxx=abs(thegma(p)./(abs(fft(b).A2);Pxx=Pxx(1:N/2);Pxx=10*log10(Pxx);n=0:N-1;n=n(1:N/2);% subplot(2,2,1);plot(n/N,Pxx);xlabel(Frequency);ylabel(Power Spectrum (dB);title(Burg Data leng N=64);grid on;mHrrceDsre WBBurg Da

11、ta leng N=64302826242220181614121000.10.20.3Frequency0.40.5Burg Data leng N=6423 00.10.20.30.40.582252Frequencywn=20*randn(1,N);例3.5.1对一个含有白噪声和两个不同频率正弦信号的随机信号x(n)的功率谱进行分析。解：设两个正弦信号 x1(n)、x2(n)的频率分别为0.15和0.2HZ,幅度分别为A, B;白噪声信号 noise(n)的方差为C,采样频率fs=1HZ,则随机信号x(n)可描述为： fs=1;N=1000;t=0:1/fs:N;x1=A*sin(2*0

12、.15*pi*t);x2=B*sin(2*0.2*pi*t);noise=C*randn(size(t);x=x1+x2+noise;当参数C A B取不同值时，随机彳t号的波形如图3.5.1所示。05010015020025030035040005010015020025030035040050-5050100150200250300350400图3.5.1 不同信噪比信号图3.5.2分别给出了基于burg递推方法，用不同阶数的模型对不同信噪比信号的功率谱进 行估计的结果。A=B=1 , C=0图3.5.2功率谱估计结果用两种不同的方法（YULEAR EQ和BURG和周期图法等）获得的结果比

13、较clc;clear;A=1;B=1;C=1;fs=1;N=100;%N=1000;概测数据的点数不同，则分析结果有差异 烟而说明了各种方法的不同t=0:1/fs:N;x1=A*sin(2*0.15*pi*t);x2=B*sin(2*0.2*pi*t);noise=C*randn(size(t);x=x1+x2+noise;P=5;P=15 nfft=1024PxxB,FB=pburg(x,P,nfft,fs);figure(1);plot(FB,10*log10(PxxB);grid on;xlabel(频率(Hz);ylabel(功率谱幅度(dB);title(命令实现Burg法功率谱估计);hold onpxx1,f1=pyulear(x,P,nfft,fs); plot(f1,10*log(pxx1), r*);grid on;title( YULEAR EQ)pp=abs(fft(x);PC=pp.*pp/N;figure(2);PCC=10*log(PC);plot(PC