Python scipy.signal.lfilter() Examples

def discount(x, gamma):
  """
  Compute discounted sum of future values
  out[i] = in[i] + gamma * in[i+1] + gamma^2 * in[i+2] + ...
  """
  return signal.lfilter([1], [1, -gamma], x[::-1], axis=0)[::-1] 

def preemphasis(wav, coeff=0.97):
    """
    Emphasize high frequency range of the waveform by increasing power(squared amplitude).

    Parameters
    ----------
    wav : np.ndarray [shape=(n,)]
        Real-valued the waveform.

    coeff: float <= 1 [scalar]
        Coefficient of pre-emphasis.

    Returns
    -------
    preem_wav : np.ndarray [shape=(n,)]
        The pre-emphasized waveform.
    """
    preem_wav = signal.lfilter([1, -coeff], [1], wav)
    return preem_wav 

def inv_preemphasis(preem_wav, coeff=0.97):
    """
    Invert the pre-emphasized waveform to the original waveform.
    
    Parameters
    ----------
    preem_wav : np.ndarray [shape=(n,)]
        The pre-emphasized waveform.

    coeff: float <= 1 [scalar]
        Coefficient of pre-emphasis.

    Returns
    -------
    wav : np.ndarray [shape=(n,)]
        Real-valued the waveform.
    """
    wav = signal.lfilter([1], [1, -coeff], preem_wav)
    return wav 

def spectrogram2wav(mag):
    '''# Generate wave file from spectrogram'''
    # transpose
    mag = mag.T

    # de-noramlize
    mag = (np.clip(mag, 0, 1) * hp.max_db) - hp.max_db + hp.ref_db

    # to amplitude
    mag = np.power(10.0, mag * 0.05)

    # wav reconstruction
    wav = griffin_lim(mag)

    # de-preemphasis
    wav = signal.lfilter([1], [1, -hp.preemphasis], wav)

    # trim
    wav, _ = librosa.effects.trim(wav)

    return wav.astype(np.float32) 

def discount_cumsum(x, discount):
    """
    Forked from rllab for computing discounted cumulative sums of vectors.

    :param x (np.ndarray or tf.Tensor)
        vector of [x0, x1, x2]
    :return output:
        [x0 + discount * x1 + discount^2 * x2,  
         x1 + discount * x2,
         x2]
    """
    return lfilter(
        b=[1],
        a=[1, float(-discount)],
        x=x[::-1],
        axis=0)[::-1] 

def preemphasis(wav, coeff=0.97):
    """
    Emphasize high frequency range of the waveform by increasing power(squared amplitude).
    
    Parameters
    ----------
    wav : np.ndarray [shape=(n,)]
        Real-valued the waveform.

    coeff: float <= 1 [scalar]
        Coefficient of pre-emphasis.

    Returns
    -------
    preem_wav : np.ndarray [shape=(n,)]
        The pre-emphasized waveform.
    """
    preem_wav = signal.lfilter([1, -coeff], [1], wav)
    return preem_wav 

def compressor(x, thresh=-24, ratio=2, attackrel=0.045, sr=44100.0, dtype=np.float32):
    """
    simple compressor effect, code thanks to Eric Tarr @hackaudio
    Inputs:
       x:        the input waveform
       thresh:   threshold in dB
       ratio:    compression ratio
       attackrel:   attack & release time in seconds
       sr:       sample rate
    """
    attack = attackrel * sr  # convert to samples
    fc = 1.0/float(attack)     # this is like 1/attack time
    b, a = scipy_signal.butter(1, fc, analog=False, output='ba')
    zi = scipy_signal.lfilter_zi(b, a)

    dB = 20. * np.log10(np.abs(x) + 1e-6)
    in_env, _ = scipy_signal.lfilter(b, a, dB, zi=zi*dB[0])  # input envelope calculation
    out_env = np.copy(in_env)              # output envelope
    i = np.where(in_env >  thresh)          # compress where input env exceeds thresh
    out_env[i] = thresh + (in_env[i]-thresh)/ratio
    gain = np.power(10.0,(out_env-in_env)/20)
    y = x * gain
    return y 

