wu.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)
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)
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
a=-inv(cxx)*cxx0;

for i=1:length(f)
   sum=0;
   for k=1:p
       sum=sum+a(k)*exp(-m*2*pi*f(i)*k);
   end
  Pxzxfc(i)=delta/(abs(1+sum))^2;
end
figure
semilogy(f,Pxzxfc);
title('修正协方差算法');


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
xn= xn(:);
x=xn';
for k=0:1:sample-1
    s=0;
    for n=1:sample-k
       s=s+conj(x(1,n))*x(1,n+k); %calculate the value of rxx   
   end
   rxx(1,k+1)=(1/sample)*s;
end
 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);
            end;
            sum=0;
         for k=p+1:sample
         sum=sum+abs(ei*V(:,k))^2;
        end 
     Pmusicf(1,i)=(1/sum);
     end
   figure
semilogy(f,Pmusicf);
title('music算法');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
pp=6;
x=xn';
M=length(x);
rxx=xcorr(x,'biased');
rxx=[rxx(M:end),0]
R=toeplitz(rxx,rxx);
Rxx=R(1:M,1:M);
Rxy=R(1:M,2:end);
[U,S,V]=svd(Rxx);
U1=U(:,1:pp);
S1=S(1:pp,1:pp);
V1=V';
V1=V1(1:pp,:);
A=U1.'*Rxy*V1.';
d=eig(S1,A);
y=angle(d)/2/pi;
figure
stem(y,ones(1,length(y)));
title('esprit算法');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
xt=xn'; 
sin_num=3;
N=length(xt);
rxx=xcorr(xt,'biased');
rxx=rxx(N:(2*sin_num)+N);
%Frequencies estimation
Rxx=toeplitz(rxx);
ev=eig(Rxx);
[S i]=min(ev);
[V D]=eig(Rxx);
a=V(:,i);
rts=roots(a);
w_est=[];
for i=1:2*sin_num
   w_est(i)= angle(rts(i));
end
F=(w_est/(2*pi))';
figure
  stem(F,ones(1,length(F)),'*');
  title('Pisarenko算法');