Python scipy.signal.filtfilt() Examples

def prepData(seqLocal = seq):
    dm = DataManager()
    dm.initHelper(dsName, subType, seqLocal)
    dt = dm.dt

    pSignal = dm.accdt_gnd
    pSignal = preClamp(pSignal)

    mSignal = dm.pr_dtr_gnd
    mSignal = preClamp((mSignal))

    mCov = dm.dtr_cov_gnd

    gtSignal = preClamp(dm.gt_dtr_gnd)
    gtSignal = filtfilt(gtSignal)
    return gtSignal, dt, pSignal, mSignal, mCov 

def filter_signal(b, a, signal):
    """
    Filter a signal.

    Simple wrapper around :func:`scipy.signal.filtfilt` to apply a
    foward-backward filter to preserve phase of the input. Requires the
    numerator and denominator polynomials from
    :func:`sensormotion.signal.build_filter`.

    Parameters
    ----------
    b : ndarray
        Numerator polynomial coefficients of the filter.
    a : ndarray
        Denominator polynomial coefficients of the filter.
    signal : ndarray
        Input array to be filtered.

    Returns
    -------
    signal_filtered : ndarray
        Filtered output of the original input signal.
    """

    return filtfilt(b, a, signal) 

def timefilter(
    xr_in, steps, step_spec, timename="time", filtertype="gaussian", stdev=0.1
):
    timedim = xr_in.dims.index(timename)
    dt = np.diff(xr_in.time.data[0:2])[0]
    cut_dt = np.timedelta64(steps, step_spec)

    if filtertype == "gaussian":
        win_length = (cut_dt / dt).astype(int)
        a = [1.0]
        win = gaussian(win_length, std=(float(win_length) * stdev))
        b = win / win.sum()
        if np.nansum(win) == 0:
            raise RuntimeError("window to short for time interval")
            print("win_length", str(win_length))
            print("stddev", str(stdev))
            print("win", str(win))

    filtered = filtfilt(b, a, xr_in.data, axis=timedim, padtype=None, padlen=0)
    out = xr.DataArray(
        filtered, dims=xr_in.dims, coords=xr_in.coords, attrs=xr_in.attrs
    )
    out.attrs.update({"filterlength": (steps, step_spec), "filtertype": filtertype})
    if xr_in.name:
        out.name = xr_in.name + "_lowpassed"
    return out 

def pb2bb(x, fs, fc, fd=None, flen=127, cutoff=None, axis=-1):
    """Convert passband signal to baseband.

    The baseband conversion uses a low-pass filter after downconversion, with a
    default cutoff frequency of `0.6*fd`, if `fd` is specified, or `1.1*fc` if `fd`
    is not specified. Alternatively, the user may specify the cutoff frequency
    explicitly.

    For communication applications, one may wish to use :func:`arlpy.comms.downconvert` instead,
    as that function supports matched filtering with a pulse shape rather than a generic
    low-pass filter.

    The convention used in that exp(2j*pi*fc*t) is a positive frequency carrier.

    :param x: passband signal
    :param fs: sampling rate of passband signal in Hz
    :param fc: carrier frequency in passband in Hz
    :param fd: sampling rate of baseband signal in Hz (``None`` => same as `fs`)
    :param flen: number of taps in the low-pass FIR filter
    :param cutoff: cutoff frequency in Hz (``None`` means auto-select)
    :param axis: axis of the signal, if multiple signals specified
    :returns: complex baseband signal, sampled at `fd`
    """
    if cutoff is None:
        cutoff = 0.6*fd if fd is not None else 1.1*_np.abs(fc)
    osc = _np.sqrt(2)*_np.exp(-2j*_np.pi*fc*time(x.shape[axis],fs))
    y = x * _utils.broadcastable_to(osc, x.shape, axis)
    hb = _sig.firwin(flen, cutoff=cutoff, nyq=fs/2.0)
    y = _sig.filtfilt(hb, 1, y, axis=axis)
    if fd is not None and fd != fs:
        y = _sig.resample_poly(y, 2*fd, fs, axis=axis)
        y = _np.apply_along_axis(lambda a: a[::2], axis, y)
    return y 

