import scipy
import scipy.signal as signal
import numpy
import numpy.fft as fft
from collections import deque
from bisect import insort



__author__ = "Andrea Maffezzoli"
__version__ = "$Revision: 0.13 - 21 Jan 2011$"
__name__ = "signalprocessing"
__all__ = ["IIRfilter", "spectrum", "smooth"]


#########################################################
# We have to append here to __all__ the modules imported
#########################################################
def IIRfilter(highband, fs, poleNr = 4, ftype='bessel'):
    """Bessel IIR filter highband pass flter"""
    
    if highband is None or not(fs):
        #print 'Filtering is disabled..'
        def filter(X): return X
    else:
        normalfreq = 2.*highband/fs
        b, a = scipy.signal.iirfilter(poleNr, normalfreq, btype='highpass',ftype='bessel')
        print 'High pass band filtering at %.3fHz..' % highband
        def filter(X): return scipy.signal.lfilter(b, a, X)
    return filter


def spectrum(T): return numpy.abs(numpy.fft.rfft(T)) ** 2

def smooth(x, window_len = 11, window = 'hanning'):
    """Smooth the data using a window with requested size.
    This method is based on the convolution of a scaled window with the signal.
    The signal is prepared by introducing reflected copies of the signal 
    (with the window size) in both ends so that transient parts are minimized
    in the begining and end part of the output signal.
    input:
        x: the input signal 
        window_len: the dimension of the smoothing window; should be an odd integer
        window: the type of window from 'flat', 'hanning', 'hamming', 'bartlett', 'blackman'
            flat window will produce a moving average smoothing.
    output:
        the smoothed signal
    example:
        t=linspace(-2,2,0.1)
        x=sin(t)+randn(len(t))*0.1
        y=smooth(x)
    
    See also: 
        numpy.hanning, numpy.hamming, numpy.bartlett, numpy.blackman, numpy.convolve
        scipy.signal.lfilter
    TODO: the window parameter could be the window itself if an array instead of a string   
    """
    
    x = numpy.array(x)
    if x.ndim != 1:
        raise ValueError, "smooth only accepts 1 dimension arrays."
    if x.size < window_len:
        raise ValueError, "Input vector needs to be bigger than window size."
    if window_len<3: return x
    for w in ('flat', 'hanning', 'hamming', 'bartlett', 'blackman'):
        if window == w: break
    else: raise ValueError, "Window is on of 'flat', 'hanning', 'hamming', 'bartlett', 'blackman'"
    s=numpy.r_[2*x[0]-x[window_len:1:-1],x,2*x[-1]-x[-1:-window_len:-1]]
    if window == 'flat': # a moving average
        w = numpy.ones(window_len,'d')
    else:
        w = eval('numpy.' + window + '(window_len)')
    y=numpy.convolve(w/w.sum(),s,mode='same')
    return y[window_len-1:-window_len+1]


