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

[Walter Roberson]

Have you considered looking at the 'same' and 'valid' options of conv() ?

[Lilly]

Hello, yes I have but the graph when using 'same' doesn't look the way that it should.

Answer 1

[Matt J]

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

[Paul]

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

[Matt J]

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.

[Paul]

No, I didn't assume that Cifull starts at t=0. Doesn't that follow from the fact that both v and Cp start from t=0?

[Matt J]

Well, v is not a signal, but rather, a convolution kernel so it doesn't really "start" at a particular time. If we assume it is the impulse response of a causal system then, yes, we would normally assume everything starts at t=0.

[Lilly]

Thank you guys! It makes sense.

I appreciate the help!