1
$\begingroup$

This is a question that I've previously asked over on math.stackexchange, and I have yet to receive a useful answer. It was suggested that I post this here.

The problem itself originally comes from simulation of wave propagation, and involves large two-dimensional arrays of complex data. I'll explain it here as a simplified one-dimensional problem, however.

Suppose I have the following integral relationship involving a convolution: $$ f(x) = \int_{-\infty}^{+\infty}dx'\, g(x')\, K(x - x') $$ Assume that the functional form of the kernel $K$ is known, and suppose that the function $g$ has compact support on $[-a/2, +a/2]$ for some width $a>0$. Assume further that I have data about the function $g$ in the form of $N$ equally-spaced samples $g_n$: $$ g_n = g(x_n)\, , $$ where $$ x_n = -\frac{a}{2} + n\frac{a}{N}\, , \qquad n = 0,..., N-1 $$ My goal is to generate some $N$ evenly-spaced samples of the function $f$.

If we suppose that we want to evaluate $f$ on the same domain on which $g$ lives, and for which we have samples of $g$, then it's easy. I discretize everything in the usual way, and the whole thing becomes a discrete convolution: $$ f_n = f(x_n) = \sum_{m = 0}^{N-1} g_m K_{n - m}\, , $$ where $$ K_{n - m} = K(x_n - x_m) = K\left(\frac{a}{N}(n-m)\right)\, . $$ Everything is simple since I can use the discrete convolution theorem: I FFT the $\{g_n\}$ and $\{K_n\}$ sequences, multiply them, do an inverse FFT, and I'm done. I have my $\{f_n\}$ data. The whole thing will be an $O(N\log N)$ process.

Suppose instead that I have to evaluate $f$ on a different-sized domain. For example, perhaps the $g$ data is narrow, but the kernel $K$ is very wide. Let's say I have to generate $N$ samples of $f$ on the domain $[-L/2, +L/2]$, for some domain width $L > a$. $L$ could even be much larger than a.

If I try to discretize everything in a similar way, I'd have to first define a sequence of points $\{X_n\}$ in the wider domain: $$ X_n = -\frac{L}{2} + n\frac{L}{N}\, , \qquad n = 0,..., N-1 $$ Then the data I seek is $\{f_n\}$ where $f_n = f(X_n)$. The original integral is still over the smaller domain, so we have: \begin{align} f_n &= \int_{-a/2}^{+a/2}dx'\, g(x')\, K(X_n - x') \\ &\approx \sum_{m = 0}^{N-1} g_m\, K(X_n - x_m) \qquad \text{(discretize)}\\ &= \sum_{m = 0}^{N-1} g_m\, K\left(\frac{L-a}{2} + \frac{a}{N}(\alpha\, n - m)\right)\, , \end{align} where $\alpha = L/a > 1$ is the ratio of the two domain sizes.

This is no longer a discrete convolution, since the second term is no longer a function of $n-m$, but instead of $\alpha\, n - m$. So I can no longer use the convolution theorem. I can imagine evaluating the sum above by brute force, but that would be an $O(N^2)$ process, which is prohibitive for my problem.

Is there any way of efficiently calculating the discrete version of this convolution integral between two different domain sizes? Thanks.

$\endgroup$
1
$\begingroup$

Essentially, to my understanding this problems boils down to a fast-convolution problem:

You have a function $g(x)$, which has a small domain and you have a function $K(\tau)$, which has a big domain and you want to calculate

$$f(x) = \int_{-\infty}^{\infty}g(x')K(x-x')dx'$$

in a discretized fashion, where you want to calculate as efficient as possible. Especially, you are interested in $f(-L/2+nL/N)$ for $n=0,\dots,N-1$.

When considering the discretization you need to take care of the bandwidth of $g$ and $K$ and you need to sample with the according Nyquist frequency at least. Let's assume that $g$ can be sufficiently sampled with $g(-\frac{a}{2}+n\frac{a}{N})$ (i.e. the sampling is fast enough regarding Nyquist frequency). Let's also assume, that $K(\tau)$ has the same or smaller bandwith, i.e. the sampling rate for $g$ is also sufficient for $K$.

Then, let's further assume $L=qa, q\in\mathbb{N}$, i.e. L is a multiple of $a$. Define the following functions:

$$g[n]=g(-\frac{L}{2}+n\frac{L}{qN}), \quad n=0,\dots,qN-1$$ and $$K[n]=K(-\frac{L}{2}+n\frac{L}{qN}), \quad n=0,\dots,qN-1$$

Note, that $g[n]$ is zero for most $n$, but that is fine here.

Now, you have

$$f[n]=\sum_n'g[n']*K[n-n'], \quad n=0,\dots,qN-1$$

Which is a higher-resolution version of your original $f_n$ (It has qN instead of only N samples within the interval $\pm\frac{L}{2}$). So, what you actually want to have is a downsampled version of this $f[n]$ given by

$f_n=f[qn]$ (where $f_n$ is the $f_n$ from your post, i.e. with $n=0,\dots,N-1$.

So, this is also a discrete convolution problem, and it can be solved via FFT (note that the FFT actually considers circular convolution, if you don't zero-pad). However, since your $g[n]$ has a much smaller domain that $K[n]$, you can resort to fast-convolution algorithms (Overlap-save or Overlap-add).

To summarize: The trick is actually just to sample both signals with the same sampling frequency (accepting that $g[n]$ contains a lot of zeros) and then understanding this problem as a simple discrete convolution.

$\endgroup$
  • $\begingroup$ Unfortunately, I can't just sample my kernel $K$ with finer sampling frequency, it would generate too much data to handle. My original problem is 2D, with an array of complex numbers about 1000 x 1000 in size. For my problem, $L/a = q = \alpha \sim 100$, which means that the over-sampled $K$ would contain about $10^{10}$ complex numbers. $\endgroup$ – John Barber Jan 25 '17 at 20:11
  • $\begingroup$ OK, I see, it's a memory problem. This would have been an important message in the question. Then, can you sample $g$ with the frequency you want to sample $K$ (i.e. is the bandwidth of $g$ small enough for that)? This would easily solve your problem. If not, you can low-pass filter $g$ and then do the sampling with the frequency you want for $K$. You can do this, assuming that there is no aliasing in sampling $K$ (i.e. after convolution, the higher parts of $g$ would be gone anyway since $K$ is much narrower in frequency domain) $\endgroup$ – Maximilian Matthé Jan 25 '17 at 21:24
  • $\begingroup$ Unfortunately, $g$ is already sampled at the appropriate sampling frequency for that function. I can't low-pass filter it without loosing too much information. My motivation behind this question is the fact that I can calculate the convolution in only a second or two when the domains are the same size ($L = a$). Changing the size of the final domain appears to leave only the brute force approach, which takes one or two days. I was amazed that that one change could make such a drastic difference, and I hoped there would be some way of doing it faster which would be well-known to dsp people. $\endgroup$ – John Barber Jan 25 '17 at 21:38
  • $\begingroup$ Consider $g$ is of bandwidth $B_g$ and $K$ is of bandwidth $B_K$ with $B_K<B_g$. Then, the convolution of both signals (in continuous time) would have bandwidth $min(B_K,B_g)=B_K$, since both spectra are multiplied in frequency domain. So, I think you can low-pass filter $g$ before the convolution to the bandwidth of $K$, and then you can apply the lower sampling rate. $\endgroup$ – Maximilian Matthé Jan 26 '17 at 7:33

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.