Discrete convolution in time/Laplace domain

Hello,

I'm dealing with a problem where I need to calculate a convolution in the time domain again and again, so efficiency is a big issue.

I need to evaluate the following formula

numerically. I know as a analytical function (which btw can also be transformed analytically in the time domain); but I have only as discrete values . That means I need to evaluate in discrete form to obtain .

I have used FFT to solve some convolution in the Fourier space, but this Laplace transforms seem to be much more challenging numerically. Is maybe z-transforms the way to go? I would be very thankfull for any guidance on the topic.

Thank you.

[Paul]

Hi Peter

Are the q_i that you have uniformly sampled in time? Do you have their associated time tags? Do you have any information or assumptions about q(t) for t less than the time tag associated with q_1 or t greater than the time tag associated with q_n? Any assumptions on the form of q(t) between the samples of q_i?

[Peter Uwsen]

Hi Paul,

Thanks for your answer!

Yes, I have q_i uniformly sampled in time, and I have the assiociated time tags. for t < 0, so before q_1, q is zero. For the time after q_n, I dont have any information. I think between the q_i q(t) can simply be linear.

Answer 1

[Paul]

Simplest and most efficient approach would be to compute samples r_i = r(i*T) where T is the sampling period and then use conv to approximate the convolution integral, i.e.,

u = T*conv(q,r);

Make sure to take enough samples r_i such that r(t) is equal to, or well appoximated by, zero outside the interval covered by the samples. This approach is simpler if r(t) = 0 for t < 0, but can still be implemented even if r(t) ~=0 for t < 0.

If q and/or r have lots of samples, the consider using cconv, which uses and FFT based approach that might be faster.

If you're willing to assume that q(t) is linear between its samples, then you can use interp1 to interpolate to an effectively higher sampling period before sampling r(t) and doing the convolution.

If you'd like, feel free to post an example and your code. The data in q can be saved in .mat file and uploaded using the Paperclip icon on the Insert menu, unless you can define q by an equation as an example. If you have R(s) as a sym, then you can include that in the .mat file as well. Otherwise, a picture (use the Image icon on the Insert menu) should be fine if R(s) isn't too complicated.

[Peter Uwsen]

Thank you again for you answer!

As I understand it would be more convienient/efficient to perform the operation in the time domain, rather then doing Laplace transforms and a multiplication in the Laplace domain?

I should then give two further information about the problem:

1. q(x,t) is actually an array in two dimensions: time and space; r is a function of one variable that is composed of x and t, i.e. r(x^2*t). The equation I want to compute is therefore actually and I need to do it for several discrete point
2. I can analytically calculate the systems response to a step input : . This is I think advantageous because in this case I dont need to take into account the full course of : I can simply interpret q as a series of steps. With a discretization in time with spacing T and index i the discrete formula would be: . And If I understand correctly this is exactly what conv() does, if q and r have the same length.

Heres a code exmaple, which should give the desired result:

N = 512;

M = 512;

x = linspace(1,2,N).'; % space discretization

t = linspace(0,1,M); % time discretization

r = 1+erf(x.^2*t); % examplary response of system to step

q = sqrt(x.^2+t.^2); % example input

u = zeros(N, M);

tic

for i = 1:N

steps = [q(i,1) diff(q(i,:))];

convolution = conv(steps, r(i,:));

u(i,:) = convolution(1:N);

end

toc

My question now is: Can this be made more efficient? Can I get rid of the loop?

[Paul]

With respect to item 1, I assume that there is a typo in the integrand and that the correct term is r(x^2*tau), not r(x^2*t). If that's correct, then it seems that u(x,t) is straighforward convolution, assuming the integral exists (and that q(x,t) = r(x,t) = 0 for t < 0.

With respect to item 2, I'm confused regarding the introduction of the step response. What is r_s(x^2*t) and how does it relate to r(x^2*t) in item (1)? Similarly, what is q_s(x) and how is it related to q(x,t) in item (1)? Does the 'dot' symbol in u_s(x,t) = q_s(x) dot r_s(x^t) indicate multiplication or convolution?

I'm sure that further clarification on item (2) will help illuminate what that code snippet is attempting to do. Having said that, it seems that code has an issue. x varies down the columns of r and q while t varies across rows. But the loop is over 1 to N, which is the number of points in t points, but should be the number of points in x. Similarly for the u(i,:) line, seems like that should be assigning from convolution(1:M).

[Peter Uwsen]

Yes Im sorry, two typos actually, the upper limit is also wrong: it should be .

And the introduction of the step response was actually completely unnecessary as I realise now. I got confused myself sorry. (r is simply the derivative of r_s). This means that r is actually a function of x AND t I realized.

Concerning your last question x is in my opinion a row vector with N points and therefore I iterate over the rows of q and u respectively. But I agree that it should be u(i,:) = convolution(1:M). Here's my code again:

N = 512;

M = 512;

x = linspace(1,2,N).'; % space discretization

t = linspace(0,5,M); % time discretization

