def butter_bandpass(lowcut, highcut, samplingrate, order=4):
    nyq = 0.5 * samplingrate
    low = lowcut / nyq
    high = highcut / nyq
    print high, low
    if high >=1. and low == 0.:
        b = np.array([1.])
        a = np.array([1.])

    elif high < 0.95 and low > 0. :
        wp = [1.05*low,high-0.05]
        ws = [0.95*low,high+0.05]
 
        print wp,ws
        order,wn = buttord(wp,ws,0., 30.)
        b, a = butter(order, wn, btype='band')
    
    elif high>= 0.95:
        print 'highpass',low,1.2*low,0.8*low
        order,wn = buttord( 15*low,0.05*low,gpass=0.0, gstop=10.0)
        print order,wn
        b, a = butter(order, wn, btype='high')
    elif low <= 0.05:
        print 'lowpass',high
        order,wn = buttord( high-0.05,high+0.05,gpass=0.0, gstop=10.0)
        b, a = butter(order, wn, btype='low')

    return b, a

def prepare_audio_filters():
    tf_rangel = 100000
    tf_rangeh = 170000

    # audio filters
    tf = SysParams["audio_lfreq"]
    N, Wn = sps.buttord(
        [(tf - tf_rangel) / (freq_hz / 2.0), (tf + tf_rangel) / (freq_hz / 2.0)],
        [(tf - tf_rangeh) / (freq_hz / 2.0), (tf + tf_rangeh) / (freq_hz / 2.0)],
        5,
        15,
    )
    Faudl = filtfft(sps.butter(N, Wn, btype="bandpass"))

    tf = SysParams["audio_rfreq"]
    N, Wn = sps.buttord(
        [(tf - tf_rangel) / (freq_hz / 2.0), (tf + tf_rangel) / (freq_hz / 2.0)],
        [(tf - tf_rangeh) / (freq_hz / 2.0), (tf + tf_rangeh) / (freq_hz / 2.0)],
        5,
        15,
    )
    Faudr = filtfft(sps.butter(N, Wn, btype="bandpass"))

    N, Wn = sps.buttord(0.016 / (afreq / 2.0), 0.024 / (afreq / 2.0), 5, 15)
    FiltAPost = filtfft(sps.butter(N, Wn))

    N, Wn = sps.buttord(3.1 / (freq / 2.0), 3.5 / (freq / 2.0), 1, 20)
    SysParams["fft_audiorf_lpf"] = Faudrf = filtfft(sps.butter(N, Wn, btype="lowpass"))

    SysParams["fft_audiolpf"] = FiltAPost  # * FiltAPost * FiltAPost

    SysParams["fft_audio_left"] = Faudrf * Faudl * fft_hilbert
    SysParams["fft_audio_right"] = Faudrf * Faudr * fft_hilbert

