Fast way to calculate a few points in convolution?

Say I have two 1000-point sequences. I can compute the entire convolution fairly efficiently using FFTs (keeping in mind that zeroes have to be padded on to avoid doing a circular convolution).

However, I'm only interested in about 500 points in this convolution. To be more specific, say the indices of the convolution are 0 to 1998; I'm only interested in 1000 - 1499 inclusive.

Is there some way to calculate only these points, without calcualting the entire convolution and throwing most of it away?

• You can use time aliasing in circular convolution for your benefit to reduce the length of DFT size. Let me put an answer for that. – Fat32 Aug 25 '17 at 21:22
• btw 500 among 2000 is far from few points... – Fat32 Aug 25 '17 at 22:50
• @Fat32 Actually, in the particular project I'm working on we only need 10 of the points. However, out of curiosity, I wanted to manufacture a situation where a brute force evaluation of the convolution sum for a "few points" is still slower than an FFT-based method on the whole sequence – Mahkoe Aug 26 '17 at 12:52
• Then that's a great mistake of yours! As you must know, if you need 10 out of 2000 time samples for the convolution then you should consider time domain processing, which will be both faster and simpler than DFT implementation. If you are loking for the decision boundary between time and frequency domain implementations, you should have asked it that way. – Fat32 Aug 26 '17 at 13:25

I'll write names of time domain sequences in lower case and their frequency-domain counterparts in upper case, for example $X = \operatorname{DTFT}(x)$, where $\operatorname{DTFT}$ stands for discrete-time Fourier transform. I chose this transform, because it is unaffected by zero-padding of the time-domain sequences, allowing equations that mix frequency-domain counterparts of sequences of different lengths. With the upper-case frequency-domain symbols, time-domain convolution is denoted by multiplication, and a time-domain delay of $N$ samples is denoted by $z^{-N},$ where $z = e^{i\omega},$ with $\omega$ as continuous frequency. No frequency-domain thinking is required here. You can mentally translate the equations while reading to time-domain convolutions, delays, and additions. Addition in frequency domain is addition in time domain.

The frequency-domain counterparts $X$ of $Y$ of time-domain sequences $x$ and $x$ each of length $2N,$ with $N=500,$ can be written in partitioned form $X = X_0 + X_1z^{-N}$ and $Y = Y_0 + Y_1z^{-N},$ where $\{X_0, X_1, Y_0, Y_1\}$ are frequency-domain counterparts of time-domain partitions $\{x_0, x_1, y_0, y_1\}$ each of length $N.$ The partitions are non-overlapping sub-sequences of $x$ and $y$, stripped of delay which is then reintroduced in the equations. Convolution of $x$ and $y$ can be written in frequency domain as:

\begin{align}&XY\\ = &(X_0 + X_1z^{-N})(Y_0 + Y_1z^{-N})\\ = &(X_0 + X_1z^{-N})Y_0 + (X_0 + X_1z^{-N})Y_1z^{-N}\\ = &X_0Y_0 + X_1Y_0z^{-N} + X_0Y_1z^{-N} + X_1Y_1z^{-N}z^{-N}\\ = &X_0Y_0 + (X_1Y_0 + X_0Y_1)z^{-N} + X_1Y_1z^{-2N}\end{align}\tag{1}

Now all convolutions between sequences are convolutions between partitions of length $N$. Each result of such a convolution is of length $2N-1$. Denoting the frequency-domain counterpart of the first $N$ samples of the convolution result by $H()$ as in "head" and the frequency-domain counterpart of the delay-stripped sub-sequence formed by the remaining $N-1$ samples by $T()$ as in "tail":

