python类butter()的实例源码

def homomorphic_envelope(x, fs=1000, f_LPF=8, order=3):
    """
    Computes the homomorphic envelope of x

    Args:
        x : array
        fs : float
            Sampling frequency. Defaults to 1000 Hz
        f_LPF : float
            Lowpass frequency, has to be f_LPF < fs/2. Defaults to 8 Hz
    Returns:
        time : numpy array
    """
    b, a = butter(order, 2 * f_LPF / fs, 'low')
    he = np.exp(filtfilt(b, a, np.log(np.abs(hilbert(x)))))
    return he

def butter_highpass(highcut, fs, order):
    nyq = 0.5 * fs
    high = highcut / nyq
    b, a = butter(order, high, btype="highpass")
    return b, a


# Sources for Batch Iterators
#
# These classes load training and test data and perform some basic preprocessing on it.
# They are then passed to factory functions that create the net. There they are used
# as data sources for the batch iterators that feed data to the net.
# All classes band pass or low pass filter their data based on min / max freq using
# a causal filter (lfilter) when the data is first loaded.
# * TrainSource loads a several series of EEG data and events, splices them together into
#   one long stream, then normalizes the EEG data to zero mean and unit standard deviation.
# * TestSource is like TrainSource except that it uses the mean and standard deviation
#   computed for the associated training source to normalize the EEG data.
# * SubmitSource is like TestSource except that it does not load and event data.

def highpass(data, BUTTER_ORDER=3, sampling_rate=10000, cut_off=500.0):
    Wn = (float(cut_off) / (float(sampling_rate) / 2.0), 0.95)
    b, a = signal.butter(BUTTER_ORDER, Wn, 'pass')
    return signal.filtfilt(b, a, data)

def __init__(self, threshold_freq, order=2, design='cheby1'):
        """
        :param threshold_freq: Threshold frequency to filter out, as a fraction of sampling freq.  E.g. 0.1
        """
        self.b, self.a = \
            butter(N=order, Wn=threshold_freq, btype='low') if design == 'butter' else \
            cheby1(N=order, rp=0.1, Wn=threshold_freq, btype='low') if design == 'cheby1' else \
            bad_value(design)
        self.filter_state = lfilter_zi(b=self.b, a=self.a)

def lpf(x, cutoff, fs, order=5):
    """
    low pass filters signal with Butterworth digital
    filter according to cutoff frequency

    filter uses Gustafsson’s method to make sure
    forward-backward filt == backward-forward filt

    Note that edge effects are expected

    Args:
        x      (array): signal data (numpy array)
        cutoff (float): cutoff frequency (Hz)
        fs       (int): sample rate (Hz)
        order    (int): order of filter (default 5)

    Returns:
        filtered (array): low pass filtered data
    """
    nyquist = fs / 2
    b, a = butter(order, cutoff / nyquist)
    if not np.all(np.abs(np.roots(a)) < 1):
        raise PsolaError('Filter with cutoff at {} Hz is unstable given '
                         'sample frequency {} Hz'.format(cutoff, fs))
    filtered = filtfilt(b, a, x, method='gust')
    return filtered

def butter_bandpass(lowcut, highcut, order=5):
    nyq = 0.5 * 250
    low = lowcut / nyq
    high = highcut / nyq
    b, a = butter(order, [low, high], btype='band')
    return b, a

def __init__(self,window_size, low, high):
        fs_Hz = 250;
        fn = fs_Hz/2
        self.filtered_data = np.array((window_size,1))


        #######################################
        # Filter Creation
        # -------------------------------------
        #
        # Create a filter using the scipy module,
        # based on specifications suggested by
        # Pan-Tompkins (bandpass from 5-15Hz)
        #
        #
        # 1) Establish constants:
        #       a) filter_order = 2
        #       b) high pass cutoff = 15Hz
        #       c) low pass cutoff = 5Hz
        # 2) Calculate the coefficients, store in variables

        filter_order = 2
        f_high = high
        f_low = low
        self.high_pass_coefficients = signal.butter(filter_order,f_low/fn, 'high')
        self.low_pass_coefficients = signal.butter(filter_order,f_high/fn, 'low')


    #######################################
    # Bandpass filter
    # -------------------------------------
    # Filter the data, using a bandpass of 
    # 5-15Hz.
    # 
    # Input: 
    #           the data buffer from Data_Buffer class
    # Output:
    #           filtered data as a numpy array

