%DEMDOA  HOSA Toolbox Demo of Direction of Arrival Estimation (DOA) 
%        Array Processing: Estimation of source bearings
%
echo off

% demos of DOA  (doagen)

% A. Swami April 15, 1993
% Copyright (c) 1991-2001 by United Signals & Systems, Inc. 
%       $Revision: 1.7 $

%     RESTRICTED RIGHTS LEGEND
% Use, duplication, or disclosure by the Government is subject to
% restrictions as set forth in subparagraph (c) (1) (ii) of the 
% Rights in Technical Data and Computer Software clause of DFARS
% 252.227-7013. 
% Manufacturer: United Signals & Systems, Inc., P.O. Box 2374, 
% Culver City, California 90231. 
%
%  This material may be reproduced by or for the U.S. Government pursuant 
%  to the copyright license under the clause at DFARS 252.227-7013. 

clear, clc, 
echo on

%	     Direction of Arrival Estimation (DOA) 
%        Array Processing: Estimation of source bearings

%  The Direction of Arrival (DOA) problem is frequently encountered
%  in sonar.  Here the signal consists of several narrow-band plane
%  waves impinging on an  (usually uniform linear) array of sensors.  
%  The sensor measurements may be corrupted by spatially correlated
%  noise.  If the spatial noise is strongly directional, second-order
%  statistics will classify it as a source.  If the spatial noise
%  is Gaussian, and the source signals are non-Gaussian, we can use
%  cumulants to estimate the bearings of the source signals. 
%  We can then use second-order statistics to characterize the spatially 
%  correlated noise. 
 
% Hit any key to continue
pause
clc

% The HOSA Toolbox offers three routines: 
% 1. DOAGEN generates synthetics for the DOA problem.  The user is 
%    prompted for all relevant parameters (number of sensors, sensor
%    spacing, source bearings, source signal and noise characteristics, etc.) 
%    DOA assumes an uniformly spaced linear array. 
%    
% 2. DOA processes the complex sensor data to estimate the direction of 
%    arrival, using the Eigenvector, MUSIC, Pisarenko, ML (Capon), AR,
%    Minimum-norm, Beamformer and ESPRIT methods, based either on 
%         the spatial fourth-order cross cumulants,  
%      or the spatial cross-covariance matrix. 
%    In the former case, DOA is insensitive to spatially and temporally 
%    correlated Gaussian noise sources. 
%
% Hit any key to continue
pause
clc

load doa1 
% Matrix ymat, in the matfile doa1.mat was generated by DOAGEN. 
% The data correspond to eight (8) sensors, spaced half a wavelength apart, 
% with two (2) sources at bearings of -15 and -25 degrees. 
% Note that the angular source separation is 10 degrees, and that the 
% beam-width of the conventional beamformer is pi/8 = 22.5 degrees. 
% The source signals are Laplace distributed, and the data are contaminated 
% with spatially correlated Gaussian noise.  The spatial correlation is such 
% that the angular spectrum of the noise has  strong peaks at +/-30 degrees. 
% The signal length is 4096 samples (snapshots). 
% The sensor and source configuration is depicted on the graphics screen. 
% The source signals and the noise have unit variance. 

echo off 
demdoat         % displays text and figure describing DOA problem

echo on 
% Hit any key to continue
pause 
clc

% We will now use DOA based on the spatial cross-correlation matrix,
%    to estimate the source bearings. 
% The sensor spacing (half wavelength) is assumed know;  we will compute the
% angular spectrum at a spacing of 1 degree;  and we will use an array 
% displacement of 1 element for ESPRIT.   Finally, we will assume that there 
% are three sources.  
% The data consists of two source signals; however, the sensor noise is
% spatially correlated, and effectively acts as a third source signal. 
% In the noiseless case, the rank of the spatial correlation matrix
% should be the number of sources;  however, spatially correlated noise
% sources will also appear as virtual sources. 
% We use the SVD to determine the effective rank:  If the number of sources 
% is p, we expect the singular values to show a sharp drop at p+1. 
%
% Hit any key to continue
pause 