def filtfilt(data):
    y = np.zeros_like(data)
    b, a = signal.butter(8, 0.1)
    for i in range(0, 3):
        y[:, i] = signal.filtfilt(b, a, data[:, i], padlen=100)
    return y 

def add_disturbance(self, input):
        if self.options['input_disturbance'] is not None:
            fc = self.options['input_disturbance']['fc']
            stdev = self.options['input_disturbance']['stdev']
            if 'mean' in self.options['input_disturbance']:
                mean = self.options['input_disturbance']['mean']
            else:
                mean = np.zeros(stdev.shape)
            n_sign = input.shape[0]
            n_samp = input.shape[1]
            disturbance = np.zeros((n_sign, n_samp))
            filt = butter(3, fc, 'low')
            for k in range(n_sign):
                disturbance[k, :] = filtfilt(filt[0], filt[1],
                                             normal(mean[k], stdev[k], n_samp))
            return input + disturbance
        else:
            return input 

def butter_lowpass_filter(data, cutoff_freq, nyq_freq, order=4):
    # Build and apply filter to data (signal)
    b, a = butter_lowpass(cutoff_freq, nyq_freq, order=order)
    y = signal.filtfilt(b, a, data)
    return y 

def to_envelopes(path,num_bands,freq_lims,downsample=True):
    """Generate amplitude envelopes from a full path to a .wav, following
    Lewandowski (2012).

    Parameters
    ----------
    filename : str
        Full path to .wav file to process.
    num_bands : int
        Number of frequency bands to use.
    freq_lims : tuple
        Minimum and maximum frequencies in Hertz to use.
    downsample : bool
        If true, envelopes will be downsampled to 120 Hz

    Returns
    -------
    2D array
        Amplitude envelopes over time
    """
    sr, proc = preproc(path,alpha=0.97)

    #proc = proc / 32768 #hack!! for 16-bit pcm
    proc = proc/sqrt(mean(proc**2))*0.03;
    bandLo = [ freq_lims[0]*exp(log(freq_lims[1]/freq_lims[0])/num_bands)**x for x in range(num_bands)]
    bandHi = [ freq_lims[0]*exp(log(freq_lims[1]/freq_lims[0])/num_bands)**(x+1) for x in range(num_bands)]

    envelopes = []
    for i in range(num_bands):
        b, a = butter(2,(bandLo[i]/(sr/2),bandHi[i]/(sr/2)), btype = 'bandpass')
        env = filtfilt(b,a,proc)
        env = abs(hilbert(env))
        if downsample:
            env = resample(env,int(ceil(len(env)/int(ceil(sr/120)))))
            #env = decimate(env,int(ceil(sr/120)))
        envelopes.append(env)
    return array(envelopes).T 

def _butter_bandpass_filter(data, low_cut, high_cut, fs, axis = 0, order=5):
    '''Apply a bandpass butterworth filter with zero-phase filtering

    Args:
        data: (np.array)
        low_cut: (float) lower bound cutoff for high pass filter
        high_cut: (float) upper bound cutoff for low pass filter
        fs: (float) sampling frequency in Hz
        axis: (int) axis to perform filtering.
        order: (int) filter order for butterworth bandpass
    
    Returns:
        bandpass filtered data.
    '''
    nyq = 0.5 * fs
    b, a = butter(order, [low_cut/nyq, high_cut/nyq], btype='band')
    return filtfilt(b, a, data, axis=axis) 

def applyLowPassFilter(y, lowcut=12.0, order=9.0):
    if len(y) - 12 <= lowcut:
        lowcut = 3
    nyq = 0.5 * len(y)
    low = lowcut / nyq
    # high = highcut / nyq
    b, a = butter(order, low)
    m = min([len(y), len(a), len(b)])
    padlen = 30 if m >= 30 else m
    fl = filtfilt(b, a, y, padlen=padlen)
    return fl 

