www.gusucode.com > DOA定位算法源码程序 > DOA定位算法源码程序/DOA.m

    % Simulation of MUSIC, ESPRIT, MVDR, Min-Norm and Classical DOA
% algorithms for a uniform linear array.
clear all;
clc ;

doas=[-30 -5 40 70]*pi/180; %DOA抯 of signals in rad.
P=[1 1 1 1]; %Power of incoming signals
N=10; %Number of array elements
K=1024; %Number of data snapshots
d=0.5; %Distance between elements in wavelengths
noise_var=1; %Variance of noise
r=length(doas); %Total number of signals
% Steering vector matrix. Columns will contain the steering vectors
% of the r signals
A=exp(-i*2*pi*d*(0:N-1)'*sin([doas(:).']));
% Signal and noise generation
sig=round(rand(r,K))*2-1; % Generate random BPSK symbols for each of the


% r signals
noise=sqrt(noise_var/2)*(randn(N,K)+i*randn(N,K)); %Uncorrelated noise
X=A*diag(sqrt(P))*sig+noise; %Generate data matrix
R=X*X'/K; %Spatial covariance matrix
[Q ,D]=eig(R); %Compute eigendecomposition of covariance matrix
[D,I]=sort(diag(D),1,'descend'); %Find r largest eigenvalues
Q=Q (:,I); %Sort the eigenvectors to put signal eigenvectors first
Qs=Q (:,1:r); %Get the signal eigenvectors
Qn=Q(:,r+1:N); %Get the noise eigenvectors
% MUSIC algorithm
% Define angles at which MUSIC 搒pectrum?will be computed
angles=(-90:0.1:90);
%Compute steering vectors corresponding values in angles
a1=exp(-i*2*pi*d*(0:N-1)'*sin([angles(:).']*pi/180));
for k=1:length(angles)
%Compute MUSIC 搒pectrum?
music_spectrum(k)=(a1(:,k)'*a1(:,k))/(a1(:,k)'*Qn*Qn'*a1(:,k));
end
figure(1)
plot(angles,abs(music_spectrum))
title('MUSIC Spectrum')
xlabel('Angle in degrees')
%ESPRIT Algorithm
phi= linsolve(Qs(1:N-1,:),Qs(2:N,:));
ESPRIT_doas=asin(-angle(eig(phi))/(2*pi*d))*180/pi;
%MVDR
IR=inv(R); %Inverse of covariance matrix
for k=1:length(angles)
mvdr(k)=1/(a1(:,k)'*IR*a1(:,k));
end
figure(gcf+1)
plot(angles,abs(mvdr))
xlabel('Angle in degrees')
title('MVDR')
%Min norm method
alpha=Qs(1,:);
Shat=Qs(2:N,:);
ghat=-Shat*alpha'/(1-alpha*alpha');
g=[1;ghat];
for k=1:length(angles)
minnorm_spectrum(k)=1/(abs(a1(:,k)'*g));
end
figure(gcf+1)
plot(angles,abs(minnorm_spectrum))
xlabel('Angle in degrees')
title('Min-Norm')
%Estimate DOA抯 using the classical beamformer
for k=1:length(angles)
Classical(k)=(a1(:,k)'*R*a1(:,k));
end
figure(gcf+1)
plot(angles,abs(Classical))
xlabel('Angle in degrees')
title('Classical Beamformer')