couple_estimate.m

存在互耦的平面阵的一种互耦参数估计方法。

%  the  data of signal
clear;
clc;
format long;
c=3*10.^8;%光速
L=10;%阵元数
sam=500;%取样点
N=500;
dx=1/2; %阵元间距
dy=1/2;
w1=0.8;%第一个信号的角频率
w2=0.9;%第二个信号的角频率
w3=1;
w4=1.1;
phase1=0;%初始相位
phase2=0;
phase3=0;
phase4=0;
singnum=4;
snr=0;
Amp=sqrt(10^(snr/10));%信号幅度
cx=0.3527+0.4854j;
cy=0.3527+0.4854j; %互相关系数
cxy=0.0927-0.2853j;
theta=[74,35;40,20;54,66;28,41];
%lamta=c/w1;
%len=lamta/2;


for t=1:sam
    s1(t)=Amp*(exp(-j*(2*pi*w1*0.1*t+phase1)));
    s2(t)=Amp*(exp(-j*(2*pi*w2*0.1*t+phase2)));
    s3(t)=Amp*(exp(-j*(2*pi*w3*0.1*t+phase2)));
    s4(t)=Amp*(exp(-j*(2*pi*w4*0.1*t+phase2)));
end
tt=1:sam;
%plot(tt,s1,'r--',tt,s2,'b--');

  s=[s1(1:N);s2(1:N);s3(1:N);s4(1:N)];

  %n=1:1800;
  %aa=(0.1*n)*pi/180;
  
  i=1:L;
  ax1=exp(-j*2*pi*(dx*(i-1)*cos(theta(1,1)*pi/180)*sin(theta(1,2)*pi/180)))';%阵列上关于信号的方向向量
  ax2=exp(-j*2*pi*(dx*(i-1)*cos(theta(2,1)*pi/180)*sin(theta(2,2)*pi/180)))';
  ax3=exp(-j*2*pi*(dx*(i-1)*cos(theta(3,1)*pi/180)*sin(theta(3,2)*pi/180)))';
  ax4=exp(-j*2*pi*(dx*(i-1)*cos(theta(4,1)*pi/180)*sin(theta(4,2)*pi/180)))';
  ay1=exp(-j*2*pi*(dy*(i-1)*sin(theta(1,1)*pi/180)*sin(theta(1,2)*pi/180)))';
  ay2=exp(-j*2*pi*(dy*(i-1)*sin(theta(2,1)*pi/180)*sin(theta(2,2)*pi/180)))';
  ay3=exp(-j*2*pi*(dy*(i-1)*sin(theta(3,1)*pi/180)*sin(theta(3,2)*pi/180)))';
  ay4=exp(-j*2*pi*(dy*(i-1)*sin(theta(4,1)*pi/180)*sin(theta(4,2)*pi/180)))';
  a=[kron(ay1,ax1),kron(ay2,ax2),kron(ay3,ax3),kron(ay4,ax4)];
  
  c1=toeplitz([1,cx,zeros(1,L-2)]);
  c2=toeplitz([cx,cxy,zeros(1,L-2)]);
  C1=kron(eye(L),c1);
  C2=kron((tril(ones(L),1)-tril(ones(L),-2)-eye(L)),c2);
  C=C1+C2;
  %a=[a1.';a2.'];
  J=[zeros(L-2,1),eye(L-2),zeros(L-2,1)];
  P=kron(J,J);
%   (L-2)*(L-2) L*L
 noise=randn( L*L,N)+j*randn(L*L,N);%零均值、方差为1的白噪声，且与信号源不相关
%P*  +noise ;%总的阵列输出向量
 z=C*a*s+noise;
 Rz=(z*z')/N;%取z的自相关矩阵
 
 [e,v]=eig(Rz);%对Rz进行特征值分解

 es=e(:,L*L-3:L*L);%信号子空间
 en=e(:,1:L*L-4);%噪声子空间
Tx11=zeros(L,2);
Tx12=zeros(L,2);
Tx13=zeros(L,2);
Tx14=zeros(L,2);
Tx21=zeros(L,2);
Tx22=zeros(L,2);
Tx23=zeros(L,2);
Tx24=zeros(L,2);
for i=1:L
    for j=1:2
        if(i+j<=L+1)
            Tx11(i,j)=ax1(i+j-1);
            Tx12(i,j)=ax2(i+j-1);
            Tx13(i,j)=ax3(i+j-1);
            Tx14(i,j)=ax4(i+j-1);
        end
        if((i>=j)&(j>=2))
          Tx21(i,j)=ax1(i-j+1);
          Tx22(i,j)=ax2(i-j+1);
          Tx23(i,j)=ax3(i-j+1);
          Tx24(i,j)=ax4(i-j+1);  
        end
    end
end
Ty11=zeros(L,2);
Ty12=zeros(L,2);
Ty13=zeros(L,2);
Ty14=zeros(L,2);
Ty21=zeros(L,2);
Ty22=zeros(L,2);
Ty23=zeros(L,2);
Ty24=zeros(L,2);
for i=1:L
    for j=1:2
        if(i+j<=L+1)
            Ty11(i,j)=ay1(i+j-1);
            Ty12(i,j)=ay2(i+j-1);
            Ty13(i,j)=ay3(i+j-1);
            Ty14(i,j)=ay4(i+j-1);
        end
        if((i>=j)&(j>=2))
          Ty21(i,j)=ay1(i-j+1);
          Ty22(i,j)=ay2(i-j+1);
          Ty23(i,j)=ay3(i-j+1);
          Ty24(i,j)=ay4(i-j+1);  
    end
    end
end
TX1=Tx11+Tx21;
TX2=Tx12+Tx22;
TX3=Tx13+Tx23;
TX4=Tx14+Tx24;
TY1=Ty11+Ty21;
TY2=Ty12+Ty22;
TY3=Ty13+Ty23;
TY4=Ty14+Ty24;
T1=kron(TY1,TX1);
T2=kron(TY2,TX2);
T3=kron(TY3,TX3);
T4=kron(TY4,TX4);
Q=[en'*T1;en'*T2;en'*T3;en'*T4];
esti_c=-[Q(:,2)+Q(:,3),Q(:,4)]\Q(:,1)