l3n2.m

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

[h0_coeff, h1_coeff] = analysis_filters3(A1c, A2c, A3c, 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);

ini_step = 1;

for ini_iter = ini_step:3
    delta = delta_g(3);
    beta = 100*norm(10^(-ini_iter-2)*ones(length(Vr(1,:)),1),2);
    
    levelc = level;    
    xx(:,iter) = x_min;
    phi(:, iter) = Vr'*(xx(:, iter) - xs);  
    
%     analysis filter coefficients
    [h0_coeff, h1_coeff] = analysis_filters3(A1c, A2c, A3c, 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;
        
        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);
        
%         compute the parameters          
        [Hd_hat, Sd_hat, Cd_hat]= L3n2_parameters(N, M1, M2, ws, wp, Vr_hat, x1, y1, x2, y2, xs_hat, phi(:, iter), H1w(1,1));
        diff_H1(1) = Sd_hat'*Sd_hat + Cd_hat;
        diff_H1_2 = Sd_hat'*Sd_hat + Cd_hat;
        
        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_filters3(A1c, A2c, A3c, 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_filters3(A1c, A2c, A3c, 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);
        
        % linear constraint on delta_phi
        A_phi = []; b_phi = [];
        A_phi(1,:) = phi(:, iter)'*(H1+H1')+S1';
        A_phi(2,:) = phi(:, iter)'*(H2+H2')+S2';
        
        b_phi(1) = phi(:, iter)'*H1*phi(:, iter) + phi(:, iter)'*S1 + C1;
        b_phi(2) = phi(:, iter)'*H2*phi(:, iter) + phi(:, iter)'*S2 + C2;
        b_phi = b_phi(:);
        
        [A_more, b_more] = VM_3lifting_2(Np, Nd, x1, y1, x2, y2, x3, y3, phi(:,iter), xs, Vr);
        A_phi = [A_phi; A_more]; 
        b_phi = [b_phi; b_more];
        
        % 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;
        phi(:, iter+1) = phi(:, iter) + delta_phi;
        xx(:, iter+1) = xs + Vr*phi(:, iter+1);
        a1n(:, iter+1) = xx(1:x1*y1, iter+1);
        
        min = g'*delta_phi;
        
        [h0_coeff, h1_coeff] = analysis_filters3(A1c, A2c, A3c, 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 = [];
    A_phi(1,:) = phi(:, iter)'*(H1+H1')+S1';
    A_phi(2,:) = phi(:, iter)'*(H2+H2')+S2';
    
    b_phi(1) = phi(:, iter)'*H1*phi(:, iter) + phi(:, iter)'*S1 + C1;
    b_phi(2) = phi(:, iter)'*H2*phi(:, iter) + phi(:, iter)'*S2 + C2;
    b_phi = b_phi(:);
    
    [A_more, b_more] = VM_3lifting_2(Np, Nd, x1, y1, x2, y2, x3, y3, phi(:,iter), xs, Vr);
    A_phi = [A_phi; A_more];
    b_phi = [b_phi; b_more];
    
    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_filters3(A1c, A2c, A3c, 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(1+x1*y1:x1*y1 + x2*y2)') x_min(1+x1*y1:x1*y1+x2*y2)'], [x2 2*y2])
        A3_min = reshape([fliplr(x_min(1 + x1*y1 + x2*y2:x1*y1 + x2*y2 + x3*y3)') x_min(1 + x1*y1 + x2*y2:x1*y1 + x2*y2 + x3*y3)'], [x3 2*y3])
        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
Total_iterations = iter - 1
save(outname)