Python scipy.signal 模块，lfilter_zi() 实例源码

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 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 rasta_filt(x):
    number = np.arange(-2., 3., 1.)
    number = -1. * number / np.sum(number*number)
    denom = np.array([1., -0.94])
    zi = lfilter_zi(number, 1)
    zi = zi.reshape(1, len(zi))
    zi = np.repeat(zi, np.size(x, 0), 0)
    y, zf = lfilter(number, 1, x[:,0:4], axis=1, zi=zi)
    y, zf = lfilter(number, denom, x, axis=1, zi=zf)
    return y

def rasta_filt(x):
    number = np.arange(-2., 3., 1.)
    number = -1. * number / np.sum(number*number)
    denom = np.array([1., -0.94])
    zi = lfilter_zi(number, 1)
    zi = zi.reshape(1, len(zi))
    zi = np.repeat(zi, np.size(x, 0), 0)
    y, zf = lfilter(number, 1, x[:,0:4], axis=1, zi=zi)
    y, zf = lfilter(number, denom, x, axis=1, zi=zf)
    return y

def initZi(self):
        # lfilter_zi does not preserve dtype of arguments, bug that should  be reported XXX - idfah
        # if above was fixed, use don't need astype below
        self.zi = spsig.lfilter_zi(self.numCoef, self.denomCoef).astype(self.dtype, copy=False)

def filter(self, s, axis=0):
        """Filter new data.
        """
        if self.bandType == 'allpass':
            return s
        if self.bandType == 'allstop':
            return np.zeros_like(s)

        """ Should be very close to filtfilt, padding? XXX - idfah
        if self.zeroPhase:
            rev = [slice(None),]*s.ndim
            rev[axis] = slice(None,None,-1)

            #ziScaled = self.scaleZi(s[rev], axis)

            y, newZi = spsig.lfilter(self.numCoef, self.denomCoef, s[rev], axis=axis, zi=newZi)
            y = y[rev]
        """

        # if zeroPhase and signal is shorter than padlen (default in filtfilt function)
        if self.zeroPhase and \
            (3*max(len(self.numCoef), len(self.denomCoef))) < s.shape[axis]:

            # need astype below since filtfilt calls lfilter_zi, which does not preserve dtype XXX - idfah
            return spsig.filtfilt(self.numCoef, self.denomCoef,
                        s, axis=axis, padtype='even').astype(self.dtype, copy=False)

        else:
            ziScaled = self.scaleZi(s, axis)

            # even padding to help reduce edge effects
            nPad = 3*max(len(self.numCoef), len(self.denomCoef))
            sPad = np.apply_along_axis(np.pad, axis, s, pad_width=nPad, mode='reflect') # edge for constant padding
            slc = [slice(nPad,-nPad) if i == axis else slice(None) for i in range(s.ndim)]

            y, newZi = spsig.lfilter(self.numCoef, self.denomCoef, sPad, axis=axis, zi=ziScaled)

            return y[slc]