Python scipy.signal.lfilter_zi() Examples

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 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 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 find_peak_vextors(price, return_ref=False, offest=0):
    """
    采用巴特沃斯信号滤波器，自适应寻找最佳极值点，决定平均周期分段数量。
    使用 scipy.Gaussian 机器学习统计算法进行第二次分析
    If you meet a Warning message, To slove this need upgrade scipy=>1.2. 
    but QUANTAXIS incompatible scipy=>1.2

    Parameters
    ----------
    price : (N,) array_like
        传入需要查找极值点的价格-时间序列。
        The numerator coefficient vector of the filter.
    return_ref : bool or None, optional
        返回作为参照的平滑曲线，平滑曲线的目的是减少锯齿抖动，减少计算的极值点。
        Return the smoothed line for reference.
    offest : int or None, optional
        传递参数时可能会被 .dropna() 或者 price[t:0] 等切片手段移除 nparray 头部
        的 np.nan 元素，因为此函数返回的是向量节点的数组索引，为了跟原始参数对应，调用者
        可以指定一个补偿偏移量，在返回的最大最小值中每个索引都会追加这个偏移量。
        The number of elements index offest, for jump np.nan in price's head.

    Returns
    -------
    x_tp_min, x_tp_max : ndarray
        包含最大值和最少值索引的数组
        The min/max peakpoint's index in array.

    """
    xn = price

    # Create an order 3 lowpass butterworth filter.
    b, a = butter(3, 0.05)

    # Apply the filter to xn.  Use lfilter_zi to choose the initial condition
    # of the filter.
    zi = lfilter_zi(b, a)
    z, _ = lfilter(b, a, xn, zi=zi * xn[0])

    # Apply the filter again, to have a result filtered at an order
    # the same as filtfilt.
    z2, _ = lfilter(b, a, z, zi=zi * z[0])

    # Use filtfilt to apply the filter.  If you meet a Warning need upgrade to
    # scipy=>1.2 but QUANTAXIS incompatible scipy=>1.2
    y = filtfilt(b, a, xn)

    # pass 1
    x_tp_min, x_tp_max = signal.argrelextrema(y, np.less)[0], signal.argrelextrema(y, np.greater)[0]
    n = int(len(price) / (len(x_tp_min) + len(x_tp_max))) * 2

    # peakutils 似乎一根筋只能查最大极值，通过曲线反相的方式查找极小点
    mirrors = (price * -1) + np.mean(price) * 2

    # pass 2 使用 peakutils 查找
    x_tp_max = peakutils.indexes(price, thres=0.01 / max(price), min_dist=n)
    x_tp_min = peakutils.indexes(mirrors, thres=0.01 / max(price), min_dist=n)

    if (return_ref):
        return x_tp_min + offest, x_tp_max + offest, y
    else:
        return x_tp_min + offest, x_tp_max + offest