function [p_min, u_min] = L2n(s_size, Np, Nu, levelc, model, ini_point, x00, sf, dm, ws, wp, outname)% optimize the coding gain subject to the error of h1 is less than dh1% with second-order cone programming% input:  filename -- size%               Np -- number of dual vanishing moments%               Nu -- number of primal vanishing moments%           levelc -- number of decomposition levels%            model -- isotropic (0) or separable (1) model%        ini_point -- initial point, 0 - Kovacevic's, 1 - random, 2 - fixed point (x00)%               sf -- scaling factor for the freq. response error function, can be 0.9, 1, or 10(no constraint)%               dm -- type of downsampling matrix, 0 - regular; 1 - twin dragon %               ws -- weight for the stopband%               wp -- weight for the passband%          outname -- name of the output file (workspace)% output:    p_min -- lifting filter coefficients for the first lifting filter%            u_min -- lifting filter coefficients for the second lifting filter% for example: %     [p, u] = L2n([2 2 2 2], 2, 2, 4, 0, 0, [], 1, 0, 1, 0, 'data/test')% Copyright (c) 2006 Yi Chenformat shortsyms a01 a02 a03 a04 a05 a06 a07 a08 a09 a10 a11 a12 a13 a14 a15 a16 a17 a18 a19 a20 a21 a22 a23 a24 a25 a26 a27 a28syms a29 a30 a31 a32 a33 a34 a35 a36 a37 a38 a39 a40 a41 a42 a43 a44 a45 a46 a47 a48 a49 a50 a51 a52 a53 a54 a55 a56 a57 a58 a59 a60 a60  a61 a62 a63 a64xa = [a01 a02 a03 a04 a05 a06 a07 a08 a09 a10 a11 a12 a13 a14 a15 a16 a17 a18 a19 a20 a21 a22 a23 a24 a25 a26 a27 a28...        a29 a30 a31 a32 a33 a34 a35 a36 a37 a38 a39 a40 a41 a42 a43 a44 a45 a46 a47 a48 a49 a50 a51 a52 a53 a54 a55 a56 a57 a58 a59 a60 a60  a61 a62 a63 a64];px = s_size(1);py = s_size(2)/2;ux = s_size(3);uy = s_size(4)/2;a1 = xa(1:px*py);a2 = xa(px*py+1:px*py+ux*uy);    x = [a1 a2]; A1c = reshape([fliplr(a1) a1], [px 2*py]);A2c = reshape([fliplr(a2) a2], [ux 2*uy]);lfilter = A1c;sizef = size(lfilter);df = - [sizef(1)-2, sizef(2)-2] / 2;f(1) = matrix2poly(lfilter,df);lfilter = A2c;sizef = size(lfilter);df = - [sizef(1), sizef(2)] / 2;f(2) = matrix2poly(lfilter,df);% compute the analysis filters (in polynomial form)[h0,h1] = lifting_quincunx(f, dm);% convert to matrix form and generate the frequency response[dh0, h0_coeff] = poly2matrix_sym(h0);[dh1, h1_coeff] = poly2matrix_sym(h1);delta_g = [0.001, 0.0001, 0.00001];delta_h1 = 1.5;beta = 0.001;if dm == 0    Mquin = [1 1; 1 -1]elseif dm == 1    Mquin = [1 -1; 1 1]endsyms w0 w1for nx = 1:px    for ny = 1:py        v_tilde(nx + px*(ny-1)) = 2*cos((nx-ny-px/2)*w0+(nx+ny-px/2-1)*w1);    endendv_tilde = v_tilde(:);% linear constraints on the lifting filters with diamond support[xs, Vr, Aeq, beq] = linear_constr2(f,Np,Nu);if rank([Aeq beq]) > rank(Aeq)    disp('not enough filter coefficients')    returnendI_hat = [eye(px*py) zeros(px*py, ux*uy)];xs_hat = I_hat*xs;Vr_hat = I_hat*Vr;N = 64;M1 = 15;M2 = 15;% compute the optimization parameters[H, S, C, H_hat, S_hat, C_hat]= L2n_parameters(N, M1, M2, ws, wp, Vr_hat, v_tilde, xs_hat);if ini_point == 0 % using Kovacevic's filter bank as the initial point        if s_size == [2 2 2 2]        xx(1,:) =  ([-2 -2 1 1]'/8)';  % 2/2    elseif s_size == [4 4 2 2]        xx(1,:) =  ([1 -10 -10 1 0 1 1 0 4 4]'/32)';  % 4/2    elseif s_size == [4 4 4 4]        xx(1,:) =  ([2 -20 -20 2 0 2 2 0 -1 10 10 -1 0 -1 -1 0]'/64)';  % 4/4    elseif s_size == [6 6 2 2]        xx(1,:) = [-3 27 -174 -174 27 -3 0 -2 27 27 -2 0 0 0 -3 -3 0 0 64 64]/2^9;  % 6/2       elseif s_size == [6 6 4 4]        xx(1,:) = [[-3 27 -174 -174 27 -3 0 -2 27 27 -2 0 0 0 -3 -3 0 0]/2^9 [-1 10 10 -1 0 -1 -1 0]/64];  % 6/4     elseif s_size == [6 6 6 6]        xx(1,:) =  [-3*2 27*2 -174*2 -174*2 27*2 -3*2 0 -2*2 27*2 27*2 -2*2 0 0 0 -3*2 -3*2 0 0 ...                3 -27 174 174 -27 3 0 2 -27 -27 2 0 0 0 3 3 0 0]/2^10;  % 6/6    elseif s_size == [8 8 2 2]        xx(1,1:32) = -[-80 850 -4470 23300 23300 -4470 850 -80 0 -75 625 -4470 -4470 625 -75 0 0 9 -75 850 850 -75 9 0 0 0 0 -80 -80 0 0 0]/2^16;        xx(1,33:34) = [1 1]/8;    elseif s_size == [8 8 4 4]        xx(1,1:32) = -[-80 850 -4470 23300 23300 -4470 850 -80 0 -75 625 -4470 -4470 625 -75 0 0 9 -75 850 850 -75 9 0 0 0 0 -80 -80 0 0 0]/2^16;        xx(1,33:40) = [-1 10 10 -1 0 -1 -1 0]/64;    elseif s_size == [8 8 6 6]        xx(1,1:32) = -[-80 850 -4470 23300 23300 -4470 850 -80 0 -75 625 -4470 -4470 625 -75 0 0 9 -75 850 850 -75 9 0 0 0 0 -80 -80 0 0 0]/2^16;        xx(1,33:50) = [3 -27 174 174 -27 3 0 2 -27 -27 2 0 0 0 3 3 0 0]/2^10;    elseif s_size == [8 8 8 8]        xx(1,1:32) = -[-80 850 -4470 23300 23300 -4470 850 -80 0 -75 625 -4470 -4470 625 -75 0 0 9 -75 850 850 -75 9 0 0 0 0 -80 -80 0 0 0]/2^16;        xx(1,33:64) = -xx(1,1:32)/2;    else        xx(1, :) = zeros(1, length(x));        ini_point = 1;    end    phi(1,:) = (Vr'*((xx(1,:))' - xs))';endif ini_point == 1 % random initial points    rand('state',sum(100*clock))    phi(1,:) = rand(1, length(Vr(1,:)))-0.5;    xx(1, :) = (xs + Vr*phi(1,:)')';elseif ini_point == 2    x00 = x00(:)';    xx(1,:) = x00;    phi(1, :) = (Vr'*((xx(1,:))' - xs))';elseif ini_point == 3    rand('state',sum(100*clock))    [A_temp, b_temp] = linear_constr2_qs(Aeq, beq, px, py, ux, uy);    r_temp = rank(A_temp);    [U_temp, S_temp, V_temp] = svd(A_temp);    Vr_temp = V_temp(:, r_temp+1:length(V_temp(1,:)));    xs_temp = pinv(A_temp)*b_temp;    phi_temp = rand(1, length(Vr_temp(1,:)))-0.5;    xx(1, :) = (xs_temp + Vr_temp*phi_temp')';endInitialPoint = xx(1,:);x_min = xx(1, :)';t = cputime;iter = 1;if ini_point == 0    ini_step = 3;else     ini_step = 1;endfor ini_iter = ini_step:3    t = cputime;    delta = delta_g(3);    beta = 50*norm(10^(-ini_iter-2)*ones(length(Vr(1,:)),1),2);      level = levelc;        xx(iter,:) =  x_min';    phi(iter,:) = (Vr'*((xx(iter,:))'- xs))';        [h0_coeff, h1_coeff] = analysis_filters2(A1c, A2c, x, xx(iter, :), Mquin);    Gmax(iter, :) = CodingGain_num(h0_coeff, h1_coeff, 0.95, level, model, Mquin);    diff_H1(iter) = phi(iter,:)*H*(phi(iter,:))'+phi(iter,:)*S+C;    lengthp = length(Vr(1,:));    Gt(iter) = -Gmax(iter,level);        f_min = 100;    diff_Gmax = 1;    diff_Gmax2 = 1;    diff_Gmax3 = 0;        while (abs(diff_Gmax) > 10^(-ini_iter-1)) && (abs(diff_Gmax2) > 10^(-ini_iter-1) && (diff_Gmax3 < 0.01))        t1 = cputime;                % compute the gradient        for i = 1:lengthp            phi_temp = zeros(lengthp,1);            phi_temp(i) = delta;            x_temp = (xx(iter,:))' + Vr*phi_temp;            [h0_coeff, h1_coeff] = analysis_filters2(A1c, A2c, x, x_temp, Mquin);            G_temp = -CodingGain_num(h0_coeff, h1_coeff, 0.95, level, model, Mquin);            G_temp = G_temp(level);            phi_temp = zeros(lengthp,1);            phi_temp(i) = -delta;            x_temp = (xx(iter, :))' + Vr*phi_temp;            [h0_coeff, h1_coeff] = analysis_filters2(A1c, A2c, x, x_temp, Mquin);            G_temp2 = -CodingGain_num(h0_coeff, h1_coeff, 0.95, level, model, Mquin);            G_temp2 = G_temp2(level);            g(i) = (G_temp - G_temp2)/2/delta;        end        % gradient around xx0        g = g(:);                % variable bound for the error function of H1        delta_h1 = sf*diff_H1(iter);        S_hat0 = H_hat*phi(iter,:)' + S_hat;        C_hat0 = phi(iter, :)*H*phi(iter, :)' + phi(iter, :)*S + C - S_hat0'*S_hat0;                % SOCP        b = g;                A1 = H_hat';        c1 = S_hat0;        b1 = zeros(lengthp,1);        d1 = sqrt(delta_h1 - C_hat0);                A2 = eye(lengthp);        c2 = zeros(lengthp,1);        b2 = zeros(lengthp,1);        d2 = beta;        %         with constraints on frequency response        At1 = -[b1 A1];        At2 = -[b2 A2];        At = [At1 At2];        bt = -b;        ct1 = [d1; c1];        ct2 = [d2; c2];        ct = [ct1; ct2];        K.q = [size(At1,2) size(At2,2)];        pars.fid = 0;        [xsp, delta_phi, info] = sedumi(At, bt, ct, K, pars);        info;                phi(iter+1, :) = phi(iter,:) + delta_phi';        xx(iter+1, :) = (xs + Vr*phi(iter+1,:)')';                min = g'*delta_phi; % + 0.5*delta_phi'*hessian*delta_phi;                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%        [h0_coeff, h1_coeff] = analysis_filters2(A1c, A2c, x, xx(iter+1, :), Mquin);        Gmax(iter+1, :) = CodingGain_num(h0_coeff, h1_coeff, 0.95, level, model, Mquin);        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%           Gt(iter+1) = -Gmax(iter+1,level);        diff_Gmax = Gt(iter+1) - Gt(iter);        if iter > 1            diff_Gmax2 = Gt(iter+1) - Gt(iter-1);        end                diff_H1(iter+1) = phi(iter+1, :)*H*phi(iter+1, :)'+phi(iter+1, :)*S+C;        (H_hat*phi(iter+1, :)'+S_hat)'*(H_hat*phi(iter+1, :)'+S_hat) + C_hat;        (H_hat*delta_phi+S_hat0)'*(H_hat*delta_phi+S_hat0) + C_hat0;                if Gt(iter+1) < f_min            f_min = Gt(iter+1);            x_min = xx(iter+1, :);            G_max = -Gt(iter+1);        end        t2 = cputime - t1;                if mod(iter, 10) == 1            iter            format long            x_min = x_min(:);            format short            save(outname)        end                iter = iter + 1;        Gopt = Gmax(iter,level);        diff_Gmax3 = G_max - Gopt;        Gmax_5iters = Gmax(max([1,iter-4]):iter, :);        Gmax_iter = Gmax(iter, :)        save(outname)    end        t = cputime - t;        ws;    x_min = x_min(:);    format short    level;    %     compute the synthesis filters    g0_coeff = h1_coeff;    for i = 1:length(g0_coeff(:,1))        j = 2-mod(i,2):2:length(g0_coeff(1,:));        g0_coeff(i,j) = -g0_coeff(i,j);    end    g0_coeff = -g0_coeff;        g1_coeff = h0_coeff;    for i = 1:length(g1_coeff(:,1))        j = 2-mod(i,2):2:length(g1_coeff(1,:));        g1_coeff(i,j) = -g1_coeff(i,j);    end    g1_coeff = -g1_coeff;end% check the vanishing momentsN_primal = Vanishing_num(double(subs(g1_coeff, x, x_min)))if N_primal == 0    sum(sum(g1_coeff));endN_dual = Vanishing_num(double(subs(h1_coeff, x, x_min)))if N_dual == 0    sum(sum(h1_coeff));endformat longp_min = reshape([fliplr(x_min(1:px*py)') x_min(1:px*py)'], [px 2*py])u_min = reshape([fliplr(x_min(px*py+1:px*py+ux*uy)') x_min(px*py+1:px*py+ux*uy)'], [ux 2*uy])format short[h0_coeff, h1_coeff] = analysis_filters2(A1c, A2c, x, x_min, Mquin);Gmax_sep = CodingGain_num(h0_coeff, h1_coeff, 0.95, 6, 1, Mquin)Gmax_iso = CodingGain_num(h0_coeff, h1_coeff, 0.95, 6, 0, Mquin)[size_h0, td] = reduce_size(h0_coeff, [0, 0], 1e-10);[size_h1, td] = reduce_size(h1_coeff, [0, 0], 1e-10);H0 = freqz2(h0_coeff, 1024);DC_gain = H0(513,513)H1 = freqz2(h1_coeff, 1024);Nyquist_gain = H1(1,1)clear H0 H1TotalIter = iter - 1save(outname)