How to graph a convolution with same time length as inputs
显示 更早的评论
Hello,
So far this is my code. What I'm trying to do is plot Ci and Cp over time, but Ci is a convolution of Cp and v. Both Cp and v are of length 27, so Ci is of length 53 since it is a convolution. The issue I'm having here is that I don't know how to get an accurate representation of Ci when graphing since time is given by specific points (see t). These specific time points are also directly related to Cp. What's the best way to do this?
Thank you in advance!
%v denotes the exponential function
v=myfunction(t);
%conce is the function for Cp(t)
conce=Cp(t);
%CC is convolution of Cp(t)*the exponential function
CC=conv(conce,v);
%Ci is the convolution having the same size of the two input vectors
Ci=CC(1:27);
plot(t,Ci,t,conce);
title("Concentration over Time");
xlabel('Time (mins)');
ylabel('Concentration')
legend("Ci(t)","Cp(t)");
%Based on plot, it appears that the highest signal strength in the brain is
%around 39.93 minutes, which should be the best time to run the PET scan.
function conce=Cp(a)
t=[0,1.08,1.78,2.30,2.75,3.30,3.82,4.32,4.80,5.28,5.95,6.32,6.98,9.83,16.30,20.25,29.67,39.93,58,74,94,100,200,300,400,500,591];
conc=[0,84.9,230,233,220,236.4,245.1,230.0,227.8,261.9,311.7,321,316.6,220.7,231.7,199.4,211.1,190.8,155.2,140.1,144.2,139.737,111.3006,82.8375,54.3744,25.9113,0.0098];
conce=conc(t==a);
end
function v=myfunction(t);
k1=0.102;
k2=0.130;
k3=0.062;
k4=0.0068;
alpha1=(k2+k3+k4+(sqrt((k2+k3+k4)^2-(4*k2*k4))))/2;
alpha2=(k2+k3+k4-(sqrt((k2+k3+k4)^2-(4*k2*k4))))/2;
A=(k1*(k4+k3-alpha1))/(alpha2-alpha1);
B=(k1*(alpha2-k4-k3))/(alpha2-alpha1);
v=A*exp(-alpha1*t)+B*exp(-alpha2*t);
%Ci=Cp(t).*(A.*exp(-1*alpha1*t)+B.*exp(-1*alpha2*t));
end
2 个评论
Walter Roberson
2023-10-11
Have you considered looking at the 'same' and 'valid' options of conv() ?
Lilly
2023-10-11
采纳的回答
Convolution doesn't make sense if your time points are not equi-spaced. You should interpolate everything so that they are:
t=[0,1.08,1.78,2.30,2.75,3.30,3.82,4.32,4.80,5.28,5.95,6.32,6.98,9.83,16.30,20.25,29.67,39.93,58,74,94,100,200,300,400,500,591];
%v denotes the exponential function
v=myfunction(t);
%conce is the function for Cp(t)
cp=Cp(t);
dt=min(diff(t));
tnew=min(t):dt:max(t);
vnew=interp1(t,v,tnew);
cpnew=interp1(t,cp,tnew);
%CC is convolution of Cp(t)*the exponential function
Ci=conv(cpnew,vnew,'same')*dt;
plot(tnew,cpnew,tnew,Ci);
title("Concentration over Time");
xlabel('Time (mins)');
ylabel('Concentration')
legend("Cp(t)","Ci(t)");
%Based on plot, it appears that the highest signal strength in the brain is
%around 39.93 minutes, which should be the best time to run the PET scan.
function conce=Cp(a)
t=[0,1.08,1.78,2.30,2.75,3.30,3.82,4.32,4.80,5.28,5.95,6.32,6.98,9.83,16.30,20.25,29.67,39.93,58,74,94,100,200,300,400,500,591];
conc=[0,84.9,230,233,220,236.4,245.1,230.0,227.8,261.9,311.7,321,316.6,220.7,231.7,199.4,211.1,190.8,155.2,140.1,144.2,139.737,111.3006,82.8375,54.3744,25.9113,0.0098];
conce=conc(t==a);
end
function v=myfunction(t);
k1=0.102;
k2=0.130;
k3=0.062;
k4=0.0068;
alpha1=(k2+k3+k4+(sqrt((k2+k3+k4)^2-(4*k2*k4))))/2;
alpha2=(k2+k3+k4-(sqrt((k2+k3+k4)^2-(4*k2*k4))))/2;
A=(k1*(k4+k3-alpha1))/(alpha2-alpha1);
B=(k1*(alpha2-k4-k3))/(alpha2-alpha1);
v=A*exp(-alpha1*t)+B*exp(-alpha2*t);
%Ci=Cp(t).*(A.*exp(-1*alpha1*t)+B.*exp(-1*alpha2*t));
end
5 个评论
But when using 'same', the values in Ci don't correspond to the values in tnew. In the plot below, Cisame has to be shifted to the right by 295 seconds to be consisten with time reference for the signals being convolved.
t=[0,1.08,1.78,2.30,2.75,3.30,3.82,4.32,4.80,5.28,5.95,6.32,6.98,9.83,16.30,20.25,29.67,39.93,58,74,94,100,200,300,400,500,591];
%v denotes the exponential function
v=myfunction(t);
%conce is the function for Cp(t)
cp=Cp(t);
dt=min(diff(t));
tnew=min(t):dt:max(t);
vnew=interp1(t,v,tnew);
cpnew=interp1(t,cp,tnew);
%CC is convolution of Cp(t)*the exponential function
Cisame = conv(cpnew,vnew,'same')*dt;
Cifull = conv(cpnew,vnew,'full')*dt;
plot(tnew,cpnew,tnew,Cisame,(0:numel(Cifull)-1)*dt,Cifull);
title("Concentration over Time");
xlabel('Time (mins)');
ylabel('Concentration')
legend("Cp(t)","Cisame(t)","Cifull(t)");
%Based on plot, it appears that the highest signal strength in the brain is
%around 39.93 minutes, which should be the best time to run the PET scan.
function conce=Cp(a)
t=[0,1.08,1.78,2.30,2.75,3.30,3.82,4.32,4.80,5.28,5.95,6.32,6.98,9.83,16.30,20.25,29.67,39.93,58,74,94,100,200,300,400,500,591];
conc=[0,84.9,230,233,220,236.4,245.1,230.0,227.8,261.9,311.7,321,316.6,220.7,231.7,199.4,211.1,190.8,155.2,140.1,144.2,139.737,111.3006,82.8375,54.3744,25.9113,0.0098];
conce=conc(t==a);
end
function v=myfunction(t);
k1=0.102;
k2=0.130;
k3=0.062;
k4=0.0068;
alpha1=(k2+k3+k4+(sqrt((k2+k3+k4)^2-(4*k2*k4))))/2;
alpha2=(k2+k3+k4-(sqrt((k2+k3+k4)^2-(4*k2*k4))))/2;
A=(k1*(k4+k3-alpha1))/(alpha2-alpha1);
B=(k1*(alpha2-k4-k3))/(alpha2-alpha1);
v=A*exp(-alpha1*t)+B*exp(-alpha2*t);
%Ci=Cp(t).*(A.*exp(-1*alpha1*t)+B.*exp(-1*alpha2*t));
end
In the plot below, Cisame has to be shifted to the right by 295 seconds to be consisten with time reference for the signals being convolved.
Only because you've assumed, arbitrarily, that Cifull starts at t=0, rather than t=-295. But yes, if we assume that v represents the response of a causal system, then that would be the appropriate starting time.
Lilly
2023-10-12
更多回答（0 个）
类别
在 帮助中心 和 File Exchange 中查找有关 Volume Visualization 的更多信息
Translated by 