\begin{align}= &H(X_0Y_0) + T(X_0Y_0)z^{-N}\\ + &\big(H(X_1Y_0) + T(X_1Y_0)z^{-N}\big)z^{-N} + \big(H(X_0Y_1) + T(X_0Y_1)z^{-N}\big)z^{-N}\\ + &\big(H(X_1Y_1) + T(X_1Y_1)z^{-N}\big)z^{-2N}\\ = &H(X_0Y_0)\\ + &\big(T(X_0Y_0) + H(X_1Y_0) + H(X_0Y_1)\big)z^{-N}\\ + &\big(T(X_1Y_0) + T(X_0Y_1) + H(X_1Y_1)\big)z^{-2N}\\ + &T(X_1Y_1)z^{-3N}\end{align}\tag{2}

Now we have neatly split the full convolution result into length $N$ partitions. You were interested in the partition of length $N$ beginning at index $2N,$ or in frequency domain:

$$T(X_1Y_0) + T(X_0Y_1) + H(X_1Y_1)$$

The time-domain counterparts of these three terms can can be calculated by three FFT convolutions of length $2N$ compared to your single FFT convolution of length $4N.$ I don't know if that is of any use. You'd still be throwing away a total of $3N$ samples of convolution results, as originally.

But wait! The time-domain counterpart of $T(X_0Y_1) + H(X_1Y_1)$ can be calculated by a circular FFT convolution of length $2N$ between $x$ and $y_1.$ First, writing that as linear (not circular) convolution:

\begin{align}&XY_1\\ = &X_0Y_1 + X_1Y_1z^{-N}\\ = &H(X_0Y_1) + T(X_0Y_1)z^{-N} + H(X_1Y_1)z^{-N} + T(X_1Y_1)z^{-2N}\\ = &H(X_0Y_1) + \big(T(X_0Y_1) + H(X_1Y_1)\big)z^{-N} + T(X_1Y_1)z^{-2N}\end{align}

Wrapping that around to make it circular, $z^{-2N}$ wraps to $z^{0} = 1:$

$$H(X_0Y_1) + T(X_1Y_1) + \big(T(X_0Y_1) + H(X_1Y_1)\big)z^{-N}$$

The second half of that has what we wanted, $T(X_0Y_1) + H(X_1Y_1).$ We still need to calculate the time-domain counterpart of $T(X1, Y0)$ by another length $2N$ FFT convolution so we are down to two FFT convolutions of length $2N,$ with only a total of $2N$ samples of the results wasted.

We can also confirm that user Fat32's suggestion of using a length $3N$ FFT convolution works, by wrapping the latter part of Eq. 2 around whereby $z^{-3N}$ becomes $z^{0} = 1$:

\begin{align}&H(X_0Y_0) + T(X_1Y_1)\\ + &\big(T(X_0Y_0) + H(X_1Y_0) + H(X_0Y_1)\big)z^{-N}\\ + &\big(T(X_1Y_0) + T(X_0Y_1) + H(X_1Y_1)\big)z^{-2N}\end{align}

The middle line has what we wanted. So a single length $3N$ FFT convolution is also an option, with only a total of $2N$ result samples wasted.

• Very interesting effort! Eventhough I couldn't understand some of its details (but have the feeling). Nevertheless computing 500 middle points among 2000 in time samples apparently do not lend into any efficient computation result at least using basic techniques. By the way, x is a time signal $x[n]$ or transform samples $X(z)$, $X(w)$ or $X[k]$ cause $z^{-N} x$ makes sense only for $X(z)$? It would be better if you declared it explicitly. – Fat32 Aug 26 '17 at 9:55
• Hmm so $x_0$ is a sequence of 1000 samples whose last 500 is zero and $x_1$ is a sequence of 1000 samples whose first 500 is zero so that $x[n] = x_0[n] + x_1[n]$ but indeed you consider $x_1[n]$ also to be resided at the $n=0$ position so you would say $x[n] = x_0[n] + x_1[n-500]$ whose transform domain samples would have $X = X_0 + z^{-N} X_1$ – Fat32 Aug 26 '17 at 10:17
• $x_0$ and $x_1$ are of length 500. – Olli Niemitalo Aug 26 '17 at 10:19
• oh yes, they are the nonzero portions of those extended verisons I defined... – Fat32 Aug 26 '17 at 10:36
• @Fat32 Now the notation is in frequency domain. – Olli Niemitalo Aug 28 '17 at 18:14

