Successive convolutions

Hello,

I have a question that I am hopeful you will help me with...

I need to numerically compute a number (K) of successive convolution
integral operators [R] to a function [f], i.e.

Rf(x) = integral[f(y) g(x-y) dy ]

My problem is: in the first step, function [f] is explicit, hence I can
obtain f(y) at any y points easily. However, in the following steps we
don't know [f] explicitly, i.e. all we have is Rf(x) from the first step.
What I mean, in the second step, we want

R(Rf(x)) = integral[Rf(y) g(x-y) dy] ; where Rf(y) is not known (especially
since the y grid is usually wider than the x grid in practice)

I know that using the convolution theorem is the first answer that one
thinks of, i.e. Rf(x) = F^-1 [Ff(u) . fg(u)] ; however, the functions I am
dealing with have Fourier transforms that decay extremely slowly, and the
convolution theorem solution returns very large errors.

I prefer as such to compute the convolutions directly using a quadrature.
The same explicit kernel is to be used in all steps, hence we have no
problems with g(x-y). 


Hopefully, I am in the right place to ask this question... If you can help
me, or point me to any papers or code that deal with a similar problem, I
will be extremely grateful..

Yours,
Tim

tamzam600 wrote:
> Hello,
> 
> I have a question that I am hopeful you will help me with...
> 
> I need to numerically compute a number (K) of successive convolution
> integral operators [R] to a function [f], i.e.
> 
> Rf(x) = integral[f(y) g(x-y) dy ]
> 
> My problem is: in the first step, function [f] is explicit, hence I can
> obtain f(y) at any y points easily. However, in the following steps we
> don't know [f] explicitly, i.e. all we have is Rf(x) from the first step.
> What I mean, in the second step, we want
> 
> R(Rf(x)) = integral[Rf(y) g(x-y) dy] ; where Rf(y) is not known (especially
> since the y grid is usually wider than the x grid in practice)
> 
> I know that using the convolution theorem is the first answer that one
> thinks of, i.e. Rf(x) = F^-1 [Ff(u) . fg(u)] ; however, the functions I am
> dealing with have Fourier transforms that decay extremely slowly, and the
> convolution theorem solution returns very large errors.
> 
> I prefer as such to compute the convolutions directly using a quadrature.
> The same explicit kernel is to be used in all steps, hence we have no
> problems with g(x-y). 
> 
> 
> Hopefully, I am in the right place to ask this question... If you can help
> me, or point me to any papers or code that deal with a similar problem, I
> will be extremely grateful..

I hope you're aware that (a * b) * c is equal to a * (b *c), where * 
indicates convolution. If you are, what's the problem?

Jerry
-- 
Everyone is entitled to an opinion; not all opinions deserve respect.
���������������������������������������������������������������������

>tamzam600 wrote:
>> Hello,
>> 
>> I have a question that I am hopeful you will help me with...
>> 
>> I need to numerically compute a number (K) of successive convolution
>> integral operators [R] to a function [f], i.e.
>> 
>> Rf(x) = integral[f(y) g(x-y) dy ]
>> 
>> My problem is: in the first step, function [f] is explicit, hence I can
>> obtain f(y) at any y points easily. However, in the following steps we
>> don't know [f] explicitly, i.e. all we have is Rf(x) from the first
step.
>> What I mean, in the second step, we want
>> 
>> R(Rf(x)) = integral[Rf(y) g(x-y) dy] ; where Rf(y) is not known
(especially
>> since the y grid is usually wider than the x grid in practice)
>> 
>> I know that using the convolution theorem is the first answer that one
>> thinks of, i.e. Rf(x) = F^-1 [Ff(u) . fg(u)] ; however, the functions I
am
>> dealing with have Fourier transforms that decay extremely slowly, and
the
>> convolution theorem solution returns very large errors.
>> 
>> I prefer as such to compute the convolutions directly using a
quadrature.
>> The same explicit kernel is to be used in all steps, hence we have no
>> problems with g(x-y). 
>> 
>> 
>> Hopefully, I am in the right place to ask this question... If you can
help
>> me, or point me to any papers or code that deal with a similar problem,
I
>> will be extremely grateful..
>
>I hope you're aware that (a * b) * c is equal to a * (b *c), where * 
>indicates convolution. If you are, what's the problem?
>
>Jerry
>-- 
>Everyone is entitled to an opinion; not all opinions deserve respect.
>���������������������������������������������������������������������
>

Hi, many thanks for your reply. 
I am aware of the associative law of convolutions...
The problem is: how can you compute the convolution operator [R] on f(x):

Rf(x) = integral [f(y) g(x-y) dy]

where f(x) is not given by an explicit function, i.e. it is quantized or
sampled at grid x. If you would like to compute the convolution by
quadrature, you will notice that we need f(y), where grid y is typically
wider than grid x, hence we are missing certain values of function f at
points y that are not included in grid x. There is no problem with the
explicit kernel g(x-y), I am only having a problem with f(y).

If I had to compute the convolution operator R only once, then I would not
have a problem since function f(.) is explicit in the first step (or first
application of the convolution operator), however, in the second step, we
wish to compute

R(Rf(x)) = integral [Rf(y) g(x-y) dy]
The problem is obvious: we have sampled numeric values of Rf(x) from the
first step, but we don't have Rf(y)...

I know it is a tough problem... Is there a solution that does not involve
the convolution theorem {f*g} (x) = F^-1[F{f}(u).F{g}(u)] (x)? 


Many thanks!
TIM





 

tamzam600 wrote:
>I need to numerically compute a number (K) of successive convolution
>integral operators [R] to a function [f], i.e.
>
>Rf(x) = integral[f(y) g(x-y) dy ]
>
>My problem is: in the first step, function [f] is explicit, hence I can
>obtain f(y) at any y points easily. However, in the following steps we
>don't know [f] explicitly, i.e. all we have is Rf(x) from the first step.
>What I mean, in the second step, we want
>
>R(Rf(x)) = integral[Rf(y) g(x-y) dy] ; where Rf(y) is not known
(especially
>since the y grid is usually wider than the x grid in practice)


Why not use the brute-force solution of recursively computing the next
integral you need at the point you need?  You mention knowing how to do
quadrature, so perhaps you can optimize that process by figuring out how to
save past values in a table, though this optimization may or may not be
difficult with matlab calling conventions (pass-by-value).  If your
quadrature uses a constant stepsize (not generally recommended), you can do
it ahead of time, easing bookkeeping.

If you don't like this option (I'm honestly surprised it hasn't occurred to
you, given the detail you gave elsewhere, so maybe I'm misunderstanding or
this is computationally prohibitive?), maybe if you said why you want this,
and what the form of your functions is (constant-frequency?  FM?  hamonics?
 noise?  etc?), someone might be able to offer an alternative approach.




>I prefer as such to compute the convolutions directly using a quadrature.
>The same explicit kernel is to be used in all steps, hence we have no
>problems with g(x-y).