l4n.m

MATLAB Code for Optimal Quincunx Filter Bank Design Yi Chen July 17, 2006 This file introduces t

    x_temp = [a1n(:,1); a2n(:,1); zeros(x3*y3, 1); zeros(x4*y4, 1)];
    [h0_coeff, h1_coeff] = analysis_filters4(A1c, A2c, A3c, A4c, x, x_temp, Mquin);
    [x10, x20, Np, Nd] = VM_two_steps(a3, a4, A3c, A4c, x_min, h0_coeff, h1_coeff, Np, Nd, a3n(:,1), a4n(:,1));
    x_min = [a1n(:,1); a2n(:,1); x10; x20];
    % examine the adjustment results
    [h0_coeff, h1_coeff] = analysis_filters4(A1c, A2c, A3c, A4c, x, x_min, Mquin);
    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
    N_primal = Vanishing_num(g1_coeff);
    N_dual = Vanishing_num(h1_coeff);
end

Initial_point = x_min;

xx(:,1) = x_min;
phi(:, 1) = Vr'*(xx(:, 1) - xs);  

I1 = zeros(x1*y1+x2*y2+x3*y3+x4*y4, 1);
I1(1:x1*y1) = ones(x1*y1, 1);
I2 = zeros(x1*y1+x2*y2+x3*y3+x4*y4, 1);
I2(1+x1*y1:x1*y1+x2*y2) = ones(x2*y2, 1);
I3 = zeros(x1*y1+x2*y2+x3*y3+x4*y4, 1);
I3(1+x1*y1+x2*y2:x1*y1+x2*y2+x3*y3) = ones(x3*y3, 1);
I4 = zeros(x1*y1+x2*y2+x3*y3+x4*y4, 1);
I4(1+x1*y1+x2*y2+x3*y3:x1*y1+x2*y2+x3*y3+x4*y4) = ones(x4*y4, 1);

iter = 1;
iter_t = 1;
t = cputime;

f_min = 100;
G_max = 100;
diff_Gmax = 1;
diff_Gmax2 = 1;
diff_Gmax3 = 0;

if ini_point == 0
    ini_step = 2;
else 
    ini_step = 1;
end

for ini_iter = ini_step:3
    delta = delta_g(3);
    beta = 50*norm(10^(-ini_iter-2)*ones(length(Vr(1,:)),1),2);
    
    level = levelc;

%     analysis filter coefficients
    [h0_coeff, h1_coeff] = analysis_filters4(A1c, A2c, A3c, A4c, x, xx(:,iter), Mquin);
    
    Gmax(iter,:) = CodingGain_num(h0_coeff, h1_coeff, 0.95, levelc, model, Mquin);
    Gt(iter) = -Gmax(iter,level);
    
    if Gt(iter) < f_min
        f_min = Gt(iter);
        x_min = xx(:, iter);
        G_max = -Gt(iter);
    end    
    
    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;
        save(outname)
        
        weight = weight_sep(N, M1, M2, wp, ws);
        [H1w,fx,fy] = freqz2(h1_coeff, N);
        H_ideal = H1w(1,1)*(1-diamond(N+1,N+1));        
