%%%%%%%%%%%%%%%%%%% 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);