Unexpected Angle of Arrival Detection Result for Sonar Application

Hi guys,

I am learning to use the phased array toolbox for sonar application. When trying the example "Locating an Acoustic Beacon with a Passive Sonar System", I noticed that when I changed the position, velocity, and/or orientation of the ULA array I might get wrong angle estimation results. Specifically, when running the below codes (I copied exact the same example from MATLAB webiste and changed the parameters of the beacon and array platforms, and the array orientation), I have detection at -45 degree (the true is 45 degree). Examine the angle power spectrum (at the last pulse), I see that there are 2 peaks both at 45 and -45 degree. I just wonder why it happens like that. It also would be great if you guys can point me to some references in the literature regarding this issue.

propSpeed = 1520;
channelDepth = 200;
OperatingFrequency = 37.5e3; 

isopaths = phased.IsoSpeedUnderwaterPaths('ChannelDepth',channelDepth,...
  'NumPathsSource','Property','NumPaths',10,'PropagationSpeed',propSpeed);

channel = phased.MultipathChannel('OperatingFrequency',OperatingFrequency);

prf = 1;
pulseWidth = 10e-3;
pulseBandwidth = 1/pulseWidth;
fs = 2*pulseBandwidth;
wav = phased.RectangularWaveform('PRF',prf,'PulseWidth',pulseWidth,...
  'SampleRate',fs);
channel.SampleRate = fs;

projector = phased.IsotropicProjector('VoltageResponse',120);

projRadiator = phased.Radiator('Sensor',projector,...
  'PropagationSpeed',propSpeed,'OperatingFrequency',OperatingFrequency);

beaconPlat = phased.Platform('InitialPosition',[1000; 1000; -2],...
  'Velocity',[0; 0.1; 0]);

hydrophone = phased.IsotropicHydrophone('VoltageSensitivity',-150);
array = phased.ULA('Element',hydrophone,...
  'NumElements',5,'ElementSpacing',propSpeed/OperatingFrequency/2,...
  'ArrayAxis','x');

arrayCollector = phased.Collector('Sensor',array,...
  'PropagationSpeed',propSpeed,'OperatingFrequency',OperatingFrequency);

arrayPlat = phased.Platform('InitialPosition',[0; 0; -100],...
  'Velocity',[0; 0.1; 0]);

rx = phased.ReceiverPreamp(...
    'Gain',20,...
    'NoiseFigure',10,...
    'SampleRate',fs,...
    'SeedSource','Property',...
    'Seed',2007);

x = wav();
numTransmits = 10;
rxsig = zeros(size(x,1),5,numTransmits);
for i = 1:numTransmits

  % Update array and acoustic beacon positions
  [pos_tx,vel_tx] = beaconPlat(1/prf);
  [pos_rx,vel_rx] = arrayPlat(1/prf);

  % Compute paths between the acoustic beacon and array
  [paths,dop,aloss,rcvang,srcang] = ...
      isopaths(pos_tx,pos_rx,vel_tx,vel_rx,1/prf);

  % Propagate the acoustic beacon waveform
  tsig = projRadiator(x,srcang);
  rsig = channel(tsig,paths,dop,aloss);
  
  % Collect the propagated signal
  rsig = arrayCollector(rsig,rcvang);
  
  % Store the received pulses
  rxsig(:,:,i) = abs(rx(rsig));
 
end

t = (0:length(x)-1)'/fs;
figure
plot(t,rxsig(:,end))
xlabel('Time (s)');
ylabel('Signal Amplitude (V)')

musicspatialspect = phased.MUSICEstimator('SensorArray',array,...
        'PropagationSpeed',propSpeed,'OperatingFrequency',...
        OperatingFrequency,'ScanAngles',-90:0.1:90,'DOAOutputPort',true,...
        'NumSignalsSource','Property','NumSignals',1);
    
numTransmits = 500;
angPassive = zeros(numTransmits,1);
angAct = zeros(numTransmits,1);

for i = 1:numTransmits

  % Update array and acoustic beacon positions
  [pos_tx,vel_tx] = beaconPlat(1/prf);
  [pos_rx,vel_rx] = arrayPlat(1/prf);

  % Compute paths between acoustic beacon and the array
  [paths,dop,aloss,rcvang,srcang] = ...
      isopaths(pos_tx,pos_rx,vel_tx,vel_rx,1/prf);
  angAct(i) = rcvang(1,1);
  
  % Propagate the acoustic beacon waveform
  tsig = projRadiator(x,srcang);
  rsig = channel(tsig,paths,dop,aloss);
  
  % Collect the propagated signal
  rsig = arrayCollector(rsig,rcvang);
  
  rxsig = rx(rsig);
 
  % Estimate the direction of arrival
  [~,angPassive(i)] = musicspatialspect(rxsig);
  
end

figure
plot([angPassive angAct])
xlabel('Pulse Number')
ylabel('Arrival angle (degrees)')
legend('Estimated DOA','Actual DOA')

figure
plotSpectrum(musicspatialspect)

Answer 1

[Aashray]

Hello Dat!

I tried executing the code in the latest version of MATLAB (R2024b) and was unable to reproduce the error (I got a single peak, corresponding to the positive azimuthal angle).

But fortunately, I also have MATLAB R2022b installed on my older laptop, and I tried to reproduce the problem and encountered the same issue you described (getting two peaks).

I believe MathWorks has addressed this issue already, and a simple version update should solve the issue.