
% This plots beam patterns -- the gain as a function of angle produced 
% by a delay-and-sum beamformer.  (It's not actually the gain that is
% plotted, but rather the improvement in signal-to-noise ratio.)  You
% specify the hydrophone or microphone positions, the beam angles you
% want plotted, the frequency or frequencies to plot, and so on.  This
% script figures out the response of a beamformer for your configuration.
%
% If you specify more than one frequency, this averages the beampatterns
% over the set of frequencies you give.  If you specify more than one 
% angle, this plots several beampatterns at once.

%%%%%%%%%%%%%%%%%%%%%%%%%%%% Configuration %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

ang       = ishBeamAngles(:);   % angles to make beams on, radians
freqs     = ishBeamFreqs(:).';    % frequencies to average over, Hz
weighting = ishBeamWeighting;     % use 1 or 0 for weighting/no
zeroIsY   = ishBeamOrient;        % is 0 deg in polar plot x-axis (0) or y-axis (1)?
p         = ishPhonePositions.';  % [x y] phone positions in columns, m
sRate     = ishSRate;             % sampling rate, samples/s
c         = ishSoundSpeed;        % speed of sound, m/s

if (size(p, 1) < 1), return; end

%%%%%%%%%%%%%%%%%%%%%%%%% End of  configuration %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if (weighting), weights = hamming(size(p, 1)).';  % Hamming weights
else            weights = ones(1, size(p,1));     % equal weights (i.e., none)
end

g = (0:1:360) * pi/180;               % angles to calculate beam response

ss = zeros(size(g,2), size(ang,1));   % for averaging over freqs
for f = freqs
  
  dx = p(:,1) - p(1,1);               % distances from phone 1
  dy = p(:,2) - p(1,2);               % distances from phone 2
  
  df = dx * sin(ang.') + dy * cos(ang.');
  delays = df / c * sRate;            % delay in samples
  delays = round(delays);             % ...and integer samples
  cycles = delays / sRate * f;
  
  gdelays = (dx * sin(g) + dy * cos(g)) / c * sRate;
  
  clf
  s = zeros(length(g), length(ang));
  for i = 1 : length(ang)
    ddiff = gdelays - repmat(delays(:,i), 1, size(gdelays,2));
    tdiff = ddiff / sRate;
    theta = tdiff * f * 2 * pi;
    s(:,i) = abs(weights * exp(j * theta)).';
  end
  ss = ss + s;
end
ss = ss / length(f);                  % average

% Divide by amount of noise (which is sqrt(nphones)), to get SNR.
ss = ss / sqrt(size(p,1));

% Scale so maximum gain is 1.
ss = ss / max(max(ss));

% Convert to decibels.
sdB = 10 * log10(ss);    % 10 log10(amplitude) = 20 log10(intensity)

% Add points at 0 and minval dB to have plot always scale itself to that 
% range. minval depends on the number of phones.
minval = -20 * log10(sqrt(size(p,1))) - 5;
sdB = [sdB ones(length(g),1) * [minval 0]];
g1 = [g.' * ones(1,length(ang)), zeros(length(g),2)];

% Plot it!
if (sscanf(version, '%d') >= 6)
  set(polar(g1, max(minval, sdB), 'auto', zeroIsY), 'LineWidth', 1.5);
else
  set(polar(g1, max(minval, sdB) - minval), 'LineWidth', 1.5);
end

% Title.
if (length(freqs) == 1), fstr = sprintf('freq %g Hz', freqs);
else fstr = sprintf('freqs %g-%g Hz', min(freqs), max(freqs));
end
title(['Angle(s)' sprintf(' %g', ang * 180/pi) ', ' fstr])
set(gcf, 'Name', 'Beam responses (dB)', 'NumberTitle', 'off');
