fourth_order_music1.m

阵列信号处理_空间谱估计原理的一些基本算法——MusicEspritMp

  %  the  data of signal
clear;
format long;
c=3*10.^8;
L=6;
sam=128*2;
N=128*2;
f1=1000;
f2=1000;
snr=0;
Amp=sqrt(2*10^(snr/10));
lamta=c/f1;
len=lamta/2;
%***************************** doa *************************************
P1=40*pi/180;
P2=50*pi/180;
%***************************** doa *************************************

%P;
pesai1=2*pi*rand(1);
pesai2=2*pi*rand(1);
for t=1:sam,
    s1(t)=Amp*(exp(j*2*pi*f1*0.1*t+pesai1));
    s2(t)=Amp*(exp(j*2*pi*f2*0.1*t+pesai2));
end
  s=[s1(1:N);s2(1:N)];
%t=1:sam ; 
%plot(t,s1)

  %n=1:1800;
  %aa=(0.1*n)*pi/180;
  
  i=1:L;
  x1=exp(j*2*pi*f1*(len*(i-1)*sin(P1))/c);
  x2=exp(j*2*pi*f2*(len*(i-1)*sin(P2))/c);

  a1=x1.';
  a2=x2.';
  
  a=[a1,a2];
  
  %a=[a1.';a2.'];
   
 noise=randn(L,N)+j*randn(L,N);
 
 z=a*s+noise;
x=s1;   
M=1;  
for k1=1:M;
     for k2=1:M;
        for k3=1:M;
            for k4=1:M;
                c1=0;
                c2=0;
                c3=0;
                c4=0;
                c5=0;
                c6=0;
                c7=0;
                for l=1:N;
                    c1=c1+x(k1,l)*conj(x(k2,l))*x(k3,l)*conj(x(k4,l));
                    c2=c2+x(k1,l)*conj(x(k2,l));
                    c3=c3+x(k3,l)*conj(x(k4,l)) ;
                    c4=c4+x(k1,l)*x(k3,l);
                    c5=c5+conj(x(k2,l))*conj(x(k4,l));
                    c6=c6+x(k1,l)*conj(x(k4,l));
                    c7=c7+conj(x(k2,l))*x(k3,l);
                end
                c=c1/N;
                c2=c2/N;
                c3=c3/N;
                c4=c4/N;
                c5=c5/N;
                c6=c6/N;
                c7=c7/N;
               R1((k1-1)*M+k2,(k3-1)*M+k4)=c1-c2*c3-c4*c5-c6*c7; 
            end
        end
    end
end

x=s2;   
M=1;  
for k1=1:M;
     for k2=1:M;
        for k3=1:M;
            for k4=1:M;
                c1=0;
                c2=0;
                c3=0;
                c4=0;
                c5=0;
                c6=0;
                c7=0;
                for l=1:N;
                    c1=c1+x(k1,l)*conj(x(k2,l))*x(k3,l)*conj(x(k4,l));
                    c2=c2+x(k1,l)*conj(x(k2,l));
                    c3=c3+x(k3,l)*conj(x(k4,l)) ;
                    c4=c4+x(k1,l)*x(k3,l);
                    c5=c5+conj(x(k2,l))*conj(x(k4,l));
                    c6=c6+x(k1,l)*conj(x(k4,l));
                    c7=c7+conj(x(k2,l))*x(k3,l);
                end
                c=c1/N;
                c2=c2/N;
                c3=c3/N;
                c4=c4/N;
                c5=c5/N;
                c6=c6/N;
                c7=c7/N;
               R2((k1-1)*M+k2,(k3-1)*M+k4)=c1-c2*c3-c4*c5-c6*c7; 
            end
        end
    end
end

T=zeros(2);
T(1,1)=R1;
T(2,2)=R2;
b1=kron(a1,conj(a1));
b2=kron(a2,conj(a2));
b=[b1,b2];
C=b*T*b';
[e,v]=eig(C);

 es=e(:,1:2);
 en=e(:,(3:L.^2));
itr=9000;
aaaa=zeros(L,itr);
 for k=1:L
     for h=1:itr
         aaaa(k,h)=exp(j*pi*(k-1)*sin(0.01*h*pi/180));
         %aaaa(k,h)=exp(j*2*pi*w1*len*(k-1)*sin(0.1*h*pi/180)/c);
     end
 end
 b=zeros(L.^2,itr);
for h=1:itr
     b(:,h)=kron(aaaa(:,h),aaaa(:,h)'.');
     aac=b(:,h);
     pp=norm(en'*aac).^(-2);
     %pp=1/((aac)'*c*aac);
     p(h)=real(pp);
     ppp(h)=pp;
 end
 p1=p(1:itr/2);
 p2=p(itr/2+1:itr);
[m1,n1]=max(p1);
doa1=n1/100
[m2,n2]=max(p2);
doa2=n2/100+45
u=0:0.01:90-0.01;
ppp=10*log(ppp);
plot(u,ppp);