psd.m

现代谱估计实验

clear all;
x=load('BBB.TXT');
N=length(x);
K=zeros(1,3);
e=zeros(3,N);
b=zeros(3,N);
aa=zeros(1,3);
a=zeros(3,3);
sigma=zeros(1,3);

%去除均值
ee=0;
for i=1:N
 e=x(i)+ee;
end
en=ee/N;
for i=1:N
x(i)=x(i)-en;
end
e0=zeros(1,N);
b0=zeros(1,N);

%初始化
% e0,b0
for i=1:N
e0(i)=x(i);
b0(i)=x(i);
end
% sigma0
sigma=0;
for i=1:N
sigma=sigma+x(i)^2;
end
sigma0=sigma/N;

%求参数
pp=13;%阶数设置为3
for p=1:pp
   m=0;
   n=0;
   if p==1  %e0,b0的情况
     for j=p+1:N
       m=e0(j)*b0(j-1)+m;
       n=e0(j)^2+b0(j-1)^2+n;
      end
    else  %其它情况
      for j=p+1:N
        m=e(p-1,j)*b(p-1,j-1)+m;
        n=e(p-1,j)^2+b(p-1,j-1)^2+n;
      end
     end
   K(p)=-2*m/n;
   if p==1
      e(p,1)=e0(1);
      b(p,1)=K(p)*e0(1);
      for i=2:N
        e(p,i)=e0(i)+K(p)*b0(i-1);
        b(p,i)=b0(i-1)+K(p)*e0(i);
      end
    else
      e(p,1)=e(p-1,1);
      b(p,1)=K(p)*e(p-1,1);
      for i=2:N
        e(p,i)=e(p-1,i)+K(p)*b(p-1,i-1);
        b(p,i)=b(p-1,i-1)+K(p)*e(p-1,i);
      end
    end
    aa(p)=K(p);
    a(p,p)=aa(p);
    if p>1 %当p>1时才运行此步
     for i=1:p-1
       a(p,i)=a(p-1,i)+aa(p)*a(p-1,p-i);
     end
    end
    % 求sigma
    if p==1
     sigma(p)=(1-abs(aa(p))^2)*sigma0;
    else
     sigma(p)=(1-abs(aa(p))^2)*sigma(p-1);
    end
end
clear j;
%求功率
wn=2*pi/N;
wx=0:wn:(2*pi-wn);
i=1;
S=zeros(1,length(wx));
for w=0:wn:(2*pi-wn)
  sm=0;
  for k=1:pp
   sm=a(pp,k)*(cos(k*w)-sin(k*w)*j)+sm;
  end
 S(i)=sigma(pp)/(abs((1+sm))^2); 
 i=i+1;
end
%归一化处理
smax=max(S);
SS=S/smax;
%xpsd=10*log10(SS);

%画图
figure(1);
subplot(211);
plot(x);
title('原始波形图');
xlabel('时间 t');
ylabel('幅值 A');
subplot(212);
plot(wx,SS);
title('信号的功率谱图');
xlabel('角频率 w');
ylabel('功率谱 S');