l4n.m

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

function [A1_min, A2_min, A3_min, A4_min] = L4n(s_size, levelc, model, ini_point, xx0, sf, Np, Nd, dm, ws, wp, outname)
% Design of quincunx filter banks with four lifting steps with SOCP algorithm
% input: 
%    s_size -- support sizes of the lifing filters
%    levelc -- levels of decomposition
%     model -- isotropic (0) or separable model (1)
% ini_point -- Kovacevic's filter bank (0) or random initial point (1) or fixed point xx0 (2)
%        sf -- constraints on frequency response
%        Np -- number of primal vanishing moments (g1_coeff)
%        Nd -- number of dual vanishing moments (h1_coeff)
%        dm -- 0 - regular sampling and 1 - twin-dragon-type
%        ws -- weight for the stopband
%        wp -- weight for the passband
%   outname -- name of the output workspace
% output: filter coefficients of the four lifting filters
% for example: [A1, A2, A3, A4] = L4n(2*ones(1,8), 4, 0, 0, [], 2, 2, 2, 0, 1, 0, 'data/test')
% Copyright (c) 2006 Yi Chen

format short

if dm == 0
    Mquin = [1 1; 1 -1]
elseif dm == 1
    Mquin = [1 -1; 1 1]
end

delta = 0.00001;
delta_g = [0.001, 0.0001, 0.00001];
delta_h1 = 1.5;
beta = 0.001;

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
syms a51 a52 a53 a54 a55 a56 a57 a58 a59 a60 a61 a62 a63 a64 a65 a66 a67 a68 a69 a70 a71 a72

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 a61 a62 a63 a64 a65 a66 a67 a68 a69 a70 a71 a72];

x1 = s_size(1);
y1 = s_size(2)/2;
x2 = s_size(3);
y2 = s_size(4)/2;
x3 = s_size(5);
y3 = s_size(6)/2;
x4 = s_size(7);
y4 = s_size(8)/2;
a1 = xa(1:x1*y1);
a2 = xa(x1*y1+1:x1*y1+x2*y2);
a3 = xa(x1*y1+x2*y2+1:x1*y1+x2*y2+x3*y3);
a4 = xa(x1*y1+x2*y2+x3*y3+1:x1*y1+x2*y2+x3*y3+x4*y4);  

A1c = reshape([fliplr(a1) a1], [x1 2*y1]);
A2c = reshape([fliplr(a2) a2], [x2 2*y2]);
A3c = reshape([fliplr(a3) a3], [x3 2*y3]);
A4c = reshape([fliplr(a4) a4], [x4 2*y4]);

x = [a1 a2 a3 a4]; 
x = x(:);

% initial point
if ini_point == 2
    xx0 = xx0(:);
    a1n(:,1) = xx0(1:x1*y1);
    a2n(:,1) = xx0(x1*y1+1:x1*y1+x2*y2);
    a3n(:,1) = xx0(x1*y1+x2*y2+1:x1*y1+x2*y2+x3*y3);
    a4n(:,1) = xx0(x1*y1+x2*y2+x3*y3+1:x1*y1+x2*y2+x3*y3+x4*y4);   
elseif ini_point == 0
%     Kovacevic's filter bank    
    if length(a1) == 2
        a1n(:, 1) = [-1 -1]'/4;
    elseif length(a1) == 8
        a1n(:, 1) = [2 -20 -20 2 0 2 2 0]'/64;
    elseif length(a1) == 18
        a1n(:, 1) = [-3 27 -174 -174 27 -3 0 -2 27 27 -2 0 0 0 -3 -3 0 0]'/2^9;
    else
        a1n(:, 1) = zeros(length(a1), 1);
        ini_point = 1;
    end    
    if length(a2) == 2
        a2n(:,1) = [1 1]'/8;
    elseif length(a2) == 8
        a2n(:,1) = [-1 10 10 -1 0 -1 -1 0]'/64;
    elseif length(a2) == 18
        a2n(:, 1) = -[-3 27 -174 -174 27 -3 0 -2 27 27 -2 0 0 0 -3 -3 0 0]'/2^10;
    else
        a2n(:, 1) = zeros(length(a2), 1);
        ini_point = 1;
    end
    sum3 = -(1+2*sum(a1n))/(4*sum(a2n)*(1+2*sum(a1n))+2);
    sum4 = -(1+4*sum(a1n)*sum(a2n)-2*sum(a2n))/(4*sum(a1n)+4*sum3+16*sum(a1n)*sum(a2n)*sum3-2-8*sum(a2n)*sum3);
    if length(a3) == 8
        a3n(:,1) = ones(length(a3),1)*sum3/6;
        a3n(5) = 0; a3n(8) = 0;
    elseif length(a3) == 18
        a3n(:,1) = ones(length(a3),1)*sum3/12;
        a3n(7) = 0; a3n(12) = 0; a3n(13) = 0; a3n(14) = 0; a3n(17) = 0; a3n(18) = 0;     
    else
        a3n(:,1) = ones(length(a3),1)*sum3/length(a3);
    end
    if length(a4) == 8
        a4n(:,1) = ones(8,1)*sum4/6;
        a4n(5) = 0; a4n(8) = 0;
    elseif length(a4) == 18
        a4n(:,1) = ones(18,1)*sum4/12;
        a4n(7) = 0; a4n(12) = 0; a4n(13) = 0; a4n(14) = 0; a4n(17) = 0; a4n(18) = 0;     
    else
        a4n(:,1) = ones(length(a4),1)*sum4/length(a4);
    end