def prepare_audio_filters():
    forder = 256
    forderd = 0

    tf_rangel = 100000
    tf_rangeh = 170000
    # audio filters
    tf = SP["audio_lfreq"]
    N, Wn = sps.buttord(
        [(tf - tf_rangel) / (freq_hz / 2.0), (tf + tf_rangel) / (freq_hz / 2.0)],
        [(tf - tf_rangeh) / (freq_hz / 2.0), (tf + tf_rangeh) / (freq_hz / 2.0)],
        1,
        15,
    )
    Baudl, Aaudl = sps.butter(N, Wn, btype="bandpass")

    tf = SP["audio_rfreq"]
    N, Wn = sps.buttord(
        [(tf - tf_rangel) / (freq_hz / 2.0), (tf + tf_rangel) / (freq_hz / 2.0)],
        [(tf - tf_rangeh) / (freq_hz / 2.0), (tf + tf_rangeh) / (freq_hz / 2.0)],
        1,
        15,
    )
    N, Wn = sps.buttord(
        [(tf - tf_rangel) / (freq_hz / 2.0), (tf + tf_rangel) / (freq_hz / 2.0)],
        [(tf - tf_rangeh) / (freq_hz / 2.0), (tf + tf_rangeh) / (freq_hz / 2.0)],
        5,
        15,
    )
    Baudr, Aaudr = sps.butter(N, Wn, btype="bandpass")

    N, Wn = sps.buttord(0.016 / (afreq / 2.0), 0.024 / (afreq / 2.0), 2, 15)
    audiolp_filter_b, audiolp_filter_a = sps.butter(N, Wn)

    # USE FIR
    audiolp_filter_b = sps.firwin(257, 0.020 / (afreq / 2.0))
    audiolp_filter_a = [1.0]

    N, Wn = sps.buttord(3.1 / (freq / 2.0), 3.5 / (freq / 2.0), 1, 20)
    audiorf_filter_b, audiorf_filter_a = sps.butter(N, Wn, btype="lowpass")

    [Baudrf_FDLS, Aaudrf_FDLS] = fdls.FDLS_fromfilt(audiorf_filter_b, audiorf_filter_a, forder, forderd, 0)
    SP["fft_audiorf_lpf"] = Faudrf = np.fft.fft(Baudrf_FDLS, blocklen)

    [Baudiolp_FDLS, Aaudiolp_FDLS] = fdls.FDLS_fromfilt(audiolp_filter_b, audiolp_filter_a, forder, forderd, 0)

    FiltAPost = np.fft.fft(Baudiolp_FDLS, blocklen)
    SP["fft_audiolpf"] = FiltAPost * FiltAPost  # * FiltAPost

    [Baudl_FDLS, Aaudl_FDLS] = fdls.FDLS_fromfilt(Baudl, Aaudl, forder, forderd, 0)
    [Baudr_FDLS, Aaudr_FDLS] = fdls.FDLS_fromfilt(Baudr, Aaudr, forder, forderd, 0)

    Faudl = np.fft.fft(Baudl_FDLS, blocklen)
    Faudr = np.fft.fft(Baudr_FDLS, blocklen)

    SP["fft_audio_left"] = Faudrf * Faudl * fft_hilbert
    SP["fft_audio_right"] = Faudrf * Faudr * fft_hilbert

    def ButterworthBandpass(self,data,sampleRate,low=0,high=None,minFreq=0,maxFreq=None,order=10,band=50):
        """ Basic IIR bandpass filter.
        Identifies order of filter, max 10.

        """
        if data is None:
            data = self.data
        if sampleRate is None:
            sampleRate = self.sampleRate
        if high is None:
            high = sampleRate/2
        if maxFreq is None:
            maxFreq = sampleRate/2
        low = max(low,0,minFreq)
        high = min(high,maxFreq,sampleRate/2)

        if low == minFreq and high == maxFreq:
            print("No filter needed!")
            return data

        nyquist = sampleRate/2

        if low == minFreq:
            # Low pass
            cut1 = high/nyquist
            cut2 = (high+band)/nyquist
            # calculate the best order
            order,wN = signal.buttord(cut1, cut2, 3, band)
            if order>10:
                order=10
            b, a = signal.butter(order,wN,btype='lowpass')
        elif high == maxFreq:
            # High pass
            cut1 = low/nyquist
            cut2 = (low-band)/nyquist
            # calculate the best order
            order,wN = signal.buttord(cut1, cut2, 3, band)
            if order>10:
                order=10
            b, a = signal.butter(order,wN, btype='highpass')
        else:
            # Band pass
            lowPass = low/nyquist
            highPass = high/nyquist
            lowStop = (low-band)/nyquist
            highStop = (high+band)/nyquist
            # calculate the best order
            order,wN = signal.buttord([lowPass, highPass], [lowStop, highStop], 3, band)
            if order>10:
                order=10
            b, a = signal.butter(order,wN, btype='bandpass')
            #b, a = signal.butter(order,[lowPass, highPass], btype='bandpass')

        return signal.filtfilt(b, a, data)

    def data_hpass(self, x, Wp, srate):
        ''' High-pass filter '''
        Wp = float(Wp*2/srate)
        Ws = Wp*float(self.lineEdit_19.text())
        Rp = float(self.lineEdit_17.text())
        Rs = float(self.lineEdit_18.text())

        tempstring = self.lineEdit_16.text()
        if tempstring == 'auto':
            if self.comboBox_2.currentIndex() == 0:
                (norder, Wn) = buttord(Wp, Ws, Rp, Rs)
            elif self.comboBox_2.currentIndex() == 1:
                (norder, Wn) = ellipord(Wp, Ws, Rp, Rs)
            else:
                (norder, Wn) = cheb1ord(Wp, Ws, Rp, Rs)
        else:
            norder = float(tempstring)
            Wn = Wp

        if self.comboBox_2.currentIndex() == 0:
            (b, a)  =  butter(norder, Wn, btype = 'high')
        elif self.comboBox_2.currentIndex() == 1:
            (b, a)  =  ellip(norder, Rp, Rs, Wn)
        else:
            (b, a)  =  cheby1(norder, Rp, Wn)


        y  =  filtfilt(b, a, x)

        return(y)

 def HPmin(self, fil_dict):
     self._get_params(fil_dict)
     self.N, self.F_PBC = buttord(self.F_PB,self.F_SB, self.A_PB,self.A_SB)
     if not self._test_N():
         return -1
     self._save(fil_dict, sig.bessel(self.N, self.F_PBC,
                         btype='highpass', analog=False, output=self.FRMT))

    def _build_lpfilter(self, fs):
        """
        builds low-pass filter with a cutoff frequency of 3/7th the resample
        frequency. The filter should be down 40 dB at 1.5 times the cutoff
        frequency (6/7th) the resample frequency.

        Parameters
        ----------
        fs : the base sampling rate

        Returns
        -------
        b, a : array_like
            Numerator (b) and denominator (a) polynomials of the IIR filter. 
        """
        nyq = fs/2. # nyquist frequency
        cutoff = (3./7.)*self.resample_fs # cutoff freq defined by Boer
        wp = cutoff * nyq # pass edge freq (pi radians / sample)
        ws = wp*2.        # pass edge freq (pi radians / sample)
        gpass = 1.5       # The maximum loss in the passband (dB)
        gstop = 40        # The minimum attenuation in the stopband (dB)
        n, wn = buttord(wp, ws, gpass, gstop)
        #print('n =',n,'wn =',wn)
        b, a = butter(n, wn, analog=True)

        return b, a

    def signal(self, fs, atten, caldb, calv):
        npts = self._duration*fs
        # start with full spectrum white noise and band-pass to get desired 
        # frequency range
        signal = self._noise[:npts]
        
        # band frequency cutoffs
        delta = 10**(3./(10.*(2*self._width)))
        low_freq = self._center_frequency / delta
        high_freq = self._center_frequency * delta
        # scipy butter function wants frequencies normalized between 0. and 1.
        nyquist = fs/2.
        low_normed = low_freq / nyquist
        high_normed = high_freq / nyquist

        order, wn = buttord([low_normed, high_normed], [low_normed-0.05, high_normed+0.05], 1, 40)

        # print 'CUTOFFS', low_freq, high_freq
        # print 'ORDER WN', order, wn, low_normed, high_normed

        b, a = butter(order, wn, btype='band')
        signal = lfilter(b, a, signal)

        if self._risefall > 0:
            rf_npts = int(self._risefall * fs) / 2
            wnd = hann(rf_npts*2) # cosine taper
            signal[:rf_npts] = signal[:rf_npts] * wnd[:rf_npts]
            signal[-rf_npts:] = signal[-rf_npts:] * wnd[rf_npts:]

        return signal

