Блог

У меня было два набора данных: функция вывода системы (временной ряд длиной 1292 записи) и функция передачи (аналогичная гауссовой с длиной 681 запись). Я хотел бы вычислить входную функцию(неизвестную) с помощью деконволюции.

Я попытался вычислить БПФ для каждой функции, а затем обратную величину деления двух БПФ:

Cin = F^-1 [F(Cout)]/[F(E)]

однако БПФ передаточной функции, по-видимому, везде равна нулю, и результирующий восстановленный сигнал имеет много шума. Я также попробовал деконволюцию Винера, но еще раз, если я вычислю свертку восстановленного сигнала с помощью передаточной функции, я не получу выходную функцию.

• Для сигналов, которые я использую, какой метод наиболее подходит? (БПФ, сверток СциПи, Сосиска)
• Нужно ли мне «обнулять» передаточную функцию, чтобы она соответствовала длине выходной функции?
• Поскольку процесс деконволюции восприимчив к шуму, нужно ли фильтровать сигналы (до/после)?
• есть подобный код в matlab или python, на который я мог бы сослаться?

Здесь вывод моего кода с помощью метода Винера (вывод синим цветом, передача красным, восстановленный ввод зеленым цветом и свертка восстановленного ввода и передачи (фиолетовый):

и код для деконволюции Винера:

 #import data
data = pandas.read_csv("data/hetre20mod.csv")
x = data.Time

#normalize data
norm = (data.CO-data.CO.min()) / (data.CO.max()-data.CO.min())

# transfer function
length = 681
pe = 27
ttau = 203
ao = rtdpy.AD_oo(tau=ttau, peclet=pe, dt=1, time_end=length)
rtd = ao.exitage

#zero pad transfer function
kernel = np.hstack((rtd, np.zeros(len(norm) - len(rtd))))

#wiener deconvolution
lambd = 0.00035
H = np.fft.fft(kernel)
deconv = np.fft.ifft(np.fft.ifft(np.fft.fft(norm)*np.conj(H)/(H*np.conj(H)   lambd**2)))


recovered = np.convolve(deconv, kernel, mode='same') 

fig, ax = plt.subplots(nrows=2, ncols=2, figsize=(4, 4))
ax[0][0].plot(x, norm, 'b', label="experimental", lw=3)
ax[0][1].plot(x, kernel, 'r', label="transfer", lw=3)
ax[1][0].plot(x, deconv, 'g', label="treated", lw=3)
ax[1][1].plot(x, recovered, 'm', label="both", lw=3)
ax[1][1].plot(x, norm, 'b',label="both", lw=3)
plt.show()