def butter_bandpass(lowcut, highpass, fs, order=4):
       nyq = 0.5 * fs
#       low = lowcut / nyq
       high = highpass / nyq
       b, a = butter(order, high, btype='highpass')
       return b, a

def butter_bandpass(self, lowcut, highcut, order=5):
        nyq = 0.5 * self.sample_rate
        low = lowcut / nyq
        high = highcut / nyq
        b, a = butter(order, [low, high], btype='band')
        return b, a

def bfilt(cop_dat, cutoff, order):
    # Butterworth filter
    b,a = signal.butter(order, cutoff)
    # filter coronal
    cor_lp = signal.filtfilt(b, a, cop_dat[:,0])
    # filter sagittal
    sag_lp = signal.filtfilt(b, a, cop_dat[:,1])
    return np.concatenate((to_sing(cor_lp),to_sing(sag_lp),to_sing(cop_dat[:,2])),axis=1)

def Filter(self, LowCorner, HighCorner, Order=3):
    """
    Butterworth bandpass filter
    """

    FS = 1./self.HDR['TSMP']

    if HighCorner >= FS/2.:
      print 'Warning: High corner must be < {0:.2f} Hz'.format(FS/2.)
      return

    if LowCorner < 0.:
      print 'Warning: Low corner must be > 0 Hz'.format(FS/2.)
      return

    # Corner frequencies
    Corners = [2.*LowCorner/FS, 2.*HighCorner/FS]

    # Butterworth filter
    b, a = _sig.butter(Order, Corners, btype='band')

    # Filtering records
    for I,S in enumerate(self.CHN):
      # self.CHN[I] = _sig.lfilter(b, a, S)
      zi = _sig.lfilter_zi(b, a);
      self.CHN[I],_ = _sig.lfilter(b, a, S, zi=zi*S[0])

#-----------------------------------------------------------------------------------------

def _analog_filter(self,
                       ftype,
                       highpass=None,
                       lowpass=None,
                       order=3,
                       zerophase=True):
        ' Use a classic analog filter on the data, currently butter or bessel'
        from scipy.signal import butter, bessel
        filter_types = {'butter': butter, 'bessel': bessel}
        afilter = filter_types[ftype]
        if highpass is None and lowpass is not None:
            b, a = afilter(order, lowpass / (self.sr / 2), btype='lowpass')
        elif highpass is not None and lowpass is None:
            b, a = afilter(order, highpass / (self.sr / 2), btype='highpass')
        elif highpass is not None and lowpass is not None:
            if highpass < lowpass:
                b, a = afilter(order,
                               (highpass / (self.sr / 2),
                                lowpass / (self.sr / 2)),
                               btype='bandpass')
            else:
                b, a = afilter(order,
                               (lowpass / (self.sr / 2),
                                highpass / (self.sr / 2)),
                               btype='bandstop')
        if zerophase:
            return self.filtfilt(b, a)
        else:
            return self.lfilter(b, a)

def butter(self, highpass=None, lowpass=None, order=3, zerophase=True):
        ' Buttworth filter the data'
        return self._analog_filter('butter', highpass, lowpass, order,
                                   zerophase)

def abs_and_smooth(x, sr, lp=100):
    abs_x = np.abs(x)
    if len(x.shape) > 1:
        abs_x = np.sum(abs_x,
                       axis=-1)  # sum over last dimension eg sum over channels
    b, a = butter(3, lp / 2 / sr, btype="low")
    filtered_x = filtfilt(b, a, abs_x, axis=0)
    return filtered_x

def test_filtfilt():
    sr = 1000
    freq = 100
    b, a = butter(3, freq / (sr / 2), 'high')
    for data in (data2, data3, data4):
        x = filtfilt(b, a, data, axis=0)
        y = Stream(data, chunksize=211, sr=sr).filtfilt(b, a).call()
        assert eq(x, y)

def low_pass_filter(self, sig, low_cut, high_cut, nyq):

        """ 
        Apply a low pass filter to the data to take the relevant frequencies
        and to smooth the drop-out regions

        No resample necessary because all files have the same sampling frequency
        """

        b, a = signal.butter(5, np.array([low_cut, high_cut]) / nyq, btype='band')
        sig_filt = signal.lfilter(b, a, sig, axis=0)

        return(np.float32(sig_filt))

