function [A, b] = VM_4lifting_2(Np, Nd, x1, y1, x2, y2, x3, y3, x4, y4, phik, xs, Vr)
% compute the approximated linear system for the desired numbers of
% vanishing moments for the optimization problem L4n
% Np primal and Nd dual vanishing moments
% xi, yi -- size of the ith lifting filter
% Copyright (c) 2006 Yi Chen

[Ny, Nx] = meshgrid(1:y1, 1:x1);
n10 = reshape(Nx+Ny-x1/2-1, [x1*y1,1]);
n11 = reshape(Nx-Ny-x1/2, [x1*y1,1]);

[Ny, Nx] = meshgrid(1:y2, 1:x2);
n20 = reshape(Nx+Ny-x2/2-1, [x2*y2,1]);
n21 = reshape(Nx-Ny-x2/2, [x2*y2,1]);

[Ny, Nx] = meshgrid(1:y3, 1:x3);
n30 = reshape(Nx+Ny-x3/2-1, [x3*y3,1]);
n31 = reshape(Nx-Ny-x3/2, [x3*y3,1]);

[Ny, Nx] = meshgrid(1:y4, 1:x4);
n40 = reshape(Nx+Ny-x4/2-1, [x4*y4,1]);
n41 = reshape(Nx-Ny-x4/2, [x4*y4,1]);

i1 = ones(x1*y1,1);
i2 = ones(x2*y2,1);
i3 = ones(x3*y3,1);
i4 = ones(x4*y4,1);

I1 = [eye(x1*y1) zeros(x1*y1, x2*y2+x3*y3+x4*y4)];
I2 = [zeros(x2*y2, x1*y1) eye(x2*y2) zeros(x2*y2, x3*y3+x4*y4)];
I3 = [zeros(x3*y3, x1*y1+x2*y2) eye(x3*y3) zeros(x3*y3, x4*y4)];
I4 = [zeros(x4*y4, x1*y1+x2*y2+x3*y3) eye(x4*y4)];

a1_hat = I1*(xs+Vr*phik);
a2_hat = I2*(xs+Vr*phik);
a3_hat = I3*(xs+Vr*phik);
a4_hat = I4*(xs+Vr*phik);

V1_hat = I1*Vr;
V2_hat = I2*Vr;
V3_hat = I3*Vr;
V4_hat = I4*Vr;

b(1) = 1 + 4*a1_hat'*i1*a2_hat'*i2 + 4*a1_hat'*i1*a4_hat'*i4 + 4*a3_hat'*i3*a4_hat'*i4...
    + 16*a1_hat'*i1*a2_hat'*i2*a3_hat'*i3*a4_hat'*i4 - 2*a2_hat'*i2 - 2*a4_hat'*i4 - 8*a2_hat'*i2*a3_hat'*i3*a4_hat'*i4;
A(1, 1:length(phik)) = 4*a1_hat'*i1*i4'*V4_hat + 4*a4_hat'*i4*i1'*V1_hat + 4*a1_hat'*i1*i2'*V2_hat + 4*a2_hat'*i2*i1'*V1_hat + 4*a3_hat'*i3*i4'*V4_hat + 4*a4_hat'*i4*i3'*V3_hat...
    - 2*i2'*V2_hat - 2*i4'*V4_hat - 8*a4_hat'*i4*a2_hat'*i2*i3'*V3_hat - 8*a4_hat'*i4*a3_hat'*i3*i2'*V2_hat - 8*a2_hat'*i2*a3_hat'*i3*i4'*V4_hat...
    + 16*i1'*a1_hat*i2'*a2_hat*i3'*a3_hat*i4'*V4_hat + 16*i2'*a2_hat*i3'*a3_hat*i4'*a4_hat*i1'*V1_hat...
    + 16*i3'*a3_hat*i4'*a4_hat*i1'*a1_hat*i2'*V2_hat + 16*i4'*a4_hat*i1'*a1_hat*i2'*a2_hat*i3'*V3_hat;
