%%%%%%%%%%%%%%%%%%% Equation39.m  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% A Matlab script file which expresses the symmetric cone programme in
% Equation (39) of Wilbur, Davidson and Reilly, ``Efficient Design of 
% Oversampled NPR GDFT Filter Banks'' (accepted for publication, 
% subject to minor revisions, in the IEEE Transactions on Signal 
% Processing, July 2003) in the input format required by
% the Matlab-based general purpose symmetric cone programme solver, SeDuMi.
%
% Update: Equation (39) of the technical report appears as equation (38)
% of the paper in the IEEE Transactions on Signal Processing, 
% 52(7):1947-1963, July 2004.

%%%%%%%%%%%%%%%%%%% 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.

% To solve the optimization problem, you will need to have the
% SeDuMi solver installed on your system. To obtain a copy, and
% for more details 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 

% At the end of the file, there are a few simple commands for
% extracting the optimal autocorrelation and the optimized
% objective value from the output provided by SeDuMi.

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

% Copyright (C) July 2003.
% 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... %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


L=48;    % Filter length
M=8;     % Number of subbands 
Kds=6;   % Downsampling factor

% Now set up the contraints. The default constraints generate
% the point denoted by the triangle in Figure 3.

passband_ripple_dB=1;       % Allowable pass-band ripple
                            % See the discussion in Example 1.

rel_stopband_levels_dB=-30; % The relative stop-band level constraint
                            % This can be a vector if desired

dist_bound_factors=1e-06;   % The normalized distortion coefficient
                            % This can be a vector if desired

normalized_TBEs=1;          % The normalized transition band energy 
                            % constraint. 
                            % This can be a vector if desired

                            % If this constraint is set to 1 then it
                            % will be inactive. 

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

% Calculate the maximum spectral component constraint

passband_level_up_dB=10*log10(Kds) + passband_ripple_dB;

% Calculate the stop-band edge

f_s=1/(2*Kds);
f_s_for_mask=f_s;
f_s_for_energy=f_s;

% Calculate the number of elements of r_p which contribute to the
% distortion coefficient

N_dist=floor((L-1)/M);

% Set up some indicies which simplify the expressions for the matrix A

no_diagsX=L;
no_diagsZ=L-1;

no_elemsX=no_diagsX^2;
no_elemsZ=no_diagsZ^2;

startRcone=1+3+1 +1; 
start_sumoffdiagsXup_pb=startRcone+L;
start_sumoffdiagsXup_sb=start_sumoffdiagsXup_pb+L;
start_sumoffdiagsZup_sb=start_sumoffdiagsXup_sb+L;
start_dist_cone=start_sumoffdiagsZup_sb+L-1;

start_Xpr= start_dist_cone+N_dist+1;
start_Xup_pb=start_Xpr+no_elemsX;
start_Xup_sb=start_Xup_pb+no_elemsX;
start_Zup_sb=start_Xup_sb+no_elemsX;
no_vars=start_Zup_sb+no_elemsZ-1;

%Set up some convenient matrices

trace_L = vectransT(0,L-1);
trace_Lmin1 = vectransT(0,L-2);

all_sumoffdiags_L =[];
for ii=1:L-1
  all_sumoffdiags_L=[all_sumoffdiags_L;vectransT(ii,L-1)];
end;

all_sumoffdiags_Lmin1 =[];
for ii=1:L-2
  all_sumoffdiags_Lmin1=[all_sumoffdiags_Lmin1;vectransT(ii,L-2)];
end;

eye_Lmin1=eye(L-1);
dist_selector=eye_Lmin1(M*(1:N_dist),:);

% Now loop over the number of relative stop-band level constraints,
% the number of transition band energy constraints, and the
% number of normalized distortion constraints 

for count_SBLs=1:length(rel_stopband_levels_dB);

rel_stopband_level_dB=rel_stopband_levels_dB(count_SBLs);
stopband_level_dB=passband_level_up_dB +rel_stopband_level_dB; 

for count_TBEs=1:length(normalized_TBEs),

for count_dist_factors=1:length(dist_bound_factors);