def butter_bandpass_filter(data, lowcut, highcut, fs, order):
    from scipy.signal import butter, lfilter

    nyq = 0.5 * fs
    low = lowcut / nyq
    high = highcut / nyq

    b, a = butter(order, [low, high], btype='bandpass')
    y = lfilter(b, a, data)
    return y

def butter_bandstop_filter(data, lowcut, highcut, fs, order):
    from scipy.signal import butter, lfilter

    nyq = 0.5 * fs
    low = lowcut / nyq
    high = highcut / nyq

    b, a = butter(order, [low, high], btype='bandstop')
    y = lfilter(b, a, data)
    return y

def butter_lowpass_filter(data, cut, fs, order, zero_phase=False):
    from scipy.signal import butter, lfilter, filtfilt

    nyq = 0.5 * fs
    cut = cut / nyq

    b, a = butter(order, cut, btype='low')
    y = (filtfilt if zero_phase else lfilter)(b, a, data)
    return y

def apply_butter_filter(self, data, order, cutoff, btype):
        cutoff /= self.nyq
        b, a = butter(N=order, Wn=cutoff, btype=btype)
        return filtfilt(b, a, data, axis=1)

def _getFiltDesign(sf, f, npts, filtname, cycle, order, axis):
    """Get the designed filter
    sf : sample frequency
    f : frequency vector/list [ex : f = [2,4]]
    npts : number of points
    - 'fir1'
    - 'butter'
    - 'bessel'
    """

    if type(f) != np.ndarray:
        f = np.array(f)

    # fir1 filter :
    if filtname == 'fir1':
        fOrder = fir_order(sf, npts, f[0], cycle=cycle)
        b, a = fir1(fOrder, f/(sf / 2))

    # butterworth filter :
    elif filtname == 'butter':
        b, a = butter(order, [(2*f[0])/sf, (2*f[1])/sf], btype='bandpass')
        fOrder = None

    # bessel filter :
    elif filtname == 'bessel':
        b, a = bessel(order, [(2*f[0])/sf, (2*f[1])/sf], btype='bandpass')
        fOrder = None

    def filtSignal(x):
        return filtfilt(b, a, x, padlen=fOrder, axis=axis)

    return filtSignal

def butter_bandpass(lowcut, highcut, fs, order=5):
    """Returns a bandpass butterworth filter."""
    nyq = 0.5 * fs
    low = lowcut / nyq
    high = highcut / nyq
    b, a = butter(order, [low, high], btype='band')
    return b, a

def butter_lowpass(cutoff, fs, order=5):
    """Returns a lowpass butterworth filter."""
    nyq = 0.5 * fs
    normal_cutoff = cutoff / nyq
    b, a = butter(order, normal_cutoff, btype='low', analog=False)
    return b, a

def butterfilt(finname, foutname, fmt, fs, fl=5.0, fh=100.0, gpass=1.0, gstop=30.0, ftype='butter', buffer_len=100000, overlap_len=100, max_len=-1):
    """Given sampling frequency, low and high pass frequencies design a butterworth filter and filter our data with it."""
    fso2 = fs/2.0
    wp = [fl/fso2, fh/fso2]
    ws = [0.8*fl/fso2,1.4*fh/fso2]
    import pdb; pdb.set_trace()
    b, a = iirdesign(wp, ws, gpass=gpass, gstop=gstop, ftype=ftype, output='ba')
    y = filtfiltlong(finname, foutname, fmt, b, a, buffer_len, overlap_len, max_len)
    return y, b, a

def periodic_derivative(x,y,max_periods):
    plot = False
    Ns =len(x)
    b,a = signalbutter(8,2.0*max_periods/Ns)
    ymid =interp(x+0.5*(x[1]-x[0]),x,y,period=2*np_pi)
    yder = diff(ymid)/diff(x) 
    #yder = Ns/(max(x)-min(x))*fftpackdiff(y,1,Ns)
    yder_filt = deepcopy(yder)
    x_filt = deepcopy(x)
    x_filt = npappend(x_filt,x_filt[-1]+x_filt[1]-x_filt[0])
    yder_filt = signalfiltfilt(b,a,npappend(yder_filt,yder_filt[0]))
    if plot:
        plt.figure(1)
        plt.subplot(311)
        plt.plot(x, y)

        plt.subplot(312)
        plt.plot(x[0:-1],yder)

        plt.subplot(313)
        plt.plot(x_filt[0:-1],yder_filt)
        plt.show()
    return yder_filt


