sfastdoa1.m

快速doa估计方法

function SFastDoa1()
% 快速Doa方法 ----线阵   单次  
clear;
                                        % t = clock;
%信源
%--------------------------------------------
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');

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

%全局变量定义
%-------------------------------------------- 
ii=sqrt(-1);
K=3;
M=10; %阵元数 
th=[-30;0;30];%角度间隔小于10度，不能分辨
nn=100;%采样数
sn=[10;10;10];
degrad=pi/180;
%--------------------------------------------

%生成离散信号
%--------------------------------------------

tt=1:nn;
S=[s1(tt).';s2(tt).';s3(tt).'];
nr=randn(M,nn);
ni=randn(M,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];
tmp1 = [0:M-K-1]';
x = length(th2);
tmp = [0:M-1]'; %阵元位置
A = exp(ii*pi*tmp*sin(th'*degrad)); %方向矩阵(M*N)
X=A*S2+Z;%阵列接收信号  (M*nn)     
                                               
Rld = X(K+1:M,:)*X(1:K,:)'/nn;%阵列数据的协方差矩阵(M*M)
U = orth(Rld);
W = eye(M-K)-U*U';
                                           
for k = 1:x
    a = exp(ii*tmp1*pi*sin(th2(k)*degrad));
    doa(k) = 1/(a'*W*a);
end
                                                 
doa = abs(doa);

res = SearchP1(doa,K);
result = th2(res')

%做图
% --------------------------------------------
% figure;
%subplot(1,2,2)
% plot(th2,abs(doa));
semilogy(th2,doa);
title('FastDOA 1-D');
xlabel('Azimuth(deg)');
ylabel('Spectrum');
% axis([-40 40 0.1 1e7]);
grid on;
%hold on;
                                             %  etime(clock,t) 
%--------------------------------------------