b(2) = 1 + 4*a2_hat'*i2*a3_hat'*i3 + 2*a1_hat'*i1 + 2*a3_hat'*i3 + 8*a1_hat'*i1*a2_hat'*i2*a3_hat'*i3;
A(2, 1:length(phik)) = 2*i1'*V1_hat + 2*i3'*V3_hat + 8*a1_hat'*i1*a2_hat'*i2*i3'*V3_hat + 8*a1_hat'*i1*a3_hat'*i3*i2'*V2_hat + 8*a2_hat'*i2*a3_hat'*i3*i1'*V1_hat...
    + 4*a2_hat'*i2*i3'*V3_hat + 4*a3_hat'*i3*i2'*V2_hat; 
b = b(:);
    
if Nd == 4
bd(1) = a1_hat'*(n10.*n11.*2.*i1) + a3_hat'*(n30.*n31.*2.*i3)...
    + a2_hat'*((n20.*n21.*2.*i2)*2*i3'+2*i2*(n30.*n31.*2.*i3)')*a3_hat...
    + a1_hat'*(n10.*n11.*2.*i1)*a2_hat'*2*i2*a3_hat'*2*i3...
    + a1_hat'*2*i1*a2_hat'*(n20.*n21.*2.*i2)*a3_hat'*2*i3...
    + a1_hat'*2*i1*a2_hat'*2*i2*a3_hat'*(n30.*n31.*2.*i3);

Ad(1, 1:length(phik)) = (n10.*n11.*2.*i1)'*V1_hat + (n30.*n31.*2.*i3)'*V3_hat...
    + a3_hat'*((n20.*n21.*2.*i2)*2*i3'+2*i2*(n30.*n31.*2.*i3)')'*V2_hat...
    + a2_hat'*((n20.*n21.*2.*i2)*2*i3'+2*i2*(n30.*n31.*2.*i3)')*V3_hat...
    + a1_hat'*(n10.*n11.*2.*i1)*a2_hat'*2*i2*(2*i3)'*V3_hat...
    + a1_hat'*(n10.*n11.*2.*i1)*a3_hat'*2*i3*(2*i2)'*V2_hat...
    + a2_hat'*2*i2*a3_hat'*2*i3*(n10.*n11.*2.*i1)'*V1_hat...
    + a2_hat'*(n20.*n21.*2.*i2)*a3_hat'*2*i3*(2*i1)'*V1_hat...
    + a2_hat'*(n20.*n21.*2.*i2)*a1_hat'*2*i1*(2*i3)'*V3_hat...
    + a1_hat'*2*i1*a3_hat'*2*i3*(n20.*n21.*2.*i2)'*V2_hat...
    + a3_hat'*(n30.*n31.*2.*i3)*a2_hat'*2*i2*(2*i1)'*V1_hat...
    + a3_hat'*(n30.*n31.*2.*i3)*a1_hat'*2*i1*(2*i2)'*V2_hat...
    + a2_hat'*2*i2*a1_hat'*2*i1*(n30.*n31.*2.*i3)'*V3_hat; 

bd(2) = a1_hat'*(n10.*n10.*2.*i1) + a3_hat'*(n30.*n30.*2.*i3)...
    + a2_hat'*((n20.*n20.*2.*i2)*2*i3'+2*i2*(n30.*n30.*2.*i3)')*a3_hat...
    + a1_hat'*(n10.*n10.*2.*i1)*a2_hat'*2*i2*a3_hat'*2*i3...
    + a1_hat'*2*i1*a2_hat'*(n20.*n20.*2.*i2)*a3_hat'*2*i3...
    + a1_hat'*2*i1*a2_hat'*2*i2*a3_hat'*(n30.*n30.*2.*i3);

Ad(2, 1:length(phik)) = (n10.*n10.*2.*i1)'*V1_hat + (n30.*n30.*2.*i3)'*V3_hat...
    + a3_hat'*((n20.*n20.*2.*i2)*2*i3'+2*i2*(n30.*n30.*2.*i3)')'*V2_hat...
    + a2_hat'*((n20.*n20.*2.*i2)*2*i3'+2*i2*(n30.*n30.*2.*i3)')*V3_hat...
    + a1_hat'*(n10.*n10.*2.*i1)*a2_hat'*2*i2*(2*i3)'*V3_hat...
    + a1_hat'*(n10.*n10.*2.*i1)*a3_hat'*2*i3*(2*i2)'*V2_hat...
    + a2_hat'*2*i2*a3_hat'*2*i3*(n10.*n10.*2.*i1)'*V1_hat...
    + a2_hat'*(n20.*n20.*2.*i2)*a3_hat'*2*i3*(2*i1)'*V1_hat...
    + a2_hat'*(n20.*n20.*2.*i2)*a1_hat'*2*i1*(2*i3)'*V3_hat...
    + a1_hat'*2*i1*a3_hat'*2*i3*(n20.*n20.*2.*i2)'*V2_hat...
    + a3_hat'*(n30.*n30.*2.*i3)*a2_hat'*2*i2*(2*i1)'*V1_hat...
    + a3_hat'*(n30.*n30.*2.*i3)*a1_hat'*2*i1*(2*i2)'*V2_hat...
    + a2_hat'*2*i2*a1_hat'*2*i1*(n30.*n30.*2.*i3)'*V3_hat; 

