l3n2.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] = L3n2(s_size, level, model, ini_point, x00, sf, Np, Nd, dm, ws, wp, outname)
% Optimize the coding gain subject to the error of h1 is less than dh1
% with second-order cone programming
% for a three-step lifting structure
% the first lifting step adds channel 1 to channel 0
% input: s_size -- support sizes of the lifing filters
%         level -- level of decomposition
%         model -- isotropic (0) or separable (1) model
%     ini_point -- 0: some good one; 1: random; and 2: a fixed point x00
%           x00 -- initial point
%            sf -- constraints on frequency response
%            Np -- number of primal vanishing moments (g1_coeff)
%            Nu -- number of dual vanishing moments (h1_coeff)
%            dm -- type of downsampling matrix, 0 - regular; 1 - twin dragon 
%            ws --  weight for the stopband
%            wp --  weight for the passband
%       outname -- name of the workspace
% output: lifting filter coefficients
% for example
% [A1_min, A2_min, A3_min] = L3n2([4 4 4 4 4 4], 1, 0, 0, [], 2, 4, 4, 0, 1, 0, 'data/test')
% [A1_min, A2_min, A3_min] = L3n2([2 2 2 2 2 2], 1, 0, 1, [], 1, 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 a51 a52 a53 a54
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];

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;

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);

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

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

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

% initial point
if ini_point == 0
%     some "good" initial point
    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
    
    sum2 = -1/(2+4*sum(a1n(:,1)));
    sum3 = -(2*sum(a1n(:,1))-1)*(2+4*sum(a1n(:,1)))/8;
    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
    x_min = [a1n(:,1); a2n(:,1); a3n(:,1)];
end
if ini_point == 1
%     random initial point
    [Aeq0, beq0] = linear_const34(A1c, A2c, Np, Nd);
    beq0 = -beq0;
    rand('state',sum(100*clock))
    if rank([Aeq0 beq0]) > rank(Aeq0)
        if length(a1) == 8
            a1n(:,1) = rand(8,1)/6;
            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
        
        sum2 = -1/(2+4*sum(a1n(:,1)));
        if length(a2) == 8
            a2n(:,1) = rand(length(a2),1);
            a2n(:,1) = a2n(:,1)*sum2/(sum(a2n(1:4,1))+sum(a2n(6:7,1)));
            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);
            a2n(:,1) = a2n(:,1)*sum2/sum(a2n(:,1));
        end
        
        sum3 = -(2*sum(a1n(:,1))-1)*(2+4*sum(a1n(:,1)))/8;
        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
    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);
    end
    x_min = [a1n(:,1); a2n(:,1); a3n(:,1)];
    
elseif ini_point == 2
    x00 = x00(:);
    a1n(:,1) = x00(1:x1*y1);
    a2n(:,1) = x00(1+x1*y1:x1*y1+x2*y2);
    a3n(:,1) = x00(1+x1*y1+x2*y2:x1*y1+x2*y2+x3*y3);
    x_min = x00;
end

line = 1;
Aeq = zeros(1, x1*y1 + x2*y2 + x3*y3);
beq = 0;

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 line > 1
    beq = zeros(line-1, 1);
end

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

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

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

H1 = 4*Vr'*I2*I1'*Vr;
S1 = 4*Vr'*(I2*I1'+I1*I2')*xs + Vr'*2*I2;
C1 = 4*xs'*I2*I1'*xs + xs'*2*I2 + 1;

H2 = 8*Vr'*I2*I3'*Vr;
S2 = 8*Vr'*(I2*I3'+I3*I2')*xs - Vr'*2*I1;
C2 = 8*xs'*I2*I3'*xs - xs'*2*I1 + 1;

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

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

Initial_point = x_min;
xx(:,1) = x_min;