pbarpcf_huzi.m

本源码是用于现代谱估计的一种方法：巴特利特主分量估计算法

clear all
m=sqrt(-1);
delta=0.101043;
a1=-0.850848;
sample=32;      %采样点
p=10;           %number of sample spot in coef method;
%f1=0.05; f2=0.40; f3=0.42;
f1=0.05; f2=0.30; f3=0.42;
fstep=0.01;
fstart=-0.5;
fend=0.5;  
f=fstart:fstep:fend;
nfft=(fend-fstart)/fstep+1;

%构造噪声：高斯白噪声，正态分布，复数un=urn+juin
urn= normrnd(0,delta/2,1,sample);
uin= normrnd(0,delta/2,1,sample);
un=urn+m*uin;

%计算 zn
for n=1:sample-1
    zn(1)=un(1);
    zn(n+1)=-a1*zn(n)+un(n+1);
end
    
%求xn
    for n=1:sample
        xn(n)=2*cos(2*pi*f1*(n-1))+2*cos(2*pi*f2*(n-1))+2*cos(2*pi*f3*(n-1))+sqrt(2)*real(zn(n));
    end

  x=xn;
M=31;
for k=0:1:sample-1
   s=0;
   for n=1:sample-k,
      s=s+conj(x(1,n))*x(1,n+k);   
  end
  rxx(1,k+1)=(1/sample)*s;   %用有偏估计        
end
 Rxx=toeplitz(rxx(1,1:32));        %变矩阵
 [U,S,V]=svd(Rxx);  

 for i=1:length(f)
     for j=1:sample
            ei(1,j)=exp(-2*pi*(j-1)*f(i)*m);     %设定E矩阵
     end;
     sum=0;
         for k=1:6
         sum=sum+S(k,k)*abs(ei*V(:,k))^2;
         end 
     Pbarpcf(1,i)=sum/M;
 end
figure(4)
plot(f,Pbarpcf);
title('巴特利特主分量估计算法f2＝0.30时')