dist_bound= sqrt(dist_bound_factors(count_dist_factors)*Kds/M/2);

% Now construct the matrix A, the vectors b and c, and the
% cone descriptor K required by SeDuMi

% Here's the cone descriptor.  The second to fifth quadratic
% cones are actually redundant. However, using them has the
% benefit that it reduces the iteration of the semidefine 
% matrix variables with the rest of the problem. This helps
% to keep the problem sparse. SeDuMi is quite good at exploiting
% sparsity in order to reduce the computational cost of solving
% the problem. 

  K.f=0;
  K.l=1+3+1;
  K.q=[L,L,L,L-1,N_dist+1];
  K.s=[L,L,L,L-1];

%Now we set up A and b to enforce equality constraints.

% Set r_p[0]=K/M;
A = sparse([1,zeros(1,no_vars-1)]);
b=Kds/M;

% Enforce the equality constraints between r_p and the X
% which enforces positive definiteness

A=[A;sparse([1,zeros(1,start_Xpr-2),-trace_L,zeros(1,no_vars-start_Xup_pb+1)])];
b=[b;0];

A=[A;sparse([zeros(L-1,startRcone),eye(L-1),zeros(L-1,start_Xpr -start_sumoffdiagsXup_pb),-all_sumoffdiags_L,zeros(L-1,no_vars-start_Xup_pb+1)])];
b=[b;zeros(L-1,1)];

% Do the same for the X which enforces the maximum spectral component
% constraint and the corresponding cone

temp_offset=1;

A=[A;sparse([zeros(1,temp_offset),1,zeros(1,start_Xup_pb-temp_offset-2),-trace_L,zeros(1,no_vars-start_Xup_sb+1)])];
b=[b;0];


A=[A;sparse([zeros(L-1,start_sumoffdiagsXup_pb),eye(L-1),zeros(L-1,start_Xup_pb-start_sumoffdiagsXup_sb),-all_sumoffdiags_L,zeros(L-1,no_vars-start_Xup_sb+1)])];
b=[b;zeros(L-1,1)];


% Do the same for the X and Z which enforce the maximum stop-band level 
% constraint 

% For the X part
temp_offset=2;
A=[A;sparse([zeros(1,temp_offset),1,zeros(1,start_Xup_sb-temp_offset-2),-trace_L,zeros(1,no_vars-start_Zup_sb+1)])];
b=[b;0];

A=[A;sparse([zeros(L-1,start_sumoffdiagsXup_sb),eye(L-1),zeros(L-1,start_Xup_sb-start_sumoffdiagsZup_sb),-all_sumoffdiags_L,zeros(L-1,(L-1)^2)])];
b=[b;zeros(L-1,1)];

% For the Z part
temp_offset=3;
A=[A;sparse([zeros(1,temp_offset),1,zeros(1,start_Zup_sb-temp_offset-2),-trace_Lmin1])];
b=[b;0];

A=[A;sparse([zeros(L-2,start_sumoffdiagsZup_sb),eye(L-2),zeros(L-2,start_Zup_sb-start_dist_cone),-all_sumoffdiags_Lmin1])];
b=[b;zeros(L-2,1)];


% Now enforce the relationships between the redundant quadratic cones
% and the cone which holds r_p

% For the maximum spectral component 
temp_offset=0;

A=[A;sparse([-1,zeros(1,temp_offset),-1,zeros(1,no_vars-temp_offset-2)])];
b=[b;-10^(passband_level_up_dB/10)];

A=[A;sparse([zeros(L-1,startRcone),-eye(L-1),zeros(L-1,start_sumoffdiagsXup_pb-start_sumoffdiagsXup_pb+1),-eye(L-1),zeros(L-1,no_vars-start_sumoffdiagsXup_sb+1)])];
b=[b;zeros(L-1,1)];

% For the stop-band level constraint

alpha=2*pi*f_s_for_mask;
beta=2*pi*(1-f_s_for_mask);

aa=sin(alpha)/(1-cos(alpha));
bb=sin(beta)/(1-cos(beta));
c0 = -(1+aa*bb);
c1= aa*bb -1 -(aa+bb)*sqrt(-1);
c1=real(c1);