def _LP_Butterworth(interval, sampling_rate, cutoff, order=5):

    nyq = sampling_rate * 0.5

    stopfreq = float(cutoff)
    cornerfreq = 0.5 * stopfreq
    Ws = stopfreq / nyq
    Wp = cornerfreq / nyq
    N, Wn = buttord(Wp, Ws, 3, 20)  # (?)
    print "The oder of LPF is: %f" % N

    """
    Wp = 2 * np.pi * 100
    Ws = 2 * np.pi * 20
    Rp = 1.5
    Rs = 20
    N, Wn = buttord(Wp, Ws, Rp, Rs)   # (?)
    """

    # for hardcoded order:
    # N = order

    b, a = butter(N, Wn, btype="low")  # should 'high' be here for bandpass?
    # b, a = butter(9, float(20.0/nyq) , btype='high')   # should 'high' be here for bandpass?
    sf = lfilter(b, a, interval)
    return sf, b, a

def psd(y):
    # Number of samplepoints
    N = 128
    # sample spacing
    T = 1.0 / 128.0
    # From 0 to N, N*T, 2 points.
    #x = np.linspace(0.0, 1.0, N)
    #y = 1*np.sin(10.0 * 2.0*np.pi*x) + 9*np.sin(20.0 * 2.0*np.pi*x)


    fs = 128.0
    fso2 = fs/2
    Nd,wn = buttord(wp=[9/fso2,11/fso2], ws=[8/fso2,12/fso2],
       gpass=3.0, gstop=40.0)

    b,a = butter(Nd,wn,'band')
    y = filtfilt(b,a,y)


    yf = fft(y)
    #xf = np.linspace(0.0, 1.0/(2.0*T), N/2)
    #import matplotlib.pyplot as plt
    #plt.plot(xf, 2.0/N * np.abs(yf[0:N/2]))
    #plt.axis((0,60,0,1))
    #plt.grid()
    #plt.show()

    return np.sum(np.abs(yf[0:N/2]))

