What's The Frequency, Joseph?

| 2 Comments

As crucial as fractals undoubtedly were for the t-shirt industry[1], the complex numbers have many important mathematical implications, the Fourier transform being one of the most significant. For a function \(f\) a common definition is
\[ F(k) = \int_{-\infty}^{+\infty} e^{-2\pi\,i\,k\,x} f(x) \, \mathrm{d}x \]
An important property of the Fourier transform is that we can recover the original function with
\[ f(x) = \int_{-\infty}^{+\infty} e^{2\pi\,i\,x\,k} F(k) \, \mathrm{d}k \]
To get a feel for what this inverse transform does it is helpful to express the exponential as a complex number
\[ \begin{align*} f(x) &= \int_{-\infty}^{+\infty} \left(\cos(2\pi x k) + i\,\sin(2\pi x k)\right) F(k) \, \mathrm{d}k\\ &= \int_{-\infty}^{+\infty} F(k) \cos(2\pi x k) \mathrm{d}k + i \int_{-\infty}^{+\infty} F(k) \sin(2\pi x k) \mathrm{d}k \end{align*} \]
Derivation 1: The Fourier Transform Of A Convolution
Given the convolution
\[ h(x) = \int_{-\infty}^{+\infty}\,f(t)\,g(x-t)\,\mathrm{d}t \]
its Fourier transform is given by
\[ \begin{align*} H(k) &= \int_{-\infty}^{+\infty} e^{-2\pi\,i\,k\,x} h(x) \, \mathrm{d}x\\ &= \int_{-\infty}^{+\infty} e^{-2\pi\,i\,k\,x} \int_{-\infty}^{+\infty}\,f(t)\,g(x-t)\,\mathrm{d}t \, \mathrm{d}x \end{align*} \]
If, for a given \(t\), we define \(y = x-t\) we have \(x=y+t\) and \(\mathrm{d}y=\mathrm{d}x\) and consequently
\[ \begin{align*} H(k) &= \int_{-\infty}^{+\infty} e^{-2\pi\,i\,k\,(y+t)} \int_{-\infty}^{+\infty}\,f(t)\,g(y)\,\mathrm{d}t \,\mathrm{d}y\\ &= \int_{-\infty}^{+\infty} e^{-2\pi\,i\,k\,t}\,f(t)\,\mathrm{d}t\,\int_{-\infty}^{+\infty}\,e^{-2\pi\,i\,k\,y}\,g(y)\,\mathrm{d}y\\ &= F(k) \times G(k) \end{align*} \]
Just as the derivative is the limit of a difference, its inverse, the integral, is the limit of a sum and so this can be considered analogous to a Taylor series representation of a function, but in terms of trigonometric functions of the argument rather than powers of it. If we think of the function \(f\) as a signal of some sort, say an audio stream, then the Fourier transform is consequently related to the frequency spectrum of that signal.

The Fourier transform has many useful properties that can simplify calculations. For example, given the convolution, \(h\), of a pair of functions \(f\) and \(g\)
\[ \begin{align*} h(x) &= (f \ast g)(x)\\ & = \int_{-\infty}^{+\infty} f(t) \, g(x-t) \, \mathrm{d}t \end{align*} \]
and their Fourier transforms \(H\), \(F\) and \(G\), we have
\[ H(k) = F(k) \times G(k) \]
as proven in derivation 1.

Derivation 2: The Fourier Transform Of The Derivative
The Fourier transform of the derivative of a function \(f\) is
\[ \int_{-\infty}^{+\infty} e^{-2\pi\,i\,k\,x} f^\prime(x) \, \mathrm{d}x \]
By integration by parts
\[ \int u\,v^\prime = [uv] - \int u^\prime v \]
If we define
\[ \begin{align*} u(x) &= e^{-2\pi\,i\,k\,x}\\ v^\prime(x) &= f^\prime(x) \end{align*} \]
then, noting that \(f(\pm \infty)=0\), we have
\[ \int_{-\infty}^{+\infty} e^{-2\pi\,i\,k\,x} f^\prime(x) \, \mathrm{d}x\\ =\left[e^{-2\pi\,i\,k\,x} f(x)\right]_{-\infty}^{+\infty} - \int_{-\infty}^{+\infty} -2\pi\,i\,k\,e^{-2\pi\,i\,k\,x} f(x) \, \mathrm{d}x\\ =2\pi\,i\,k \int_{-\infty}^{+\infty} e^{-2\pi\,i\,k\,x}\,f(x)\,\mathrm{d}x\\ =2\pi\,i\,k\,F(k) \]
For another example, given a function \(f\) for which
\[ \begin{align*} f(-\infty) &= 0\\ f(+\infty) &= 0 \end{align*} \]
the Fourier transform of the derivative of \(f\) is given by
\[ 2 \pi \, i \, k \, F(k) \]
as proven in derivation 2.

