function rob_reject_constr()
% Matlab script for the paper "Array Signal Processing with
% Robust Rejection Constraints via Second-Order Cone Programming",
% by A. Mutapcic, S.-J. Kim, and S. Boyd.
% In Proceedings of the 40th Asilomar Conference on Signals, Systems,
% and Computers, Oct. 29 - Nov. 1, 2006, Pacific Grove, CA.
% 
% This script requires CVX software, which is available at
% http://www.stanford.edu/~boyd/cvx/

% set random seed to repeat experiments
rand('state',1);
randn('state',0);

% angle sampling 
theta  = [1:1:360]';

%********************************************************************
% antenna specs
%********************************************************************
lambda = 1;           % wavelength
theta_tar = 45;       % desired (target) direction

% interference directions
theta_rej = [10 65 120 200];
% rejection level
EPS = 0.05; % 20*log10(0.05) = -26 dB

% random array of n antenna elements
% (uniformly distributed on [0,L]-by-[0,L] square)
n = 20;
L = 5;
loc = L*rand(n,2);

%********************************************************************
% construct optimization data
%********************************************************************
% build matrix A that relates w and y(theta), ie, y = A*w
A = kron(cos(pi*theta/180), loc(:,1)') + kron(sin(pi*theta/180), loc(:,2)');
A = exp(2*pi*i/lambda*A);

% desired angle constraint
Atar = A(theta_tar,:);

% stopband constraint matrix
Areject = [];
Arej = {};
Prej = {};

for k = 1:length(theta_rej)
  Areject = [Areject; A(theta_rej(k),:)];
  Arej{k} = A(theta_rej(k),:)';
  Prej{k} = eye(n);
end


%********************************************************************
% nominal optimal solution
%********************************************************************
cvx_begin
  variable wnom(n) complex
  minimize( norm(wnom) )
  subject to
    real( Atar*wnom ) >= 1;
    for k = 1:length(Arej)
      abs( Arej{k}'*wnom ) <= EPS;
    end
cvx_end

% check if problem was successfully solved
disp(['Problem is ' cvx_status])
if ~strcmp(cvx_status,'Solved')
  error('Check problem solution...')
end

fprintf(1,'The nominal max interference angle gain is %3.2f dB.\n\n',...
          20*log10( max(abs( Areject*wnom )) ) );


wc_rej_array = [];
wc_rej_array_nom = [];
obj_value_array = [];

rho = linspace(0,0.2,20);
for r = 1:length(rho)
  fprintf(1,'processing rho = %f\n',rho(r))
  %********************************************************************
  % robust optimal solution
  %********************************************************************
  cvx_begin
    variable w(n) complex
    minimize( norm(w) )
    subject to
      real( Atar*w ) >= 1;

      for k = 1:length(Arej)
        sup_ellip_abs( w, Arej{k}, rho(r)*Prej{k} ) <= EPS;
      end

  cvx_end

  % check if problem was successfully solved
  disp(['Problem is ' cvx_status])
  if ~strcmp(cvx_status,'Solved') && ~strcmp(cvx_status,'Inaccurate/Solved')
    disp('Check problem solution...')
    continue
  end

  obj_value_array = [obj_value_array norm(w)];

  fprintf(1,'The max interference angle gain is %3.2f dB.\n\n',...
            20*log10( max(abs( Areject*w )) ) );

  %********************************************************************
  % check the worst-case performance
  %********************************************************************
  wc_rej     = [];
  wc_rej_nom = [];

  cvx_quiet(true);
  for k = 1:length(Arej)
    wc_rej(k,1)     = sup_ellip_abs( w, Arej{k}, rho(r)*Prej{k} );
    wc_rej_nom(k,1) = sup_ellip_abs( wnom, Arej{k}, rho(r)*Prej{k} );
  end
  cvx_quiet(false);

  wc_rej_array     = [wc_rej_array max(wc_rej)];
  wc_rej_array_nom = [wc_rej_array_nom max(wc_rej_nom)];

end % for loop

%********************************************************************
% plots
%********************************************************************
% worst-case rejection level versus rho
figure(1);
plot(rho,obj_value_array,rho,norm(wnom)*ones(length(rho),1),'r--');
xlabel('rho'), ylabel('Pw')
set(gca,'FontSize',14);
axis([0 0.2 0.235 0.25]), grid on,

% objective value versus rho
figure(2);
plot(rho,wc_rej_array,rho,wc_rej_array_nom,'r--');
xlabel('rho'), ylabel('wc-null')
grid on
set(gca,'FontSize',14);

%********************************************************************
% helper function for the robust gain constraint
%********************************************************************
function cvx_optval = sup_ellip_abs(w,a_nom,P)
% SUP_ELLIP_ABS robust absolute value function over an ellipsoid.
%
% Gives an epigraph formulation of the following convex function g(w):
%
%   g(w) = sup_{ a \in E } |w*a| 
%
% where ellipsoid E = { Pu + a_nom | ||u|| <= 1 }.

cvx_optval = abs(a_nom'*w) + norm(P'*w);