l2n.m

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

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 Chen

format short
syms 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
syms 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

xa = [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]
end

syms w0 w1

for 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);
    end
end
v_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')
    return
end

I_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))';
end
if 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')';
end

InitialPoint = xx(1,:);
x_min = xx(1, :)';
t = cputime;
iter = 1;

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

for 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 moments
N_primal = Vanishing_num(double(subs(g1_coeff, x, x_min)))
if N_primal == 0
    sum(sum(g1_coeff));
end
N_dual = Vanishing_num(double(subs(h1_coeff, x, x_min)))
if N_dual == 0
    sum(sum(h1_coeff));
end

format long
p_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 H1
TotalIter = iter - 1
save(outname)