| 查看: 990 | 回复: 0 | ||
jbsxh铜虫 (小有名气)
|
[求助]
Matlab 求解微信方程组，遇到反卷积的问题
|
|
数据处理的最后一步，是解一个3个微分方程联立的方程组的多个系数。 每一个独立方程中的Y都是实验采集到的随时间变化的信号CA,CB, CC，该信号隐含了仪器响应函数IRF. 如果忽略IRF，直接使用CA,CB, CC去接微分方程组，得到的系数与报道值差别非常大。请问这种情况下，应该如何处理，才能获得正确的微分方程系数？ dCAdt = -k0 * CA; dCBdt = k1 *CA - k2*CB+k3CC; dCCdt = k2 *CA + k2*CB- k3CC; k1+k2=2.3; k3+k4=5 另一种思路是对CA,CB, CC先进行去卷积的操作。目前又被去卷积的问题困住了。 求大神指点，谢谢！ 反卷积的问题如下： 问题简化如下： 步骤1. 成功 用一个指数衰减函数E与一个高斯函数G进行卷积，得到目标函数 F1=conv(E,G). 然后反卷积【q,r】=deconv(F1,G)也可以顺利得到指数函数E； 步骤2. 失败 目标函数F2采用解析方式获得： 同样是用E和G进行卷积。但是奇怪的是，F2和F1在x值较小的位置并不重合。这个感觉很奇怪！原因也不知道。 然后进行反卷积【q,r】=deconv(F2,G)，得到的q竟然与原来的F2非常相似，仅是在X轴上进行了移动。 步骤3.失败 对实验数据进行去卷积操作deconv，发现均是无法实现获得正常的结果。 代码如下： clc clear all dx=0.001; xx=-0.5:dx:10; xt=0.0:dx:7; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% exponential decay k1=1/0.179; PY= (exp(-k1.*(xx))) ; size(PY) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% gaussian c=0.0857/2.355; xc=+0.0233; a=1/sqrt(2*pi)/c; GS_Y = a.*exp(-1.*(xt-xc).*(xt-xc)/2/c/c); size(GS_Y); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% convolution with E decay and gaussian x_GS = (xx(1)+xt(1)):dx: ((length(xx)+length(xt)-2)*dx+(xx(1)+xt(1))); conv_Y = conv(PY,GS_Y); Max_Y = max(abs(conv_Y)); conv_Y = (conv_Y) / Max_Y; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% convolution with equation, X_data = -0.5:dx:10; KIT_parameters_PRE =[ 1; 0.0857; 0; xc; 0.179]; for II=1:length(X_data) YY(II) = KIT_parameters_PRE(3)+KIT_parameters_PRE(1)*(1-erf(-(1.6651/KIT_parameters_PRE(2)).*(X_data(II)-KIT_parameters_PRE(4))+(1.0/KIT_parameters_PRE(5))/(2*(1.6651/KIT_parameters_PRE(2))))).*exp(-(1.0/KIT_parameters_PRE(5)).*(X_data(II)-KIT_parameters_PRE(4))+(1.0/KIT_parameters_PRE(5)).*(1.0/KIT_parameters_PRE(5))./(4*(1.6651/KIT_parameters_PRE(2)).*(1.6651/KIT_parameters_PRE(2)))); end MAX_YY = max(YY); YY=YY/MAX_YY; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% plot %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% [q,r] = deconv(YY, GS_Y); X_data_GS = (xt(1)+X_data(1)): dx: ((length(X_data)-length(xt))*dx+(xt(1)+X_data(1))); subplot(3,1,1) plot(X_data_GS,q/max(abs(q)),'ro',X_data,YY,'ko' ) axis([-1 2 -0.1 1]) %by curve fitting [q,r] = deconv(conv_Y, GS_Y); X_data_GS = (xt(1)+x_GS(1)): dx: ((length(x_GS)-length(xt))*dx+(xt(1)+x_GS(1))); subplot(3,1,2) plot(X_data_GS,q/max(abs(q)),'ro',x_GS,conv_Y,'ko' ) axis([-1 2 -0.1 1]) %compare subplot(3,1,3) plot(X_data,YY,'k',x_GS,conv_Y,'r' ) axis([-1 2 -0.1 1]) |
回复此楼» 猜你喜欢
瑞典林雪平大学博士后招聘｜PEC Water Splitting 方向
已经有130人回复
-
PRB投稿，Acknowledgment sent to author状态十天了
已经有3人回复
-
物理学I论文润色/翻译怎么收费?
已经有238人回复
-
PVA溶解
已经有2人回复
-
基金申请
已经有44人回复
-
CSC与新西兰维多利亚大学PhD奖学金项目
已经有0人回复
-
新西兰Robinson研究所 招聘CSC公派访问人员
已经有0人回复
-
帮我的英语口语老师找学生
已经有0人回复
-
什么时候开奖？
已经有6人回复
-
散金币祈福
已经有94人回复
40