burg2.m

信号的谱估计的相关问题例子

clear;
clc;
%Burg Data leng N=64
clear;
clc;
N=1024;
n=[0:N-1];
randn('state',0);
wn=randn(1,N);
xn=sqrt(20)*sin(2*pi*0.1*n)+sqrt(20)*sin(2*pi*0.4*n)+wn;
ef(1,:)=xn;
eb(1,:)=xn;
tempsum=0;
for i=1:1:N
  tempsum=tempsum+xn(i).*xn(i);
end
thegma(1)=tempsum/N;

%p=1
sum1=0;
sum2=0;
for n=2:N,
    sum1=sum1+ef(1,n)*eb(1,n-1);
    sum2=sum2+abs(ef(1,n))^2+abs(eb(1,n-1))^2;
end
a(1,1)=-2*sum1/sum2;
thegma(2)=(1-abs(a(1,1))^2)*thegma(1);
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=1;
IP=14;
while p<IP
    p=p+1;
    sum1=0;
    sum2=0;
    for i=p+1:N
       sum1=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;
    end
     thegma(p+1)=[1-(abs(a(p,p)))^2]*thegma(p);
   for i=p+2:N
       ef(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);
   end
   for k=1:1:p-1
        a(p,k)=a(p-1,k)+a(p,p)*a(p-1,p-k);
    end
end
b=[1,a(p,:),zeros(1,N-p-1)];
Pxx=abs(thegma(p))./((abs(fft(b))).^2);
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;
%   
% %Burg Data leng N=128
% clear ;
% clc;
% N=128;
% n=[0:N-1];
% randn('state',0);
% wn=randn(1,N);
% xn=sqrt(20)*sin(2*pi*0.2*n)+sqrt(20)*sin(2*pi*0.3*n)+wn;
% ef(1,:)=xn;
% eb(1,:)=xn;
% tempsum=0;
% for i=1:1:N
%   tempsum=tempsum+xn(i).*xn(i);
% end
% thegma(1)=tempsum/N;
% p=1;
% IP=14;
% while p<IP
%     p=p+1;
%     sum1=0;
%     sum2=0;
%     for i=p:1:N-1
%        sum1=sum1+ ef(p-1,i)*eb(p-1,i-1);
%        sum2=sum2+ef(p-1,i)*ef(p-1,i)+eb(p-1,i-1)*eb(p-1,i-1);
%        a(p,p)=-2*sum1/sum2;
%     end
%     for k=1:1:p-1
%         a(p,k)=a(p-1,k)+a(p,p)*a(p-1,p-k);
%     end
%     thegma(p)=[1-(abs(a(p,p)))^2]*thegma(p-1);
%    for i=p+1:1:N-1
%        ef(p,i)=ef(p-1,i)+a(p,p)*eb(p-1,i-1);
%        eb(p,i)=eb(p-1,i-1)+a(p,p)*ef(p-1,i);
%    end
% end
% b=[1,a(IP,:),zeros(1,N-IP-1)];
% Pxx=abs(thegma(IP))./((abs(fft(b))).^2);
% Pxx=Pxx(1:N/2);
% Pxx=10*log10(Pxx);
% n=n(1:N/2);
% subplot(2,2,2);
% plot(n/N,Pxx);
% xlabel('Frequency');
% ylabel('Power Spectrum (dB)');
% title('Burg Data leng N=128');
% grid on;
% 
% %Burg Data leng N=512
% clear;
% clc;
% N=512;
% n=[0:N-1];
% randn('state',0);
% wn=randn(1,N);
% xn=sqrt(20)*sin(2*pi*0.2*n)+sqrt(20)*sin(2*pi*0.3*n)+wn;
% ef(1,:)=xn;
% eb(1,:)=xn;
% tempsum=0;
% for i=1:1:N
%   tempsum=tempsum+xn(i).*xn(i);
% end
% thegma(1)=tempsum/N;
% p=1;
% IP=14;
% while p<IP
%     p=p+1;
%     sum1=0;
%     sum2=0;
%     for i=p:1:N-1
%        sum1=sum1+ ef(p-1,i)*eb(p-1,i-1);
%        sum2=sum2+ef(p-1,i)*ef(p-1,i)+eb(p-1,i-1)*eb(p-1,i-1);
%        a(p,p)=-2*sum1/sum2;
%     end
%     for k=1:1:p-1
%         a(p,k)=a(p-1,k)+a(p,p)*a(p-1,p-k);
%     end
%     thegma(p)=[1-(abs(a(p,p)))^2]*thegma(p-1);
%    for i=p+1:1:N-1
%        ef(p,i)=ef(p-1,i)+a(p,p)*eb(p-1,i-1);
%        eb(p,i)=eb(p-1,i-1)+a(p,p)*ef(p-1,i);
%    end
% end
% b=[1,a(IP,:),zeros(1,N-IP-1)];
% Pxx=abs(thegma(IP))./((abs(fft(b))).^2);
% Pxx=Pxx(1:N/2);
% Pxx=10*log10(Pxx);
% n=n(1:N/2);
% subplot(2,2,3);
% plot(n/N,Pxx);
% xlabel('Frequency');
% ylabel('Power Spectrum (dB)');
% title('Burg Data leng N=512');
% grid on;
% 
% %Burg Data leng N=1024
% clear;
% clc;
% N=1024;
% n=[0:N-1];
% randn('state',0);
% wn=randn(1,N);
% xn=sqrt(20)*sin(2*pi*0.2*n)+sqrt(20)*sin(2*pi*0.3*n)+wn;
% ef(1,:)=xn;
% eb(1,:)=xn;
% tempsum=0;
% for i=1:1:N
%   tempsum=tempsum+xn(i).*xn(i);
% end
% thegma(1)=tempsum/N;
% p=1;
% IP=14;
% while p<IP
%     p=p+1;
%     sum1=0;
%     sum2=0;
%     for i=p:N-1
%        sum1=sum1+ ef(p-1,i)*eb(p-1,i-1);
%        sum2=sum2+ef(p-1,i)*ef(p-1,i)+eb(p-1,i-1)*eb(p-1,i-1);
%        a(p,p)=-2*sum1/sum2;
%     end
%     for k=1:1:p-1
%         a(p,k)=a(p-1,k)+a(p,p)*a(p-1,p-k);
%     end
%     thegma(p)=[1-(abs(a(p,p)))^2]*thegma(p-1);
%    for i=p+1:1:N-1
%        ef(p,i)=ef(p-1,i)+a(p,p)*eb(p-1,i-1);
%        eb(p,i)=eb(p-1,i-1)+a(p,p)*ef(p-1,i);
%    end
% end
% b=[1,a(IP,:),zeros(1,N-IP-1)];
% Pxx=abs(thegma(IP))./((abs(fft(b))).^2);
% Pxx=Pxx(1:N/2);
% Pxx=10*log10(Pxx);
% n=n(1:N/2);
% subplot(2,2,4);
% plot(n/N,Pxx);
% xlabel('Frequency');
% ylabel('Power Spectrum (dB)');
% title('Burg Data leng N=1024');
% grid on;