def load_series(self, subject, series):
        min_freq = self.min_freq
        max_freq = self.max_freq
        key = (subject, series, min_freq, max_freq) 
        if key not in self._series_cache:
            while len(self._series_cache) > self.MAX_CACHE_SIZE:
                # Randomly throw away an item
                self._series_cache.popitem()    
            print("Loading", subject, series)
            data = read_csv(path(subject, series, "data"))
            # Filter here since it's slow and we don't want to filter multiple
            # times. `lfilter` is CAUSAL and thus doesn't violate the ban on future data.
            if (self.min_freq is None) or (self.min_freq == 0):
                print("Low pass filtering, f_h =", max_freq)
                b, a = butter_lowpass(max_freq, SAMPLE_RATE, FILTER_N)
            else:
                print("Band pass filtering, f_l =", min_freq, "f_h =", max_freq)
                b, a = butter_bandpass(min_freq, max_freq, SAMPLE_RATE, FILTER_N)                
            self._series_cache[key] = lfilter(b, a, data, axis=0)
        return self._series_cache[key] 

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 geterrors(self, params):
        #copied from sandbox.tsa.arima.ARIMA
        p, q = self.nar, self.nma
        ar = np.concatenate(([1], -params[:p]))
        ma = np.concatenate(([1], params[p:p+q]))

        #lfilter_zi requires same length for ar and ma
        maxlag = 1+max(p,q)
        armax = np.zeros(maxlag)
        armax[:p+1] = ar
        mamax = np.zeros(maxlag)
        mamax[:q+1] = ma
        #remove zi again to match better with Skipper's version
        #zi = signal.lfilter_zi(armax, mamax)
        #errorsest = signal.lfilter(rhoy, rhoe, self.endog, zi=zi)[0] #zi is also returned
        errorsest = signal.lfilter(ar, ma, self.endog)
        return errorsest 

def inv_preemphasis(preem_wav, coeff=0.97):
    """
    Invert the pre-emphasized waveform to the original waveform.

    Parameters
    ----------
    preem_wav : np.ndarray [shape=(n,)]
        The pre-emphasized waveform.

    coeff: float <= 1 [scalar]
        Coefficient of pre-emphasis.

    Returns
    -------
    wav : np.ndarray [shape=(n,)]
        Real-valued the waveform.
    """
    wav = signal.lfilter([1], [1, -coeff], preem_wav)
    return wav 

def generate_garch(nobs, ar, ma, mu=1., scale=0.1):
    '''simulate standard garch

    scale : float
       scale/standard deviation of innovation process in GARCH process
    '''

    eta = scale*np.random.randn(nobs)
    # copied from armageneratesample
    h = signal.lfilter(ma, ar, eta**2)

    #
    #h = (mu+h)**2
    #h = np.abs(h)
    #h = np.exp(h)
    #err = np.sqrt(h)*np.random.randn(nobs)
    err = np.sqrt(h)*eta #np.random.standard_t(8, size=nobs)
    return err, h 

def butter_filter(sig, fs, ftype="low", low_cut=50, high_cut=2000, order=5):
    """
    Apply filter to signal.

    Args:
        - sig      (array) : the signal array to filter.
        - fs       (float) : the sampling rate.
        - ftype      (str) : the filter type, by default defined to a low pass filter
        - low_cut  (float) : the low cutoff frequency, by default defined to  50Hz
        - high_cut (float) : the high cutoff frequency, by default defined to 2000Hz.
        - order      (int) : order of the filter, by default defined to 5.

    Returns:
        array of the filtered signal.
    """
    if   ftype == "band" : b, a = butter_bandpass(low_cut, high_cut, fs, order)
    elif ftype == "high" : b, a = butter_highpass(high_cut, fs, order)
    else                 : b, a = butter_lowpass(low_cut,  fs, order)

    # filter signal
    y = lfilter(b, a, sig)
    return y 