end
if ini_point == 1
%     random initial point
    [Aeq0, beq0] = linear_const34(A1c, A2c, Np, Nd);
    rand('state',sum(100*clock))
    % if the first two lifting steps are not able to provide Np and Nd VMs
    if rank([Aeq0 beq0]) > rank(Aeq0)
        rand('state',sum(100*clock))
        if length(a1) == 8
            a1n(:,1) = -rand(8,1)/8;
            a1n(5) = 0; a1n(8) = 0;
        elseif length(a1) == 18
            a1n(:,1) = -rand(18,1)/12;
            a1n(7) = 0; a1n(12) = 0; a1n(13) = 0; a1n(14) = 0; a1n(17) = 0; a1n(18) = 0; 
        else
            a1n(:,1) = rand(length(a1),1)/length(a1);
        end
        if length(a2) == 8
            a2n(:,1) = -rand(8,1)/16;
            a2n(5) = 0; a2n(8) = 0;
        elseif length(a2) == 18
            a2n(:,1) = -rand(18,1)/12;
            a2n(7) = 0; a2n(12) = 0; a2n(13) = 0; a2n(14) = 0; a2n(17) = 0; a2n(18) = 0;     
        else
            a2n(:,1) = rand(length(a2),1)/length(a2);
        end
        sum3 = -(1+2*sum(a1n))/(4*sum(a2n)*(1+2*sum(a1n))+2);
        if length(a3) == 8
            a3n(:,1) = rand(length(a3),1);
            a3n(:,1) = a3n(:,1)*sum3/(sum(a3n(1:4,1))+sum(a3n(6:7,1)));
            a3n(5) = 0; a3n(8) = 0;
        elseif length(a3) == 18
            a3n(:,1) = rand(length(a3),1);
            a3n(:,1) = a3n(:,1)*sum3/(sum(a3n(1:6,1))+sum(a3n(8:11,1))+sum(a3n(15:16,1)));
            a3n(7) = 0; a3n(12) = 0; a3n(13) = 0; a3n(14) = 0; a3n(17) = 0; a3n(18) = 0;     
        else
            a3n(:,1) = rand(length(a3),1);
            a3n(:,1) = a3n(:,1)*sum3/sum(a3n(:,1));
        end
        sum4 = -(1+4*sum(a1n)*sum(a2n)-2*sum(a2n))/(4*sum(a1n)+4*sum3+16*sum(a1n)*sum(a2n)*sum3-2-8*sum(a2n)*sum3);
        if length(a4) == 8
            a4n(:,1) = rand(length(a4),1);
            a4n(:,1) = a4n(:,1)*sum4/(sum(a4n(1:4,1))+sum(a4n(6:7,1)));
            a4n(5) = 0; a4n(8) = 0;
        elseif length(a4) == 18
            a4n(:,1) = rand(length(a4),1);
            a4n(:,1) = a4n(:,1)*sum4/(sum(a4n(1:6,1))+sum(a4n(8:11,1))+sum(a4n(15:16,1)));
            a4n(7) = 0; a4n(12) = 0; a4n(13) = 0; a4n(14) = 0; a4n(17) = 0; a4n(18) = 0;     
        else
            a4n(:,1) = rand(length(a4),1);
            a4n(:,1) = a4n(:,1)*sum4/sum(a4n(:,1));
        end
    else
        r = rank(Aeq0);
        [U0,S0,V0] = svd(Aeq0);
        % last n-r columns of V
        Vr0 = V0(:, r+1:length(V0(1,:)));
        xs0 = pinv(Aeq0)*beq0;
        phi0 = rand(length(Vr0(1,:)), 1)-0.5;
        x0 = xs0 + Vr0*phi0;
        a1n(:, 1) = x0(1:x1*y1);
        a2n(:, 1) = x0(x1*y1+1:x1*y1+x2*y2);
        a3n(:, 1) = zeros(length(a3),1);
        a4n(:, 1) = zeros(length(a4),1);
    end
