vm_two_steps.m

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

function [x10, x20, Np, Nd] = VM_two_steps(a3, a4, A3c, A4c, x_min, h0_coeff, h1_coeff, Np, Nd, a3n, a4n)
% add two lifting steps to increase the number of vanishing moments for
% four-lifting steps
 % Copyright (c) 2006 Yi Chen

x0 = x_min(:)';

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 a29 a30 a31 a32;
syms z0 z1 z00 z01

% transfer functions of h0, h1
sizeh = size(h0_coeff);
n0 = -(sizeh(1)-1)/2:(sizeh(1)-1)/2;
n1 = -(sizeh(2)-1)/2:(sizeh(2)-1)/2;
[n1,n0] = meshgrid(n1,n0);
h0 = sum(sum(h0_coeff.*(z0.^(-n0)).*(z1.^(-n1))));

sizeh = size(h1_coeff);
n0 = -(sizeh(1)-1)/2-1:(sizeh(1)-1)/2-1;
n1 = -(sizeh(2)-1)/2:(sizeh(2)-1)/2;
[n1,n0] = meshgrid(n1,n0);
h1 = sum(sum(h1_coeff.*(z0.^(-n0)).*(z1.^(-n1))));

g1 = subs(h0,{z0,z1},{-z0, -z1});
g0 = -subs(h1,{z0,z1},{-z0, -z1});

x1 = a3;
x2 = a4;
length1 = length(a3);
length2 = length(a4);

f1_coeff = A3c;
f2_coeff = A4c;


sizef1 = size(f1_coeff);
df1 = - [sizef1(1)-2, sizef1(2)-2] / 2;

sizef2 = size(f2_coeff);
df2 = - [sizef2(1), sizef2(2)] / 2;

f1 = matrix2poly(f1_coeff, df1);
f2 = matrix2poly(f2_coeff, df2);

f1_up = subs(subs(f1,{z0,z1},{z00*z01,z00*z01^(-1)}),{z00,z01},{z0,z1});
f2_up = subs(subs(f2,{z0,z1},{z00*z01,z00*z01^(-1)}),{z00,z01},{z0,z1});


% the first lifting step (from channel 0 to channel 1)
h1f = expand(h1 + f1_up*h0);
g0f = expand(g0 - f1_up*g1);

[dh1f, h1f_coeff] = poly2matrix_sym(h1f);

sizeh1f = size(h1f_coeff);

n0 = -(sizeh1f(1)-1)/2:(sizeh1f(1)-1)/2;
n1 = -(sizeh1f(2)-1)/2:(sizeh1f(2)-1)/2;

[n1,n0] = meshgrid(n1,n0);


ind = 1;
for j = 0:2:Nd-1
    for k = 0:j;
        matrix_h = (n0.^k).*(n1.^(j-k));
        hx(ind) = expand(sum(sum(matrix_h.*h1f_coeff)));
        ind = ind + 1;
    end
end

% find the linear constraints on the lifting filter coefficients for Nd
% dual vanishing moments
for ind = 1:length(hx)
    for j = 1:length1
        Aeq1(ind, j) = double(maple('coeff', hx(ind), x1(j)));
    end
    beq1(ind) = -double(subs((hx(ind) - Aeq1(ind, :)*transpose(x1)), x1, zeros(1, length1)));
end

beq1 = beq1(:);

if length(a3) < size(Aeq1, 1)
    disp('not enough filter coefficients')
elseif rank([Aeq1, beq1]) > rank(Aeq1)
    disp('change to two primal vanishing moments')
    Aeq1 = Aeq1(1, :);
    beq1 = beq1(1);
    Np = 2;
end

if length1 == 8
    lb = length(beq1);
    Aeq1(lb+1, :) = zeros(1, 8); Aeq1(lb+1, 5) = 1;
    Aeq1(lb+2, :) = zeros(1, 8); Aeq1(lb+2, 8) = 1;
    beq1(lb+1) = 0; beq1(lb+2) = 0;
end

beq1 = beq1(:);
x10 = pinv(Aeq1)*beq1;

hx1 = double(subs(hx, x1, x10));

g0f = subs(g0f, x1, x10);
[dg0f, g0f_coeff] = poly2matrix(g0f);
h1f = subs(h1f, x1, x10);
f1 = subs(f1, x1, x10);
g0f_coeff = double(subs(g0f_coeff, x1, x10));
h1f_coeff = double(subs(h1f_coeff, x1, x10));

% the second lifting step (from channel 0 to channel 1)
h0f = expand(h0 + f2_up*h1f);
g1f = expand(g1 - f2_up*g0f);
[dh0f, h0f_coeff] = poly2matrix_sym(h0f);
[dg1f, g1f_coeff] = poly2matrix_sym(g1f);

sizeg1f = size(g1f_coeff);
n0 = [0:-1:(1-sizeg1f(1))] - dg1f(1);
n1 = [0:-1:(1-sizeg1f(2))] - dg1f(2);

[n1,n0] = meshgrid(n1,n0);


% compute the linear equations to be solved for the Np primal vanishing moments
% the number of equations should be N(N-1)/2
% double(sum(sum(matrix_g.*g1f_coeff))) is used, this part can be used to
% verify the correctness of the optimization results with a certain nubmer
% of vanishing moments (i.e. gx(ind) should be all zero)
ind = 1;
for j = 0:2:Np-1
    for k = 0:j;
        matrix_g = (n0.^k).*(n1.^(j-k));
        gx(ind) = expand(sum(sum(matrix_g.*g1f_coeff)));
        ind = ind + 1;
    end
end

% find the linear constraints on the lifting filter coefficients for Np
% primal vanishing moments
for ind = 1:length(gx)
    for j = 1:length2
        Aeq2(ind, j) = double(maple('coeff', gx(ind), x2(j)));
    end
    beq2(ind) = -double(subs((gx(ind) - Aeq2(ind, :)*transpose(x2)), x2, zeros(1, length2)));
end

if length(a4) < size(Aeq2, 1)
    disp('not enough filter coefficients')
elseif rank([Aeq2, beq2(:)]) > rank(Aeq2)
    disp('change to two dual vanishing moments')
    Aeq2 = Aeq2(1, :);
    beq2 = beq2(1);
    Nd = 2;
end

if length2 == 8
    lb = length(beq2);
    Aeq2(lb+1, :) = zeros(1, 8); Aeq2(lb+1, 5) = 1;
    Aeq2(lb+2, :) = zeros(1, 8); Aeq2(lb+2, 8) = 1;
    beq2(lb+1) = 0; beq2(lb+2) = 0;
end

beq2 = beq2(:);
x20 = pinv(Aeq2)*beq2(:);

gx2 = double(subs(gx, x2, x20));