Given two discrete time sequences $x[n]$ and $h[n]$ which have lengths of $1000$ points each, their linear convolution $z[n] = x[n] \star h[n]$ will have a length of $L_z = 1000 + 1000 -1 = 1999$ pts.

The theory of DTFT (discrete-time Fourier transform) provides a relationship between the linear convolution of two sequences and multiplication of their respective DTFTs as;

$$x[n] \star h[n] \longleftrightarrow X(e^{jw}) H(e^{jw})$$ hence $$z[n] = \mathcal{F}^{-1} ( X(e^{jw}) H(e^{jw}) )$$

In practice the above relationship is implemented using the DFT (discrete Fourier transform) instead of DTFT which is more of a theoretical tool. Therefore, the relationship is practically implemented such that the resulting sequence is not $z[n]$ but a periodic extention of the theoretical signal $z[n]$ given by the relationship: $$\tilde{z} [n] = \sum_k z[n - k N]$$ where $N$ is the number of points used in computing the inverse (and as well as the forward) DFTs in implementing the algorithm.

Now if $N$ is chosen such that $N \geq 1999$ samples then the first period of the periodic signal $\tilde{z}[n]$ will be equal to the theoretical result of the linear convolution;i.e., $\tilde{z}[n] = z[n] ~\text{ for} ~~ n=0,1,...,1998$ iff $N \geq 1999$. Now if $N$ is chosen to be less than $1999$ then there will be aliasing in the time samples of the periodic extension $\tilde{z}[n]$; i.e. shifted versions of $z[n]$ will overlapp in time. This overlapp is undesirable for most applicaitons, but in your case since you do not want a full linear convolution but samples of $z[n]$ between $n=1000$ and $n=1499$ , you may use the following DFT length: $$M = 1500$$ will create time aliased $z[n]$ such that samples from $n=0$ to $n=499$ and samples from $n=1500$ to $n=1998$ will be aliased and those from $n=500$ to $n=1499$ will be clean.

Note that when you perform the $1500$ point inverse DFT, your time sequence will involve $1500$ time samples from which you will discard the first $1000$ and retain only the last $500$ which yield your desired set of samples.

Also note that you are still computing 3x more samples than (500) you want to and thus waste 2x of them, but that's all you can do with the method I presented here and therefore I'm not sure whether you find this a big saving or not... Btw the best you can achieve with this method would happen had you selected a sample range of $500$ sampes in between (about) $n=750$ and $n=1250$ for which a DFT size about $N = 1250$ would suffice.

I'ld like to provide the same (or similar) answer posted by @OlliNiemitalo using a slightly different notation.

Consider two discrete-time sequences $x[n]$ and $y[n]$ of length $L=1000$ points each. We're interested in their linear convolution $z[n] = x[n] \star y[n]$ of length $L_z = 1999$ points, for not all but only those samples from $n=1000$ to $n=1499$ using a partitioning of the input signal $x[n]$ and $y[n]$. The method also relies on the DFT based implementation of the discrete time linear convolutions, which states:

$$z[n] = x[n] \star y[n] \longleftrightarrow Z(w) = X(w)Y(w) \tag{1}$$ therefore $z[n] = \mathcal{F}^{-1} \{ X(w) Y(w) \}$ , in practice we shall use DFT instead of DTFT, however.

Consider the half-length partitioning of $x[n]$ and $y[n]$ into $x_0[n], x_1[n], y_0[n], y_1[n]$ of length $N=500$ each, which are defined as:

\begin{align} x_0[n] &= x[n] &, \text{ for } &n = 0,1,...,N-1 \\ x_1[n] &= x[n+N] &, \text{ for } &n = 0,1,...,N-1 \\ y_0[n] &= y[n] &, \text{ for } &n = 0,1,...,N-1 \\ y_1[n] &= y[n+N] &, \text{ for } &n = 0,1,...,N-1 \\ \end{align}

