sfastdoa2.m

快速doa估计方法

function SFastDoa2()
% 快速doa方法  ---平面阵   单次
clear;
                                         
%信源
%--------------------------------------------
Pd=2000;
Fd=1;% Fd is the sampling rate of the message signal.
Fs=4*Fd;
R=0.5;
Delay=5;
No=1;
M=4;
x1=randint(Pd,1,M);
x2=randint(Pd,1,M);
x3=randint(Pd,1,M);
y1=modmap(x1,Fd,Fs,'qask',M);
y2=modmap(x2,Fd,Fs,'qask',M);
y3=modmap(x3,Fd,Fs,'qask',M);
[rcv_a1,ti]=rcosflt(y1,Fd,Fs,'fir/sqrt/Fs',R,Delay);
[rcv_a2,ti]=rcosflt(y2,Fd,Fs,'fir/sqrt/Fs',R,Delay);
[rcv_a3,ti]=rcosflt(y3,Fd,Fs,'fir/sqrt/Fs',R,Delay);
s1=amodce(rcv_a1,10,'qam');
s2=amodce(rcv_a2,10,'qam');
s3=amodce(rcv_a3,10,'qam');

% x4=randint(Pd,1,M);                  
% x5=randint(Pd,1,M);                 
% x6=randint(Pd,1,M);                 
% y4=modmap(x4,Fd,Fs,'qask',M);      
% y5=modmap(x5,Fd,Fs,'qask',M);        
% y6=modmap(x6,Fd,Fs,'qask',M);       
% [rcv_a4,ti]=rcosflt(y4,Fd,Fs,'fir/sqrt/Fs',R,Delay);              
% [rcv_a5,ti]=rcosflt(y5,Fd,Fs,'fir/sqrt/Fs',R,Delay);              
% [rcv_a6,ti]=rcosflt(y6,Fd,Fs,'fir/sqrt/Fs',R,Delay);            
% s4=amodce(rcv_a4,10,'qam');                          
% s5=amodce(rcv_a5,10,'qam');                        
% s6=amodce(rcv_a6,10,'qam');  
% 
% x7=randint(Pd,1,M);                  
% x8=randint(Pd,1,M);                              
% y7=modmap(x7,Fd,Fs,'qask',M);      
% y8=modmap(x8,Fd,Fs,'qask',M);              
% [rcv_a7,ti]=rcosflt(y7,Fd,Fs,'fir/sqrt/Fs',R,Delay);              
% [rcv_a8,ti]=rcosflt(y8,Fd,Fs,'fir/sqrt/Fs',R,Delay);                        
% s7=amodce(rcv_a7,10,'qam');                          
% s8=amodce(rcv_a8,10,'qam'); 
% --------------------------------------------

%全局变量定义
%-------------------------------------------- 
                                             
ii=sqrt(-1);        
dx=0.5;
dy=0.5;
K=3; % 信号源数目
M=K; %阵元数
N=10; 
th=[-50,0,50];
fi=[150,100,50];
% th=[-50,0,50,-70,-30,30,65,80];
% fi=[150,100,50,170,120,70,30,15];
nn=250;%采样数
sn=[10;10;10];
% sn=[10;10;10;10;10;10;10;10];
degrad=pi/180; 
%--------------------------------------------

%生成离散信号
%--------------------------------------------
tt=1:nn;
S=[s1(tt).';s2(tt).';s3(tt).'];
% S=[s1(tt).';s2(tt).';s3(tt).';s4(tt).';s5(tt).';s6(tt).';s7(tt).';s8(tt).'];
nr=randn(M*N,nn);
ni=randn(M*N,nn);
Z=nr+ii*ni; %每路噪声方差为2 ？.
Ps=S*S'/nn;%  (N*N)  
ps=diag(Ps);%平均功率(p*1)
refp=2*10.^(sn/10);%10*log10(S/N)=SNR
tmp=sqrt(refp./ps);
S2=diag(tmp)*S;%构造信号(N*nn)，信噪比为sn

%--------------------------------------------

%协方差矩阵
%-------------------------------------------
th2=[-90:90]';
fi2=[0:180]';
x1 = length(th2);
x2 = length(fi2);                           %t = clock;
                                            
for i=1:K
    A(:,i) = FastOriVec(M,N,th(i),fi(i),dx,dy);
end
X = A*S2+Z;
% d = floor(K/N)+1;                        
                                             
Rld = X(K+1:M*N,:)*X(1:K,:)'/nn;  % (MN-K)*K

                                   

U = orth(Rld);   

                                   
% e = zeros(M*N-K,1);      %(1)
% e(1) = 1;                %(1)
% V = e*e';                %(1)   
% V = eye(M*N-K);          %(2)
% V = zeros(M*N-K);           %(3) 
% n = K+2;                    %(3)
% e = ones(n,1);            %(3)
% for i = 1:M*N-K-n+1           %(3)
%     V(i:i+n-1,i) = e/n;   %(3)
% end                         %(3)

W = eye(M*N-K)-U*U';   % (MN-K)*(MN-K)
                                            % tic
for i = 1:x1
    for j = 1:x2
       r = FastOriVec(M,N,th2(i),fi2(j),dx,dy); 
       a = r(1:M*N-K);   % (M*N-K)*1
       doa(i,j) = a'*a/(a'*W*a);
%        doa(i,j) = 1/(a'*W*V*W*a);   % 加权V  (MN-K)*(MN-K)
    end
end
                                        
doa = abs(doa);
     
res = SearchP2(doa,K);
result = [th2(res(:,1)),fi2(res(:,2))];            
%                                                    toc 
%                                                    etime(clock,t)

%--------------------------------------------

% subplot(1,2,2)
a = th2*ones(1,length(fi2));
b = ones(length(th2),1)*fi2';
% doa=reshape(doa,length(th2),length(fi2));
figure
mesh(a,b,doa);
% title('FastDOA 2-D')
title('2-D DOA Nonuniform MUSIC ');
xlabel('Azimuth(deg)');
ylabel('Elevation(deg)');
zlabel('Spectrum');
grid on;
                                           
%--------------------------------------------