def attach_to_controller(self, controller):
        b, a = butter(3, 150.0 / controller.model.sample_rate * 2.0, 'high')

        @controller.raw_data_filter.add_filter
        def high_pass(arr, axis=0):
            return filtfilt(b, a, arr, axis=axis) 

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 low_pass_filter(data, cutoff, fs, n_taps=255):
    """Apply low-pass filter

    Parameters
    ----------
    data : array, shape (`T`, `dim`)
        Array of sequence.
    cutoff : int,
        Cutoff frequency
    fs : int,
        Sampling frequency
    n_taps : int, optional
        Tap number

    Returns
    -------
    modified data: array, shape (`T`, `dim`)
        Array of modified sequence.
    """
    if data.shape[0] < n_taps * 3:
        raise ValueError(
            'Length of data should be three times longer than n_taps.')

    fil = firwin(n_taps, cutoff, pass_zero=True, nyq=fs//2)
    modified_data = filtfilt(fil, 1, data, axis=0)
    return modified_data 

def demodulate(x, Fs, freq):
    """return decimated and demodulated audio signal envelope at a known CW frequency """
    t = np.arange(len(x))/ float(Fs)
    mixed =  x*((1 + np.sin(2*np.pi*freq*t))/2 )

    #calculate envelope and low pass filter this demodulated signal
    #filter bandwidth impacts decoding accuracy significantly 
    #for high SNR signals 40 Hz is better, for low SNR 20Hz is better
    # 25Hz is a compromise - could this be made an adaptive value?
    low_cutoff = 25. # 25 Hz cut-off for lowpass
    wn = low_cutoff/ (Fs/2.)    
    b, a = butter(3, wn)  # 3rd order butterworth filter
    z = filtfilt(b, a, abs(mixed))
    
    decimate = int(Fs/64) # 8000 Hz / 64 = 125 Hz => 8 msec / sample 
    Ts = 1000.*decimate/float(Fs)
    o = z[0::decimate]/max(z)
    return o 

def upscale_log(log, freq=20):
    """
    downscale a well log with a lowpass butterworth filter
    """
    depth = np.array(log.depth)
    data = np.array(log.data)
    mask = np.isfinite(data)
    func = interp1d(depth[mask], data[mask])
    interp_data = func(depth[log.start_idx: log.stop_idx])
    nyq = 10000 / 2
    dw = freq / nyq
    b, a = butter(4, dw, btype='low', analog=False)
    filtered = filtfilt(b, a, interp_data, method='gust')
    downscale_data = np.array(data)
    downscale_data[log.start_idx: log.stop_idx] = filtered
    log_downscale = Log()
    log_downscale.name = log.name + "_downscale_" + str(freq)
    log_downscale.units = log.units
    log_downscale.descr = log.descr
    log_downscale.depth = log.depth
    log_downscale.data = downscale_data
    return log_downscale 

def detrend(x, dt=0.02, stop_hz=0.01, order=5):
    orig_shape = x.shape
    x = np.atleast_2d(x)
    nyquist = 0.5 / dt
    stop = stop_hz / nyquist
    b, a = signal.butter(order, Wn=stop, btype='lowpass')
    y = signal.filtfilt(b, a, x, axis=1)
    return (x - y).reshape(orig_shape) 

def apply(data):
        return scipy.signal.filtfilt(self.b, self.a, data, axis=-1, padtype='odd', padlen=None) 

def apply_filter(x, filter=None, axis=0):
    """Apply a filter to an array."""
    if isinstance(x, list):
        x = np.asarray(x)
    if x.shape[axis] == 0:
        return x
    b, a = filter
    return signal.filtfilt(b, a, x[:], axis=axis) 

def _decimate(x, q):
    """
    Downsample the signal after low-pass filtering to avoid aliasing.
    An order 16 Chebyshev type I filter is used.

    Parameters
    ----------
    x : ndarray
        The signal to be downsampled, as an N-dimensional array.
    q : int
        The downsampling factor.

    Returns
    -------
    y : ndarray
        The down-sampled signal.
    """
    if not isinstance(q, int):
        raise TypeError("q must be an integer")

    b, a = signal.filter_design.cheby1(16, 0.025, 0.98 / q)

    y = signal.filtfilt(b, a, x, axis=-1)

    sl = [slice(None)] * y.ndim
    sl[-1] = slice(None, None, q)
    return y[sl] 

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 _butterworth(ts, low_frequency, high_factor, order, sampling_frequency):
    """Butterworth filter

    Parameters
    ----------
    ts: np.array
        T  numpy array, where T is the number of time samples
    low_frequency: int
        Low pass frequency (Hz)
    high_factor: float
        High pass factor (proportion of sampling rate)
    order: int
        Order of Butterworth filter
    sampling_frequency: int
        Sampling frequency (Hz)

    Notes
    -----
    This function can only be applied to a one dimensional array, to apply
    it to multiple channels use butterworth

    Raises
    ------
    NotImplementedError
        If a multidmensional array is passed
    """

    low = float(low_frequency) / sampling_frequency * 2
    high = float(high_factor) * 2
    b, a = butter(order, low, btype='high', analog=False)

    if ts.ndim == 1:
        return filtfilt(b, a, ts)
    else:
        T, C = ts.shape
        output = np.zeros((T, C), 'float32')
        for c in range(C):
            output[:, c] = filtfilt(b, a, ts[:, c])

        return output 

def test_basic(self):
        zpk = tf2zpk([1, 2, 3], [1, 2, 3])
        out = self.filtfilt(zpk, np.arange(12))
        assert_allclose(out, arange(12), atol=1e-11) 

def _high_frequency_completion(self, x, transformed):
        """
        Please see Sect. 3.2 and 3.3 in the following paper to know why we complete the
        unvoiced synthesized voice of the original voice into high frequency range
        of F0 transformed voice.

        - K. Kobayashi et al., "F0 transformation techniques for statistical voice
        conversion with direct waveform modification with spectral differential,"
        Proc. IEEE SLT 2016, pp. 693-700. 2016.
        """
        # construct feature extractor and synthesis
        feat = FeatureExtractor(fs=self.fs)
        f0, spc, ap = feat.analyze(x)
        uf0 = np.zeros(len(f0))

        # synthesis
        synth = Synthesizer(fs=self.fs)
        unvoice_anasyn = synth.synthesis_spc(uf0, spc, ap)

        # HPF for synthesized speech
        fil = firwin(255, self.f0rate, pass_zero=False)
        HPFed_unvoice_anasyn = filtfilt(fil, 1, unvoice_anasyn)

        if len(HPFed_unvoice_anasyn) > len(transformed):
            return transformed + HPFed_unvoice_anasyn[:len(transformed)]
        else:
            transformed[:len(HPFed_unvoice_anasyn)] += HPFed_unvoice_anasyn
            return transformed 

def high_pass_filter(data, cutoff, fs, n_taps=255):
    """Apply high-pass filter

    Parameters
    ----------
    data : array, shape (`T`, `dim`)
        Array of sequence.
    cutoff : int,
        Cutoff frequency
    fs : int,
        Sampling frequency
    n_taps : int, optional
        Tap number

    Returns
    -------
    modified data: array, shape (`T`, `dim`)
        Array of modified sequence.
    """
    if data.shape[0] < n_taps * 3:
        raise ValueError(
            'Length of data should be three times longer than n_taps.')

    fil = firwin(n_taps, cutoff, pass_zero=False, nyq=fs//2)
    modified_data = filtfilt(fil, 1, data, axis=0)
    return modified_data 

def _filter_obliquity(OBL, F, Kx, vel, critical, ntaper, Ky=0):
    """Apply masking of ``OBL`` based on critical angle and tapering at edges

    Parameters
    ----------
    OBL : :obj:`np.ndarray`
        Obliquity factor
    F : :obj:`np.ndarray`
        Frequency grid
    Kx : :obj:`np.ndarray`
        Horizonal wavenumber grid
    vel : :obj:`float`
        Velocity along the receiver array (must be constant)
    critical : :obj:`float`, optional
        Percentage of angles to retain in obliquity factor
    ntaper : :obj:`float`, optional
        Number of samples of taper applied to obliquity factor around critical
        angle
    Ky : :obj:`np.ndarray`, optional
        Second horizonal wavenumber grid

    Returns
    -------
    OBL : :obj:`np.ndarray`
        Filtered obliquity factor

    """
    critical /= 100.
    mask = np.sqrt(Kx**2 + Ky**2) < critical * np.abs(F) / vel
    OBL *= mask
    OBL = filtfilt(np.ones(ntaper) / float(ntaper), 1, OBL, axis=0)
    OBL = filtfilt(np.ones(ntaper) / float(ntaper), 1, OBL, axis=1)
    if isinstance(Ky, np.ndarray):
        OBL = filtfilt(np.ones(ntaper) / float(ntaper), 1, OBL, axis=2)
    return OBL 

def test_basic(self):
        out = signal.filtfilt([1, 2, 3], [1, 2, 3], np.arange(12))
        assert_equal(out, arange(12)) 

def test_sine(self):
        rate = 2000
        t = np.linspace(0, 1.0, rate + 1)
        # A signal with low frequency and a high frequency.
        xlow = np.sin(5 * 2 * np.pi * t)
        xhigh = np.sin(250 * 2 * np.pi * t)
        x = xlow + xhigh

        b, a = butter(8, 0.125)
        z, p, k = tf2zpk(b, a)
        # r is the magnitude of the largest pole.
        r = np.abs(p).max()
        eps = 1e-5
        # n estimates the number of steps for the
        # transient to decay by a factor of eps.
        n = int(np.ceil(np.log(eps) / np.log(r)))

        # High order lowpass filter...
        y = filtfilt(b, a, x, padlen=n)
        # Result should be just xlow.
        err = np.abs(y - xlow).max()
        assert_(err < 1e-4)

        # A 2D case.
        x2d = np.vstack([xlow, xlow + xhigh])
        y2d = filtfilt(b, a, x2d, padlen=n, axis=1)
        assert_equal(y2d.shape, x2d.shape)
        err = np.abs(y2d - xlow).max()
        assert_(err < 1e-4)

        # Use the previous result to check the use of the axis keyword.
        # (Regression test for ticket #1620)
        y2dt = filtfilt(b, a, x2d.T, padlen=n, axis=0)
        assert_equal(y2d, y2dt.T) 

def test_axis(self):
        # Test the 'axis' keyword on a 3D array.
        x = np.arange(10.0 * 11.0 * 12.0).reshape(10, 11, 12)
        b, a = butter(3, 0.125)
        y0 = filtfilt(b, a, x, padlen=0, axis=0)
        y1 = filtfilt(b, a, np.swapaxes(x, 0, 1), padlen=0, axis=1)
        assert_array_equal(y0, np.swapaxes(y1, 0, 1))
        y2 = filtfilt(b, a, np.swapaxes(x, 0, 2), padlen=0, axis=2)
        assert_array_equal(y0, np.swapaxes(y2, 0, 2)) 

def test_poly_vs_filtfilt(self):
        # Check that up=1.0 gives same answer as filtfilt + slicing
        random_state = np.random.RandomState(17)
        try_types = (int, np.float32, np.complex64, float, complex)
        size = 10000
        down_factors = [2, 11, 79]

        for dtype in try_types:
            x = random_state.randn(size).astype(dtype)
            if dtype in (np.complex64, np.complex128):
                x += 1j * random_state.randn(size)

            # resample_poly assumes zeros outside of signl, whereas filtfilt
            # can only constant-pad. Make them equivalent:
            x[0] = 0
            x[-1] = 0

            for down in down_factors:
                h = signal.firwin(31, 1. / down, window='hamming')
                yf = filtfilt(h, 1.0, x, padtype='constant')[::down]

                # Need to pass convolved version of filter to resample_poly,
                # since filtfilt does forward and backward, but resample_poly
                # only goes forward
                hc = convolve(h, h[::-1])
                y = signal.resample_poly(x, 1, down, window=hc)
                assert_allclose(yf, y, atol=1e-7, rtol=1e-7) 

def filtfilt(self, zpk, x, axis=-1, padtype='odd', padlen=None,
                 method='pad', irlen=None):
        if self.filtfilt_kind == 'tf':
            b, a = zpk2tf(*zpk)
            return filtfilt(b, a, x, axis, padtype, padlen, method, irlen)
        elif self.filtfilt_kind == 'sos':
            sos = zpk2sos(*zpk)
            return sosfiltfilt(sos, x, axis, padtype, padlen)