aepmethod.m

空时自适应信号处理

clc;clear all;close all;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
N= 8;                      % number of sensors
M= 8;                      % number of pulse
lamda=0.03;               % wavelength
V_p=90;                    % platform velocity
d=0.015;                   % spacing of sensors
PRF=12000;                 % pulse repetition frequency
T=1/PRF;                   % pulse repetition interval
H=500;                     % platform height    
R=1000;                    % range
theta=asin(H/R);           % depression angle
Bc=0;                      % clutter bandwidth
mm=0:1:N-1;mm=mm';
nn=0:1:M-1;nn=nn';
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%V_d=20;
%fai_darget=pi/2;
%Wd_d=4*pi*V_d*cos(fai_darget)/(lamda*PRF);         %目标多普勒角频率
%Ws_d=2*pi*d*cos(fai_darget)/lamda;                 %目标空间角频率
%bM_d=exp(j*Wd_d*mm);                                %M×1目标维时间导向矢量
%aN_d=exp(j*Ws_d*nn);                                %N×1目标维空间导向矢量
%S_d=kron(bM_d,aN_d); 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% %Clutter Covariance Matrix% %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for m=1:1:M
    for p=1:1:M
        for i=1:1:N
            for k=1:1:N
                l=(m-1)*N+i;
                n=(p-1)*N+k;
                phi=linspace(0,2*pi,61);
                D=0.5*(1+cos(2*(phi-pi/2)));
                D=D.^2;                
                                
                %D=1;% 有无加权影响很大！
                G=1;                    
                t_phase=exp(j*2*pi/lamda*2*V_p*(m-p)*T*cos(phi));
                s_phase=exp(j*2*pi/lamda*(i-k)      *d*cos(phi));
                integral=D.*t_phase.*s_phase.*G;
                
                Q(l,n)=sum(integral)/60;
               % Q(l,n)=Q(l,n)*exp(-Bc*Bc*(m-p)^2/8);
            end
        end
    end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
CNR=trace(Q)/trace(0.00001*eye(N*M,N*M))
Q=Q+0.001*eye(N*M,N*M);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fd_d=0;
AA=exp(j*2*pi*nn);
BB=exp(j*2*pi*mm*fd_d);
S_d=kron(AA,BB);         %search channel
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%特征分解AEP
Q_tr=trace(Q);           %
Q_inv=inv(Q);            %
Q_d=eig(Q);
%d=flipud(d);
%[V,D] = eig(Q,'nobalance')
[V,D]=eig(Q);            %对杂波协方差矩阵进行特征分解
D_inv=inv(D);
T_c=V(:,1:50);          %auxiliary channel
%T_c=T_c';
TT=T_c;
%TT=[S_d,T_c];           %transform matrix
%QT=TT'*Q*TT;
K=50;  
for i=1:K-1
    QT(i,i)=0.001;
end
QT(K,K)=Q_d(K);
QT_inv=inv(QT);

%% IF %%
ii=1;
for fd=-0.5:0.01:0.5
    
    AAi=exp(j*2*pi*nn);
    BBi=exp(j*2*pi*mm*fd);
    S=kron(BBi,AAi);
    ST=TT'*S;
    %Wopt=Q_inv*S;
    %WT=QT_inv*ST;
    US1=zeros(M*N,1);
    US2=zeros(K,1);
  for i1=1:M*N
     US1(i1,1)=V(:,i1)'*S/Q_d(i1);
  end
    % Wopt=V*(D_inv*(V'*S));
    Wopt=V*US1; 
  for i=1:K
     US2(i,1)=V(:,i)'*S/Q_d(i);
  end
     WT=US2;
     %WT=QT_inv*V(:,1:50)'*S;
    SINR1(ii)=abs(Wopt'*S);
    SINR2(ii)=abs(WT'*ST);
    ii=ii+1;
end
SINR1=10*log10(SINR1/max(SINR1));
SINR2=10*log10(SINR2/max(SINR2));    
    
figure(1);
mesh(abs(Q));title('Modulus');

figure(2);
fd=-0.5:0.01:0.5;
plot(fd,SINR1,'-b');
hold on;
plot(fd,SINR2,'r');