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

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


# ---------- 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-2 Balancing Conditions`);
k := 'k':
for k from 0 to 10 do
	lprint(coeff(R2,z,k),`,`);
od;


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

lprint(`// 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)):

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

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


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


lprint(`// Normalization`);

nor0 := subs(z=1,H0) - 1:
nor1 := subs(z=1,H1) - 1:

lprint(nor0,`,`);
lprint(nor1,`;`);

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

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