使用Python进行反向过滤

Question

给定的脉冲响应h和输出y(均为一维阵列),我试图找到一种方法来计算逆滤波器x,使得h * x = y,其中*表示卷积乘积.

例如,假设脉冲响应h是[1,0.5],并且输出是阶跃函数(即,由所有1s 组成).可以看出第一个系数应该是[1, 0.5, 0.75],产生一个输出[1, 1, 1, 0.375].最后一个术语包含错误,但这不是一个问题,因为我只关心输出到一定的最大时间.

我想将这种反向滤波自动化并"扩展"为更长,更复杂的脉冲响应函数.到目前为止,我已经想出获得系数的唯一方法是使用sympy生成Z变换的级数展开.(注意,阶跃函数的Z变换是1 /(1-z)).

但是,我注意到,在计算系数时,同情心很慢:即使是一个简单的简短示例,也需要0.8秒,如下面的脚本所示:

import numpy as np
from scipy import signal
from sympy import Symbol, series
import time

h = np.array([1,0.5])           # Impulse response function
y = np.array([1,1,1,1,1])       # Ideal output is a step function

my_x = np.array([1,0.5,0.75])   # Inverse filter response (calculated manually)

my_y = signal.convolve(h,my_x)  # Actual output (should be close to ideal output)

print(my_y)

start = time.time()
z = Symbol('z')
print(series(1/((1-z)*(1+z/2)),z))
end = time.time()
print("The power series expansion took "+str(end - start)+" seconds to complete.")

这会产生输出:

[ 1.     1.     1.     0.375]
1 + z/2 + 3*z**2/4 + 5*z**3/8 + 11*z**4/16 + 21*z**5/32 + O(z**6)
The power series expansion took 0.798881053925 seconds to complete.

简而言之,功率扩展的系数与所需的滤波器响应相匹配,但使用sympy来做这件事似乎很麻烦.有没有更好的方法来计算Python中的逆滤波器系数？

Answer 1

这称为解卷积:scipy.signal.deconvolve将为您完成此操作.您知道原始输入信号的示例x:

import numpy as np
import scipy.signal as signal

# We want to try and recover x from y, where y is made like this:
x = np.array([0.5, 2.2, -1.8, -0.1])
h = np.array([1.0, 0.5])
y = signal.convolve(x, h)

# This is called deconvolution:
xRecovered, xRemainder = signal.deconvolve(y, h)
# >>> xRecovered
# array([ 0.5,  2.2, -1.8, -0.1])
# >>> xRemainder
# array([  0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
#          0.00000000e+00,  -1.38777878e-17])

# Check the answer
assert(np.allclose(xRecovered, x))
assert(np.allclose(xRemainder, 0))

这听起来像你不知道原来的信号,所以你xRemainder会不会为0到机床精度,稍后即代表了记录信号中的噪声y.