Trapezoidal Rule of Convolution with Non-Uniform Intervals
20 次查看(过去 30 天)
显示 更早的评论
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
0 个评论
采纳的回答
Catalytic
2023-10-12
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.
0 个评论
更多回答(1 个)
Matt J
2023-10-12
编辑:Matt J
2023-10-12
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 个评论
Paul
2023-10-13
编辑:Paul
2023-10-13
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)
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
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Matrices and Arrays 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!