31.2 Fourier analysis
In homage to “The Hobbit” let us call this section “Fourier Analysis, or there and back again”.
The convolution of two vectors and is defined as the vector
that is, the th member of is formed by multiplying all the members of and whose indices sum to .
Note that we define for all indices which are not part of the vector.
This is quite a useful operation, for example you may recognise it from the multiplication of polynomials where we have
Multiplication of polynomials
which should feel natural because every term in the product is formed by multiplying one term in with one term in , that is each term in the result is formed from
and our sum just collects all the and which sum to give .
Another example comes from probability, for example consider the probability that the sum of two sided dice is .
Probability of the sum of dice
Convolution is a very useful operation – it turns out that it is also an operation that can be computed efficiently. We will see that the periodic convolution (which we can write the general convolution just introduced in this section) can be calculated very quickly using the Fourier transform (in a later section).
31.2.2 Periodic convolution
Sometimes we encounter “periodic” sequences. For example the sequence
is periodic because in general
that is, it repeats every terms. The nice thing about periodic sequences is that although they are infinite, we only need to look at a finite portion of them to derive results which apply to the rest of the sequence.
Let be -periodic integer sequences. Then the periodic convolution of and is an -periodic integer sequence given by
We can also define the periodic convolution of two vectors (and not just infinite sequences) with (almost) the same formula
where refers to modulo arithmetic (i.e. , ) which is explored in more detail in Chapter 22.2.
We are interested in discrete convolution because it is useful in computing interesting things, and we are interested in periodic convolution because of these two attributes (which will we will show later)
the discrete convolution of two sequences can be computed using only the periodic convolution
the periodic convolution of two vectors can be computed very quickly using the Fourier transform
The first one we will show right now; let be an dimensional vector and let be an dimensional vector. We will then seek to find some modified vectors and a natural number such that
We know that is an dimensional vector. Therefore our periodic convolution must be a periodic sequence with period at least . It turns out that if we guess we will get lucky.
and we want to obey the formula
that is, the th element of will be
and will just be with zeroes padded on the end.
31.2.3 The Fourier series
A Taylor series is an infinite power series such that
and one way to think about this is that as we add more and more terms to our series we get a better and better approximation for . In the discrete (i.e. finite) case we can fit an -degree polynomial to points. In the infinite case we can perfectly represent a function (provided it is infinitely differentiable).
Another way to interpret this is as the projection of a function onto the set of polynomials.
But why restrict ourselves to polynomials? We can also represent functions in terms of trigonometric functions. For example, we can write a function as
We are essentially writing as a polynomial in and . It helps to remember that we can write and as th degree polynomials in and respectively (see Example 14.7.2 for the case for if you’re not sure).
To find the coefficients we can integrate with respect to or . For example to find the coefficient we can integrate both sides with respect to between and , i.e.
and it turns out that and are orthogonal when considering integration as the inner product. That is to say, most of the integrals are zero.
Let us now turn to the discrete case with equally spaced points. Instead of a function , let us assume that we have a series of points
for which we want to find a set of points
We use in lieu of a purely real number formulation because working over the complex numbers makes the algebra much easier.
If we consider Equation 31.86 we can form an coefficient matrix and try to solve the linear system
If like me you find this a bit confusing I would recommend trying out the cases for some small values of and .
In the general case we will use
as the value for the th equation. Note that is the principal th root of unity, a fact that will come in very handy later.
This means that the coefficient at in will be . For example the matrix will look like
Of course, on the face of it this isn’t very helpful because computing the solution to takes which isn’t great. It turns out that the inverse of is very easy to compute, and we can shave this down to , and then using a clever divide-and-conquer algorithm we will be able to compute the discrete Fourier transform in which is really fast!
31.2.5 …and back again
It turns out that we can go back the other way using the matrix
and we can compute in . This is actually a very slow algorithm, and in Section 31.2.7 it turns out that it is possible to cleverly compute this in time (which is very fast and essentially almost linear).
is the inverse of
Proof let be natural numbers, then
in the case where we therefore find that and thus
Now we turn to the case where , in this case
Note that the sum of the th roots of unity is , so we obtain here. A proof of this can be found in Theorem 14.8.1. Another nice proof uses the sum for the geometric series
which we can apply to to get
31.2.6 The convolution theorem
The convolution theorem
The convolution theorem tells us that in the “frequency domain” we can compute the Fourier transform using element-wise multiplication
where is the element-wise multiplication operator.
Consider which we can write as
Here is a slightly tedious proof, but in Section 31.2.10 there is an alternative proof of this using “circulant” matrices.
Let us consider
and now we can re-order the sums giving us
which is not fun, but if you enjoy this kind of thing bully for you!
31.2.7 “Geht es besser?” (The Fast Fourier Transform)
This is a great German phrase which seems to be the basis (pun unintended) for a lot of discoveries in algorithms. It means “can we do better”, and the question is usually used rhetorically when introducing a faster solution to an algorithmic problem.
It turns out that we are carrying out a lot of extra work when computing the Fourier transform. Currently we are carrying out a matrix-vector computation to compute it, which we can do in time (very slow). This section outlines the Cooley–Tukey FFT (Fast Fourier Transform) algorithm which achieves a much better time complexity of .
Let us consider the -dimensional case. Here we want to compute such that
The matrix is written in terms of the principal th root of unity, which has the important attribute that for
We will use this when computing
by splitting the sum
The first part of the sum satisfies
and the second part of the sum satisfies
Therefore we can write in terms of
This algebra takes a moment to digest, so I would read it a few times to get the idea.
We can now analyse the time complexity of this algorithmic approach. Let denote the number of operations to compute the DFT for dimension , then at each stage we
compute of all the odd coefficients, which takes time
compute of all the even coefficients, which also takes time
combine the results by adding them together and doing a bit of multiplication which we can do in time
therefore in aggregate our recurrence looks like
where is some cutoff and is a constant. Iterating the recurrence, we find that
which is a very fast algorithm. Hence it is called the “Fast Fourier Transform”.
31.2.8 Filtering with Fourier
31.2.9 Circulant matrices
A circulant matrix is a matrix which has this structure
in words; we “rotate” the first by one to obtain each of the subsequence rows.
Let us consider the eigenvalues of this special matrix. When we compute eigenvalues we want quite a lot of things to stay the same; we know that every row contains the same values (if in a different order) and addition is commutative, so if we just take all the values (i.e. use the vector ) we will find that
if we let then we can write
and therefore is an eigenvalue.
We can find the rest of the eigenvalues without too much trouble by thinking about rotations. We know that in general the th roots of unity have the rather amazing property that for any we essentially rotate them when multiplying through by
this rotation hopefully feels quite a bit like a circulant matrix. In fact consider what happens when we multiply a circulant matrix by the vector composed of the th roots of unity;
Although it looks like we are stuck in algebra soup, we are on the hunt for eigenvalues, so we must try to pull out our vector. It turns out that because the roots of unity rotate so nicely, we are in luck and we find that
provided that . I am getting tired of typing so if you would like to verify this it is possible using Equation 31.130 and the above.