gusucode.com > ​均匀线阵和均匀圆阵的DOA估计matlab源码程序 > code/ULCAMUS.m

    %-------JIN Wei
%----ULA and UCA estimation  based on the MUSIC method
%Remark 1: The UCA estimation has a worse performance than ULA
%estimation under the MUSIC ,espatially in the low snr,though it has the same resolution in
%360 degree. The reason is that the vector of the UCA convariance matrix is
%not orthogonal to each other.
%Remark 2: The program is based on the assumption of narrowband signal
%model.
%Remark 3:If we increase the element number or snr, we will find the the
%merit of UCA array.
clc;
clear;
clear all;

msens = 16; 
nsamp = 256; 

fc = 10000;%the carrier frequency
f  = [123,183,263,333];%the signal frequency
ts = 1/2/fc;%the sampling time inteval
v=3*10^8;%the speed of electromagnetic wave
lambda=v/fc;%
msource=4;
bearing  =  [20 40 65 85]*pi/180; 


%--------------ULA,uniform linear array DOA estimation-------------
dspace = 0.5; 
for i=1:length(bearing)
    smat_ULA(i,:)=hilbert(cos(2*pi*f(i)*(0:nsamp-1)*ts).*cos(2*pi*fc*(0:nsamp-1)*ts));%signal matrix of the sources
end
%-----------------generate data--------------------------------
amat_ULA = exp(sqrt(-1) * [0:msens-1]' * (2*pi*dspace*sin(bearing)));%array manifold matrix
smat_ULA = amat_ULA*smat_ULA;%signal matrix of the array elements
theta=0:pi/100:pi/2;
for run=1:30
%-----------signal + noise-------------------------------------
snr=20;
for i=1:msens,
    	noise=	hilbert(randn(1,nsamp));
    	confi=(norm(smat_ULA(i,:))^2)/(10^(snr/10))/(norm(noise)^2);
    	smat_ULA(i,:)=smat_ULA(i,:)+noise.*sqrt(confi);
end
R=smat_ULA*smat_ULA'/nsamp;
[V,eigvalue]=eig(R);
[Y,i]=sort(diag(eigvalue));
amat1 = exp(sqrt(-1) * [0:msens-1]' * (2*pi*dspace*sin(theta)));
P=zeros(size(theta'));
for k=1:msens-msource,
    P=P+abs(amat1'*V(:,i(k))).^2;
end	
P=1./P;
% P=P/max(P);
P=10*log10(P);
figure(1)
plot(theta*180/pi,P,'b')
hold on
end





%--------------UCA,uniform circular array DOA estimation-------------
R = lambda/4/sin(pi/msens); 
beta=2*pi*R/lambda;
for i=1:length(bearing)
    smat_UCA(i,:)=hilbert(cos(2*pi*f(i)*(0:nsamp-1)*ts).*cos(2*pi*fc*(0:nsamp-1)*ts));%signal matrix of the sources
end
%-----------------generate data--------------------------------
for m=1:length(bearing)
    amat_UCA(:,m) = exp(sqrt(-1)*(beta*cos(bearing(m)-[0:msens-1]'*2*pi/msens)));%array manifold matrix
end
a=size(amat_UCA)
b=size(smat_UCA)
smat_UCA = amat_UCA*smat_UCA;%signal matrix of the array elements
theta=0:pi/100:pi/2;
for run=1:30
%-----------signal + noise-------------------------------------
snr=20;
for i=1:msens,
    	noise=	hilbert(randn(1,nsamp));
    	confi=(norm(smat_UCA(i,:))^2)/(10^(snr/10))/(norm(noise)^2);
    	smat_UCA(i,:)=smat_UCA(i,:)+noise.*sqrt(confi);
end
R=smat_UCA*smat_UCA'/nsamp;
[V,eigvalue]=eig(R);
[Y,i]=sort(diag(eigvalue));
for m=1:length(theta)
    amat1(:,m) = exp(sqrt(-1)*(beta*cos(theta(m)-[0:msens-1]'*2*pi/msens)));
end
P=zeros(size(theta'));
for k=1:msens-msource,
    P=P+abs(amat1'*V(:,i(k))).^2;
end	
P=1./P;
% P=P/max(P);
P=10*log10(P);
figure(2)
plot(theta*180/pi,P,'r')
hold on
end