def lpf(m):
    cutoff = 9000
    norm_pass = cutoff/(m.fs/2)
    (N, Wn) = signal.buttord(wp=norm_pass, ws=1.5*norm_pass, gpass=2, gstop=50, analog=0)
    (b, a) = signal.butter(N, Wn, btype='lowpass', analog=0, output='ba')
    m.msg = signal.lfilter(b, a, m.msg)
    m.msg *= 1.0/np.abs(m.msg).max()
    return m

    def set_lowPass_filter(self, wp=20., ws=40., gpass=1., gstop=10.):
        Nq = self.fs/2.
        wp, ws = float(wp)/Nq, float(ws)/Nq
        gpass, gstop = float(gpass), float(gstop)
        N_filtr, Wn_filtr = buttord(wp, ws, gpass, gstop)
        self.b_L, self.a_L = butter(N_filtr, Wn_filtr, btype='low')

        self.N_L, self.Wn_L = N_filtr, Wn_filtr

 def BSmin(self, fil_dict):
     self._get_params(fil_dict)
     self.N, self.F_PBC = buttord([self.F_PB, self.F_PB2],
         [self.F_SB, self.F_SB2], self.A_PB,self.A_SB, analog = self.analog)
     if not self._test_N():
         return -1
     self._save(fil_dict, sig.butter(self.N, self.F_PBC, btype='bandstop',
                                    analog=self.analog, output=self.FRMT))

    def data_lpass(self, x, Wp, srate):
        ''' Low-pass filter using various filter type '''
        tempstring = self.lineEdit_16.text()
        
        if tempstring == 'auto':
            Wp = float(Wp*2/srate)
            Ws = Wp*float(self.lineEdit_19.text())
            Rp = float(self.lineEdit_17.text())
            Rs = float(self.lineEdit_18.text())
            
            if self.comboBox_2.currentIndex() == 0:
                (norder, Wn) = buttord(Wp, Ws, Rp, Rs)
            elif self.comboBox_2.currentIndex() == 1:
                (norder, Wn) = ellipord(Wp, Ws, Rp, Rs)
            else:
                (norder, Wn) = cheb1ord(Wp, Ws, Rp, Rs)

        else:
            norder = float(tempstring)
            Wp = float(Wp*2/srate)
            Ws = Wp*2
            self.lineEdit_19.setText(str(Ws/Wp))
            Rp = 3
            self.lineEdit_17.setText(str(Rp))
            Rs = 0.3*norder*20
            self.lineEdit_18.setText(str(Rs))
            
            if self.comboBox_2.currentIndex() == 0:
                (norder, Wn) = buttord(Wp, Ws, Rp, Rs)
            elif self.comboBox_2.currentIndex() == 1:
                (norder, Wn) = ellipord(Wp, Ws, Rp, Rs)
            else:
                (norder, Wn) = cheb1ord(Wp, Ws, Rp, Rs)
            
        if self.comboBox_2.currentIndex() == 0:
            (b, a)  =  butter(norder, Wn)
        elif self.comboBox_2.currentIndex() == 1:
            (b, a)  =  ellip(norder, Rp, Rs, Wn)
        else:
            (b, a)  =  cheby1(norder, Rp, Wn)
            
        
        y  =  filtfilt(b, a, x)
        
        return(y)

    def compute_parameters(self, target='stopband'):
        """ This function computes the order and the -3 dB-frequency
            of the filter for the specific parameters.

            Arguments:

                target: The optimization goal for the filter computation.
                Choices are:
                    - stopband: optimize to the stopband (like MATLAB)
                    - passband: optimize to the passband
        """

        if target not in ['passband', 'stopband', None]:
            raise ValueError("Target must be one of passband or stopband, \
                             or not given if filter is not Butterworth.")
        else:
            self.filter_target = target

        if True: # Change here to be more verbose.
            print("Ws = ", self.Ws)
            print("Wp = ", self.Wp)
            print("Rp = ", self.passband_attenuation)
            print("Rs = ", self.stopband_attenuation)

        if self.filter_class == 'butterworth':
            if target == 'passband':
                self.N, self.Wn = signal.buttord(self.Wp, self.Ws,
                                                 self.passband_attenuation,
                                                 self.stopband_attenuation,
                                                 analog=True)
            elif target == 'stopband':
                self.N, self.Wn = custom.custom_buttord(self.Wp, self.Ws,
                                                        self.passband_attenuation,
                                                        self.stopband_attenuation,
                                                        analog=True)
            else:
                raise ValueError("Butterworth filters must match either the \
                                 passband or the stopband.")
        elif self.filter_class == 'chebyshev_1':
            self.N, self.Wn = signal.cheb1ord(self.Wp, self.Ws,
                                              self.passband_attenuation,
                                              self.stopband_attenuation,
                                              analog=True)
        elif self.filter_class == 'chebyshev_2':
            self.N, self.Wn = signal.cheb2ord(self.Wp, self.Ws,
                                              self.passband_attenuation,
                                              self.stopband_attenuation,
                                              analog=True)
        elif self.filter_class == 'elliptical':
            self.N, self.Wn = signal.ellipord(self.Wp, self.Ws,
                                              self.passband_attenuation,
                                              self.stopband_attenuation,
                                              analog=True)
        else:
            raise NotImplementedError(
                "Filter family {} not yet implemented".format(self.filter_class))
        pass

