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

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

Например, предположим, что импульсный отклик h равен [1,0.5] , а выход – ступенчатой ​​(т. [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?

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