def filter(self, x):
        '''
        filter a timeseries with the ARMA filter

        padding with zero is missing, in example I needed the padding to get
        initial conditions identical to direct filter

        Initial filtered observations differ from filter2 and signal.lfilter, but
        at end they are the same.

        See Also
        --------
        tsa.filters.fftconvolve

        '''
        n = x.shape[0]
        if n == self.fftarma:
            fftarma = self.fftarma
        else:
            fftarma = self.fftma(n) / self.fftar(n)
        tmpfft = fftarma * fft.fft(x)
        return fft.ifft(tmpfft) 

def filter2(self, x, pad=0):
        '''filter a time series using fftconvolve3 with ARMA filter

        padding of x currently works only if x is 1d
        in example it produces same observations at beginning as lfilter even
        without padding.

        TODO: this returns 1 additional observation at the end
        '''
        from statsmodels.tsa.filters import fftconvolve3
        if not pad:
            pass
        elif pad == 'auto':
            #just guessing how much padding
            x = self.padarr(x, x.shape[0] + 2*(self.nma+self.nar), atend=False)
        else:
            x = self.padarr(x, x.shape[0] + int(pad), atend=False)

        return fftconvolve3(x, self.ma, self.ar) 

def spectrogram2wav(mag):
    '''# Generate wave file from spectrogram'''
    # transpose
    mag = mag.T

    # de-noramlize
    mag = (np.clip(mag, 0, 1) * hp.max_db) - hp.max_db + hp.ref_db

    # to amplitude
    mag = np.power(10.0, mag * 0.05)

    # wav reconstruction
    wav = griffin_lim(mag)

    # de-preemphasis
    wav = signal.lfilter([1], [1, -hp.preemphasis], wav)

    # trim
    wav, _ = librosa.effects.trim(wav)

    return wav.astype(np.float32) 

def exactprocess(self, xzero, nobs, ddt=1., nrepl=2):
        '''ddt : discrete delta t



        should be the same as an AR(1)
        not tested yet
        '''
        t = np.linspace(ddt, nobs*ddt, nobs)
        #expnt = np.exp(-self.lambd * t)
        expddt = np.exp(-self.lambd * ddt)
        normrvs = np.random.normal(size=(nrepl,nobs))
        #do I need lfilter here AR(1) ? if mean reverting lag-coeff<1
        #lfilter doesn't handle 2d arrays, it does?
        inc = self._exactconst(expddt) + self._exactstd(expddt) * normrvs
        return signal.lfilter([1.], [1.,-expddt], inc) 

def exactprocess(self, xzero, nobs, ddt=1., nrepl=2):
        '''ddt : discrete delta t

        should be the same as an AR(1)
        not tested yet
        # after writing this I saw the same use of lfilter in sitmo
        '''
        t = np.linspace(ddt, nobs*ddt, nobs)
        expnt = np.exp(-self.lambd * t)
        expddt = np.exp(-self.lambd * ddt)
        normrvs = np.random.normal(size=(nrepl,nobs))
        #do I need lfilter here AR(1) ? lfilter doesn't handle 2d arrays, it does?
        from scipy import signal
        #xzero * expnt
        inc = ( self.mu * (1-expddt) +
                self.sigma * np.sqrt((1-expddt*expddt)/2./self.lambd) * normrvs )

        return signal.lfilter([1.], [1.,-expddt], inc) 

def impz(b,a):
    """Pseudo implementation of the impz method of MATLAB"""
#% Compute time vector
# M = 0;  NN = [];
# if isempty(N)
#   % if not specified, determine the length
#   if isTF
#     N = impzlength(b,a,.00005);
#   else
#     N  = impzlength(b,.00005);
#   end
    p = np.roots(a)
    N = stableNmarginal_length(p, 0.00005, 0)
    N = len(b) * len(b) * len(b) # MATLAB AUTOFINDS THE SIZE HERE... 
    #TODO: Implement some way of finding the autosieze of this... I used a couple of examples... matlab gave 43 as length we give 64
    x = zeros(N)
    x[0] = 1
    h = lfilter(b,a, x)
    return h 

