def fm_afskDemod(sig, TBW=4, N=74, fs=48000.0):
    #TODO: add docstring
    #  non-coherent demodulation of afsk1200
    # function returns the NRZI (without rectifying it)
    baud = 1200.0
    bandwidth = 2*500.0 + baud

    #TODO fix this
    M = int(2/(bandwidth/fs))
    h = signal.firwin(74.0, bandwidth, nyq=fs/2)

    t = r_[0.0:len(h)]/fs
    fc = 1700.0
    h = h*np.exp(1j*2*np.pi*fc*t)
    output = signal.fftconvolve(h, sig)

    temp = output*np.conjugate(np.roll(output, 1))
    NRZ_fm = np.angle(temp)/3

    h2 = signal.firwin(74.0, 1200, nyq=fs/2)
    NRZ_fm = signal.fftconvolve(NRZ_fm, h2)

    NRZ_fm = (NRZ_fm*fs/(2.0*np.pi)-550.0)/500.0

    return NRZ_fm

def nc_afsk1200Demod(sig, fs=48000.0, TBW=2.0):
    #  non-coherent demodulation of afsk1200
    # function returns the NRZ (without rectifying it)
    # 
    # sig  - signal
    # baud - The bitrate. Default 1200
    # fs   - sampling rate in Hz
    # TBW  - TBW product of the filters
    #
    # Returns:
    #     NRZ  
    # your code here
    taps = fs/1200-1
    bandpass = signal.firwin(taps, 1200, nyq=fs/2)
    spacepass = bandpass * np.exp(1j*2*np.pi*1200*np.r_[0.0:taps]/fs)
    markpass = bandpass * np.exp(1j*2*np.pi*3600*np.r_[0.0:taps]/fs)
    spaces = signal.fftconvolve(sig, spacepass, mode='same')
    marks = signal.fftconvolve(sig, markpass, mode='same')

    analog = np.abs(spaces)-np.abs(marks)
    lowpass = signal.firwin(taps, 2400*1.2, nyq=fs/2)
    filtered = signal.fftconvolve(analog, lowpass, mode='same')
    NRZ = filtered
    
    return NRZ

    def fir_filter(self, fir_ac=None, fir_dc=None, f_ac=None, f_dc=None,
                   a_ac=10, a_dc=10, alpha=None, filter_name=None, **kwargs):
        """Apply filters to generate the lock-in and dc components of phi"""

        if filter_name == 'bessel_matched':
            N_pts = kwargs.get('N_pts', int(self.ks / self.k0_dc * 6))
            dec = kwargs.get('dec', 32)
            n_pts_eval_fir = kwargs.get('n_pts_eval_fir', 2**16)
            window = kwargs.get('window', 'hann')

            fir_ac, fir_dc = _matched_filters(self.ks, self.x_m, N_pts, dec, window,
                                              n_pts_eval_fir)

            self.fir_ac = fir_ac
            self.fir_dc = fir_dc
        else:
            if fir_ac is None:
                if f_ac is None and alpha is None:
                    f_ac = self.fx * 0.5
                elif alpha is not None:
                    f_ac = self.v_tip/self.x_m * alpha
                self.fir_ac = signal.firwin(self.fs / (f_ac) * a_ac,
                                            f_ac, nyq=0.5 * self.fs,
                                            window='blackman')
            else:
                self.fir_ac = fir_ac

            if fir_dc is None:
                if f_dc is None and alpha is None:
                    f_dc = self.fx * 0.5
                elif alpha is not None:
                    f_dc = self.v_tip/self.x_m * alpha
                self.fir_dc = signal.firwin(self.fs/(f_dc) * a_dc,
                                            f_dc, nyq=0.5*self.fs,
                                            window='blackman')
            else:
                self.fir_dc = fir_dc

        indices = np.arange(self.phi.size)
        fir_ac_size = self.fir_ac.size
        fir_dc_size = self.fir_dc.size

        fir_max_size = max(fir_ac_size, fir_dc_size)

        self.m = indices[fir_max_size//2: -fir_max_size//2]
        self.tm = self.t[self.m]

        self._lock = np.exp(np.pi * 2j * self.fx * self.t)

        self.phi_lock = signal.fftconvolve(self.phi * self._lock * 2,
                                           self.fir_ac,
                                           mode='same')

        self.V_lock = self.phi_lock

        self.phi_lock_a = np.abs(self.phi_lock)
        self.phi_lock_phase = np.angle(self.phi_lock)

        self.phi_dc = signal.fftconvolve(self.phi, self.fir_dc, mode='same')
        self.V_dc = self.phi_dc

    def demodulate2(self,samples):
        # DEMODULATION CODE

        # LIMITER goes here

        # low pass & down sampling
        h = signal.firwin(128,80000,nyq=1.2e5)
        lp_samples = signal.fftconvolve(samples, h)

    # polar discriminator

        A = lp_samples[1:lp_samples.size]
        B = lp_samples[0:lp_samples.size-1]

        dphase = ( A * np.conj(B) ) / np.pi

        dphase.resize(dphase.size+1)
        dphase[dphase.size-1] = dphase[dphase.size-2]

        h = signal.firwin(128,16000,nyq=1.2e5)
        rebuilt = signal.fftconvolve(dphase,h)

        output = rebuilt[::self.decim_r2]

        output = self.lowpass(output, self.audioFilterSize)

        return np.real(output)

def designLinearBandpass(fa, fb, s_step, n=1001, show=False):
    f_sampling = 1/s_step    
    nyq = f_sampling / 2
    
    fa_ny = fa / nyq
    fb_ny = fb / nyq
    
    a = signal.firwin(n, cutoff = fa_ny, window = 'blackmanharris')
    #Highpass filter with spectral inversion
    b = - signal.firwin(n, cutoff = fb_ny, window = 'blackmanharris')
    # b[n/2] = b[n/2] + 1
    #Combine into a bandpass filter
    # d = - (a+b); d[n/2] = d[n/2] + 1


    b[nyq/2] = b[nyq/2] + 1
    #Combine into a bandpass filter
    d = - (a+b)
    d[nyq/2] = d[nyq/2] + 1


    #Frequency response
    if show:
        mfreqz(d)
        pylab.show()
    
    return d

def test_fir():
    print "Testing dsp-fir"

    #ref = ones(512)
    ref = (2.0 * random.rand(512)) - 1.0

    #test short mono fir
    writeaudio(ref)
    h = signal.firwin(21, 0.4)
    savetxt("test_coeffs.txt", h)
    expected = signal.lfilter(h, 1, ref)
    writeaudio(expected, 'expected.wav')
    os.system("../file-qdsp -n 64 -i test_in.wav -o test_out.wav -p fir,h=test_coeffs.txt")
    compareaudio(expected, readaudio(), 1e-6)

    #test long mono fir
    writeaudio(ref)
    h = signal.firwin(312, 0.4)
    savetxt("test_coeffs.txt", h)
    expected = signal.lfilter(h, 1, ref)
    os.system("../file-qdsp -n 64 -i test_in.wav -o test_out.wav -p fir,h=test_coeffs.txt")
    compareaudio(expected, readaudio(), 1e-6)

    #test short stereo fir, mono coeffs
    writeaudio(transpose([ref,-ref]))
    h = signal.firwin(21, 0.4)
    savetxt("test_coeffs.txt", h)
    expected = signal.lfilter(h, 1, ref)
    os.system("../file-qdsp -n 64 -i test_in.wav -o test_out.wav -p fir,h=test_coeffs.txt")
    compareaudio(transpose([expected, -expected]), readaudio(), 1e-6)

    #test long stereo fir, mono coeffs
    writeaudio(transpose([ref,-ref]))
    h = signal.firwin(312, 0.4)
    savetxt("test_coeffs.txt", h)
    expected = signal.lfilter(h, 1, ref)
    os.system("../file-qdsp -n 64 -i test_in.wav -o test_out.wav -p fir,h=test_coeffs.txt")
    compareaudio(transpose([expected, -expected]), readaudio(), 1e-6)

    #test asymmetric mono fir
    writeaudio(ref)
    impulse = concatenate(([1], zeros(499)))
    b, a = signal.butter(2, 500.0/24000, 'low')
    h = signal.lfilter(b, a, impulse)
    savetxt("test_coeffs.txt", h)
    expected = signal.lfilter(h, 1, ref)
    os.system("../file-qdsp -n 64 -i test_in.wav -o test_out.wav -p fir,h=test_coeffs.txt")
    compareaudio(expected, readaudio(), 1e-6)

    #test asymmetric stereo fir
    writeaudio(transpose([ref,-ref]))
    impulse = concatenate(([1], zeros(499)))
    b, a = signal.butter(2, 500.0/24000, 'low')
    h = signal.lfilter(b, a, impulse)
    savetxt("test_coeffs.txt", h)
    expected = signal.lfilter(h, 1, ref)
    os.system("../file-qdsp -n 64 -i test_in.wav -o test_out.wav -p fir,h=test_coeffs.txt")
    compareaudio(transpose([expected, -expected]), readaudio(), 1e-6)

    os.remove('test_coeffs.txt')

 def test_scaling(self):
     """
     For one lowpass, bandpass, and highpass example filter, this test
     checks two things:
       - the mean squared error over the frequency domain of the unscaled
         filter is smaller than the scaled filter (true for rectangular
         window)
       - the response of the scaled filter is exactly unity at the center
         of the first passband
     """
     N = 11
     cases = [
         ([.5], True, (0, 1)),
         ([0.2, .6], False, (.4, 1)),
         ([.5], False, (1, 1)),
     ]
     for cutoff, pass_zero, expected_response in cases:
         h = firwin(N, cutoff, scale=False, pass_zero=pass_zero, window='ones')
         hs = firwin(N, cutoff, scale=True, pass_zero=pass_zero, window='ones')
         if len(cutoff) == 1:
             if pass_zero:
                 cutoff = [0] + cutoff
             else:
                 cutoff = cutoff + [1]
         assert_(self.mse(h, [cutoff]) < self.mse(hs, [cutoff]),
             'least squares violation')
         self.check_response(hs, [expected_response], 1e-12)

    def __init__(self, fs = 48000.0, Abuffer = 1024, Nchunks=43):
        
        #  Implementation of an afsk1200 TNC. 
        #
        #  The TNC processes a `Abuffer` long buffers, till `Nchunks` number of buffers are collected into a large one.
        #  This is because python is able to more efficiently process larger buffers than smaller ones.
        #  Then, the resulting large buffer is demodulated, sampled and packets extracted.
        #
        # Inputs:
        #    fs  - sampling rate
        #   TBW  -  TBW of the demodulator filters
        #   Abuffer - Input audio buffers from Pyaudio
        #   Nchunks - Number of audio buffers to collect before processing
        #   plla    - agressivness parameter of the PLL
        
        
        ## compute sizes based on inputs
        self.TBW = 2.0   # TBW for the demod filters
        self.N = (int(fs/1200*self.TBW)//2)*2+1   # length of the filters for demod
        self.fs = fs     # sampling rate   
        self.BW = self.TBW/(1.0*self.N/fs)      # BW of filter based on TBW
        self.Abuffer = Abuffer             # size of audio buffer
        self.Nchunks = Nchunks             # number of audio buffers to collect
        self.Nbuffer = Abuffer*Nchunks+self.N*3-3         # length of the large buffer for processing
        self.Ns = 1.0*fs/1200 # samples per symbol
        
        ## state variables for the modulator
        self.prev_ph = 0  # previous phase to maintain continuous phase when recalling the function
        
        ##  Generate Filters for the demodulator
        self.h_lp = signal.firwin(self.N,self.BW/fs*1.0,window='hanning')
        self.h_lpp = signal.firwin(self.N,self.BW*2*1.2/fs,window='hanning')
        self.h_space = self.h_lp*exp(1j*2*pi*(2200)*r_[-self.N/2:self.N/2]/fs)
        self.h_mark = self.h_lp*exp(1j*2*pi*(1200)*r_[-self.N/2:self.N/2]/fs)
        self.h_bp = signal.firwin(self.N,self.BW/fs*2.2,window='hanning')*exp(1j*2*pi*1700*r_[-self.N/2:self.N/2]/fs)


        ## PLL state variables  -- so conntinuity between buffers is preserved
        self.dpll = np.round(2.0**32 / self.Ns).astype(int32)    # PLL step
        self.pll =  0                # PLL counter
        self.ppll = -self.dpll       # PLL counter previous value -- to detect overflow
        self.plla = 0.74             # PLL agressivness (small more agressive)
        

        ## state variable to NRZI2NRZ
        self.NRZIprevBit = True  
        
        ## State variables for findPackets
        self.state='search'   # state variable:  'search' or 'pkt'
        self.pktcounter = 0   # counts the length of a packet
        self.packet = bitarray.bitarray([0,1,1,1,1,1,1,0])   # current packet being collected
        self.bitpointer = 0   # poiter to advance the search beyond what was already searched in the previous buffer

        ## State variables for processBuffer
        self.buff = zeros(self.Nbuffer)   # large overlapp-save buffer
        self.chunk_count = 0              # chunk counter
        self.oldbits = bitarray.bitarray([0,0,0,0,0,0,0])    # bits from end of prev buffer to be copied to beginning of new
        self.Npackets = 0                 # packet counter

def fir_filt_prep(length, start_ind, stop_ind, numtaps=1000, window="blackmanharris"):
    """preps the FIR filter."""
    if start_ind==0:
        return firwin(numtaps, stop_ind, pass_zero=True,
                  nyq=length/2.0,
                  window=window)
    return firwin(numtaps, [start_ind, stop_ind], pass_zero=False,
                  nyq=length/2.0,
                  window=window)

 def init(self):                             
        self.count0 = 0
        self.count1 = 0
        
        # Open the wav music file
        MusicFile = open(self.musicFileName, 'r')
        
        # Read the music samples into an array
        self.MusicFileArray = MusicFile.read()
        self.MusicFileArray = self.MusicFileArray[:(len(self.MusicFileArray) - 1)]
        self.MusicFileArray = self.MusicFileArray.split("\n")
        
        # Filter0
        coeffFilter0 = signal.firwin(self.FilterDict["Filter0"]['N'], [self.FilterDict["Filter0"]["fc0"], self.FilterDict["Filter0"]["fc1"]], nyq = (0.5 * self.FilterDict["Filter0"]["fs"]), window = 'hamming')
        ## Calculate coefficients in C
        print "C implementation"
        set_trace()
        ###############
        for i in range(len(coeffFilter0)):
               coeffFilter0[i] = example.hamming_win(self.FilterDict["Filter0"]['N'], i, self.FilterDict["Filter0"]["fc1"], self.FilterDict["Filter0"]["fs"]) 
        
        
        set_trace()
        ##################
        
        
        coeffFilter0Fp = []      
        for i in coeffFilter0:
               coeffFilter0Fp.append(int(fp_sign_to_bin(i, self.fpNotation), 2))
       
        self.FilterDict["Filter0"]["Coefficients"] = coeffFilter0Fp
        
        # Filter1
        coeffFilter1 = signal.firwin(self.FilterDict["Filter1"]['N'], [self.FilterDict["Filter1"]["fc0"], self.FilterDict["Filter1"]["fc1"]], nyq = (0.5 * self.FilterDict["Filter1"]["fs"]), window = 'hamming')
        coeffFilter1Fp = []      
        for i in coeffFilter1:
               coeffFilter1Fp.append(int(fp_sign_to_bin(i, self.fpNotation), 2))
       
        self.FilterDict["Filter1"]["Coefficients"] = coeffFilter1Fp              
        
        # Filter2
        coeffFilter2 = signal.firwin(self.FilterDict["Filter2"]['N'], [self.FilterDict["Filter2"]["fc0"], self.FilterDict["Filter2"]["fc1"]], nyq = (0.5 * self.FilterDict["Filter2"]["fs"]), window = 'hamming')
        coeffFilter2Fp = []      
        for i in coeffFilter2:
               coeffFilter2Fp.append(int(fp_sign_to_bin(i, self.fpNotation), 2))
       
        self.FilterDict["Filter2"]["Coefficients"] = coeffFilter2Fp              
        
        # Filter3
        coeffFilter3 = signal.firwin(self.FilterDict["Filter3"]['N'], [self.FilterDict["Filter3"]["fc0"], self.FilterDict["Filter3"]["fc1"]], nyq = (0.5 * self.FilterDict["Filter3"]["fs"]), window = 'hamming')
        coeffFilter3Fp = []      
        for i in coeffFilter3:
               coeffFilter3Fp.append(int(fp_sign_to_bin(i, self.fpNotation), 2))
       
        self.FilterDict["Filter3"]["Coefficients"] = coeffFilter3Fp

    def test_response(self):
        N = 51
        f = .5
        # increase length just to try even/odd
        h = firwin(N, f)  # low-pass from 0 to f
        self.check_response(h, [(.25,1), (.75,0)])

        h = firwin(N+1, f, window='nuttall')  # specific window
        self.check_response(h, [(.25,1), (.75,0)])

        h = firwin(N+2, f, pass_zero=False)  # stop from 0 to f --> high-pass
        self.check_response(h, [(.25,0), (.75,1)])

        f1, f2, f3, f4 = .2, .4, .6, .8
        h = firwin(N+3, [f1, f2], pass_zero=False)  # band-pass filter
        self.check_response(h, [(.1,0), (.3,1), (.5,0)])

        h = firwin(N+4, [f1, f2])  # band-stop filter
        self.check_response(h, [(.1,1), (.3,0), (.5,1)])

        h = firwin(N+5, [f1, f2, f3, f4], pass_zero=False, scale=False)
        self.check_response(h, [(.1,0), (.3,1), (.5,0), (.7,1), (.9,0)])

        h = firwin(N+6, [f1, f2, f3, f4])  # multiband filter
        self.check_response(h, [(.1,1), (.3,0), (.5,1), (.7,0), (.9,1)])

        h = firwin(N+7, 0.1, width=.03)  # low-pass
        self.check_response(h, [(.05,1), (.75,0)])

        h = firwin(N+8, 0.1, pass_zero=False)  # high-pass
        self.check_response(h, [(.05,0), (.75,1)])

def get_h_parameters(NFIR, fcut):
    """NFIR : length of FIR filter
    fcut: fraction of Nyquist for filter"""
    h = si.firwin(NFIR+1,fcut) * np.exp(2*1j*np.pi*np.arange(NFIR+1)*0.125)
    n = np.arange(2,NFIR+2)
    g = h[(1-n)%NFIR]*(-1)**(1-n)
    NFIR = np.fix((3./2.*NFIR))
    h1 = si.firwin(NFIR+1,2./3*fcut)*np.exp(2j*np.pi*np.arange(NFIR+1)*0.25/3.)
    h2 = h1*np.exp(2j*np.pi*np.arange(NFIR+1)/6.)
    h3 = h1*np.exp(2j*np.pi*np.arange(NFIR+1)/3.)  
    return (h, g, h1, h2, h3)

    def process(self,samples):

        #samples = sdr.read_samples(2.56e6)

        h = signal.firwin(256,80000,nyq=1.2e5)
        output = signal.fftconvolve(samples,h)
        output[:h.size/2] += self.prevConv1[h.size/2:]   #add the latter half of tail end of the previous convolution
        outputa = np.append(self.prevConv1[:h.size/2], output) # also delayed by half size of h so append the first half
        self.prevConv1 = output[output.size-h.size:]   # set the tail for next iteration
        lp_samples = outputa[:output.size-h.size]  # chop off the tail and decimate

        #lp_samples = output[::5]

        dmod = np.zeros(lp_samples.size)

        A = lp_samples[1:]
        B = lp_samples[:lp_samples.size-1]

        dmod[1:] = np.real(np.angle(A * np.conj(B))) / (np.pi)
        dmod[0] = np.real(np.angle(lp_samples[0] * np.conj(self.prevB))) / (np.pi)
        self.prevB = lp_samples[lp_samples.size-1]



        h = signal.firwin(256,1.6e4,nyq=1.2e5)
        output = signal.fftconvolve(dmod,h)
        output[:h.size/2] += self.prevConv2[h.size/2:]   #add the latter half of tail end of the previous convolution
        outputa = np.append(self.prevConv2[:h.size/2], output) # also delayed by half size of h so append the first half
        self.prevConv2 = output[output.size-h.size:]   # set the tail for next iteration
        audible = outputa[:output.size-h.size:5]  # chop off the tail and decimate
    

        #h = signal.firwin(128,1.6e4,nyq=24000)
        #output = signal.fftconvolve(audible,h)
        #output[:h.size/2] += prevConv3[h.size/2:]   #add the latter half of tail end of the previous convolution
        #outputa = np.append(prevConv3[:h.size/2], output) # also delayed by half size of h so append the first half
        #prevConvo3 = output[output.size-h.size:]   # set the tail for next iteration
        #audible = outputa[:output.size-h.size:5]  # chop off the tail and decimate
    
        #print audible.size
        #spec = gen_spec(audible,256)
        #show_image(spec)
   
        self.spec = np.roll(self.spec,26,axis=1)

        self.spec[:,:26] = gen_spec(np.real(audible),512)   ##np.abs(np.fft.fft(audible)[:audible.size/2:-4])
        spec = cv2.GaussianBlur(self.spec,(5,5),1,.75)
        spectsc = cv2.convertScaleAbs(self.spec,alpha=255/np.max(spec))
        spect = cv2.applyColorMap(spectsc,cv2.COLORMAP_JET)
        
        cv2.imshow('Spectrum',spect)
        cv2.waitKey(1)

        return np.real(.5*audible)

def createFilters():
    fs = 500 
    f_nyquist = fs/2 # Calling nyquist fs/2 because this corresponds to pi radians
    numtaps = 71 # must be less than sig_length / 3

    h_lp_ecg = signal.firwin(numtaps, 150./f_nyquist)
    h_lp_icg = signal.firwin(numtaps, 1./f_nyquist)
    h_lp_ppg = signal.firwin(numtaps, 30./f_nyquist) # Can be used for ICG - cardiac
    h_60 = signal.firwin(numtaps, [55./f_nyquist, 65./f_nyquist])
    h_hp = signal.firwin(numtaps, 5./f_nyquist, pass_zero=False)

    return(h_lp_ecg, h_lp_icg, h_lp_ppg, h_60, h_hp)

 def bandreject(self, S, low=700, high=800):
     # http://www.scipy.org/Cookbook/FIRFilter
     # http://mpastell.com/2010/01/18/fir-with-scipy/
     nyq_rate = self.spl_rate / 2.0
     width = 50.0 / nyq_rate
     ripple_db = 60.0
     N, beta = signal.kaiserord(ripple_db, width)
     tapsL = signal.firwin(N, low / nyq_rate, window=("kaiser", beta))
     tapsH = signal.firwin(N, high / nyq_rate, window=("kaiser", beta))
     tapsB = -(tapsL + tapsH)
     tapsB[N / 2] = tapsB[N / 2] + 1
     return signal.lfilter(tapsB, 1.0, S)

def nc_mafsk1200Demod(sig, fs=48000.0, baud=1200, TBW=2.0, fc = 2700, fd = 1000):
    #  non-coherent demodulation of afsk1200
    # function returns the NRZ (without rectifying it)
    # 
    # sig  - signal
    # baud - The bitrate. Default 1200
    # fs   - sampling rate in Hz
    # TBW  - TBW product of the filters
    #
    # Returns:
    #     NRZ 
    
    N = (int(fs/baud*TBW)//2)*2+1
    f11 = fc - 2*fd
    f10 = fc - fd
    f01 = fc + fd
    f00 = fc + 2*fd
    
    # your code here
    taps = TBW*fs/1200-1
    taps = N
    filt = signal.firwin(taps, baud/2, window='hanning', nyq=fs/2)
    #plt.plot(np.fft.fft(filt))
    #plt.plot(filt)
    #f1 = 1200
    #f2 = 2200
    t2 = (r_[0:fs]/fs)[:taps]
    filt11 = filt* np.exp(t2*1j*f11*-2*np.pi)
    filt10 = filt* np.exp(t2*1j*f10*-2*np.pi)
    sig11 = signal.fftconvolve(sig, filt11, mode="same")
    sig10 = signal.fftconvolve(sig, filt10, mode="same")
    filt01 = filt* np.exp(t2*1j*f01*-2*np.pi)
    filt00 = filt* np.exp(t2*1j*f00*-2*np.pi)
    sig01 = signal.fftconvolve(sig, filt01, mode="same")
    sig00 = signal.fftconvolve(sig, filt00, mode="same")
    midsig = 0#(max(sig00)+min(sig00)+max(sig11)+min(sig11))/4
    
    sig11r = sig11 - midsig
    sig10r = sig10 - midsig
    sig01r = sig01 - midsig
    sig00r = sig00 - midsig
    
    return sig11r, sig10r, sig01r, sig00r
    
    diff = np.abs(sig12k)-np.abs(sig22k)
    return diff
    opt = signal.firwin(taps, baud*1.2, window='hanning', nyq=fs/2)
    ana = signal.fftconvolve(diff, opt, mode="same")
    #sign = np.sign(ana)

    NRZ = ana
    return NRZ

def initialize_online_filter(fsample, highpass, lowpass, order, x, axis=-1):
    # boxcar, triang, blackman, hamming, hann, bartlett, flattop, parzen, bohman, blackmanharris, nuttall, barthann
    filtwin = 'nuttall'
    nyquist = fsample / 2.
    ndim = len(x.shape)
    axis = axis % ndim

    if highpass != None:
        highpass = highpass/nyquist
        if highpass < 0.01:
            highpass = None
        elif highpass > 0.99:
            highpass = None

    if lowpass != None:
        lowpass = lowpass/nyquist
        if lowpass < 0.01:
            lowpass = None
        elif lowpass > 0.99:
            lowpass = None

    if not(highpass is None) and not(lowpass is None) and highpass>=lowpass:
        # totally blocking all signal
        print 'using NULL filter', [highpass, lowpass]
        b = np.zeros(window)
        a = np.ones(1)
    elif not(lowpass is None) and (highpass is None):
        print 'using lowpass filter', [highpass, lowpass]
        b = firwin(order, cutoff = lowpass, window = filtwin, pass_zero = True)
        a = np.ones(1)
    elif not(highpass is None) and (lowpass is None):
        print 'using highpass filter', [highpass, lowpass]
        b = firwin(order, cutoff = highpass, window = filtwin, pass_zero = False)
        a = np.ones(1)
    elif not(highpass is None) and not(lowpass is None):
        print 'using bandpass filter', [highpass, lowpass]
        b = firwin(order, cutoff = [highpass, lowpass], window = filtwin, pass_zero = False)
        a = np.ones(1)
    else:
        # no filtering at all
        print 'using IDENTITY filter', [highpass, lowpass]
        b = np.ones(1)
        a = np.ones(1)

    # initialize the state for the filtering based on the previous data
    if ndim == 1:
        zi = zi = lfiltic(b, a, x, x)
    elif ndim == 2:
        f = lambda x : lfiltic(b, a, x, x)
        zi = np.apply_along_axis(f, axis, x)

    return b, a, zi

    def __init__(self, fs = 48000, baud = 2400, fc = 1800):
        self.fs = fs
        self.fc = fc
        self.baud = baud
        self.f_mark = self.fc - baud / 4
        self.f_space = self.fc + baud / 4
        self.delta_f = (self.f_space - self.f_mark) / 2
        self.spb = self.fs / baud

        taps = 2 * int(fs * 2 / baud) + 1
        h = signal.firwin(taps, self.baud / 2., nyq = fs / 2., window = 'hanning')
        self.bp_mark = np.exp(2j * np.pi * self.f_mark * np.r_[0:taps] / float(fs)) * h
        self.bp_space = np.exp(2j * np.pi * self.f_space * np.r_[0:taps] / float(fs)) * h
        self.h_nrz = signal.firwin(taps, self.baud * 1.2, nyq = fs / 2., window = 'hanning')

def bandpass(data, sampling, fmin, fmax, ripple_db=50, width=2.0,\
             return_filter=False, verbose=False):

  """
    This function will bandpass filter data in the given [fmin,fmax] band
    using a kaiser window.

    Arguments:

      data : numpy.ndarray
        array of data points
      sampling : int
        number of data points per second
      fmin : float
        frequency of lowpass
      fmax : float
        frequency of highpass

    Keyword arguments:

      ripple_db : int
        Attenuation in the stop band, in dB
      width : float
        Desired width of the transition from pass to stop, in Hz
      return_filter: boolean
        Return filter
      verbose : boolean
  """

  # construct filter
  order, beta = signal.kaiserord(ripple_db, width*2/sampling)

  lowpass = signal.firwin(order, fmin*2/sampling, window=('kaiser', beta))
  highpass = - signal.firwin(order, fmax*2/sampling, window=('kaiser', beta))
  highpass[order//2] = highpass[order//2] + 1

  bandpass = -(lowpass + highpass); bandpass[order//2] = bandpass[order//2] + 1

  # filter data forward then backward
  data = signal.lfilter(bandpass,1.0,data)
  data = data[::-1]
  data = signal.lfilter(bandpass,1.0,data)
  data = data[::-1]

  if verbose: sys.stdout.write("Bandpass filter applied to data.\n")

  if return_filter:
    return data, bandpass
  else:
    return data

def testfir(st,cb,ct,n):
    from scipy import signal
    import numpy as np
    
    sr = st[0].stats.sampling_rate/2.
    xx = np.empty([st.count(),n],)
    a = signal.firwin(n, cutoff = cb/sr, window = 'hamming')
    b = - signal.firwin(n, cutoff = ct/sr, window = 'hamming'); b[n/2] = b[n/2] + 1
    d = - (a+b); d[n/2] = d[n/2] + 1
    fft1 = np.abs(np.fft.fft(d))
    for i in range(st.count()):
        fft = np.fft.fft(st[i][:n])*fft1
        xx[i] = np.fft.ifft(fft)
    return xx