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 parse(self, current):
        """
        Apply the filter-derivative method to filter the ionic current.
        """

        # Filter the current using a first order Bessel filter twice, one in
        # both directions to preserve phase
        from scipy import signal

        nyquist = self.sampling_freq / 2.0
        b, a = signal.bessel(1, self.cutoff_freq / nyquist, btype="low", analog=0, output="ba")
        filtered_current = signal.filtfilt(b, a, np.array(current).copy())

        # Take the derivative
        deriv = np.abs(np.diff(filtered_current))

        # Find the edges of the blocks which fulfill pass the lower threshold
        blocks = np.where(deriv > self.low_threshold, 1, 0)
        block_edges = np.abs(np.diff(blocks))
        tics = np.where(block_edges == 1)[0] + 1

        # Split points are points in the each block which pass the high
        # threshold, with a maximum of one per block
        split_points = [0]

        for start, end in it.izip(tics[:-1:2], tics[1::2]):  # For all pairs of edges for a block..
            segment = deriv[start:end]  # Save all derivatives in that block to a segment
            if (
                np.argmax(segment) > self.high_threshold
            ):  # If the maximum derivative in that block is above a threshold..
                split_points = np.concatenate((split_points, [start, end]))  # Save the edges of the segment
                # Now you have the edges of all transitions saved, and so the states are the current between these transitions
        tics = np.concatenate((split_points, [current.shape[0]]))
        tics = map(int, tics)
        return [Segment(current=current[tics[i] : tics[i + 1]], start=tics[i]) for i in xrange(0, len(tics) - 1, 2)]

	def filterData(self, icurr, Fs):
		"""
			Denoise an ionic current time-series and store it in self.eventData

			:Parameters:
				- `icurr` :	ionic current in pA
				- `Fs` :	original sampling frequency in Hz
		"""
		self.eventData=icurr
		self.Fs=Fs

		#pad the data with 10x the transient time at both ends to manually eliminate edge effects of the filter
		#for some reason I can't get good results using the pad method in filtfilt so manual it is
		#this means there may be some numerical artefacts but they should be well below the level of noise
		
		padding = int(10 * self.Fs/float(self.filterCutoff))
		paddedsignal = np.pad(self.eventData,pad_width=padding,mode='edge')
		
		b, a=sig.bessel(
							N=self.filterOrder, 
							Wn=(self.filterCutoff/(float(self.Fs)/2.0)), 
							btype='lowpass', 
							analog=False, 
							output='ba',
							norm='mag'
						)

		self.eventData=sig.filtfilt(b, a, paddedsignal, padtype=None, method='pad')[padding:-padding]

def calculate_dvdt(v, t, filter=None):
    """Low-pass filters (if requested) and differentiates voltage by time.

    Parameters
    ----------
    v : numpy array of voltage time series in mV
    t : numpy array of times in seconds
    filter : cutoff frequency for 4-pole low-pass Bessel filter in kHz (optional, default None)

    Returns
    -------
    dvdt : numpy array of time-derivative of voltage (V/s = mV/ms)
    """

    if has_fixed_dt(t) and filter:
        delta_t = t[1] - t[0]
        sample_freq = 1. / delta_t
        filt_coeff = (filter * 1e3) / (sample_freq / 2.) # filter kHz -> Hz, then get fraction of Nyquist frequency
        if filt_coeff < 0 or filt_coeff > 1:
            raise FeatureError("bessel coeff (%f) is outside of valid range [0,1]. cannot compute features." % filt_coeff)
        b, a = signal.bessel(4, filt_coeff, "low")
        v_filt = signal.filtfilt(b, a, v, axis=0)
        dv = np.diff(v_filt)
    else:
        dv = np.diff(v)

    dt = np.diff(t)
    dvdt = 1e-3 * dv / dt # in V/s = mV/ms

    # Remove nan values (in case any dt values == 0)
    dvdt = dvdt[~np.isnan(dvdt)]

    return dvdt

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 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)
     if not self._test_N():
         return -1
     self._save(fil_dict, sig.bessel(self.N, self.F_PBC,
                         btype='bandstop', analog=False, output=self.FRMT))

 def filter_data(self):
     cutoff = float(self.cutoff_entry.get())
     order = int(self.order_entry.get())
     Wn = 2.0 * cutoff/float(self.samplerate)
     b, a = bessel(order,Wn,'low')
     padding = 1000
     padded = np.pad(self.filtered_data, pad_width=padding, mode='median')
     self.filtered_data = filtfilt(b, a, padded, padtype=None)[padding:-padding]