So that $x[n]$ and $y[n]$ are now defined in terms of their partitions as: $$x[n]=x_0[n] + x_1[n-N] \longleftrightarrow X(w) = X_0(w) + e^{-jwN} X_1(w) \tag{2}$$ $$y[n]=y_0[n] + y_1[n-N] \longleftrightarrow Y(w) = Y_0(w) + e^{-jwN} Y_1(w)$$

Writing the linear convolution $z[n] = x[n] \star y[n]$ in terms of the partitions becomes: $$z[n] = x[n] \star y[n] = (x_0[n] + x_1[n-N]) \star (y_0[n] + y_1[n-N]) \tag{3}$$

$$z[n] = x_0[n]\star y_0[n] + x_0[n]\star y_1[n-N] + x_1[n-N] \star y_0[n] + x_1[n-N] \star y_1[n-N]$$

Defining the sub-convolutions as: \begin{align} z_0[n] =& ~x_0[n] \star y_0[n]\\ z_1[n-N] =& ~x_0[n] \star y_1[n-N] + x_1[n-N] \star y_0[n] \tag{4}\\ z_2[n-2N] =& ~x_1[n-N] \star y_1[n-N]\\ \end{align}

or equivalently as:

\begin{align} z_0[n] &= x_0[n] \star y_0[n]\\ z_1[n] &= x_0[n] \star y_1[n] + x_1[n] \star y_0[n] \tag{5}\\ z_2[n] &= x_1[n] \star y_1[n]\\ \end{align}

yields the following for $z[n]$: $$z[n] = z_0[n] + z_1[n-N] + z_2[n-2N] \tag{6}$$

The following figure, reveals the time domain relationship between $z[n]$, $z_0[n]$, $z_1[n]$ and $z_2[n]$ :

As the figure shows that the region of $z[n]$ that we are interested in, between $n=1000$ to $n=1499$ is given by the sum of the second half of $z_1[n]$ and the first half of $z_2[n]$. Note that $z_0[n]$ is irrelevant.

Hence we should only compute the subsequences $z_1[n]$ and $z_2[n]$ which we choose to do so using DFTs of length $2N = 1000$ points for the forward and inverse transforms. Note that in each subsequence we are computing $L=2N=1000$ time points of which one half is wasted.

Hence we need to compute the following subsequences $z_1[n]$ and $z_2[n]$ of length $L=2N=1000$ points each using 1000 point inverse DFT's of $Z_1[k]$ and $Z_2[k]$ which are $1000$ point DFTs related to $1000$ point DFTs of partitions from eq. $(5)$ as:

\begin{align} Z_1[k] &= X_0[k] Y_1[k] + X_1[k] Y_0[k] \tag{7}\\ Z_2[k] &= X_1[k] Y_1[k]\\ \end{align}

where each multiplicative term $X_0[k]$, $X_1[k]$, $Y_0[k]$, and $Y_1[k]$ will be implemented using $L=2N=1000$ points DFTs of corresponding partition signal $x_0[n]$, $x_1[n]$, $y_0[n]$, and $y_1[n]$ of length $N=500$ each.

Hence we have: \begin{align} z_1[n] &= \mathcal{F}^{-1} \{Z_1[k] \} = \mathcal{F}^{-1} \{ X_0[k] Y_1[k] + X_1[k] Y_0[k] \} \tag{8}\\ z_2[n] &= \mathcal{F}^{-1} \{Z_2[k] \} = \mathcal{F}^{-1} \{ X_1[k] Y_1[k] \} \\ \end{align}

A simple complexity analysis yeilds that for obtaining $Z_1[k]$ and $Z_2[k]$ we need to compute a total of $4$ , $2N=1000$ point FFTs and $3$ , $2N=1000$ point complex multiplications, which makes about $$4 \times (2N/2) \log_2(2N) \times 4 + 3 \times 2N \times 4 \tag{9}$$ real MACs, giving $(160 + 24)N$ real MACs.