#x = numpy.array(range(100))
#y = numpy.sin(twopi*x/100)+numpy.sin(twopi*x/10)
#periodic_derivative(x,y,4)

def butterplot(fs,fc):
    """
    https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.butter.html
    """
    b, a = signal.butter(4, 100, 'low', analog=True)
    w, h = signal.freqs(b, a)
    ax = figure().gca()
    ax.semilogx(fs*0.5/np.pi*w, 20*np.log10(abs(h)))
    ax.set_title('Butterworth filter frequency response')
    ax.set_xlabel('Frequency [Hz]')
    ax.set_ylabel('Amplitude [dB]')
    ax.grid(which='both', axis='both')
    ax.axvline(fc, color='green') # cutoff frequency
    ax.set_ylim(-50,0)

def butter_lowpass(cutoff, fs, order=2):
    """
     This function calculates the butterworth filter coefficients for a given cutoff frequency and order.

     :param cutoff: Cutoff frequency as a proportion of the Nyquist frequency (Freq (Hz) / (Sample rate / 2))
     :param fs:     Sample rate of signal to be filtered
     :param order:  Filter order, defaults to second order filter, as required by timbral_metallic
     :return: returns the coefficients of a Butterworth filter
    """
    nyq = 0.5 * fs
    normal_cutoff = cutoff / nyq
    b, a = butter(order, normal_cutoff, btype='low', analog=False)
    return b, a

def butter_highpass(cutoff, fs, order=5):
    nyq = 0.5 * fs
    normal_cutoff = cutoff / nyq
    b, a = butter(order, normal_cutoff, btype='highpass', analog=False)
    return b, a

def butter_highpass(cutoff, fs, order=5):
        nyq = 0.5 * fs
        normal_cutoff = cutoff / nyq
        b, a = signal.butter(order, normal_cutoff, btype='high', analog=False)
        return b, a

def filter(self):
        path = os.getcwd()+'/trialGraspEventDetection_dataFiles'
        self.Fgr = np.sum(self.values[:, 9:15], axis=1)  # SAI
        self.Fgl = np.sum(self.values[:, 0:7], axis=1)  # SAI

        # can use this to plot in matlab graspeventdetection_plot.m
        np.savetxt(path+'/SAI_Fgr.txt', self.Fgr)
        # can use this to plot in matlab
        np.savetxt(path+'/SAI_Fgl.txt', self.Fgl)

        # 0.55*pi rad/samples
        b1, a1 = signal.butter(1, 0.55, 'high', analog=False)
        self.f_acc_x = signal.lfilter(b1, a1, self.acc_x, axis=-1, zi=None)
        self.f_acc_y = signal.lfilter(b1, a1, self.acc_y, axis=-1, zi=None)
        self.f_acc_z = signal.lfilter(b1, a1, self.acc_z, axis=-1, zi=None)
        # self.f_eff = signal.lfilter(b1, a1, self.eff, axis=-1, zi=None)
        # type(eff)
        self.FAII = np.sqrt(np.square(self.f_acc_x) +
                            np.square(self.f_acc_y) +
                            np.square(self.f_acc_z))
        # can use this to plot in matlab
        np.savetxt(path+'/FAII.txt', self.FAII)

        # subtract base values from the values array
        self.values1 = self.values - self.values.min(axis=0)
        # pass the filter for each sensor
        self.fvalues1 = np.zeros(self.values1.shape)
        # 0.48*pi rad/samples
        b, a = signal.butter(1, 0.48, 'high', analog=False)
        for i in range(16):
            self.fvalues1[:, i] = signal.lfilter(b, a, self.values1[:, i],
                                                 axis=-1, zi=None)
        self.FAI = np.sum(self.fvalues1, axis=1)
        # can use this to plot in matlab
        np.savetxt(path+'/FAI.txt', self.FAI)