Вопрос о деконволюции сигнала с использованием Python scipy

Я пытаюсь немного изучить обработку сигналов, в частности, используя Python. Вот пример кода, который я написал.

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import deconvolve

a = np.linspace(-1,1,50)
b = np.linspace(-1,1,50)**2
c = np.convolve(a,b,mode='same')
quotient,remainder = deconvolve(c,b);
plt.plot(a/max(a),"g")
plt.plot(b/max(b),"r")
plt.plot(c/max(c),"b")
plt.plot(remainder/max(remainder),"k")
#plt.plot(quotient/max(quotient),"k")

plt.legend(['a_original','b_original','convolution_a_b','deconvolution_a_b'])

Насколько я понимаю, деконволюция запутанного массива должна возвращать точно такой же массив «a», поскольку я использую «b» в качестве фильтра. Это явно не так, как видно из графиков ниже.

Я не совсем уверен, неправильно ли мое математическое понимание деконволюции или что-то не так с кодом. Любая помощь приветствуется!

введите здесь описание изображения


person CasualPythoner    schedule 16.10.2020    source источник


Ответы (1)


Вы используете mode='same', и это кажется несовместимым с scipy deconvolve. Попробуйте с mode='full', должно работать намного лучше.

Вот исправленный код:

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import deconvolve

a = np.linspace(-1,1,50)
b = np.linspace(-1,1,50)**2
c = np.convolve(a,b,mode='full')
quotient,remainder = deconvolve(c,b)
plt.plot(a,"g")
plt.plot(b,"r")
plt.plot(c,"b")
plt.plot(quotient,"k")
plt.xlim(0,50)
plt.ylim(-6,2)

plt.legend(['a_original','b_original','convolution_a_b','deconvolution_c_b'])

person pygri    schedule 16.10.2020
comment
Ваш код отлично подходит для этого случая. Однако фактические данные, которые я получу и с которыми мне придется работать, не будут работать с режимом = 'полный'. Любые идеи, почему mode = 'same' является проблемой? Я действительно смущен этим. - person CasualPythoner; 16.10.2020
comment
это не то, что это проблема как таковая, просто она не поддерживается функцией деконволюции scipy - если это такая большая проблема для вас, вы должны развернуть свою собственную функцию свертки/деконволюции (это не так сложно, посмотрите, как для этого я думаю, что большинство реализаций основаны на быстром преобразовании Фурье (БПФ)) - person pygri; 17.10.2020