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

编程问答 2022-04-25

问题描述

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

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

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

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

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

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

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

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

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

wiener method output

和维纳反卷积的代码：

#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,'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)
plt.show()

解决方法

暂无找到可以解决该程序问题的有效方法，小编努力寻找整理中！

如果你已经找到好的解决方法，欢迎将解决方案带上本链接一起发送给小编。

小编邮箱:dio#foxmail.com (将#修改为@）

convolutiondeconvolutionmatlabmatlabpythonsignal-processing