root_mvdrmse.m

信号与信息处理——阵列信号处理DOA估计的matlab算法

clc;
clear all;
q=1;
save RootMSE.mat;
%=========== Parameters setting =================
M=50;      % sensor number
S=3;      % signal source number
J=1;
snap=200; % snapshot
%===== signal model parameters setting =====
fs=1000;      % sampling frquence (Hz)
f(1)=100;      % the desired signal frequence (Hz)
f(2)=200;      % the coherent interference frequence (Hz)
f(3)=300;      % the uncorrelated interference frequence (Hz)

Th(2)=-40;     % DOA of the coherent interference (degree)
Th(3)=40;    % DOA of the uncorrelated interference (degree)
iiii=1;
for SNR=-20:4:20
    SNR1=SNR;
    SNR2=SNR;
    SNR3=SNR;
    Np=1;         % Var. (power) of noise 
    Sp(1)=Np*10^(SNR1/10);     % Var. (power) of the desired signal
    Sp(2)=Np*10^(SNR2/10);     % Var. (power) of the coherent interference 
    Sp(3)=Np*10^(SNR3/10);     % Var. (power) of the uncorrelated interference

              %===== transform signal power to amplitude ===== 
    for ii=1:1:length(Sp)
        Ap(ii)=sqrt(Sp(ii));
    end
    An=sqrt(Np);  % Amplitude of noise  
       %==== interelement spacing ====
    for ii=1:1:M
        if ii==1
            d(ii)=0;
        else
            d(ii)=0.5;
        end	
    end
    Th(1)=0;      % DOA of the desired signal (degree)
    %===== construct true sterring vector matrix =====
    for ii=1:1:M
        for iii=1:1:S
            A(ii,iii)=exp(j*(ii-1)*2*pi*d(ii)*sin(Th(iii)*pi/180));      % the estimated
        end
    end
    eTh(1)=0*pi/180;      % estimate DOA of the desired signal (degree)
%========== Weighting ===============
    Rxx=zeros(M,M);
    for t=1:1:snap
%====== construct signal vector ======
        for ii=1:1:S
            s(ii,1)=Ap(ii)*exp(-j*2*pi*(f(ii)/fs)*(t-1));
        end
    %======= construct spatial white noise vector =======
        N=An*randn(M,1);
        X=A*s+N;
        Rxx=Rxx+X*X';
    end
    eRxx=Rxx/snap;     % estimation of correlation matrix
    AA=pinv(eRxx);
    for i=-(M-1):(M-1)
        z(i+M)=sum(diag(AA,i))
    end
    aa=roots(z);
    aa2=abs(aa);
    [cc,dd]=sort(abs(1-abs(aa)));
    for iii=1:1:S
        ccc=sum(sqrt(cc(1*iii,1)^2));
    end
    MVDR(iiii)=ccc;
    iiii=iiii+1;
end
x=real(aa);
y=imag(aa);
zz=atand(y./x);
save RootMSE.mat MVDR -append
figure(1)
snr=-20:4:20
semilogy(snr,MVDR)
xlabel('SNR')
ylabel('RMSE'),

Ed=-90:1:90;
figure(2)
plot(x ,y,'b.'),

ylabel('image'),
xlabel('real')
axis([-2 2 -2 2])
title('Root-MVDR')