%         the sum below should be equal to diff_H1(iter)
        diff_H1_1 = sum(sum((real(H1w) - H_ideal(1:N,1:N)).^2.*weight(1:N,1:N)))*((2*pi/N)^2);
        
        size1 = size(h1_coeff); 
        i = -(size1(1)-1)/2:(size1(1)-1)/2;
        j = -(size1(2)-1)/2:(size1(2)-1)/2;
        [j,i] = meshgrid(j,i);
        D = sum(sum(h1_coeff.*((-1).^(i+j))));
        H1w(1,1);
        [Hd_hat, Sd_hat, Cd_hat]= L4n_parameters(N, M1, M2, ws, wp, Vr_hat, x1, y1, x2, y2, x3, y3, xs_hat, phi(:, iter), H1w(1,1));
        diff_H1(iter) = Sd_hat'*Sd_hat + Cd_hat;
        
        % linear constraint on delta_phi
        [A_phi, b_phi] = VM_4lifting_2(Np, Nd, x1, y1, x2, y2, x3, y3, x4, y4, phi(:,iter), xs, Vr);
        
        lengthp = length(Vr(1,:));
        
        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_filters4(A1c, A2c, A3c, A4c, 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_filters4(A1c, A2c, A3c, A4c, 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;
            g1(i) = (G_temp - Gt(iter))/delta;
        end
        % gradient around xx0
        g = g(:);
        g1 = g1(:);       

        % variable bound for the error function of H1
        delta_h1 = sf*diff_H1(iter);
        
        % SOCP
        b = g;
        
        A1 = Hd_hat';
        c1 = Sd_hat;
        b1 = zeros(lengthp,1);
        d1 = sqrt(delta_h1 - Cd_hat);
        
        A2 = eye(lengthp);
        c2 = zeros(lengthp,1);
        b2 = zeros(lengthp,1);
        d2 = beta;
        
        A3 = A_phi';
        c3 = b_phi;
        b3 = zeros(lengthp,1);
        d3 = 1e-5;
        
%         with constraints on frequency response
        At1 = -[b1 A1];
        At2 = -[b2 A2];
        At3 = -[b3 A3];
        At = [At1 At2 At3];
        bt = -b;
        ct1 = [d1; c1];
        ct2 = [d2; c2];
        ct3 = [d3; c3];
        ct = [ct1; ct2; ct3];
        K.q = [size(At1,2) size(At2,2) size(At3,2)];
              
        pars.fid = 0;
        [xsp, delta_phi, info] = sedumi(At, bt, ct, K, pars);
        info;
        
%         % insert a line-search step here
%         lambda = line_search4_2(1e-4, level, model, phi(:, iter), delta_phi, A1c, A2c, A3c, A4c, x, Vr, xs);
%         delta_phi = lambda*delta_phi;
        
        phi(:, iter+1) = phi(:, iter) + delta_phi;
        xx(:, iter+1) = xs + Vr*phi(:, iter+1);
        
        min = g'*delta_phi;
        
        [h0_coeff, h1_coeff] = analysis_filters4(A1c, A2c, A3c, A4c, x, xx(:, iter+1), Mquin);
        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
        N_primal = Vanishing_num(g1_coeff);
        N_dual = Vanishing_num(h1_coeff);
        if mod(iter,10)==1
            iter
            format long
            x_min = x_min(:);
            format short
        end
        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) = (Hd_hat*delta_phi+Sd_hat)'*(Hd_hat*delta_phi+Sd_hat) + Cd_hat;
        diff_H1_3 = (Hd_hat*delta_phi+Sd_hat)'*(Hd_hat*delta_phi+Sd_hat) + Cd_hat;
        
        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;
        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
    
    iter_t = iter;
    
    t = cputime - t;
    
    % adjust the result    
    xx(:,iter) = x_min;
    phi(:, iter) = Vr'*(xx(:, iter) - xs);
    
    [A_phi, b_phi] = VM_4lifting_2(Np, Nd, x1, y1, x2, y2, x3, y3, x4, y4, phi(:,iter), xs, Vr);
    
    Hphi = zeros(length(phi(:,1)), length(phi(:,1))); Sphi = zeros(length(phi(:,1)), 1);
    for i = 1:length(b_phi)
        Hphi = A_phi(i,:)'*A_phi(i,:) + Hphi;
        Sphi = 2*b_phi(i)*A_phi(i,:)' + Sphi;
    end
    phi_min = -0.5*pinv(Hphi)*Sphi;
    
    phi(:, iter) = phi(:, iter) + phi_min;
    xx(:,iter) = Vr*phi(:,iter) + xs;
    [h0_coeff, h1_coeff] = analysis_filters4(A1c, A2c, A3c, A4c, x, xx(:, iter), Mquin);
    Gmax(iter,:) = CodingGain_num(h0_coeff, h1_coeff, 0.95, level, model, Mquin);
    Gt(iter) = -Gmax(iter,level);
    f_min = Gt(iter);
    x_min = xx(:, iter);
    G_max = -Gt(iter);

    x_min = x_min(:);
    
    if ini_iter == 3        
        format long
        A1_min = reshape([fliplr(x_min(1:x1*y1)') x_min(1:x1*y1)'], [x1 2*y1])
        A2_min = reshape([fliplr(x_min(x1*y1+1:x1*y1+x2*y2)') x_min(x1*y1+1:x1*y1+x2*y2)'], [x2 2*y2])
        A3_min = reshape([fliplr(x_min(x1*y1+x2*y2+1:x1*y1+x2*y2+x3*y3)') x_min(x1*y1+x2*y2+1:x1*y1+x2*y2+x3*y3)'], [x3 2*y3])
        A4_min = reshape([fliplr(x_min(x1*y1+x2*y2+x3*y3+1:x1*y1+x2*y2+x3*y3+x4*y4)') x_min(x1*y1+x2*y2+x3*y3+1:x1*y1+x2*y2+x3*y3+x4*y4)'], [x4 2*y4])
        format short
        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)
        
        % check the vanishing moments
        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;
        N_primal = Vanishing_num(g1_coeff)
        N_dual = Vanishing_num(h1_coeff)
    end    
    level;
end

H0 = freqz2(h0_coeff, 1024);
DC_gain = H0(513,513)
H1 = freqz2(h1_coeff, 1024);
Nyquist_gain = H1(1,1)

clear H0 H1
TotalIter = iter - 1
save(outname)