####################################################################
# from http://code.activestate.com/recipes/577059-running-median/
####################################################################
def grouper(seq, n):
    'Split sequence into n-length deques'
    return [deque(seq[i:i+n]) for i in range(0, (len(seq)+n-1)//n*n, n)]


"""
# Some test
if __name__ == 'signalprocessing':
    
    sin, pi = numpy.sin, numpy.pi
    fs = 1000.
    t = numpy.linspace(-pi, pi, int(2*pi*fs))
    #t = numpy.arange(-pi, pi, 1/fs.)
    #x = sin(2*pi*t) + sin((2*pi/10)*t) +
    x = sin(5*2*pi*t) #+2*r.rand(n)
    _filter = IIRfilter(5, fs, poleNr = 8)
    pylab.figure()
    xf = _filter(x)
    pylab.plot(t, xf, '-')
    pylab.show()
"""


"""
def FIRfilter(lowband = None, highband = None):

    n, Va = 101, 1
    if lowband and not(highband):
        #print 'Low pass filter enabled'
        b = signal.firwin(n, cutoff = 2*lowband/self.Frequency, window = "hamming")
    elif not(lowband) and highband:
        #print 'High pass filter enabled'
        b = -signal.firwin(n, cutoff = 2*highband/self.Frequency, window = "hanning")
        b[n/2] = b[n/2] + 1
    else:
        #print 'Band pass filter enabled'
        b1 = signal.firwin(n, cutoff = 2*lowband/self.Frequency, window = 'blackmanharris')
        b2 = - signal.firwin(n, cutoff = 2*highband/self.Frequency, window = 'blackmanharris')
        b2[n/2] = b2[n/2] + 1
        b = - (b1+b2); b[n/2] = b[n/2] + 1
    filter = lambda X: signal.lfilter(b,a,X)
    return filter
"""
"""
def FIRfilter(lowband = None, highband = None):

    print lowband, highband 
    if not(lowband) and not(highband):
        filter = lambda X: X
    else:
        n = 101
        a = 1
        if lowband and not(highband):
            print 'Low pass filter enabled:', 2*lowband/Frequency
            b = signal.firwin(n, cutoff = 2*lowband/Frequency, window = "hamming")
        elif not(lowband) and highband:
            print 'High pass filter enabled:', 2*highband/Frequency
            b = -signal.firwin(n, cutoff = 2*highband/Frequency, window = "hanning")
            b[n/2] = b[n/2] + 1
        else:
            print 'Band pass filter enabled:', 2*lowband/Frequency,  2*highband/Frequency
            b1 = signal.firwin(n, cutoff = 2*lowband/Frequency, window = 'blackmanharris')
            b2 = - signal.firwin(n, cutoff = 2*highband/Frequency, window = 'blackmanharris')
            b2[n/2] = b2[n/2] + 1
            b = - (b1+b2); b[n/2] = b[n/2] + 1
        filter = lambda X: signal.lfilter(b,a,X)
    return filter


def mfreqz(b,a=1):
    w,h = signal.freqz(b,a)
    h_dB = 20 * log10 (abs(h))
    subplot(211)
    plot(w/max(w),h_dB)
    ylim(-150, 5)
    ylabel('Magnitude (db)') 
    xlabel(r'Normalized Frequency (x$\pi$rad/sample)')
    title(r'Frequency response')
    subplot(212)
    h_Phase = unwrap(arctan2(imag(h),real(h)))
    plot(w/max(w),h_Phase)
    ylabel('Phase (radians)') 
    xlabel(r'Normalized Frequency (x$\pi$rad/sample)')
    title(r'Phase response')
    h_Phase = unwrap(arctan2(imag(h),real(h)))
    plot(w/max(w),h_Phase)
    ylabel('Phase (radians)') 
    xlabel(r'Normalized Frequency (x$\pi$rad/sample)')
    title(r'Phase response')
    subplots_adjust(hspace=0.5)

def impz(b,a=1):
    impulse = repeat(0.,50); impulse[0] =1.
    x = arange(0,50)
    response = signal.lfilter(b,a,impulse)
    subplot(211)
    stem(x, response)
    ylabel('Amplitude') 
    xlabel(r'n (samples)')
    title(r'Impulse response')
    subplot(212)
    step = cumsum(response)
    stem(x, step)
    ylabel('Amplitude') 
    xlabel(r'n (samples)')
    title(r'Step response')
    subplots_adjust(hspace=0.5)
"""

"""
Python code for basis FIR filter design

from pylab import *
import scipy.signal as signal

#Plot frequency and phase response
def mfreqz(b,a=1):
    w,h = signal.freqz(b,a)
    h_dB = 20 * log10 (abs(h))
    subplot(211)
    plot(w/max(w),h_dB)
    ylim(-150, 5)
    ylabel('Magnitude (db)')
    xlabel(r'Normalized Frequency (x$\pi$rad/sample)')
    title(r'Frequency response')
    subplot(212)
    h_Phase = unwrap(arctan2(imag(h),real(h)))
    plot(w/max(w),h_Phase)
    ylabel('Phase (radians)')
    xlabel(r'Normalized Frequency (x$\pi$rad/sample)')
    title(r'Phase response')
    subplots_adjust(hspace=0.5)

#Plot step and impulse response
def impz(b,a=1):
    l = len(b)
    impulse = repeat(0.,l); impulse[0] =1.
    x = arange(0,l)
    response = signal.lfilter(b,a,impulse)
    subplot(211)
    stem(x, response)
    ylabel('Amplitude')
    xlabel(r'n (samples)')
    title(r'Impulse response')
    subplot(212)
    step = cumsum(response)
    stem(x, step)
    ylabel('Amplitude')
    xlabel(r'n (samples)')
    title(r'Step response')
    subplots_adjust(hspace=0.5)


#Lowpass FIR filter 
figure(1)
n = 61
a = signal.firwin(n, cutoff = 0.3, window = "hamming")
mfreqz(a)
show()
figure(4)
impz(a)
show()


#Highpass FIR filter
figure(2)
n = 101
a = signal.firwin(n, cutoff = 0.3, window = "hanning")
a = -a
a[n/2] = a[n/2] + 1
mfreqz(a)
show()

#Bandpass FIR 0.3 - 0,5
figure(3)
n = 1001
a = signal.firwin(n, cutoff = 0.3, window = 'blackmanharris')
b = - signal.firwin(n, cutoff = 0.5, window = 'blackmanharris'); b[n/2] = b[n/2] + 1
d = - (a+b); d[n/2] = d[n/2] + 1
#Frequency response
mfreqz(d)
show()

"""
