# -*- coding: utf-8 -*-
"""
Created on Mon Nov 22 10:37:38 2010

@author: daryl

Presents a filter design based on the least squares
criterion. This is the Eigenfilter method described
by Vaidyanathan and Nguyen. It generalizes the method
for generating Prolate Spheroidal Sequences by
including some parameters that affect the passband.
This implementation uses quasi monte carlo integration
simply because I like it, though better results might
by achievable by obtaining an analytical solution to
the integrals in the paper and using them instead.

Eventually, I'll add routines for designing prolate
filters and nyquist filters as well.

If you find my code useful, I'd appreciate it if
you'd email me my email address is wasden <at>
eng <dot> utah <dot> edu.
"""

import numpy as np

def getKorobov(a,n,s):
    '''
    This generates a korobov lattice (a quasi monte carlo
    sequence used for qmc integration in order to generate
    the power matrix). This method is described in Monte
    Carlo and Quasi Monte Carlo Methods by
    Christiane Lemieux.
    
    a is the base/power
    n is the number of samples
    s is the dimension per sample 
    '''
    
    # Initialization
    z = np.zeros(s)
    z[0] = 1
    
    for idx in xrange(1,s):
        z[idx] = np.mod(a*z[idx-1],n)

    u = np.zeros((n,s),float)

    # Sequence Generation
    u[0,:] = np.mod(z/n,1)
    
    for idx in xrange(1,n):
        u[idx,:] = np.mod(u[idx-1,:] + z/n,1)
        
    # Perform Barker Transformation
    test = (u < 0.5)
    u[test] = 2*u[test]
    u[test==0] = 2*(1-u[test==0]) 
    
    test = (u == 1.0)
    u[test] = 1.0 - 1.0/n
    
    test = (u == 0.0)
    u[test] = 1.0/n
    
    return u

def eigenLS(wp,ws,N,alpha=0.5):
    '''
    Designs an optimal low-pass FIR filter using the
    eigenfilter method. This is an optimal filter in 
    the least squares sense. For more information see
    the paper by Vaidyanathan and Nguyen (1987). 
    
    wp = pass band edge (rads per sample)
    ws = stop band edge (rads per sample)
    N = filter length (must be odd)
    alpha = relative importance of stopband
    '''
    
    # Generate QMC sequences for numerical integration
    #w = getKorobov(123421,61711,5)
    w = getKorobov(123421,6171,5)
    w = w[:,2:4]    
    
    # Force the filter length to be odd.
    if N/2 == (N+1)/2:
        N = N + 1
    
    # Calculate the dimension of P
    M = (N-1)/2
    
    # Allocate memory for P
    P = np.zeros((M,M),float)
    
    # Calculate stopband error
    w1 = (np.pi - ws)*w + ws
    for n in xrange(0,M):
        for m in xrange(0,n+1):
            X = np.mean(np.cos(n*w1[:,0])*np.cos(m*w1[:,1]))
            X += np.mean(np.cos(n*w1[:,1])*np.cos(m*w1[:,0]))
            X = alpha*X/(2.0*np.pi)
            P[n,m] = X
            #P[m,n] = P[n,m]
    
    # Calculate passband error
    w2 = wp*w
    for n in xrange(0,M):
        for m in xrange(0,n+1):
            X = np.mean((1-np.cos(n*w2[:,0]))*(1-np.cos(m*w2[:,1])))
            X += np.mean((1-np.cos(n*w2[:,1]))*(1-np.cos(m*w2[:,0])))
            X = (1.0-alpha)*X/(2.0*np.pi)
            P[n,m] += X
            P[m,n] = P[n,m]
    
    # Compute eigenvalues/eigenvectors
    (q,r) = np.linalg.eig(P)
        
    # Retain the eigenvector corresponding to the smallest eigenvalue
    best = r[:,np.abs(q).argmin()].flatten()
    
    # Put the coefficients in the appropriate order
    h = np.zeros(N)
    
    h[M] = best[0]
    for idx in xrange(1,M):
        h[M-idx] = 0.5*best[idx]
        h[M+idx] = 0.5*best[idx]
    
    # Return the filter
    return h/np.sum(h)

# Begin Script
if __name__ == '__main__':
    import matplotlib.pyplot as pp
    import numpy.random as npr
    
    # Define Passband and Stopband edges 
    wp = 0.1*pi
    ws = 0.5*pi

    # Define relative importance of stopband to passband
    alpha = 0.5
    
    # Define filter length
    N = 31
    
    # Calculate filter
    h = eigenLS(wp,ws,N,alpha)
    
    # Compute Frequency Response
    waxis = 2*np.pi*np.arange(0,1.0,0.001)
    H = np.zeros(waxis.shape)
    Hang = np.zeros(waxis.shape)
    for idx,widx in enumerate(waxis):
        tmp = np.dot(h,np.exp(-1j*widx*np.arange(0,h.size)))
        H[idx] = 20*np.log10(np.abs(tmp))
        Hang[idx] = 180.0/np.pi*np.angle(tmp)
    
    # Plot Responses
    f1 = pp.figure()
    pp.subplot(2,1,1)
    pp.plot(waxis/(2.0*np.pi), H, linewidth=2.0)
    pp.xlabel('Frequency (normalized by 2 pi)', fontsize=22)
    pp.ylabel('Magnitude', fontsize=22)
    pp.subplot(2,1,2)
    pp.plot(waxis/(2.0*np.pi), Hang, linewidth=2.0)
    pp.xlabel('Frequency (normalized by 2 pi)', fontsize=22)
    pp.ylabel('Phase (degrees)', fontsize=22)
    
    # Plot Time Domain
    f2 = pp.figure()
    pp.plot(h, linewidth=2.0)
    pp.xlabel('Sample Number', fontsize=22)
    pp.ylabel('Coefficient', fontsize=22)