инверсной фильтрации с использованием Python

Question

Учитывая импульсную характеристику h и выход y (как одномерные массивы), я пытаюсь найти способ вычисления обратного фильтра x таким образом, что h * x = y, где * обозначает свертку продукт.

Например, предположим, что импульсная характеристика h равна [1,0.5], а выход представляет собой ступенчатую функцию (то есть состоит из всех 1). Можно выяснить, что первые коэффициенты должны быть [1, 0.5, 0.75], что дает выход [1, 1, 1, 0.375]. Последний термин содержит ошибку, но это не проблема, так как я занимаюсь только выходом до определенного максимального времени.

Я хотел бы автоматизировать и «масштабировать» эту обратную фильтрацию до более длинной, более сложной функции импульсного отклика. До сих пор единственный способ, по которому я выяснил, получить коэффициенты, - это генерировать серийное расширение Z-преобразования с помощью sympy. (Заметим, что 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

Это называется deconvolution: 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 до машинной точности-it'll представляет собой шум в записанном сигнал y.