def bessel_bandpass_filter(data, lowcut, highcut, fs, order=2):
    nyq = 0.5 * fs
    low = lowcut / nyq
    high = highcut / nyq

    # bessel() and lfilter() are from scipy.signal

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

 def test_bad_filter(self):
     """Regression test for #651: better handling of badly conditioned
     filter coefficients."""
     warnings.simplefilter("error", BadCoefficients)
     try:
         b, a = bessel(20, 0.1)
         z, p, k = tf2zpk(b, a)
         raise AssertionError("tf2zpk did not warn about bad "\
                              "coefficients")
     except BadCoefficients:
         pass
     finally:
         warnings.simplefilter("always", BadCoefficients)

def filter_design(sf, N, f, filtname='fir1', cycle=3, order=3, axis=0):
    if filtname == 'fir1':
        fOrder = fir_order(sf, N, f[0], cycle=cycle)
        b, a = fir1(fOrder, f/(sf / 2))
    elif filtname == 'butter':
        b, a = butter(order, [(2*f[0])/sf, (2*f[1])/sf], btype='bandpass')
        fOrder = None
    elif filtname == 'bessel':
        b, a = bessel(order, [(2*f[0])/sf, (2*f[1])/sf], btype='bandpass')
        fOrder = None

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

    return FiltDesign

def _gufunc_filtfilt_bessel(x, y, N, Wn, out=None):  # pragma: no cover
    # Pre-process
    xynm = preprocess_interp1d_nan_func(x, y, out)
    if xynm is None:  # all nan
        return
    x, y, x_even, y_even, num_nan, mask = xynm

    # filter function
    b, a = signal.bessel(N=N[0], Wn=Wn[0])
    # filter even data
    yf_even = signal.filtfilt(b, a, y_even, method='gust')

    # Post-process
    postprocess_interp1d_nan_func(x, x_even, yf_even, num_nan, mask, out)

def filt(sf, f, x, btype='bandpass', order=3, method='butterworth',
         way='filtfilt', axis=0):
    """Filt data.

    Parameters
    ----------
    sf : float
        The sampling frequency
    f : array_like
        Frequency vector (2,)
    x : array_like
        The data to filt.
    btype : {'bandpass', 'bandstop', 'highpass', 'lowpass'}
        If highpass, the first value of f will be used. If lowpass
        the second value of f will be used.
    order : int | 3
        The filter order.
    method : {'butterworth', 'bessel'}
        Filter type to use.
    way : {'filtfilt', 'lfilter'}
        Specify if the filter has to be one way ('lfilter') or two ways
        ('filtfilt').
    axis : int | 0
        The axis along which the filter is applied.

    Returns
    -------
    xfilt : array_like
        Filtered data.
    """
    # Normalize frequency vector according to btype :
    if btype in ['bandpass', 'bandstop']:
        fnorm = np.divide(f, .5 * sf)
    elif btype == 'lowpass':
        fnorm = np.array(f[-1] / (.5 * sf))
    elif btype == 'highpass':
        fnorm = np.array(f[0] / (.5 * sf))

    # Get filter coefficients :
    if method == 'butterworth':
        b, a = butter(order, fnorm, btype=btype)
    elif method == 'bessel':
        b, a = bessel(order, fnorm, btype=btype)

    # Apply filter :
    if way == 'filtfilt':
        return filtfilt(b, a, x, axis=axis)
    elif way == 'lfilter':
        return lfilter(b, a, x, axis=axis)

    def filter(self, order=1, cutoff=2000.0):
        """
        Performs a bessel filter on the selected data, normalizing the cutoff frequency by the 
        nyquist limit based on the sampling rate. 
        """

        if type(self) != Event:
            raise TypeError("Cannot filter a metaevent. Must have the current.")
        from scipy import signal

        nyquist = self.second / 2.0

        (b, a) = signal.bessel(order, cutoff / nyquist, btype="low", analog=0, output="ba")
        self.current = signal.filtfilt(b, a, self.current)
        self.filtered = True
        self.filter_order = order
        self.filter_cutoff = cutoff

def besselfilter(ECGdata):
    """ Filters the data using IIR bessel filter

        Description:
            Digital filter which returns the filtered signal using bessel
            4th order low pass design. The cutoff frequency is 0-35Hz with 100Hz
            as sampling frequency.
        Input:
            ECGdata -- list of integers (ECG data)
        Output:
            lfilter(b,a,ECGdata)-- filtered data along one-dimension with IIR
                                     bessel filter
    """        
    fs = 500.00
    f = 35.00
    N=4
    [b,a]=bessel(N,f/fs)
    return lfilter(b,a,ECGdata)

