import numpy as np
from numpy import *
import matplotlib
matplotlib.use('TkAgg')
from matplotlib import pyplot as plt
from matplotlib import animation

# low pass filter
# b0*1 + b1*1 + b2*1 + b3*1 = 0
# b0*0 + b1*1 + b2*2 + b3*3 = 0

# high pass filter
# c0*b0 + c1*b1 + c2*b2 + c3*b3 = 0
# b0 = c3, b1 = -c2, b2 = c1, b3 = -c0

dim = 12
N = 2**dim

def formCoeffMatrix(c,b,N):
	# Form the matrix of coefficients
	rowc = append(c, zeros(N-len(c)))
	rowb = append(b, zeros(N-len(b)))
	cc = zeros([N,N])
	for i in arange(0,N,2):
		cc[i] = roll(rowc, i)
		cc[i+1] = roll(rowb, i)
	return cc

def formSplititngMatrix(N):
	# Form the splitting matrix
	row1 = append([1, 0], zeros(N-2))
	row2 = append([0, 1], zeros(N-2))
	sp = zeros([N,N])
	for i in arange(0,N/2,2):
		sp[i] = roll(row1, 2*i)
		sp[i+1] = roll(row2, 2*i+1)
	j = 0
	for i in arange(N/2,N-1,2):
		sp[i] = roll(row2, 2*j)
		sp[i+1] = roll(row1, 2*j+3) 
		j = j+2
	sp[N-1] = roll(row2, 2*j)
	return sp

c = [(1+pow(3,.5))/(4*pow(2,.5)), (3+pow(3,.5))/(4*pow(2,.5)), 
	 (3-pow(3,.5))/(4*pow(2,.5)), (1-pow(3,.5))/(4*pow(2,.5))]
b = [-c[3], -c[2], c[1], -c[0]]
xw = zeros(N)
xw[4] = 1
xw[29] = 1

for i in range(2,dim+1):
	print i
	filt = formCoeffMatrix(c, b, 2**i)
	split = formSplititngMatrix(2**i)
	xw[:2**i] = dot(dot(filt.T,split.T),xw[:2**i])

fig = plt.figure()
ax = plt.axes(xlim=(0,len(xw)), ylim=(min(xw)*1.2,max(xw)*1.2))
ax.plot(xw)
plt.show()