python - Python/Matlab中系统响应的反卷积

Question

我有两组数据，系统的输出函数（长度为 1292 个条目的时间序列）和传递函数（类似于长度为 681 个条目的高斯函数）。我想使用反卷积计算输入函数（未知）。

我尝试计算每个函数的 FFT，然后计算两个 fft 的除法的倒数：

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

然而，传递函数的 FFT 似乎到处都是零，并且得到的恢复信号有很多噪声。我也尝试了 Wiener 反卷积，但再一次，如果我计算恢复信号与传递函数的卷积，我不会获得输出函数。

• 对于我使用的信号，哪种方法最合适？（FFT、SciPy 卷积、维纳）

• 我是否需要“零填充”传递函数以匹配输出函数的长度？

• 由于去卷积过程容易受到噪声的影响，我是否需要过滤信号（之前/之后）？

• 我可以参考 matlab 或 python 中的类似代码吗？

这里是我的代码使用 wiener 方法的输出（蓝色输出，红色传输，绿色恢复输入以及恢复输入和传输的卷积（紫色）：

和维纳反卷积的代码：

#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()

标签： pythonmatlabsignal-processingconvolutiondeconvolution