使用 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)）。

但是，我注意到 sympy 计算系数的速度非常慢：即使是一个简单、简短的示例也需要 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 中的噪声。 p>