2017-08-30 77 views
1

鉴于一些数据x递归:用`IIR滤波器scipy.lfilter`

from pandas_datareader.data import DataReader as dr 
x = np.squeeze(dr('DTWEXB', 'fred').dropna().values) 

我想计算的另一种载体y如下:

enter image description here

阿尔法等于0.03 , 在这种情况下。

我可以用scipy.lfilter?来做到这一点。类似的问题here,但在这种情况下,结果的起始值是0,这是抛出一些东西。

我尝试:

from scipy.signal import lfilter 
a = 0.03 
b = 1 - a 
y0 = x[0] 
y = lfilter([a], [y0, -b], x) 

结果应该是:

true_y = np.empty(len(x)) 
for k in range(0, len(true_y)): 
    if k == 0: 
     true_y[k] = x[0] 
    else: 
     true_y[k] = a*x[k] + b*true_y[k-1] 
print(true_y) 
[ 101.1818  101.176862 101.16819314 ..., 120.9813121 120.92484874 
    120.85786628] 

回答

2

传递函数系数的正确参数是[a][1, -b]

处理您所需的初始条件,则可以使用scipy.signal.lfiltic创建过滤器的正确的初始状态:

zi = lfiltic([a], [1, -b], y=[x[0]]) 

然后调用lfilterzi说法:

y, zo = lfilter([a], [1, -b], x, zi=zi) 

这里有一个xy(使用lfilterzi计算)和您的true_y

In [37]: x 
Out[37]: array([ 3., 1., 2., 0., -1., 2.]) 

In [38]: y 
Out[38]: 
array([ 3.  , 2.94  , 2.9118 , 2.824446 , 2.70971262, 
     2.68842124]) 

In [39]: true_y 
Out[39]: 
array([ 3.  , 2.94  , 2.9118 , 2.824446 , 2.70971262, 
     2.68842124]) 
1

把你的滤波器的z变换,这为您提供分子b和分母a值:

 alpha 
y(z) = ------------------ x(z) 
     1 - (1-alpha) z^-1 

所以你运行

import scipy.signal 
x[0] /= alpha 
y = scipy.signal.lfilter([alpha], [1, - 1 + alpha], x) 

其中产量

array([ 101.1818 , 101.176862 , 101.16819314, ..., 120.9813121 , 
     120.92484874, 120.85786628]) 

注意我攀登x[0]考虑到你要的初始条件。

+0

似乎是在正确的轨道上,但有些事情是关闭的。结果'y'应该从上面等于'true_y'。我也从设置'a [0] = 0',即'[0,1-alpha]'('ValueError:BUG:filter coefficient a [0] == 0 yet supported')中得到错误。 –

+1

递归关系可以写成'y [k] - (1-a)* y [k-1] = a * x [k]'。取*表示*的z变换并求解'Y(z)'得到传递函数。 –

+0

我需要刷一下我的z转换技能看起来:-D –