def bessel_bandpass_filter(data, lowcut, highcut, fs, order=2):
    """This is a wrapper for the bessel bandpass filter.
    :param: data - 1d numpy array to be filtered
    :param: lowcut - low pass frequency, in Hz
    :param: highcut - high pass frequency, in Hz
    :param: fs - sampling frequency, in samples / second (i.e.: 10000)
    :param: order - filter order
    :returns: filtered data
    """

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

    # bessel() and lfilter() are from scipy.signal

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

def bandpass(x, dt, f_lo, f_hi):
    """
    Highpass filter

    Parameters
    ----------
    x : stfio_plot.Timeseries
        Input data
    dt : float
        Sampling interval in ms
    f_c : float
        Cutoff frequency in kHz (-3 dB)

    Returns
    -------
    x convolved with a Gaussian filter kernel.
    """
    fs = 1.0/dt
    B, A = bessel(1, [f_lo / (fs / 2), f_hi / (fs / 2)], btype='bandpass') # 1st order Butterworth low-pass
    return lfilter(B, A, x, axis=0)

def lowpass(x, dt, f_c):
    """
    Lowpass filter

    Parameters
    ----------
    x : stfio_plot.Timeseries
        Input data
    dt : float
        Sampling interval in ms
    f_c : float
        Cutoff frequency in kHz (-3 dB)

    Returns
    -------
    x convolved with a Gaussian filter kernel.
    """
    fs = 1.0/dt
    cutoff = f_c
    B, A = bessel(1, cutoff / (fs / 2), btype='low') # 1st order Butterworth low-pass
    return lfilter(B, A, x, axis=0)

    def bessel(data, **kwargs):
        """
        Implement a Bessel type analog filter
        ----------
        data: xarray DataSet as obtained from a TDI call

        Returns
        -------
        Smoothed array and time derivative
        """
        time = data.time.values
        dt = (time.max() - time.min()) / (time.size - 1)
        Ny = numpy.round(
            1. / ((time.max() - time.min()) / (time.size - 1))) / 2
        # implement an appropriate Bessel analog filter
        fcutoff = kwargs.get('fcutoff', 30.)
        _Wn = fcutoff / Ny
        b, a = signal.bessel(4, _Wn)
        # create a copy of the signals
        sm = data.values.transpose().copy()
        smd = signal.filtfilt(b, a, numpy.gradient(sm, dt, axis=-1), axis=-1)
        return sm.transpose(), smd.transpose()

 def design(self, ripple=None):
     if self.filter_class == 'butterworth':
         self.B, self.A = signal.butter(self.N, self.Wn,
                                        self.filter_type, analog=True,
                                        output='ba')
     elif self.filter_class == 'chebyshev_1':
         if ripple is None or ripple <= 0:
             raise ValueError("Must give a ripple that is > 0")
         self.B, self.A = signal.cheby1(self.N, ripple, self.Wn,
                                        self.filter_type, analog=True,
                                        output='ba')
     elif self.filter_class == 'chebyshev_2':
         self.B, self.A = signal.cheby2(self.N, self.stopband_attenuation, self.Wn,
                                        self.filter_type, analog=True, output='ba')
     elif self.filter_class == 'elliptical':
         self.B, self.A = signal.ellip(self.N, self.passband_attenuation,
                                       self.stopband_attenuation, self.Wn,
                                       self.filter_type, analog=True, output='ba')
     elif self.filter_class == 'bessel':
         self.B, self.A = signal.bessel(self.N, self.Wn, self.filter_type, analog=True)
     else:
         raise NotImplementedError("Computation of {} not implemented yet.".format(self.filter_class))

    def _getFiltDesign(self, sf, f, npts):
        """Get the designed filter
        sf : sample frequency
        f : frequency vector/list [ex : f = [2,4]]
        npts : number of points
        """
        if type(f) != n.ndarray:
            f = n.array(f)
        if self.filtname == 'fir1':
            fOrder = fir_order(sf, npts, f[0], cycle=self.cycle)
            b, a = fir1(fOrder, f/(sf / 2))
        elif self.filtname == 'butter':
            b, a = butter(self.order, [(2*f[0])/sf,
                                       (2*f[1])/sf], btype='bandpass')
            fOrder = None
        elif self.filtname == 'bessel':
            b, a = bessel(self.order, [(2*f[0])/sf,
                                       (2*f[1])/sf], btype='bandpass')
            fOrder = None

        def filSignal(x): return filtfilt(b, a, x, padlen=fOrder,
                                          axis=self.axis)
        return filSignal