Also obtaining $z_1[n]$ and $z_2[n]$ requires $2$ inverse $2N=1000$ point DFTs would is about $2 \times (2N)/2 \log_2(2N) \times 4$ real MACs, which gives $80N$ real MACs. Hence we have a total of $$264 \times N \tag{10}$$ real MACs.

(note I have replaced $\log_2(2N)$ with $10$ for $2N=1000$, otherwise the total real MAC count would be $24 (1 + \log_2(2N)) N$ real MACs)

A time domain implementation of the full convolution of $x[n]$ and $y[n]$ would require about $$4 N^2 = 4N \times N \tag{11}$$ real MACS.

A time domain implementation of the partial convolutions for $z_1[n]$ and $z_2[n]$ would require about: $$3 N^2 = 3N \times N \tag{12}$$ real MACS.

A frequency domain implementation of the full convolution requires about: $$3 \times (4N/2) \log_2(4N) \times 4 + 4 N = 24 N (6 + \log_2(N)) \tag{13}$$ real MACS which for $N=500$ is about $$360 \times N \tag{14}$$ real MACs.

And finally the frequency domain implementation partial convolution to compute $z_1[n]$ and $z_2[n]$ results in:

$$24 (2 + log_2(N)) N \tag{15}$$ real MACs which for $N=500$ is about $$264 \times N$$ real MACS.

It's clearly an imprevement over time domain implmentations and the full frequency domain implementation.

Note that the in my original answer the method of time aliased computation of the linear convolution would require, with $M=1500 = 3 N$ about $$18 (4.6 + \log_2(N)) N \tag{16}$$ real MACs which makes about $$244 \times N$$ real MACS for $N=500$.

Lets conclude with the note that, instead of a time domain partial convolutions approch we could have reached the same results by following a frequency domain based approach as follows , given the orginal partitioning:

$x[n] \star y[n] \leftrightarrow X(w)Y(w) = [X_0(w) + e^{-jwN} X_1(w)] \cdot [Y_0(w) + e^{-jwN} Y_1(w)] \tag{17}$

$$z[n] \longleftrightarrow X_0 Y_0 + e^{-jwN} [X_1 Y_0 + X_0 Y_1] + e^{-jw 2N} X_1 Y_1 \tag{18}$$

At this point $X_0$, $X_1$, $Y_0$, and $Y_1$ are all the DTFTs of the half sized sequences $x_0$, $x_1$, $y_0$, and $y_1$. Let's make the re-definitions:

\begin{align} Z_0(w) &= X_0(w) Y_0(w)\\ Z_1(w) &= X_0(w) Y_1(w) + X_1(w) Y_0(w) \tag{19}\\ Z_2(w) &= X_1(w) Y_1(w)\\ \end{align}

So that the relationship becomes:

$$z[n] \longleftrightarrow Z(w) = Z_0(w) + e^{-jwN} Z_1(w) + e^{-jw 2N} Z_2(w) \tag{20}$$

defining the following signals each being $2N=1000$ points as;

\begin{align} z_0[n] &= \mathcal{F}^{-1} \{ X_0(w) Y_0(w) \} \\ z_1[n] &= \mathcal{F}^{-1} \{ X_0(w) Y_1(w) + X_1(w) Y_0(w) \} \tag{21}\\ z_2[n] &= \mathcal{F}^{-1} \{ X_1(w) Y_1(w) \} \\ \end{align}

We therefore have: $$z[n] = z_0[n] + z_1[n-N] + z_2[n-2N]$$

where the remainig parts are identical to the first approach.

• You said you're really only looking for about ten points. Do you know which ten coefficients you want? Or just the ten highest? Or something else? A few years ago an FFT algorithm came out of MIT for sparse problems; that is, problems in which only a small number of coefficients are nonzero (or above the noise floor). Although the algorithm targets problems that are sparse in the Fourier domain, it's possible it could be adapted to your problem. Perhaps. I don't have the reference on front of me now, but if you're interested I could look it up. – Rodney Price Aug 28 '17 at 14:46