clf
[sp2,th2,b2] = doa(ymat,0.5,1,3,2,1); 
set (gcf, 'Name', 'HOSA correlation-based DOA')

% The singular value plot indicates the possible presence of three sources.
% Recall that the observed signal consists of two sources at -15 and -25
% degrees, and that the signal is contaminated with spatially correlated
% noise which acts as a virtual source at a bearing of 30 degrees.
% As expected, the ML, AR, MUSIC and Min-Norm methods resolve the two sources
% and the virtual source due to the noise, whereas the beamformer does not. 
% The MUSIC and min-norm estimates are virtually identical in this example. 

% Hit any key to continue
pause

% Let us zoom in on the display
echo off 
for k=2:8
  eval(['subplot(42' int2str(k) ')' ]) 
  a = axis; axis([-45 45 a(3), a(4)])
end 
echo on 
% 

% Hit any key to continue
pause

% We can also use fourth-order cross-cumulants in the DOA problem. 
% DOA allows you to estimate the number of (non-Gaussian) source signals.
% A matrix of fourth-order cross cumulants of the sensor signals is computed.
% It is known that the rank of this matrix equals the number of non-Gaussian
% sources:  we use the Singular Value Decomposition (SVD) of the matrix
% to estimate the rank (If the matrix size is M, and the number of sources is 
% p, we expect the singular values to show a sharp drop at p+1). 
% 
% You will notice from the singular value plot that the singular values 
% are essentially zero for n > 2, indicating two non-Gaussian sources;
% justifying our assumption of 2 sources: 
% 

% Hit any key to start the DOA computations 
pause

clf
[sp4,th4,b4] = doa(ymat,0.5,1,2,4,1);
set(gcf,'Name','HOSA cumulant-based DOA')

% The angular spectra have sharp peaks at -15 and -25 degrees. 
% The angular spectrum due to the noise is barely visible 

% Hit any key to continue
pause

% Let us zoom in on the display
echo off 
for k=2:8
  eval(['subplot(42' int2str(k) ')' ]) 
  a = axis; axis([-45 45 a(3), a(4)])
end 
echo on 
% 

% Hit any key to continue
pause

% Let us look at both estimates together 
echo off
clf
set(gcf,'Name','DOA c2 and c4')
subplot(211)
semilogy(th2, sp2), grid on 
title('DOA-c2 "spectra" ')
xlabel('angle in degrees') 
hold on 
 a=axis;  
   plot([-15 -15], [a(3) a(4)], 'g-')
   plot([-25 -25], [a(3) a(4)], 'g-')
   plot([ 30  30], [a(3) a(4)], 'r-')
   plot(-15,a(4),'go'), plot(-25,a(4),'go'), plot(30,a(4),'gx')
 hold off


subplot(212) 
semilogy(th4, sp4), grid on 
title('DOA-c4 "spectra" ')
xlabel('angle in degrees') 
hold on 
 a=axis;  
   plot([-15 -15], [a(3) a(4)], 'g-')
   plot([-25 -25], [a(3) a(4)], 'g-')
   plot([ 30  30], [a(3) a(4)], 'r-')
   plot(-15,a(4),'go'), plot(-25,a(4),'go'), plot(30,a(4),'gx')
 hold off

echo on
%The true source bearings are shown by circles, and the noise source by
%a cross; the vertical lines are added for clarity.

%Hit any key to continue
pause

% Let us zoom in on the display
subplot(211), a=axis; axis([-45 0 a(3) a(4)])
subplot(212), a=axis; axis([-45 0 a(3) a(4)])


% 
% Routine TDE solves the two-sensor time-delay estimation problem.
% Choose the `Time Delay Estimation' option from the main menu to learn
%    more about HOSA's TDE routine.  
% Hit any key to return to main menu 
pause 
echo off
clc


return