Minimize sidelobe level of an FIR broadband far-field antenna array

% "Antenna array pattern synthesis via convex optimization"
% by H. Lebret and S. Boyd
% (figures are generated)
%
% Designs a broadband antenna array with the far-field wave model such that:
% - it minimizes sidelobe level outside the beamwidth of the pattern
% - it has a unit sensitivity at some target direction and for some frequencies
%
% This is a convex problem (after sampling it can be formulated as an SOCP).
%
%   minimize   max |y(theta,f)|        for theta,f outside the desired region
%       s.t.   y(theta_tar,f_tar) = 1
%
% where y is the antenna array gain pattern (complex function) and
% variables are w (antenna array weights or shading coefficients).
% Gain pattern is a linear function of w: y(theta,f) = w'*a(theta,f)
% for some a(theta,f) describing antenna array configuration and specs.
%
% Written for CVX by Almir Mutapcic 02/02/06

% select array geometry
ARRAY_GEOMETRY = '2D_UNIFORM_LATTICE';
% ARRAY_GEOMETRY = '2D_RANDOM';

%********************************************************************
% problem specs
%********************************************************************
P = 2;                % number of filter taps at each antenna element
fs = 8000;            % sampling rate = 8000 Hz
T = 1/fs;             % sampling spacing
c = 2000;             % wave speed

theta_tar = 70;       % target direction
half_beamwidth = 10;  % half beamwidth around the target direction
f_low  = 1500;        % low frequency bound for the desired band
f_high = 2000;        % high frequency bound for the desired band

%********************************************************************
% random array of n antenna elements
%********************************************************************
if strcmp( ARRAY_GEOMETRY, '2D_RANDOM' )
  % set random seed to repeat experiments
  rand('state',0);

  % uniformly distributed on [0,L]-by-[0,L] square
  n = 20;
  L = 0.45*(c/f_high)*sqrt(n);
  % loc is a column vector of x and y coordinates
  loc = L*rand(n,2);

%********************************************************************
% uniform 2D array with m-by-m element with d spacing
%********************************************************************
elseif strcmp( ARRAY_GEOMETRY, '2D_UNIFORM_LATTICE' )
  m = 6; n = m^2;
  d = 0.45*(c/f_high);

  loc = zeros(n,2);
  for x = 0:m-1
    for y = 0:m-1
      loc(m*y+x+1,:) = [x y];
    end
  end
  loc = loc*d;

else
  error('Undefined array geometry')
end

%********************************************************************
% construct optimization data
%********************************************************************
% discretized grid sampling parameters
numtheta = 180;
numfreqs = 6;

theta = linspace(1,360,numtheta)';
freqs = linspace(500,3000,numfreqs)';