end

%     linear constraints on the support lifting filter
Aeq = zeros(1, x1*y1 + x2*y2 + x3*y3 + x4*y4);
beq = 0;
line = 1;
if length(a1) == 8
    Aeq(line, 5) = 1; Aeq(line+1, 8) = 1; 
    line = line + 2;
elseif length(a1) == 18
    Aeq(line, 7) = 1;
    Aeq(line+1, 12) = 1;
    Aeq(line+2, 13) = 1;
    Aeq(line+3, 14) = 1;
    Aeq(line+4, 17) = 1;
    Aeq(line+5, 18) = 1;
    line = line + 6;
end
if length(a2) == 8
    Aeq(line, x1*y1 + 5) = 1; 
    Aeq(line+1, x1*y1 + 8) = 1;
    line = line + 2;
elseif length(a2) == 18
    Aeq(line, x1*y1 + 7) = 1;
    Aeq(line+1, x1*y1 + 12) = 1;
    Aeq(line+2, x1*y1 + 13) = 1;
    Aeq(line+3, x1*y1 + 14) = 1;
    Aeq(line+4, x1*y1 + 17) = 1;
    Aeq(line+5, x1*y1 + 18) = 1;
    line = line + 6;    
end
if length(a3) == 8
    Aeq(line, x1*y1 + x2*y2 + 5) = 1; 
    Aeq(line+1, x1*y1 + x2*y2 + 8) = 1; 
    line = line + 2;
elseif length(a3) == 18
    Aeq(line, x1*y1 + x2*y2 + 7) = 1;
    Aeq(line+1, x1*y1 + x2*y2 + 12) = 1;
    Aeq(line+2, x1*y1 + x2*y2 + 13) = 1;
    Aeq(line+3, x1*y1 + x2*y2 + 14) = 1;
    Aeq(line+4, x1*y1 + x2*y2 + 17) = 1;
    Aeq(line+5, x1*y1 + x2*y2 + 18) = 1;
    line = line + 6;    
end
if length(a4) == 8
    Aeq(line, x1*y1 + x2*y2 + x3*y3  + 5) = 1; 
    Aeq(line+1, x1*y1 + x2*y2 + x3*y3  + 8) = 1; 
    line = line + 2;
elseif length(a4) == 18
    Aeq(line, x1*y1 + x2*y2 + x3*y3 + 7) = 1;
    Aeq(line+1, x1*y1 + x2*y2 + x3*y3 + 12) = 1;
    Aeq(line+2, x1*y1 + x2*y2 + x3*y3 + 13) = 1;
    Aeq(line+3, x1*y1 + x2*y2 + x3*y3 + 14) = 1;
    Aeq(line+4, x1*y1 + x2*y2 + x3*y3 + 17) = 1;
    Aeq(line+5, x1*y1 + x2*y2 + x3*y3 + 18) = 1;
    line = line + 6;     
end
if line > 1
    beq = zeros(line-1, 1);
end

syms w0 w1

N = 64;
M1 = 15;
M2 = 15;


x_min = [a1n(:,1); a2n(:,1); a3n(:,1); a4n(:,1)];

r = rank(Aeq);
[Ua,Sa,Va] = svd(Aeq);
Vr = Va(:, r+1:length(Va(1,:)));
xs = pinv(Aeq)*beq;
xs = xs(:);

I_hat = [eye(x1*y1+x2*y2+x3*y3) zeros(x1*y1+x2*y2+x3*y3, x4*y4)];
xs_hat = I_hat*xs;
Vr_hat = I_hat*Vr;

% adjust the initial point to have enough numbers of VMs    
[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);

if (N_primal < Np) || (N_dual < Nd)