Trapezoidal Rule of Convolution with Non-Uniform Intervals

23 views (last 30 days)
Hi all,
I'm supposed to do the numerical integration of a convolution for t, which is given by specific non-uniform timepoints. For t, there is also a function Cp which is a value associated with the t value.
Here's what I have so far. It doesn't work but I hope I'm on the right track?
Thanks in advance!
t=4.80;
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);
function trape=traperule(t)
time=[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];
[~,pos] = ismembertol(t, time, 1E-8)
tp=time(1:pos)
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];
C=conc(1:pos);
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);
trape=0;
for n=1:pos
tprime=tp(n)
if n>1
told=tp(n-1)
else
told=0;
end
Cp=C(n);
fun = @(tprime) 0.0675*Cp.*exp(-alpha1*(t-tprime)+Cp.*0.0345*exp(-alpha1*(t-tprime)));
Ci=integral(fun,told,tprime);
h=tprime-told;
area = Ci * h / 2 % area of the trapezoid
trape = trape + area;
end
end

Accepted Answer

Catalytic
Catalytic on 12 Oct 2023
I don't believe there is such a thing as the "Trapezoidal Rule of Convolution". Undoubtedly, you are thinking of the trapezoidal rule of integration, however, that will not help you deal with the non-uniform time sampling. Trapezoidal integration assumes you can approximate the integrand with linear interpolation. However, you do not have samples of the convolution integrand with which to interpolate. You could linearly interpolate the two signals you are trying to convolve and then multiply them together, but that would be a different kind of integral apprxoimation.

More Answers (1)

Matt J
Matt J on 12 Oct 2023
Edited: Matt J on 12 Oct 2023
It is not clear from your post which two functions you are convolving, so I will assume here that it is the same as your other recent post here. It does not give significantly different results as the other post, however, because in both posts the signals being convolved are approximated at some stage by trapezoids.
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);
vlin=@(tq) interp1(t,v,tq,'linear',0);
cplin=@(tq) interp1(t,cp,tq,'linear',0);
dt=min(diff(t));
tnew=min(t):dt:2*max(t);
fun=@(tau) integral( @(t) vlin(t).*cplin(tau-t), min(t),max(t) );
Ci=arrayfun(fun,tnew);
%CC is convolution of Cp(t)*the exponential function
plot(t,cp,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
  8 Comments
Lilly
Lilly on 13 Oct 2023
In that case, thank you so much for your help! I appreciate it
Paul
Paul on 13 Oct 2023
Edited: Paul on 13 Oct 2023
Assuming that Cp(t) is defined by the linear interpolation between the points, and that Cp(t) = 0 for t < 0 and for t > 591, and that v(t) = 0 for t < 0 (which are the same assumptions used above), then a closed form expression for the convolution Cp(t)*v(t) can be obtained.
tconc = [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];
syms t real
v(t) = myfunction(t);
syms Cp(t)
Cp(t) = 0;
for ii = 1:numel(tconc)-1
Cp(t) = Cp(t) + trapezoidalPulse(tconc(ii),conc(ii),tconc(ii+1),conc(ii+1),t);
end
figure
fplot(Cp(t),[0 600])
hold on
syms tau real
Ci(t) = int(Cp(tau)*v(t-tau),tau,0,t);
% use vpa to make the display marginally readable
vpa(Ci(t),5)
ans = 
fplot(Ci(t),[0 1200])
In the interval 0 < t < 591, this solution for Ci(t) is very close to that developed by @Matt J with a small difference presumably due to the difference between the exact solution returned by int and the approximate solution returned by integral and that Matt J's solution used linear interpolation for v(t). However, for t > 591 the solutions diverge. I believe that's because Matt J's solution assumes that v(t) = 0 for t > 591 whereas this soution uses the definition of v(t) as given, which has infinite duration.
Matt J's solution could be modified so that the integrand would be defined as
@(t) myfunction(t).*cplin(tau-t)
and then the solutions would (presumably, I didn't try it) match up quite well. Or the solution here could (presumably, I didn't try it) be modified to window v(t) over the inteval 0 < t < 591 if that's the appropriate assumption.
function y = trapezoidalPulse(x1,y1,x2,y2,x)
y = (y2 - y1)/(x2 - x1)*(x - x1) + y1;
y = y*rectangularPulse(x1,x2,x);
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);
end

Sign in to comment.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!