part_matrix_of_Cs=zeros(L-1,L-2);
for ii=1:L-1
  if ii <=L-2
    part_matrix_of_Cs(ii,ii)=2*c0;
  end;
  if ii<= L-3
   part_matrix_of_Cs(ii,ii+1)=c1;
  end;
  if ii>=2,
    part_matrix_of_Cs(ii,ii-1)=conj(c1);
  end;
end;

temp_offset=1;

A=[A;sparse([-1,zeros(1,temp_offset),-1,-c0,zeros(1,start_sumoffdiagsZup_sb-temp_offset-3),-c1,zeros(1,L-3),zeros(1,no_vars-start_dist_cone+1)])];
b=[b;-10^(stopband_level_dB/10)];

A=[A;sparse([zeros(L-1,temp_offset+2),[-conj(c1);zeros(L-2,1)],zeros(L-1,startRcone-temp_offset-3),-2*eye(L-1),zeros(L-1,start_sumoffdiagsXup_sb-start_sumoffdiagsXup_pb+1),-2*eye(L-1),zeros(L-1,1),-part_matrix_of_Cs,zeros(L-1,no_vars-start_dist_cone+1)])];
b=[b;zeros(L-1,1)];

% Now set the bound on the cone which contains copies of the elements of r_p
% which cause distortion

A=[A;sparse([zeros(1,start_dist_cone-1),1,zeros(1,no_vars-start_dist_cone)])];
b=[b;dist_bound];

% Now ensure that equality between the elements of that cone and
% the true elements of r_p

A=[A;sparse([zeros(N_dist,startRcone),-dist_selector,zeros(N_dist,start_dist_cone-start_sumoffdiagsXup_pb+1),eye(N_dist),zeros(N_dist,no_vars-start_Xpr+1)])];
b=[b;zeros(N_dist,1)];

% Now set up the constraint on the transition band energy
% In order to do this via equality constraints, we need a "slack" variable

m_vec=1:L-1;
m_vec=m_vec(:);
one_on_m_vec=1./m_vec;

A=[A;sparse([1/Kds-1/M-normalized_TBEs(count_TBEs),zeros(1,startRcone-3),1,0,2/pi*(one_on_m_vec.').*(sin(pi*(m_vec.')/Kds)-sin(pi*(m_vec.')/M)),zeros(1,no_vars-start_sumoffdiagsXup_pb+1)])];
b=[b;0];


% Construct c so that c'*x is (half) the stop-band energy

c=[1/2-f_s_for_energy;zeros(startRcone-1,1);-1/pi*sin(2*pi*m_vec*f_s_for_energy).*one_on_m_vec;zeros(no_vars-start_sumoffdiagsXup_pb+1,1)];

% The "pars" parameter can be used to control the way in which
% SeDuMi operates. See the SeDuMi Users' Guide

pars=[];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%% 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,y,info]=sedumi(A,b,c,K,pars);

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

% If the filter specifications are feasible, as recorded in
% the info record above (see the SeDuMi Users' Guide for the
% details), then extract the solution

if info.pinf==0 & info.dinf==0 & info.numerr<2,

% The achieved normalized SBE is stored in a matrix so that trade-offs can be
% easily plotted   

  SBEvDist_tradeoff(count_SBLs,count_dist_factors,count_TBEs)= 2*c'*x/(Kds/M);

% You may also wish to store the optimal autocorrelation for each
% set of constraints.
% The following command stores the optimal r_p[0], ..., r_p[L-1] 

  optimal_autocorrelations(count_SBLs,count_dist_factors,count_TBEs).non_neg_part=[x(1);x(startRcone+1:startRcone+L-1)];


else

% SeDuMi could not (reliably) find a feasible solution

  SBEvDist_tradeoff(count_SBLs,count_dist_factors,count_TBEs)= inf;
  optimal_autocorrelations(count_SBLs,count_dist_factors,count_TBEs).non_neg_part=NaN;
end;

end; %count_Dist_factors
end; %count_TBEs
end; %count_SBLs