By repeatedly applying this identity we find that the Fourier transform of the \(n^{th}\) derivative, \(f^{(n)}\), is given by
\[ (2 \pi \, i \, k)^n \, F(k) \]
This is an extremely useful result for solving differential equations for which we wish to find functions that satisfy some relationship between them and their derivatives.

If the function \(f\) only returns real numbers, then the Fourier transform satisfies the relation
\[ F(-k) = F(k)^\ast \]
where the asterisk denotes the conjugate of a complex number, or in other words that
\[ \begin{align*} \Re(F(-k)) &= \phantom{-}\Re(F(k))\\ \Im(F(-k)) &= -\Im(F(k)) \end{align*} \]
where \(\Re\) and \(\Im\) represent the real and imaginary parts of their argument.

The Discrete Fourier Transform

The Fourier transform is not particularly well suited to calculation with computers since, even if we program in the rules of integration, many integrals do not have explicit solutions and of those that have, it is not always a trivial matter to find them. For these reasons, the usual tactic for solving integrals with computers is to approximate them numerically.
Figure 1: Approximating The Integrand

The simplest algorithm for doing this involves approximating the function we wish to integrate, known as the integrand, with one composed of adjacent regions of equal width in which it takes the value of the integrand at the start of that region, as illustrated in figure 1.

Noting that the integral is equal to the area under (or minus the area above for negative regions) the graph of the function, the integral of the approximate function is simply the sum of the areas (both positive and negative) of the rectangles defined by those regions. We then use the integral of the approximate function as an approximation of the integral of the original function
\[ \begin{align*} \int_a^b f(x) \mathrm{d}x &\approx \sum_{n=0}^{N-1} \frac{b-a}{N} \times f\left(a + n \times \frac{b-a}{N}\right)\\ &=\frac{b-a}{N} \times \sum_{n=0}^{N-1} f\left(a + n \times \frac{b-a}{N}\right)\\ &=\Delta \sum_{n=0}^{N-1} f\left(a + n\Delta\right) \end{align*} \]
where \(\sum\) is the summation sign and \(N\) is the number of regions.

Approximating the Fourier transform in the same fashion we have
\[ \begin{align*} F(k) &= \int_{-\infty}^{+\infty} e^{-2\pi\,i\,k\,x} f(x) \, \mathrm{d}x\\ &\approx \Delta \sum_{n=-\infty}^{+\infty} e^{-2\pi\,i\,k\,n\Delta} f(n\Delta) \end{align*} \]
Now an infinite sum like this isn't really any more convenient than the integral it approximates so, on the face of it, it might not look like we've made any headway. However, if we note that an exponential of a purely imaginary number is trigonometric in nature and that consequently adding an integer multiple of \(2 \pi i\) to its argument yields the same result, we have
\[ \begin{align*} e^{-2\pi\,i\,k\,\left(n + \frac{j}{\Delta}\right)\Delta} &= e^{-2\pi\,i\,k\,n\Delta -2\pi\,i\,k\,j}\\ &= e^{-2\pi\,i\,k\,n\Delta} \end{align*} \]
when \(j\) and \(k\) are integers.
If we choose
\[ \Delta = \frac{1}{N} \]
for some positive integer \(N\), we can rearrange the approximation
\[ \begin{align*} F(k) &\approx \frac{1}{N} \sum_{n=0}^{N-1} e^{-2\pi\,i\,k\,\frac{n}{N}} \sum_{j=-\infty}^{+\infty} f\left(\frac{n}{N} + j\right)\\ & = \frac{1}{N} \sum_{n=0}^{N-1} e^{-2\pi\,i\,k\,\frac{n}{N}} f_n \end{align*} \]
The sum in this approximation is known as the discrete Fourier transform
\[ F_k = \sum_{n=0}^{N-1} e^{-2\pi\,i\,k\,\frac{n}{N}} f_n \]
Now both \(f_n\) and \(F_n\) are periodic, satisfying
\[ \begin{align*} f_{n+j\,N} &= f_n\\ F_{k+j\,N} &= F_k \end{align*} \]
for any integer \(j\); the former since every integer multiple of \(N\) is included in the sum and so adding another won't change its value, and the latter since the index \(k\) only appears in the exponential term and hence the observation we exploited to simplify the approximation applies again.
It might seem a little restrictive that the points at which we evaluate \(f\) are entirely governed by the choice of \(N\), but we can always define a new function
\[ g(x) = f(a\,x + b) \]
and compute that discrete Fourier transform of that instead.

