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 := : 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 := : 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 := : 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 := : 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 := : 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 := : de := M.de: a := a + c * z[1]^de[1]*z[2]^de[2]: end do: end do: return a: end proc: