## Athor: Jordi Pascual Prado
## CDI  L512

from cdi_D6 import *

canvas = plt.figure
plot = plt.plot
close = plt.close
axis = plt.axis

'''
    L512
 By analogy with L423, choose your favorite signal f (to be useful it should be not too wild)
 and compute the arrays Ar(f,r) and Dr(f,r) by using the scaling and wavelet vectors and by 
 using the method of filters. Present the results in suitable graphical representations.  
'''

# Orthogonal projection in orthonormal basis
def proj(f,V):
    x = zeros(len(V[0]))
    for v in V:
        x = x + dot(f,v)*v
    return x  

# Projection coefficients
def proj_coeffs(f,V):
    return array([dot(f,v) for v in V])

# auxiliary functions
def up_sample(a):
    x = []
    for t in a:
        x += [t, 0]
    return array(x)
U_ = up_sample

def dual(h):
    s = 1
    hd = []
    for t in reversed(h):
        hd += [s*t]
        s = -s
    return hd
    
def filter(h,x):
    m = len(h); n = len(x)
    y = []
    for k in range(m+n-1):
        a = max(0,k-n+1); b = min(m,k+1)
        s = sum(h[j]*x[k-j] for j in range(a,b))
        y += [s]
    for i in range(n,n+m-1):
        y[i%n] += y[i]
    return y[:n]
    
def H6(x): return filter(h6, x)
def G6(x): return filter(g6, x)

def HF6(f, r=1):
    if r == 0: return array(f)
    a = D6trend(f, r)
    for _ in range(r):
        a = H6(U_(a))
    return array(a)
    
def LF6(f, r=1):
    if r == 0: return zeros(len(f))
    d = G6(U_(D6fluct(f, r)))
    for _ in range(r - 1):
        d = H6(U_(d))
    return array(d)


def A(f,r): return proj(f,V[r])
def A1(f,r): return array(HF6(f,r))
def D(f,r): return proj(f,W[r-1])
def D1(f,r): return array(LF6(f,r))

# Mildly oscillating function
n = 10; N = 2**n
F = lambda t: 1500*t**2*(1-t)**4*((t-0.35)**2+0.01)*np.cos(41*t)
t = ls(0,1,N)
f = F(t)

lo =-1; hi=15; sep=1.7 

# Examples

V, W = D6VW(N)
close('all')

canvas('1. Scaling VS High-filters')
axis([0, N, lo, hi])

plt.subplot(211)
for r in range(n+1):
    plt.xlim(0,N)
    plot(sep*r+A(f,r))
    
plt.subplot(212)
for r in range(n+1):
    plt.xlim(0,N)
    plot(sep*r+A1(f,r))
'''   
plt.subplot(313)
for r in range(n+1):
    plot(sep*r+ (A(f,r)-A1(f,r)) )'''


canvas('2. Wavelets VS Low-filters')
axis([0, N, lo, hi])
# Be careful: index 0 for W is level 1
plt.subplot(211)
for r in range(1,n+1):
    plt.xlim(0,N)
    plot(sep*r+D(f,r))

plt.subplot(212)
for r in range(1,n+1):
    plt.xlim(0,N)
    plot(sep*r+D1(f,r))
    

'''
plt.subplot(313)
for r in range(1,n+1):
    plot(sep*r+ (D(f,r)-D1(f,r)) )'''