def butter_bandpass(filt_freq,fs):
    nyq = 0.5 * fs

    wp=[filt_freq[1]/nyq, filt_freq[2]/nyq, ]
    ws=[filt_freq[0]/nyq, filt_freq[3]/nyq]
    N, wn = buttord(wp, ws, 3, 16)
    print N
    b, a = butter(N, wn, btype='band')
    return b, a

def smooth(time, vals, dt=None, gapFactor=20, T=T_M2):
    from scipy import signal
    if dt is None:
        dt = np.diff(time).mean()
    ta = timeArray.timeArray(time, 'epoch')
    # try to calculate exact dt by omitting large gaps
    gaps, ranges, t = ta.detectGaps(dt=dt, gapFactor=gapFactor)
    diff = []
    for i in range(ranges.shape[0]):
        twin = time[ranges[i, 0]:ranges[i, 1]]
        diff.append(np.diff(twin))
    diff = np.concatenate(tuple(diff), axis=0)
    dt = diff.mean()
    # filter design, low-pass butterworth
    T0 = (2 * dt)  # period of Nyquist frequency
    Tpass = 8 * T  # period of pass frequency
    Gpass = 3.0       # max dB loss in pass band
    Tstop = 1 * T  # period of stop frequency
    Gstop = 30.0     # min dB atennuation in stop band
    o, Wn = signal.buttord(T0 / Tpass, T0 / Tstop, Gpass, Gstop)
    if o < 0:
        raise Exception(
            'Cannot create tidal filter. Data sampling frequency may be too low, dt=' +
            str(dt))
    b, a = signal.butter(o, Wn, 'low')
    newvals = []
    newtime = []
    # filter each contiquous data range separately
    for i in range(ranges.shape[0]):
        twin = time[ranges[i, 0]:ranges[i, 1]]
        vwin = vals[ranges[i, 0]:ranges[i, 1]]
        if len(vwin) > 3 * len(a):
            try:
                # default forward-backward filter
                # filtered = signal.filtfilt(b, a, vwin, padtype='constant')
                # forward-backward filter with custom boundary conditions
                # pad with mean of 1/2 pass window lenght
                N_init = int(np.ceil(Tpass / dt / 2 / 4))
                # forward filter
                x_init = vwin[:N_init]
                y_init = x_init.mean() * np.ones_like(x_init)
                z_init = signal.lfiltic(b, a, y_init, x_init)
                filtered, _ = signal.lfilter(b, a, vwin, zi=z_init)
                # backward filter
                x_init = vwin[-N_init:][::-1]
                y_init = x_init.mean() * np.ones_like(x_init)
                z_init = signal.lfiltic(b, a, y_init, x_init)
                filtered, _ = signal.lfilter(b, a, filtered[::-1], zi=z_init)
                filtered = filtered[::-1]
                newvals.append(filtered)
                newtime.append(twin)
            except Exception as e:
                print a.shape, vwin.shape
                raise e
    newvals = np.concatenate(tuple(newvals), axis=0)
    newtime = np.concatenate(tuple(newtime), axis=0)
    return newtime, newvals

def highpass(s,stopBand,fs):
    ''' 
    highpass signal for noise reduction
    finch songdata should be in range 350Hz to 11kHz 
    '''
    N,wn = ss.buttord(wp=2*np.pi*(stopBand+100.)/fs,ws=2*np.pi*stopBand*1.0/fs,gpass=2.0, gstop=30.0,analog=0)
    b, a = ss.butter(N, wn,btype='high',analog=0)  
    filteredsong = ss.lfilter(b, a, s)
    filteredsong = np.round(filteredsong).astype('int16') # scipy wav needs 16 bit depth
    return filteredsong

def ex4_1():
    wp = 0.2*np.pi
    ws = 0.3*np.pi
    fs = 4000.0
    T = 1/fs
    Wp = wp/T
    Ws = ws/T
    n, x = signal.buttord(Wp, Ws, 2, 40, analog=True)
    b, a = signal.butter(n, x, analog=True)
    z, p, k = signal.tf2zpk(b, a)
    print z, p, k

 def calculate_bandpass(self):
     # Bandpass filter
     #self.width = 10000.0/self.nyquist_rate
     #self.cutoff = 2000.0/self.nyquist_rate
     self.ripple_db = 60.0
     self.high_pass = (self.tx_freq+1000.0)/self.nyquist_rate
     self.low_pass = (self.tx_freq-1000.0)/self.nyquist_rate
     self.high_cutoff = (self.tx_freq+1500.0)/self.nyquist_rate
     self.low_cutoff = (self.tx_freq-1500.0)/self.nyquist_rate
     ord, wn = signal.buttord([self.low_pass, self.high_pass],[self.low_cutoff,self.high_cutoff],1.0,30.0)
     self.b,self.a = signal.butter(ord,wn, btype="band")