如何解决Python/Matlab 中系统响应的反卷积
我有两组数据,系统的输出函数(长度为 1292 个条目的时间序列)和传递函数(类似于一个长度为 681 个条目的高斯)。我想使用反卷积计算输入函数(未知)。
我尝试计算每个函数的 FFT,然后计算两个 fft 的除法的倒数:
Cin = F^-1 [F(Cout)]/[F(E)]
然而,传递函数的 FFT 似乎处处为零,并且得到的恢复信号有很多噪声。 我也尝试过 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,'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 举报,一经查实,本站将立刻删除。