Estimating the direction of arrival of a sound signal using time delay sum beamforming

I have a recorded/simulated audio file which is basically an impulse signal when I plot. I want to use a linear microphone array of 8 elements separated by some distance 'd' to record the audio(here it is simulation only). So now considering the audio signal being recorded at all microphones, I want to estimate the direction of arrival of the signal(assuming it is far-field scenario) using time-delay-sum beamformer(Phased array system toolbox). As I am new to this, I am quite confused with how to give the input file to the microphones and how to input the angles and estimate them using the beamforming technique. Here is the sample code Iam attaching which I got it from matlab exapmles help section.

array = phased.ULA('NumElements',11,'ElementSpacing',0.04);
sensorpos = getElementPosition(array);
fs = 8000;
t = 0:1/fs:0.3;
x = chirp(t,0,1,500);
c = 340;
collector = phased.WidebandCollector('Sensor',array,...
    'PropagationSpeed',c,'SampleRate',fs,'ModulatedInput',false);
incidentAngle = [90;10];
x = step(collector,x.',incidentAngle);
sigma = 0.2;
noise = sigma*randn(size(x));
rx = x + noise;
beamformer = phased.TimeDelayBeamformer('SensorArray',array,...
    'SampleRate',fs,'PropagationSpeed',c,...
    'Direction',incidentAngle);
y = step(beamformer,rx);

Here, in 'x' place I have some other audio but what I couldn't understand is there is nothing related to estimate the angle. I just want to give some input audio file to an array of sensors with some azimuth and elevation angles, and I want to estimate that DOA using time delay sum beamforming. Could anyone please help me to get through this. Thanks in advance!!

Answer 1

[Chunru]

See the comments inside the code

array = phased.ULA('NumElements',11,'ElementSpacing',0.04);
sensorpos = getElementPosition(array);
fs = 8000;
% the chirp has higher freq 1000-3000 for better doa estimation
t = 0:1/fs:0.3;
x = chirp(t,1000,0.3,3000);  

c = 340;
collector = phased.WidebandCollector('Sensor',array,...
    'PropagationSpeed',c,'SampleRate',fs,'ModulatedInput',false);
incidentAngle = [30; 0];  % elevation =0;
x = step(collector,x.',incidentAngle);
sigma = 0.2;
noise = sigma*randn(size(x));
rx = x + noise;

% The following beamformer is looking at the direction incidentAngle
beamformer = phased.TimeDelayBeamformer('SensorArray',array,...
    'SampleRate',fs,'PropagationSpeed',c,...
    'DirectionSource','Input port');   % lood direction is specified later
% The following will produce the waveform from the beamformer wiht the
% specified look direction incidentAngle
%y = step(beamformer,rx);

% For doa estimation with broadband signal, you don't know where signal
% comes from.  You have to SEARCH over different directions.  Then you can
% determine the signal power (for example) at each looking direction.  Then
% based on the signal power at different direction, you can determine DOA.

% Let's say we are looking at different az directions (el=0)
az = [-180:3:180];
el = 0;
p = zeros(length(az), 1);       % power at different look direction
for i=1:length(az)
    y = step(beamformer,rx, [az(i); el]);       % beamformer at different direction
    p(i) = std(y);                              % power
end

plot(az, p);grid;
hold on
xline(incidentAngle(1))

It can be seen that the DOA is 30 deg (with left/right ambiguity for ULA).

[Rakesh Reddy]

Thanks a lot. This cleared my doubts to a lot of extent...

[Rakesh Reddy]

How to estimate the elevation angle in this scenario if I give some elevation angle instead of zero ??

[Chunru]

Genrally speaking, you need a 2d search over alll possible elevation and azimuth angles.

For ULA, it has ambiguity over a cone-shaped pattern.

[Rakesh Reddy]

The audio plot is having two impulse signals at times t1 and t2. Now If I perform the delay and sum beamforming on that audio, whether the estimated angle belongs to the signal1 or signal2 ? How to estimate the angle for the preferred one in the audio using this method? If we can't find using this, then is there any method to estimate that?

[Chunru]

You need to differentiate the time and doa of the signal. We are dealing with spatial-temporal signal. If a sound source from a fixe direction (spatial) and have some time events (pulse 1 and 2) (hope I understand you correctly), then you can do the delay and sum beamforming to form the time series at different direction (the beamformer output y as in the code above). If signal varying with time, you cannot just do std (power) over beamformer output. You need to do instataneous or moving power estimte. Try "movstd" with a suitable time window. Then plot the output as a function of DOA and time (in sonar, this is called bearing-time recording or BTR display). From there, you can see the signal at different time and differnet location.