r = (2*x.^2.*exp(-x.^4.*t.^2))/sqrt(pi); % examplary function r

q = sqrt(x.^2+t.^2); % example input

u = zeros(N, M);

tic

for i = 1:N

convolution = conv(q(i,:), r(i,:));

u(i,:) = convolution(1:M);

end

toc

[Paul]

Yes, I agree that iterating over the rows is to iterate over the values in x, but x is column vector.

As long as q(x,t) = r(x,t) = 0 for t <= 0, then the upper limit of the integral can be inf or t.

The convolution sum has to be multiplied by dt properly approximate the convolution integral.

N = 512;
M = 512;
x = linspace(1,2,N).'; % space discretization
t = linspace(0,5,M); % time discretization

Define dt, knowing that t(1) = 0

dt = t(2);

Defining functions for q and r. Not really necessary for this problem, but will come in handy later.

rfunc = @(x,t) (2*x.^2.*exp(-x.^4.*t.^2))/sqrt(pi); % examplary function r
qfunc = @(x,t) sqrt(x.^2+t.^2); % example input

Get the values of q and r and compute the convolution sum

r = rfunc(x,t);
q = qfunc(x,t);
u = zeros(N, M);
tic
for i = 1:N
    convolution = conv(q(i,:), r(i,:));
    u(i,:) = dt*convolution(1:M);  % multiply by dt
end
toc

Elapsed time is 0.017662 seconds.

For comparison, define u(x,t) as the integral

ufunc = @(x,t) integral(@(tau) qfunc(x,t-tau).*rfunc(x,tau),0,t);

Compare to the convolution sum for one value of x

figure
plot(t,u(100,:)),hold on
uintegral = zeros(1,M);
tic
for jj = 1:M
    uintegral(jj) = ufunc(x(100),t(jj));
end
toc

Elapsed time is 0.063278 seconds.

plot(t,uintegral,'r')
xlabel('t')

The numerical integral evaluation took about 3x as long to compute only 1/512th of the results.

Depending on the value of M, cconv might be faster than conv.

For this particular example, under the same assumption as used with the solution using conv that q(x,t) is constant between samples, we can get the closed form solution for the pulse response of r(x,t) wrt t, and then sum up all of the small pulse responses (I think you were thinking about a solution along these lines.

syms x t tau p dt Q real
assumeAlso([p dt],'positive');

syms r(x,t)
r(x,t) = (2*x.^2.*exp(-x.^4.*t.^2))/sqrt(sym(pi))

r(x, t) = 

Define a pulse function of amplitude Q over the interval p < t < p + dt

p(t) = Q*(heaviside(t-p) - heaviside(t-(p+dt)));

Convolution of p(t) and r(x,t)

z(x,t) = simplify(int(p(t-tau)*r(x,tau),tau,0,t),'IgnoreAnalyticConstraints',true)

z(x, t) = 

Function for numerical evaluation

zfunc  = matlabFunction(z(x,t))

zfunc = function_handle with value:
    @(Q,dt,p,t,x)Q.*(erf(x.^2.*(p-t)).*(heaviside(p-t)-1.0)-erf(x.^2.*(dt+p-t)).*(heaviside(dt+p-t)-1.0))

x = linspace(1,2,N).'; % space discretization
t = linspace(0,5,M); % time discretization

dt = t(2);

r = rfunc(x,t);
q = qfunc(x,t);

tic
u = squeeze(sum(zfunc(q,dt,t,reshape(t,1,1,[]),x),2));
toc

Elapsed time is 1.003539 seconds.

This solution is incredibly slow, probably because heaviside is a Symbolic function (we could overcome this) and I suspect that the evaluations of erf are costly.

But the solution does pretty well match the numerical solution of the integral.

figure
plot(t,u(100,:),t,uintegral,'r')
xlabel('t')

[Peter Uwsen]

Thank you so much for the detailed response, really appreciate it! I now understand everything better.

I'm still concerned about the computation time though: Is there maybe a way to use fft() and ifft() (its just so much faster), if I make assumptions for q outside tmax=5; for example q = const for t>tmax?

[Paul]

I've suggested above to investigate cconv, which uses the fft/ifft approach and may be a lot faster for large data sets. cconv does come with some small overhead. Interestingly, it does not pad to nextpow2 in an attempt to gain some efficiency. Search the doc for "Linear and Circular Convolution" to see how to correctly use fft/ifft to perform linear convolution if you want to roll your own.

rng(100);
N = [1e3 1e6];

for n = N
q = rand(1,n);
r = rand(1,n);

tic,c1 = conv(q,r);toc
tic,c2 = cconv(q,r);toc
norm(c1-c2,'inf')
end

Elapsed time is 0.039565 seconds.

Elapsed time is 0.029024 seconds.

ans = 5.6843e-13

Elapsed time is 8.851229 seconds.

Elapsed time is 0.181574 seconds.

ans = 2.6368e-08

[Peter Uwsen]

Ok great, thank you! I will check the doc.