# filter H0
H0 := h0(0) + h0(1)*z + h0(0)*z^2;

# filter H1
H1 := h1(0) + h1(1)*z + h1(2)*z^2 + h1(3)*z^3 + 
	h1(2)*z^4  + h1(1)*z^5 + h1(0)*z^6;


# balancing conditions

P := z^3+z^2+z+1;
V1 := 1;
V2 := (3-z^4)/2;

R1 := rem(H0+V1*H1,P^1,z);
R2 := rem(H0+V2*H1,P^2,z);

interface(screenwidth=300);
writeto(`eqs`);

lprint(`// Order-1 Balancing Conditions`);

k := 'k':
for k from 0 to 5 do
        lprint(coeff(R1,z,k),`,`);	
od;


# ---------- orthogonality conditions ----------

H0H0 := expand(collect(H0*subs(z=1/z,H0),z)):
H1H1 := expand(collect(H1*subs(z=1/z,H1),z)):
H0H1 := expand(collect(H0*subs(z=1/z,H1),z)):

lprint(`// Cross Orthogonality Conditions `);
k := 'k':
for k from -7 to 7 do
	lprint(coeff(H0H1,z,4*k),`,`);
od;

lprint(`// Self Orthogonality Conditions `);
k := 'k':
for k from 1 to 5 do
        lprint(coeff(H0H0,z,4*k),`,`);
	lprint(coeff(H1H1,z,4*k),`,`);
od;


# ---------- normalization conditions ----------

nor0 := coeff(expand(H0H0),z,0) - 2:
nor1 := coeff(expand(H1H1),z,0) - 2:
norA := subs(z=1,expand(H0)) - 2:
norB := subs(z=1,expand(H1)) - 2:

lprint(`// Normalizations`);
lprint(norA,`,`);
lprint(norB,`,`);
lprint(nor0,`,`);
lprint(nor1,`;`);

        # note: put a semicolon at the end of the last eq

writeto(terminal);
interface(screenwidth=80);