Convolutions and Fourier Transforms
Convolutions
Convolution is the process of
\[(x * y)_n = \sum_{m=0}^{N-1} x_m \cdot y_{n-m}\]x is a vector of size N, y, the kernel, can be any size, if the index is outside its range, the value would be 0. Generally n goes from 0 to N-1. I think. Using modulo arithmetic, a circular convolution can be performed, which can be restated using DFTs:
\[(x * y_N)_n = \sum_{m=0}^{N-1} x_m \cdot y_{n-m \textrm{ mod } N} = F^{-1}(X \cdot Y)_n\]For n = 0, and n = 1, the summations look like
\[ \begin{align*} & x_0 \cdot y_0 + x_1 \cdot y_{N-1} + \dots + x_{N-2} \cdot y_2 + x_{N-1} \cdot y_1 \\ & x_0 \cdot y_1 + x_1 \cdot y_0 + x_2 \cdot y_{N-1} + \dots + x_{N-1} \cdot y_2 \end{align*} \]Each time n is incremented, the y values are shifted right, or rotated, which is why the convolution is called "circular". For the final value N-1 the first y index is N-1, and it goes down to 0 at the end. Circular convolutions are useful because they can be performed as 2 forwards DFTs and one inverse, and we can use the an FFT for those. To use Fourier transforms to calculate convolutions, the convolution must first be converted to a circular convolution. So for, say n=1, the orignial convolution and it's circular convolution equivalent are:
\[ \begin{align*} & x_0 \cdot y_1 + x_1 \cdot y_0 + x_2 \cdot y_{-1} + \dots + x_{N-1} \cdot y_{-(N-2)} \\ & x_0 \cdot y_1 + x_1 \cdot y_0 + x_2 \cdot y_{N-1} + \dots + x_{N-1} \cdot y_2 \end{align*} \]The solution would appear to be to set index N-1 to the original value of index -1, and so on, however for n = N-1 the terms contain identical indices for the whole series, which implies that the y values are the same between the 2 convolutions at the same indices. Instead you can extend the circular convolution to a length M, such that for the same case the two convolutions are
\[ \begin{align*} & x_0 \cdot y_1 + x_1 \cdot y_0 + x_2 \cdot y_{-1} + \dots + x_{N-1} \cdot y_{-(N-2)} \\ & x_0 \cdot y_1 + x_1 \cdot y_0 + x_2 \cdot y_{M-1} + \dots + x_{N-1} \cdot y_{M-(N-2)} + \dots + x_{M-1} \cdot y_2 \end{align*} \]The lowest modulo index is M-N+1, which, to avoid collisions should be > N-1, giving the requirement M > sum of the sizes of the x and y series. To make the circular and linear convolutions equal, we should get rid of x indices > N-1, as they don't appear in the original, so set them to 0, and keep the rest of the x values the same. The y values from 0 to N-1 should be the same, for the cases n > size of y - 1, the result is outside the bounds of the original convolution (n iterated over the y vector) so it can be discarded. Which just leaves the negative indices, which simply have M added to the index, so e.g. -1 goes to M-1. So when we expand the vectors for the circular convolution, for indices 0 to ... L-1, if L-1 was the highest n we originally went up to, the values are the same as the linear case. For indices M-N+1 to M-1 the values are just the values of y at those indices minus M. Between L-1 and M-N+1, the values are irrelevant because they are multiplied by 0 for n less than L, and we discard the results when they are eventually rotated around to x index 0 for n > L-1. The value of M can be as large as we want, so it can be set to a power of 2 greater than N+L for use in a simple radix-2 fft. Once we pad the vectors out, the convolution is simple:
/* convolve 2 samples, assuming they have been zero padded to appropriate lengths. Results are in ax/ay. */
void fft_convolve(float *ax, float *ay, float *bx, float *by, unsigned int M) {
if(M&(M-1)) return; /* not power of 2 */
float t;
int n;
/* DFT */
fft_ctip(ax, ay, M);
fft_ctip(bx, by, M);
/* pointwise multiplication of DFTs with scaling for conjugate inverse */
for(n = 0; n < M; n++) {
t = (ax[n] * bx[n] - ay[n] * by[n]) / M;
ay[n] = - (ax[n] * by[n] + bx[n] * ay[n]) / M;
ax[n] = t;
}
/* iDFT */
BitReverse(ax, ay, M);
fft_ctip(ax, ay, M);
BitReverse(ax, ay, M);
return;
}
References
Prev: Cooley-Tukey algorithm
Up: Audio
- Uses javascript SyntaxHighlighter for syntax highlighting.
- Uses javascript MathJax for LateX rendering.
To the extent possible under law,
Craig Hughes
has waived all copyright and related or neighboring rights to
the text and source code in this example/tutorial.