bd(3) = a1_hat'*(n11.*n11.*2.*i1) + a3_hat'*(n31.*n31.*2.*i3)...
    + a2_hat'*((n21.*n21.*2.*i2)*2*i3'+2*i2*(n31.*n31.*2.*i3)')*a3_hat...
    + a1_hat'*(n11.*n11.*2.*i1)*a2_hat'*2*i2*a3_hat'*2*i3...
    + a1_hat'*2*i1*a2_hat'*(n21.*n21.*2.*i2)*a3_hat'*2*i3...
    + a1_hat'*2*i1*a2_hat'*2*i2*a3_hat'*(n31.*n31.*2.*i3);

Ad(3, 1:length(phik)) = (n11.*n11.*2.*i1)'*V1_hat + (n31.*n31.*2.*i3)'*V3_hat...
    + a3_hat'*((n21.*n21.*2.*i2)*2*i3'+2*i2*(n31.*n31.*2.*i3)')'*V2_hat...
    + a2_hat'*((n21.*n21.*2.*i2)*2*i3'+2*i2*(n31.*n31.*2.*i3)')*V3_hat...
    + a1_hat'*(n11.*n11.*2.*i1)*a2_hat'*2*i2*(2*i3)'*V3_hat...
    + a1_hat'*(n11.*n11.*2.*i1)*a3_hat'*2*i3*(2*i2)'*V2_hat...
    + a2_hat'*2*i2*a3_hat'*2*i3*(n11.*n11.*2.*i1)'*V1_hat...
    + a2_hat'*(n21.*n21.*2.*i2)*a3_hat'*2*i3*(2*i1)'*V1_hat...
    + a2_hat'*(n21.*n21.*2.*i2)*a1_hat'*2*i1*(2*i3)'*V3_hat...
    + a1_hat'*2*i1*a3_hat'*2*i3*(n21.*n21.*2.*i2)'*V2_hat...
    + a3_hat'*(n31.*n31.*2.*i3)*a2_hat'*2*i2*(2*i1)'*V1_hat...
    + a3_hat'*(n31.*n31.*2.*i3)*a1_hat'*2*i1*(2*i2)'*V2_hat...
    + a2_hat'*2*i2*a1_hat'*2*i1*(n31.*n31.*2.*i3)'*V3_hat; 
bd = bd(:);
end