clear Atotal;
for k = 1:numfreqs
  % FIR portion of the main matrix
  Afir = kron( ones(numtheta,n), -[0:P-1]/fs );

  % cos/sine part of the main matrix
  Alocx = kron( loc(:,1)', ones(1,P) );
  Alocy = kron( loc(:,2)', ones(1,P) );
  Aloc = kron( cos(pi*theta/180)/c, Alocx ) + kron( sin(pi*theta/180)/c, Alocy );

  % create the main matrix for each frequency sample
  Atotal(:,:,k) = exp(2*pi*i*freqs(k)*(Afir+Aloc));
end

% single out indices so we can make equalities and inequalities
inbandInd    = find( freqs >= f_low & freqs <= f_high );
outbandInd   = find( freqs < f_low | freqs > f_high );
thetaStopInd = find( theta > (theta_tar+half_beamwidth) | ...
                     theta < (theta_tar-half_beamwidth) );
[diffClosest, thetaTarInd] = min( abs(theta - theta_tar) );

% create target and stopband constraint matrices
Atar = []; As = [];
% inband frequencies constraints
for k = [inbandInd]'
  Atar = [Atar; Atotal(thetaTarInd,:,k)];
  As = [As; Atotal(thetaStopInd,:,k)];
end
% outband frequencies constraints
for k = [outbandInd]'
  As = [As; Atotal(:,:,k)];
end

%********************************************************************
% optimization problem
%********************************************************************
cvx_begin
  variable w(n*P) complex
  minimize( max( abs( As*w ) ) )
  subject to
    % target direction equality constraint
    Atar*w == 1;
cvx_end

% check if problem was successfully solved
disp(['Problem is ' cvx_status])
if ~strfind(cvx_status,'Solved')
  return
end

fprintf(1,'The minimum sidelobe level is %3.2f dB.\n\n',...
          20*log10(cvx_optval) );

%********************************************************************
% plots
%********************************************************************
figure(1); clf;
plot(loc(:,1),loc(:,2),'o')
title('Antenna locations')
axis('square')

% plots of array patterns (cross sections for different frequencies)
figure(2); clf;
clr = { 'r' 'r' 'b' 'b' 'r' 'r' };
linetype = {'--' '--' '-' '-' '--' '--'};
for k = 1:numfreqs
  plot(theta, 20*log10(abs(Atotal(:,:,k)*w)), [clr{k} linetype{k}]);
  hold on;
end
axis([1 360 -15 0])
title('Passband (blue solid curves) and stopband (red dashed curves)')
xlabel('look angle'), ylabel('abs(y) in dB');
hold off;

% cross section polar plots
figure(3); clf;
bw = 2*half_beamwidth;
subplot(2,2,1); polar_plot_ant(abs( Atotal(:,:,2)*w ),theta_tar,bw,'f = 1000 (stop)');
subplot(2,2,2); polar_plot_ant(abs( Atotal(:,:,3)*w ),theta_tar,bw,'f = 1500 (pass)');
subplot(2,2,3); polar_plot_ant(abs( Atotal(:,:,4)*w ),theta_tar,bw,'f = 2000 (pass)');
subplot(2,2,4); polar_plot_ant(abs( Atotal(:,:,5)*w ),theta_tar,bw,'f = 2500 (stop)');
 
Calling sedumi: 3184 variables, 145 equality constraints
   For improved efficiency, sedumi is solving the dual problem.
------------------------------------------------------------
SeDuMi 1.21 by AdvOL, 2005-2008 and Jos F. Sturm, 1998-2003.
Alg = 2: xz-corrector, Adaptive Step-Differentiation, theta = 0.250, beta = 0.500
eqs m = 145, order n = 2123, dim = 3186, blocks = 1062
nnz(A) = 304322 + 0, nnz(ADA) = 21025, nnz(L) = 10585
 it :     b*y       gap    delta  rate   t/tP*  t/tD*   feas cg cg  prec
  0 :            1.06E+03 0.000
  1 :  -1.80E+00 3.57E+02 0.000 0.3363 0.9000 0.9000   1.37  1  1  9.7E+02
  2 :  -2.44E+00 8.95E+01 0.000 0.2504 0.9000 0.9000   0.55  1  1  2.4E+02
  3 :  -8.45E-01 2.29E+01 0.000 0.2559 0.9000 0.9000   3.64  1  1  2.0E+01
  4 :  -7.47E-01 1.25E+01 0.000 0.5458 0.9000 0.9000   1.99  1  1  9.5E+00
  5 :  -6.25E-01 5.52E+00 0.000 0.4414 0.9000 0.9000   1.53  1  1  3.5E+00
  6 :  -5.62E-01 2.77E+00 0.000 0.5018 0.9000 0.9000   1.33  1  1  1.6E+00
  7 :  -5.40E-01 1.72E+00 0.000 0.6218 0.9000 0.9000   1.20  1  1  9.3E-01
  8 :  -5.18E-01 9.22E-01 0.000 0.5355 0.9000 0.9000   1.15  1  1  4.8E-01
  9 :  -5.05E-01 4.43E-01 0.000 0.4809 0.9000 0.9000   1.09  1  1  2.2E-01
 10 :  -4.97E-01 1.67E-01 0.000 0.3758 0.9000 0.9000   1.05  1  1  8.3E-02
 11 :  -4.96E-01 7.94E-02 0.000 0.4769 0.9000 0.1736   1.02  1  1  4.0E-02
 12 :  -4.93E-01 1.98E-02 0.000 0.2490 0.9201 0.9000   1.02  1  1  1.0E-02
 13 :  -4.92E-01 7.99E-03 0.000 0.4039 0.9068 0.9000   1.01  1  1  4.1E-03
 14 :  -4.92E-01 3.00E-03 0.000 0.3758 0.9051 0.9000   1.00  1  1  1.6E-03
 15 :  -4.92E-01 1.12E-03 0.000 0.3741 0.9161 0.9000   1.00  1  1  5.8E-04
 16 :  -4.92E-01 4.95E-04 0.000 0.4410 0.9000 0.9091   1.00  1  1  2.6E-04
 17 :  -4.92E-01 2.03E-04 0.000 0.4093 0.9000 0.9165   1.00  2  2  1.0E-04
 18 :  -4.92E-01 6.79E-05 0.000 0.3348 0.9000 0.9154   1.00  2  2  3.5E-05
 19 :  -4.92E-01 1.77E-05 0.000 0.2614 0.9000 0.9061   1.00  2  2  9.1E-06
 20 :  -4.92E-01 4.26E-06 0.000 0.2400 0.9000 0.9013   1.00  3  3  2.2E-06
 21 :  -4.92E-01 6.81E-07 0.000 0.1600 0.9082 0.9000   1.00  4  4  3.5E-07
 22 :  -4.92E-01 3.21E-08 0.000 0.0472 0.9901 0.9900   1.00  7  8  1.6E-08
 23 :  -4.92E-01 2.34E-09 0.000 0.0730 0.9902 0.9900   1.00 24 24  1.2E-09

iter seconds digits       c*x               b*y
 23      4.0   9.3 -4.9187930504e-01 -4.9187930526e-01
|Ax-b| =   1.2e-09, [Ay-c]_+ =   2.4E-12, |x|=  9.1e-01, |y|=  7.0e+00

Detailed timing (sec)
   Pre          IPM          Post
4.500E-01    4.020E+00    1.000E-02    
Max-norms: ||b||=1, ||c|| = 1,
Cholesky |add|=1, |skip| = 3, ||L.L|| = 7975.18.
------------------------------------------------------------
Status: Solved
Optimal value (cvx_optval): +0.491879
Problem is Solved
The minimum sidelobe level is -6.16 dB.