The discrete Fourier transform has similar properties to the Fourier transform itself. For example, it is related to a sample of the frequency spectrum and has the inverse transform
\[ f_n = \frac{1}{N}\sum_{k=0}^{N-1} e^{2\pi\,i\,\frac{n}{N}\,k}F_k \]

The Fast Fourier Transform

As it stands, computing the discrete Fourier transform is a fairly expensive \(O(N^2)\) operation; each of the \(N\) unique terms of \(F_k\) is a sum involving all \(N\) terms of \(f_n\). Fortunately, with a little further algebraic trickery we can significantly improve performance[2]. Specifically, if \(N\) is even and we split the sum into even and odd indexed terms, we have
\[ \begin{align*} F_k &= \sum_{n=0}^{\frac{1}{2}N-1} e^{-2\pi\,i\,k\,\frac{2n}{N}} f_{2n} + \sum_{n=0}^{\frac{1}{2}N-1} e^{-2\pi\,i\,k\,\frac{2n+1}{N}} f_{2n+1}\\ &= \sum_{n=0}^{\frac{1}{2}N-1} e^{-2\pi\,i\,k\,\frac{2n}{N}} f_{2n} + e^{-2\pi\,i\,k\,\frac{1}{N}} \sum_{n=0}^{\frac{1}{2}N-1} e^{-2\pi\,i\,k\,\frac{2n}{N}} f_{2n+1}\\ &= \sum_{n=0}^{\frac{1}{2}N-1} e^{-2\pi\,i\,k\,\frac{n}{N/2}} f_{2n} + e^{-2\pi\,i\,k\,\frac{1}{N}} \sum_{n=0}^{\frac{1}{2}N-1} e^{-2\pi\,i\,k\,\frac{n}{N/2}} f_{2n+1}\\ \end{align*} \]
The first sum is the discrete Fourier transform of the even indexed terms and the second that of the odd indexed terms. Denoting them as \(F^\mathrm{e}_k\) and \(F^\mathrm{o}_k\) respectively, we can rewrite \(F_k\) as
\[ F_k = F^\mathrm{e}_k + e^{-2\pi\,i\,k/N}\,F^\mathrm{o}_k \]
Because the even and odd transforms have \(\frac{1}{2}N\) terms each, their periodicity implies that
\[ \begin{align*} F^\mathrm{e}_{k+\frac{1}{2}N} &= F^\mathrm{e}_k\\ F^\mathrm{o}_{k+\frac{1}{2}N} &= F^\mathrm{o}_k \end{align*} \]
Furthermore, we have
\[ \begin{align*} e^{-2\pi\,i\,(k+\frac{1}{2}N) / N} &= e^{-2\pi\,i\,k/N\,-\,\pi\,i}\\ &= e^{-2\pi\,i\,k/N}\,e^{-\pi\,i}\\ &= -e^{-2\pi\,i\,k/N} \end{align*} \]
Putting these together yields
\[ \begin{align*} F_k &= F^\mathrm{e}_k + e^{-2\pi\,i\,k/N}\,F^\mathrm{o}_k\\ F_{k+\frac{1}{2}N} &= F^\mathrm{e}_k - e^{-2\pi\,i\,k/N}\,F^\mathrm{o}_k \end{align*} \]
So we can express the discrete Fourier transform in terms of two transforms of length \(\frac{1}{2}N\) for a computational cost of some constant multiple of
\[ 2 \times \left(\frac{1}{2}N\right)^2 = \frac{N^2}{2} \]
plus some lower order terms, effectively halving the effort required to compute it. If we choose \(N\) to be some power of two, we can repeat this process until we are calculating the discrete Fourier transforms of single values, which are trivially equal to those values, for a total cost of \(O(n \ln n)\).
This is known as the fast Fourier transform and, given how great an improvement in performance it represents, is pretty much the de facto approach for approximating the Fourier transform.

ak.fft

