with(LinearAlgebra):
with(ArrayTools):
with(combinat):
with(PolynomialTools):

CosetShift := proc(a)
    description "Shift coset with Quincunx Dilation":
    return eval(a, {z[1]=-z[1], z[2]=-z[2]}):
end proc:

# currently only works for real input
hc2D := proc(a)
    description "Hermitian Conjugate of 2D filter":
    # only works for real case now.
    
    return eval(a, {z[1]=1/z[1], z[2]=1/z[2]}):
end proc:

# currently only works for real input
getM_QCX := proc(a)
    description "Get matrix M for Quincunx Dilation":

    local a1, v, v1:

    a1:= eval(a, {z[1] = -z[1], z[2] = -z[2]}):
    v := <a, a1>:
    v1 := Transpose(eval(v, {z[1]=1/z[1], z[2]=1/z[2]})):

    return v.v1;
end proc:

getCoset_QCX := proc(a)

    local an, a0, a1, a0new, a1new, M, Minv, de, de1, de2, c, t:
    # find 2 coset sequences
    an := eval(a, {z[1]=-z[1], z[2]=-z[2]}):
    a0 := simplify(a+an)/2:
    a0 := expand(a0):
    a1 := (a - a0)/z[1]:
    a1 := expand(a1):

    a0new := 0:
    a1new := 0:

    # Quincunx Dilation
    M := Matrix([[1, 1], [1, -1]]):
    Minv := MatrixInverse(M):

    # downsample
    if type(a0, `+`) then
        for t in [op(a0)] do
            de1 := degree(t, z[1]):
            de2 := degree(t, z[2]):
            de := <de1, de2>:
            de := Minv.de:
            c := coeff(coeff(t, z[1], de1), z[2], de2):

            a0new := a0new + c * z[1]^de[1] * z[2]^de[2]:
        od:
    else
        de1 := degree(a0, z[1]):
        de2 := degree(a0, z[2]):
        de := <de1, de2>:
        de := Minv.de:
        c := coeff(coeff(a0, z[1], de1), z[2], de2):
        a0new := a0new + c * z[1]^de[1] * z[2]^de[2]:
    end if;
    
    if type(a1, `+`) then
        for t in [op(a1)] do
            de1 := degree(t, z[1]):
            de2 := degree(t, z[2]):
            de := <de1, de2>:
            de := Minv.de:
            c := coeff(coeff(t, z[1], de1), z[2], de2):

            a1new := a1new + c * z[1]^de[1] * z[2]^de[2]:
        od:
    else
        de1 := degree(a1, z[1]):
        de2 := degree(a1, z[2]):
        de := <de1, de2>:
        de := Minv.de:
        c := coeff(coeff(a1, z[1], de1), z[2], de2):
        a1new := a1new + c * z[1]^de[1] * z[2]^de[2]:
    end if;

    return a0new, a1new:

end proc:

upsample_QCX := proc(a0)
    description "Upsample 2D filter a with Quincunx Dilation Matrix":

    local a, M, ld1, d1, ld2, d2, de, di1, di2, c:
    a := 0:
    M := Matrix([[1, 1], [1, -1]]):

    ld1 := ldegree(a0, z[1]):
    d1 := degree(a0, z[1]):
    ld2 := ldegree(a0, z[2]):
    d2 := degree(a0, z[2]):
    for di1 from ld1 to d1 do
        for di2 from ld2 to d2 do
            c := coeff(coeff(a0, z[1], di1), z[2], di2):
            de := <di1, di2>:
            de := M.de:
            a := a + c * z[1]^de[1]*z[2]^de[2]:
        end do:
    end do:

    return a:
    
end proc: