Python signal.buttord函数代码示例

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