
%%%%%%%%%%%%%%%%%%% rpam_ex1.m  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% A Matlab script file which expresses the symmetric cone programme in
% Example 1 of T. N. Davidson, ``Efficient Design of Waveforms for
% Robust Pulse Amplitude Modulation'' (to appear in the IEEE Transactions
% on Signal Processing, 2001 or 2002) in the input format required by
% the Matlab-based general purpose symmetric cone programme solver, SeDuMi.


%%%%%%%%%%%%%%%%%%% Instructions for use %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% The design specifications for the filter are made at the beginning of
% the file.  These can be changed, as desired.

% Following the specifications, is some code which translates these
% specifications into the input format required by SeDuMi. You should
% not need to change this section.

% The comments at the end of the file describe the syntax for calling
% SeDuMi and the way in which the SeDuMi output can be processed to
% extract the optimal solution. For more information on SeDuMi
% please visit http://fewcal.kub.nl/sturm/software/sedumi.html
% UPDATE, August 2006: 
% The home page for SeDuMi has moved to: http://sedumi.mcmaster.ca

%%%%%%%%%%%%%%%%% Legal Statements %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Copyright (C) July 2001.
% Timothy N. Davidson,
% Dept. Elec. & Comp. Engineering,
% McMaster University,
% Hamilton, Ontario, Canada.
% davidson@mcmaster.ca

% This program is free software; you can redistribute it and/or
% modify it under the terms of the GNU General Public License
% as published by the Free Software Foundation; either version 2
% of the License, or (at your option) any later version.

% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.

% A copy of the GNU General Public License is available on the web at:
% http://www.gnu.org/copyleft/gpl.html; or by writing to the Free Software
% Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%% The m-file begins... %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%% Data from the IS95 specification


% Oversampling factor

N=4;


% Pass-band and stop-band edges of the spectral mask

f_p=590*10^3/(4*1.2288*10^6); 
f_s=740*10^3/(4*1.2288*10^6);


% Relative spectral mask bounds specified in the standard, in dB

delta_one_dB=1.5;  % pass-band ripple
delta_two_dB=-40;  % stop-band level


% Relative spectral mask bounds achieved by the IS95 filter, in dB

%delta_one_dB=.975; % pass-band ripple
%delta_two_dB=-45;  % stop-band level


% Filter length

L=48;


