wudisong.m

谱估计 (1)修正协方差法 (2)多重信号分类（MUSIC）算法 (3)ESPRIT算法 (4)皮萨论科（Pisarenko）谐波分解法 并对各算法进行分析。 解压后即可在MATLAB上运

%AR model parameter estimate 
clear
m=sqrt(-1);
delta=0.1;
a1=0.5;
sample=32;      %number of sample spot
p=10;           %number of sample spot in coef method;
f1=0.05; f2=0.20; f3=0.45;  %三个频率分量
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;

%calculate zn   噪声
for n=1:sample-1
    zn(1)=un(1);
    zn(n+1)=a1*zn(n)+un(n+1);
end
    
%calculate 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
    


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%calculate cxx(j,k)   公式3.102
for j=1:p
     for k=1:p
         s=0;
         for n=p+1:sample
            s=s+1/2/(sample-p)*(conj(xn(n-j))*xn(n-k)+xn(n-p+j)*conj(xn(n-p+k)));
         end
       cxx(j,k)=s;
   end
end
%calculate cxx0(j,1)  公式 3.103
for j=1:p
  s=0;
  for n=p+1:sample
      s=s+1/2/(sample-p)*(conj(xn(n-j))*xn(n)+xn(n-p+j)*conj(xn(n-p)));
  end
cxx0(j,1)=s; 
end
%calculate a    公式3.104
a=-inv(cxx)*cxx0;
%a(k)*exp(-j2pifk)累加
for i=1:length(f)   
   sum=0;
   for k=1:p
       sum=sum+a(k)*exp(-m*2*pi*f(i)*k);
   end
%a(k)cxx(1,k)累加 
   sun=0;     
  for k=1:p
      sun=sun+a(k)*cxx(1,k);
  end
   delta1=cxx(1,1)+sun;
 %AR谱计的形式
  XIE(i)=delta1/(abs(1+sum))^2;
end
figure
semilogy(f,XIE);
title('修正协方差算法');


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
xn= xn(:);
x=xn';
%calculate the value of rxx 
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
%calculate the value of Rx 
 Rx=zeros(sample,sample);
 Rx=toeplitz(rxx(1,1:32));
 [U,S,V]=svd(Rx);  %特征值分解
 Pmusicf=zeros(1,1/fstep+1);
 ei=zeros(1,sample); 
     for i=1:length(f)
            for j=1:sample
            ei(1,j)=exp(-2*pi*(j-1)*f(i)*m);  %calculate the value of ei
            end;
            sum=0;
         for k=p+1:sample
         sum=sum+abs(ei*V(:,k))^2;  %累加e*vi的绝对值
        end 
     Pmusicf(1,i)=(1/sum);
     end
   figure
semilogy(f,Pmusicf);
title('music算法');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
pp=6; %6 ge fen liang 
x=xn';
M=length(x);
rxx=xcorr(x,'biased'); % x de zi xiang guan
rxx=[rxx(M:end),0];     %qu hou 32 wei 
R=toeplitz(rxx,rxx);   %qiu Rxx,Rxy (公式5.158 5.159)
Rxx=R(1:M,1:M);
Rxy=R(1:M,2:end);
D=eig(Rxx);            %Rxx de te zheng zhi
[delta2 i]=min(D);     %qiu zhi xiao te zheng zhi
for i=pp+1:M           
    ey=eye(7,7);
    Z=ey(1:pp,2:end);  %qiu Z
    Cxx=Rxx(1:pp,1:pp)-delta2*ey(1:pp,1:pp); %p182 (4)
    Cxy=Rxy(1:pp,1:pp)-delta2*Z;
    D1=eig(Cxx,Cxy);         %p182 (5)          
    y=angle(D1)/2/pi;
end
figure
stem(y,ones(1,length(y)));
title('esprit算法');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
xt=xn'; 
sin_num=3;  %frequency fen liang
N=length(xt);
rxx=xcorr(xt,'biased');%x de zi xiang guang
rxx=rxx(N:(2*sin_num)+N);%rxx de qu zhi 
%Frequencies estimation
Rxx=toeplitz(rxx);  %qiu Rxx
ev=eig(Rxx);        %qiu Rxx de te zheng zhi 
[S i]=min(ev);      %qiu zhui xiao te zheng zhi 
[V D]=eig(Rxx);      %te zheng zhi fen jie
a=V(:,i);            %qiu zui xiao te zheng zhi dui ying de yi lie te zheng xiang liang 
rts=roots(a);       %qiu gen
w_est=[];
for i=1:2*sin_num
   w_est(i)= angle(rts(i));  %qiu jiao du 
end
F=(w_est/(2*pi))';       %qiu pin lv
figure
stem(F,ones(1,length(F)),'*');
title('Pisarenko算法');