def bandpass_filter(data, low, high, fs, order=5):
    """
    Does a bandpass filter over the given data.

    :param data: The data (numpy array) to be filtered.
    :param low: The low cutoff in Hz.
    :param high: The high cutoff in Hz.
    :param fs: The sample rate (in Hz) of the data.
    :param order: The order of the filter. The higher the order, the tighter the roll-off.
    :returns: Filtered data (numpy array).
    """
    if not scipy_imported:
        raise NotImplementedError("This function is unusable without Scipy")

    nyq = 0.5 * fs
    low = low / nyq
    high = high / nyq
    b, a = signal.butter(order, [low, high], btype='band')
    y = signal.lfilter(b, a, data)
    return y 

def lowpass_filter(data, cutoff, fs, order=5):
    """
    Does a lowpass filter over the given data.

    :param data: The data (numpy array) to be filtered.
    :param cutoff: The high cutoff in Hz.
    :param fs: The sample rate in Hz of the data.
    :param order: The order of the filter. The higher the order, the tighter the roll-off.
    :returns: Filtered data (numpy array).
    """
    if not scipy_imported:
        raise NotImplementedError("This function is unusable without Scipy")

    nyq = 0.5 * fs
    normal_cutoff = cutoff / nyq
    b, a = signal.butter(order, normal_cutoff, btype='low', analog=False)
    y = signal.lfilter(b, a, data)
    return y 

def __call__(self, x, axis=-1):
        """Apply filter to signal.

        Parameters
        ----------
        x : 1d array
        axis : int
        """
        return lfilter(self.taps, 1.0, x, axis=axis) 

def preemp(input, p):
    """Pre-emphasis filter."""
    return lfilter([1., -p], 1, input) 

def integrated_water_vapour_synthesis(self, lat, lon, Ns, Ts=1, n=None):
        # A Estimation of κ and λ
        kappa, lambd = self.integrated_water_vapour_coefficients(lat, lon)

        # B Low-pass filter parameter
        beta_V = 3.24e-6

        # Step C: Time series synthesis
        # Step C1: Synthesize a white Gaussian noise time series
        if n is None:
            n = np.random.normal(0, 1, (Ns * Ts + 5000000))[::Ts]
            discard_samples = True
        else:
            discard_samples = False

        # Step C3: Filter the noise time series, with two recursive low-pass
        # filters
        rho = np.exp(-beta_V * Ts)
        G_v = lfilter([np.sqrt(1 - rho**2)], [1, -rho], n, 0)
        # Step C4: Compute Compute V(kTs),
        V = lambd * (- np.log10(stats.norm.sf(G_v)))**(1 / kappa)
        # Step C5: Discard the first 5 000 000 samples from the synthesized
        if discard_samples:
            V = V[np.ceil(5000000/Ts).astype(int):]

        return V.flatten() 

def discount(x, gamma):
  """
  Compute discounted sum of future values
  out[i] = in[i] + gamma * in[i+1] + gamma^2 * in[i+2] + ...
  """
  return signal.lfilter([1], [1, -gamma], x[::-1], axis=0)[::-1] 

def inv_preemphasis(x):
	return signal.lfilter([1], [1, -hparams.preemphasis], x) 

def preemphasis(x):
	return signal.lfilter([1, -hparams.preemphasis], [1], x) 

def inv_preemphasis(x):
	return signal.lfilter([1], [1, -hparams.preemphasis], x) 

def preemphasis(x):
	return signal.lfilter([1, -hparams.preemphasis], [1], x) 

def spectrogram2wav(mag):
    '''# Generate wave file from linear magnitude spectrogram

    Args:
      mag: A numpy array of (T, 1+n_fft//2)

    Returns:
      wav: A 1-D numpy array.
    '''
    # transpose
    mag = mag.T

    # de-noramlize
    mag = (np.clip(mag, 0, 1) * max_db) - max_db + ref_db

    # to amplitude
    mag = np.power(10.0, mag * 0.05)

    # wav reconstruction
    wav = griffin_lim(mag**power)

    # de-preemphasis
    wav = signal.lfilter([1], [1, -preemphasis], wav)

    # trim
    wav, _ = librosa.effects.trim(wav)

    return wav.astype(np.float32)