% Multiplying factor for the number of frequency domain points
% at which the spectral mask will be enforced. The number of
% points is mult*L uniformly spaced points, plus the `corners'
% of the mask

mult=15;

% Half of the beta value achieved by the IS95 filter

beta_IS95_on2 = 0.24504023151123; 

% Choose a value for the bound on beta for the designed filter.

epsilon_on2=beta_IS95_on2/10; 



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% You should not need to change anything below this line. %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


% Conversion of dB spectral mask specifications 

delta_p_up=10^(delta_one_dB/10)-1;
delta_p_low=1-10^(-delta_one_dB/10);
delta_two=10^(delta_two_dB/10);


% Number of terms needed to calculate beta. 

N_isi=floor((L-1)/N);


% Number of uniformly spaced points for discretization of the spectral mask

N_f_uniform=mult*L;


% The vector f is the set of points in frequency for 
% discretization of the mask. It includes the corners of the mask
% N_f_p is the index of the pass-band edge in f
% N_f_s is the index of the stop-band edge in f

f=linspace(0,0.5,N_f_uniform+1);
f=f(1:N_f_uniform);
i_p = floor(2*N_f_uniform*f_p)+1;
i_s=ceil(2*N_f_uniform*f_s)+1;
f=[f(1:i_p),f_p,f(i_p+1:i_s-1),f_s,f(i_s:N_f_uniform)];
N_f_p=i_p+1;
N_f_s=i_s+1;
N_f=length(f);


% Now we begin to assemble the matrix A, the vectors b and c, and
% the record K.  

% The cone has the following structure: 
K.l=1+N_f_p+N_f+1+3*N_isi;
K.q=0;
K.r=L+1;
K.s=L;

% The parts of the variable x (SeDuMi notation) in the linear cone
% are (in order) zeta, slack variables for the lower bound of  the
% pass-band mask (i.e., non-negative numbers describing the amount by 
% which the mask is satisfied), slack variables for the (whole) upper 
% bound on the mask, r_g[0], the mu_i's,  slacks for the lower bound on
% r_g[Ni], and slacks for the upper bound on r_g[Ni]. 

% Note that r_g[0] could have been eliminated from the formulation
% analytically, but it has been left in for paedegogical reasons.

% The rotated second-order cone contains r_g[1],...,r_g[L-1]

% The semidefinite cone contains the elements of X.

% Set r_g[0] = 1;
A=sparse([zeros(1,1+N_f_p+N_f),1,zeros(1,3*N_isi+1),0,zeros(1,L-1+L^2)]);
b=1;


% Set the scaling of the rotated second-order cone to 2r_g[0]

A=[A;sparse([zeros(1,1+N_f_p+N_f),.5,zeros(1,3*N_isi+1),-1,zeros(1,L-1+L^2)])];
b=[b;0];


% Enforce the trace component of (9)

A=[A;sparse([zeros(1,1+N_f_p+N_f),1,zeros(1,3*N_isi+L+1),-vec(eye(L))'])];
b=[b;0];


% Enforce off-diagonal component of (9)

for k=1:L-1,
   I_offset = [zeros(L-k,k),eye(L-k);zeros(k,L)];
   A=[A;sparse([zeros(1,1+N_f_p+N_f+1+3*N_isi+2+k-1),1,zeros(1,L-1-k),-0.5*(vec(I_offset)'+vec(I_offset')')])];
end;

b=[b;zeros(L-1,1)];

% Enforce the bound on the sum on the right hand side of (11)  

A=[A;sparse([zeros(1,1+N_f_p+N_f),-epsilon_on2,ones(1,N_isi),zeros(1,2*N_isi+L+1+L^2)])];
b=[b;0];

% Now make the slack variables equal to the true slacks 

for ii=1:N_isi,
  kk=N*ii;
A=[A;sparse([zeros(1,1+N_f_p+N_f+1),zeros(1,ii-1),-1,zeros(1,N_isi-ii),zeros(1,ii-1),1,zeros(1,N_isi-ii),zeros(1,N_isi),zeros(1,2),zeros(1,kk-1),1,zeros(1,L-1-kk),zeros(1,L^2)])];

A=[A;sparse([zeros(1,1+N_f_p+N_f+1),zeros(1,ii-1),1,zeros(1,N_isi-ii),zeros(1,N_isi),zeros(1,ii-1),-1,zeros(1,N_isi-ii),zeros(1,2),zeros(1,kk-1),1,zeros(1,L-1-kk),zeros(1,L^2)])];
end;

b=[b;zeros(2*N_isi,1)];


%%% Now enforce the mask constraints

%two vectors used to calculate R(e^{j\omega})

freq_multiplier = 1:L-1;
cos_matrix_primal=cos(2*pi*f'*freq_multiplier);


% Make the slack variables for the lower bound in the  pass-band 
% equal to the true slacks

A=[A;sparse([-10^(-delta_one_dB/10)*ones(N_f_p,1),-eye(N_f_p),zeros(N_f_p,N_f),ones(N_f_p,1),zeros(N_f_p,3*N_isi+2),2*cos_matrix_primal(1:N_f_p,:),zeros(N_f_p,L^2)])];
b=[b;zeros(N_f_p,1)]; 

% Note that the lower bound in the stop band is equal to zero
% and hence is dealt with by the Positive Real Lemma and constraint (9)


% Now deal with the slacks for the upper bound. 

upper_mask_bound=[10^(delta_one_dB/10)*ones(N_f_s-1,1);10^(delta_two_dB/10)*ones(N_f- N_f_s+1,1)];

A=[A;sparse([upper_mask_bound,zeros(N_f,N_f_p),-eye(N_f),-ones(N_f,1),zeros(N_f,3*N_isi+2),-2*cos_matrix_primal,zeros(N_f,L^2)])];

b=[b;zeros(N_f,1)];


% Finally, set c for the objective.  

c=[zeros(1+N_f_p+N_f+1+3*N_isi,1);1;zeros(L+L^2,1)];




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%% Call to SeDuMi %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Now that the matrix A, the vectors b and c, and the record K have
% been set, the problem can be solved using:
%
% [x_opt,y_opt,sedumi_info]=sedumi(A,b,c,K);




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%% Extracting the optimal solution %%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% If the filter specifications are feasible, as recorded in
% the sedumi_info record above (see the SeDuMi Users' Guide for the
% details), then the components of the optimal solution can be
% extracted as follows:

% A vector  containing the optimal r_g[0], ..., r_g[L-1] is:
%
% r_g_opt = [x_opt(1+N_f_p+N_f+1);x_opt(1+N_f_p+N_f+1+3*N_isi+3:1+N_f_p+N_f+1+3*N_isi+3+L-2)];
%
% The optimal zeta is: 
% zeta_opt=x_opt(1);
%
% The optimal value of B_g is 
%
% B_g_opt = sqrt(2*c'*x_opt);