if Np == 4
    m10 = n10; m11 = n11; 
    m20 = n20; m21 = n21;
    m30 = n30; m31 = n31;
    m40 = n40; m41 = n41;
    u1 = m10.*m11.*i1; u2 = m20.*m21.*i2; u3 = m30.*m31.*i3; u4 = m40.*m41.*i4; 
    bp(1) = 4*a1_hat'*(u1*i2'+i1*u2')*a2_hat + 4*a3_hat'*(u3*i4'+i3*u4')*a4_hat...
        + 4*a1_hat'*(u1*i4'+i1*u4')*a4_hat - 2*u2'*a2_hat - 2*u4'*a4_hat...
        - 8*(a4_hat'*u4*a2_hat'*i2*a3_hat'*i3 + a4_hat'*i4*a2_hat'*u2*a3_hat'*i3 + a4_hat'*i4*a2_hat'*i2*a3_hat'*u3)...
        + 16*a1_hat'*u1*a2_hat'*i2*a3_hat'*i3*a4_hat'*i4 + 16*a2_hat'*u2*a1_hat'*i1*a3_hat'*i3*a4_hat'*i4...
        + 16*a3_hat'*u3*a2_hat'*i2*a1_hat'*i1*a4_hat'*i4 + 16*a4_hat'*u4*a2_hat'*i2*a3_hat'*i3*a1_hat'*i1; 
    Ap(1, 1:length(phik)) = 4*a1_hat'*(u1*i2'+i1*u2')*V2_hat + 4*a2_hat'*(u2*i1'+i2*u1')*V1_hat...
        + 4*a1_hat'*(u1*i4'+i1*u4')*V4_hat + 4*a4_hat'*(u4*i1'+i4*u1')*V1_hat...
        + 4*a3_hat'*(u3*i4'+i3*u4')*V4_hat + 4*a4_hat'*(u4*i3'+i4*u3')*V3_hat...
        - 2*u2'*V2_hat - 2*u4'*V4_hat...
        - 8*a4_hat'*u4*a2_hat'*i2*i3'*V3_hat - 8*a4_hat'*u4*a3_hat'*i3*i2'*V2_hat - 8*a2_hat'*i2*a3_hat'*i3*u4'*V4_hat...
        - 8*a2_hat'*u2*a3_hat'*i3*i4'*V4_hat - 8*a2_hat'*u2*a4_hat'*i4*i3'*V3_hat - 8*a4_hat'*i4*a3_hat'*i3*u2'*V2_hat...
        - 8*a3_hat'*u3*a2_hat'*i2*i4'*V4_hat - 8*a3_hat'*u3*a4_hat'*i4*i2'*V2_hat - 8*a2_hat'*i2*a4_hat'*i4*u3'*V3_hat...
        + 16*u1'*a1_hat*i2'*a2_hat*i3'*a3_hat*i4'*V4_hat + 16*i2'*a2_hat*i3'*a3_hat*i4'*a4_hat*u1'*V1_hat...
        + 16*i3'*a3_hat*i4'*a4_hat*u1'*a1_hat*i2'*V2_hat + 16*i4'*a4_hat*u1'*a1_hat*i2'*a2_hat*i3'*V3_hat...        
        + 16*i1'*a1_hat*u2'*a2_hat*i3'*a3_hat*i4'*V4_hat + 16*u2'*a2_hat*i3'*a3_hat*i4'*a4_hat*i1'*V1_hat...
        + 16*i3'*a3_hat*i4'*a4_hat*i1'*a1_hat*u2'*V2_hat + 16*i4'*a4_hat*i1'*a1_hat*u2'*a2_hat*i3'*V3_hat...        
        + 16*i1'*a1_hat*i2'*a2_hat*u3'*a3_hat*i4'*V4_hat + 16*i2'*a2_hat*u3'*a3_hat*i4'*a4_hat*i1'*V1_hat...
        + 16*u3'*a3_hat*i4'*a4_hat*i1'*a1_hat*i2'*V2_hat + 16*i4'*a4_hat*i1'*a1_hat*i2'*a2_hat*u3'*V3_hat...        
        + 16*i1'*a1_hat*i2'*a2_hat*i3'*a3_hat*u4'*V4_hat + 16*i2'*a2_hat*i3'*a3_hat*u4'*a4_hat*i1'*V1_hat...
        + 16*i3'*a3_hat*u4'*a4_hat*i1'*a1_hat*i2'*V2_hat + 16*u4'*a4_hat*i1'*a1_hat*i2'*a2_hat*i3'*V3_hat;  
    
    m10 = n10; m11 = n10; 
    m20 = n20; m21 = n20;
    m30 = n30; m31 = n30;
    m40 = n40; m41 = n40;
    u1 = m10.*m11.*i1; u2 = m20.*m21.*i2; u3 = m30.*m31.*i3; u4 = m40.*m41.*i4; 
    bp(2) = 4*a1_hat'*(u1*i2'+i1*u2')*a2_hat + 4*a3_hat'*(u3*i4'+i3*u4')*a4_hat...
        + 4*a1_hat'*(u1*i4'+i1*u4')*a4_hat - 2*u2'*a2_hat - 2*u4'*a4_hat...
        - 8*(a4_hat'*u4*a2_hat'*i2*a3_hat'*i3 + a4_hat'*i4*a2_hat'*u2*a3_hat'*i3 + a4_hat'*i4*a2_hat'*i2*a3_hat'*u3)...
        + 16*a1_hat'*u1*a2_hat'*i2*a3_hat'*i3*a4_hat'*i4 + 16*a2_hat'*u2*a1_hat'*i1*a3_hat'*i3*a4_hat'*i4...
        + 16*a3_hat'*u3*a2_hat'*i2*a1_hat'*i1*a4_hat'*i4 + 16*a4_hat'*u4*a2_hat'*i2*a3_hat'*i3*a1_hat'*i1; 
    Ap(2, 1:length(phik)) = 4*a1_hat'*(u1*i2'+i1*u2')*V2_hat + 4*a2_hat'*(u2*i1'+i2*u1')*V1_hat...
        + 4*a1_hat'*(u1*i4'+i1*u4')*V4_hat + 4*a4_hat'*(u4*i1'+i4*u1')*V1_hat...
        + 4*a3_hat'*(u3*i4'+i3*u4')*V4_hat + 4*a4_hat'*(u4*i3'+i4*u3')*V3_hat...
        - 2*u2'*V2_hat - 2*u4'*V4_hat...
        - 8*a4_hat'*u4*a2_hat'*i2*i3'*V3_hat - 8*a4_hat'*u4*a3_hat'*i3*i2'*V2_hat - 8*a2_hat'*i2*a3_hat'*i3*u4'*V4_hat...
        - 8*a2_hat'*u2*a3_hat'*i3*i4'*V4_hat - 8*a2_hat'*u2*a4_hat'*i4*i3'*V3_hat - 8*a4_hat'*i4*a3_hat'*i3*u2'*V2_hat...
        - 8*a3_hat'*u3*a2_hat'*i2*i4'*V4_hat - 8*a3_hat'*u3*a4_hat'*i4*i2'*V2_hat - 8*a2_hat'*i2*a4_hat'*i4*u3'*V3_hat...
        + 16*u1'*a1_hat*i2'*a2_hat*i3'*a3_hat*i4'*V4_hat + 16*i2'*a2_hat*i3'*a3_hat*i4'*a4_hat*u1'*V1_hat...
        + 16*i3'*a3_hat*i4'*a4_hat*u1'*a1_hat*i2'*V2_hat + 16*i4'*a4_hat*u1'*a1_hat*i2'*a2_hat*i3'*V3_hat...        
        + 16*i1'*a1_hat*u2'*a2_hat*i3'*a3_hat*i4'*V4_hat + 16*u2'*a2_hat*i3'*a3_hat*i4'*a4_hat*i1'*V1_hat...
        + 16*i3'*a3_hat*i4'*a4_hat*i1'*a1_hat*u2'*V2_hat + 16*i4'*a4_hat*i1'*a1_hat*u2'*a2_hat*i3'*V3_hat...        
        + 16*i1'*a1_hat*i2'*a2_hat*u3'*a3_hat*i4'*V4_hat + 16*i2'*a2_hat*u3'*a3_hat*i4'*a4_hat*i1'*V1_hat...
        + 16*u3'*a3_hat*i4'*a4_hat*i1'*a1_hat*i2'*V2_hat + 16*i4'*a4_hat*i1'*a1_hat*i2'*a2_hat*u3'*V3_hat...        
        + 16*i1'*a1_hat*i2'*a2_hat*i3'*a3_hat*u4'*V4_hat + 16*i2'*a2_hat*i3'*a3_hat*u4'*a4_hat*i1'*V1_hat...
        + 16*i3'*a3_hat*u4'*a4_hat*i1'*a1_hat*i2'*V2_hat + 16*u4'*a4_hat*i1'*a1_hat*i2'*a2_hat*i3'*V3_hat;
    
    m10 = n11; m11 = n11; 
    m20 = n21; m21 = n21;
    m30 = n31; m31 = n31;
    m40 = n41; m41 = n41;
    u1 = m10.*m11.*i1; u2 = m20.*m21.*i2; u3 = m30.*m31.*i3; u4 = m40.*m41.*i4; 
    bp(3) = 4*a1_hat'*(u1*i2'+i1*u2')*a2_hat + 4*a3_hat'*(u3*i4'+i3*u4')*a4_hat...
        + 4*a1_hat'*(u1*i4'+i1*u4')*a4_hat - 2*u2'*a2_hat - 2*u4'*a4_hat...
        - 8*(a4_hat'*u4*a2_hat'*i2*a3_hat'*i3 + a4_hat'*i4*a2_hat'*u2*a3_hat'*i3 + a4_hat'*i4*a2_hat'*i2*a3_hat'*u3)...
        + 16*a1_hat'*u1*a2_hat'*i2*a3_hat'*i3*a4_hat'*i4 + 16*a2_hat'*u2*a1_hat'*i1*a3_hat'*i3*a4_hat'*i4...
        + 16*a3_hat'*u3*a2_hat'*i2*a1_hat'*i1*a4_hat'*i4 + 16*a4_hat'*u4*a2_hat'*i2*a3_hat'*i3*a1_hat'*i1; 
    Ap(3, 1:length(phik)) = 4*a1_hat'*(u1*i2'+i1*u2')*V2_hat + 4*a2_hat'*(u2*i1'+i2*u1')*V1_hat...
        + 4*a1_hat'*(u1*i4'+i1*u4')*V4_hat + 4*a4_hat'*(u4*i1'+i4*u1')*V1_hat...
        + 4*a3_hat'*(u3*i4'+i3*u4')*V4_hat + 4*a4_hat'*(u4*i3'+i4*u3')*V3_hat...
        - 2*u2'*V2_hat - 2*u4'*V4_hat...
        - 8*a4_hat'*u4*a2_hat'*i2*i3'*V3_hat - 8*a4_hat'*u4*a3_hat'*i3*i2'*V2_hat - 8*a2_hat'*i2*a3_hat'*i3*u4'*V4_hat...
        - 8*a2_hat'*u2*a3_hat'*i3*i4'*V4_hat - 8*a2_hat'*u2*a4_hat'*i4*i3'*V3_hat - 8*a4_hat'*i4*a3_hat'*i3*u2'*V2_hat...
        - 8*a3_hat'*u3*a2_hat'*i2*i4'*V4_hat - 8*a3_hat'*u3*a4_hat'*i4*i2'*V2_hat - 8*a2_hat'*i2*a4_hat'*i4*u3'*V3_hat...
        + 16*u1'*a1_hat*i2'*a2_hat*i3'*a3_hat*i4'*V4_hat + 16*i2'*a2_hat*i3'*a3_hat*i4'*a4_hat*u1'*V1_hat...
        + 16*i3'*a3_hat*i4'*a4_hat*u1'*a1_hat*i2'*V2_hat + 16*i4'*a4_hat*u1'*a1_hat*i2'*a2_hat*i3'*V3_hat...        
        + 16*i1'*a1_hat*u2'*a2_hat*i3'*a3_hat*i4'*V4_hat + 16*u2'*a2_hat*i3'*a3_hat*i4'*a4_hat*i1'*V1_hat...
        + 16*i3'*a3_hat*i4'*a4_hat*i1'*a1_hat*u2'*V2_hat + 16*i4'*a4_hat*i1'*a1_hat*u2'*a2_hat*i3'*V3_hat...        
        + 16*i1'*a1_hat*i2'*a2_hat*u3'*a3_hat*i4'*V4_hat + 16*i2'*a2_hat*u3'*a3_hat*i4'*a4_hat*i1'*V1_hat...
        + 16*u3'*a3_hat*i4'*a4_hat*i1'*a1_hat*i2'*V2_hat + 16*i4'*a4_hat*i1'*a1_hat*i2'*a2_hat*u3'*V3_hat...        
        + 16*i1'*a1_hat*i2'*a2_hat*i3'*a3_hat*u4'*V4_hat + 16*i2'*a2_hat*i3'*a3_hat*u4'*a4_hat*i1'*V1_hat...
        + 16*i3'*a3_hat*u4'*a4_hat*i1'*a1_hat*i2'*V2_hat + 16*u4'*a4_hat*i1'*a1_hat*i2'*a2_hat*i3'*V3_hat;
    bp = bp(:);
end

if Np == 2 && Nd == 2
    A = A;
    b = b;
elseif Np == 4 && Nd == 2
    A = [A; Ap]; 
    b = [b; bp];
elseif Np == 2 && Nd == 4
    A = [A; Ad];
    b = [b; bd];
elseif Np == 4 && Nd == 4
    A = [A; Ap; Ad];
    b = [b; bp; bd];
end