# (c) Rory Clune 4/14/2011
from numpy import*

npts = 2**12

#w contains the wavelet coefficients that we start with
w = zeros((npts))
w[4] = 1 
w[29] = 1

#num_cycles is how many times we need to do the inverse matrix operation
num_cycles = log(npts)/log(2)-1

#The c-coefficients for the matrix M
c0 = (1+sqrt(3))/(4*sqrt(2))
c1 = (3+sqrt(3))/(4*sqrt(2))
c2 = (3-sqrt(3))/(4*sqrt(2))
c3 = (1-sqrt(3))/(4*sqrt(2))

for cycle in range(0,int(num_cycles)):
    n = 2**(cycle+1)
    # w vector is made up of x's (all together at the start) and w's (all together at the end)
    # the first part (size n) of w is the x coefficients...
    x = w[0:n]
    # ...and the next part (size n) of w is the w coefficients
    w_small = w[n:2*n]   
        
    # set up the w vector (2n*2n) that will be premultiplied by the matrix to give the new (bigger..size 2n) x. w now has alternating x's and w's
    w[range(0,2*n,2)] = x
    w[range(1,2*n,2)] = w_small
        
    # Premultiply (just use elements instead of building the whole transform matrix c0, c1, c2, c3) to give x. 
    # The matrix, if I had constructed it, would be 2nx2n, sp the new x has to be of length 2n
    x = zeros((2*n))
    x[0] = c0*w[0] + c3*w[1] + c2*w[2*n-2] + c1*w[2*n-1] #should this be npts or 2*n?
    x[1] = c1*w[0] - c2*w[1] + c3*w[2*n-2] - c0*w[2*n-1]
#    x[0] = c0*w[0] + c3*w[1] + c2*w[npts-2] + c1*w[npts-1] #should this be npts or 2*n?
#    x[1] = c1*w[0] - c2*w[1] + c3*w[npts-2] - c0*w[npts-1]
    
    for i in range(2,2*n,2):        
        x[i] = c2*w[i-2] + c1*w[i-1] + c0*w[i] + c3*w[i+1]
        x[i+1] = c3*w[i-2] - c0*w[i-1] + c1*w[i] - c2*w[i+1]   
    #put the new x vector at the start of the w, and repeat until all the w's have become x's when 2n = 2**12
    w[0:2*n] = x    
    
from pylab import*
plot(w)
show()