Rather than require the user to supply input data that is some power of two in length we shall append zeros until they are of such length, or at least the implementations of the Fourier transform and its inverse will behave as if we had. We'll consequently need to calculate the smallest power of two that is greater or equal to the length of the input data, which we do with the nextPowerOfTwo helper function given in listing 1.

Listing 1: nextPowerOfTwo
function nextPowerOfTwo(n) {
 var n2 = 1;
 while(n2<n) n2 *= 2;
 return n2;
}

If we ignore the constant factor of \(\frac{1}{N}\), the only difference between the discrete Fourier transform and its inverse is the sign of the constant in the exponential. It therefore makes sense to define another helper function that performs the recursive calculation for either sign, as provided in listing 2.

Listing 2: fft
function fft(F, F0, f, f0, n, d, sTwoPi) {
 var n2, d2, F1, f1, j, theta, twiddle, F0j, F1j;

 if(n===1) {
  F[F0] = f0<f.length ? ak.complex(f[f0]) : ak.complex(0);
 }
 else {
  n2 = n/2;
  d2 = d*2;
  F1 = F0+n2;
  f1 = f0+d;

  fft(F, F0, f, f0, n2, d2, sTwoPi);
  fft(F, F1, f, f1, n2, d2, sTwoPi);

  for(j=0;j<n2;++j) {
   theta   = sTwoPi * j/n;
   twiddle = ak.complex(Math.cos(theta), Math.sin(theta));

   F0j = F[F0+j];
   F1j = ak.mul(twiddle, F[F1+j]);

   F[F0+j] = ak.add(F0j, F1j);
   F[F1+j] = ak.sub(F0j, F1j);
  }
 }
}

We'll use the argument sTwoPi to switch between using \(-2\pi\) and \(+2\pi\) in the exponent for the transform and its inverse respectively. Note that we're explicitly constructing the exponential in terms of trigonometric functions of theta since we know that its argument has no real part. Finally, in the terminating condition when n equals one, in which we perform the transforms for a single value, we convert the input to an ak.complex and implicitly pad it with zeros.
Implementing the fast Fourier transform and its inverse in terms of these helper functions is a trivial matter, as illustrated in listing 3.

Listing 3: ak.fft And ak.fftInv
ak.fft = function(f) {
 var n = nextPowerOfTwo(f.length);
 var F = new Array(n);

 fft(F, 0, f, 0, n, 1, -2*ak.PI);
 return F;
};

ak.fftInv = function(F) {
 var n = nextPowerOfTwo(F.length);
 var f = new Array(n);

 fft(f, 0, F, 0, n, 1, 2*ak.PI);
 f.forEach(function(x, i, f){f[i] = ak.div(x, n);});
 return f;
};

Note that we divide the elements of f by n after calling fft in the inverse transform.
The implementations of these functions can be found in FFT.js.

Using ak.fft

Program 1 illustrates the use of ak.fft, plotting a graph of the input and a bar chart of the output.

Program 1: Using ak.fft

The default implementation of f is a simple sinusoid of frequency one, but the result of ak.transform has two bars rather than one. This is a consequence of the function returning only real numbers and hence having a transform that satisfies
\[ F(-k) = F(k)^\ast \]
and therefore
\[ |F(-k)| = |F(k)^\ast| = |F(k)| \]
where the horizontal bars denote the magnitude, as calculated by ak.abs. Taking this, and the periodicity of the result of the transform, into account, the rightmost bar corresponds to a sinusoid of frequency minus one.
The example I really wanted to use was to allow you to select an audio file, use an audio object to load it and then use ak.fft to calculate the frequency spectrum of a sample of it. Unfortunately the audio object doesn't work like that, so the above will have to suffice.

Note that we're using another new UI object; the bar chart. This differs from most charts in that the x-axis bounds don't represent values, but are instead used to determine the number of bars with a calculation along the lines of

ak.round(bounds[1][0]-bounds[0][0])

This is because the positions of the bars are determined only by indices of the input array which, for the same reason, contains the heights of the bars rather than [x, y] coordinates.

References

[1] Our Imaginary Friends, www.thusspakeak.com, 2014.

[2] Cooley, J.W. & Tukey, J.W., An algorithm for the machine calculation of complex Fourier series, Mathematics of Computation, Vol. 19. Pages 297-301., 1965.

2 Comments

Nice. In Derivation 2 there is a type. Define v(x) = f(x) and not f'(x).

Leave a comment