%  Array beam pattern calculation using equations from
%  Digital Signal Processing
%  DeFatta, Lucas, Hodgkiss
%  0 and 180 degrees are broadside
%  90 and -90 degrees are endfire
%
clear
c = 1500.0;
d2r = pi / 180;
num_phones = 4;

freq = [100 200 500 1000 2000 5000 10000 20000 50000];
wave_length = c ./ freq;
spacing = min(wave_length) / 2.0;
k = (2.0 * pi) ./ wave_length;

%weight = hanning(num_phones)';
weight = [1 1 1 1]';

n = 0:(num_phones - 1);
steer_angle = 90.0;

for arrival_angle = 0:180,
  pattern(arrival_angle + 1,1) = sum(weight .* exp(-1 * j * n * k(1) * ...
          spacing * (sin(steer_angle * d2r) - sin(arrival_angle * d2r))));
  pattern(arrival_angle + 1,2) = sum(weight .* exp(-1 * j * n * k(2) * ...
          spacing * (sin(steer_angle * d2r) - sin(arrival_angle * d2r))));
  pattern(arrival_angle + 1,3) = sum(weight .* exp(-1 * j * n * k(3) * ...
          spacing * (sin(steer_angle * d2r) - sin(arrival_angle * d2r))));
  pattern(arrival_angle + 1,4) = sum(weight .* exp(-1 * j * n * k(4) * ...
          spacing * (sin(steer_angle * d2r) - sin(arrival_angle * d2r))));
end  
temp = 0:180;
theta = temp .* pi / 180.0;
normalizer = max(max(abs(pattern)));
rho = 25 + (20 .* log10(abs(pattern)/normalizer));
for index = 1:181,
  for jndex = 1:4,
    if( rho(index, jndex) < 0 )
      rho(index, jndex) = 0;
    end
  end
end
polar(theta', rho(:,1), 'r')
title('8 kHz PATTERN (solid = actual dashed = desired)')
grid