load('E:\grad_project\multiTarget\2Dmusic\simulink\walk.mat')
twpi=2*pi;
derad=pi/180;
radeg=180/pi;
kelmx=2;   %横向阵元个数
kelmy=2;   %纵向阵元个数
snr=10;
% 多径参数 3路信号 方位角覆盖范围90 正负 15度，俯仰角覆盖范围 45 - 90度

iwave=2;
n=500;
dd=0.5;   %距离波长比
dx=0:dd:(kelmx-1)*dd;
dy=0:dd:(kelmy-1)*dd;

for time = 1:20
    [theta0,phi0] = walk(0);
    for i = 1:length(theta0)
        Ax0=exp(-1j*twpi*dx.'*(sin(theta0(i)*derad).*cos(phi0(i)*derad))); %横向导向矢量
        Ay0=exp(-1j*twpi*dy.'*(sin(theta0(i)*derad).*sin(phi0(i)*derad))); %纵向导向矢量
        B(:,i)=kron(Ax0,Ay0);%直积
    end
    %------------------------非相干信源----------------------------
    att = [0.4 0.8]; %头、胸、胸、胳膊、胳膊、腿、腿
    Gamma = [];
    for i = 1:iwave
        Gamma = [Gamma; att(i)*randn(1,n)]; % complex attuation
    end
    %  S=randn(iwave,n);
    X0=B*Gamma;
    X1=awgn(X0,snr,'measured'); %在信号中加入高斯白噪声
    Rxx0=X1*X1';
    issp = 1;  
    % spatial smoothing music
    if issp == 1
        Rxx = msspFor2D(Rxx0);
    elseif issp == 2
        % 仅针对2*M大小的天线阵，进行水平平滑
        Rxx = sspFor2D(Rxx0,kelmx - 1);
        dx=0:dd:(kelmx-2)*dd;
    else
        Rxx = Rxx0;
    end
    invRxx=inv(Rxx);  %没用到
    [EV,D]=eig(Rxx);
    EVA=diag((D)');
    [EVA,I]=sort(EVA);
    EVA=fliplr(EVA);
    EV=fliplr(EV(:,I));
    En=EV(:,iwave+1:end); 
    theta = 40:90;
    phi = 60:120;
    for ang1=1:length(theta)
        for ang2=1:length(phi)
            a1=exp(-1j*twpi*dx.'*sin(theta(ang1)*derad)*cos(phi(ang2)*derad));
            a2=exp(-1j*twpi*dy.'*sin(theta(ang1)*derad)*sin(phi(ang2)*derad));
            a=kron(a1,a2);%直积
            m=En*En';
            SP(ang1,ang2)=1/(a'*(En*En')*a);
        end 
    end
    SP=abs(SP);
    SPmax=max(max(SP));
    SP=SP/SPmax;
    result.image(:,:,time) = SP;
    result.theta = theta;
    result.phi = phi;
end
save walk.mat result
% save '..\data\bendKnee2.mat' result;


% h=mesh(theta,fy,SP);
% set(h,'Linewidth',2)
% xlabel('theta(degree)')
% ylabel('phi(degree)')
% zlabel('magnitude(dB)')