Table of Contents

New LaTeX commands declared in this cell. Expand to view. $$ \newcommand\abs[1]{\left|#1\right|} \newcommand\dx[2]{\frac{\partial#1}{\partial#2}} \newcommand\ddx[2]{\frac{\partial^2#1}{\partial#2^2}} \newcommand\divv[1]{\nabla\cdot\vect{#1}} \newcommand\curl[1]{\nabla\times\vect{#1}} \newcommand\grad[1]{\nabla\vect{#1}} \newcommand{\ket}[1]{\left|#1\right\rangle} \newcommand{\bra}[1]{\left\langle#1\right|} \newcommand{\braket}[1]{\langle#1\rangle} \newcommand{\aint}[0]{\int_{-\infty}^{\infty}} \newcommand{\vect}[1]{\mathbf{#1}} \newcommand{\order}[1]{\mathcal{O}{#1}} $$

Fourier Series

See discussions in:

Fourier Transform

Basics

The following discussion assumes evenly spaced steps in time $\Delta t$. The Fourier series of a function is its decomposition into individual sines and cosines:

$$ y(t) = \sum_{j=1}^n y_j \sin(2 \pi f_j t + \phi_j) \\ \text{With } _j = \frac{\omega_j}{2\pi} $$

.

If we take this to the continuum limit, summing over all frequencies,

$$ y(t) = \int_\infty^\infty Y(f) e^{-2\pi ift} df = \frac{1}{2\pi} \int_\infty^\infty Y(\frac{\omega}{2\pi}) ^{-i\omega t} d\omega \\ Y(f) = \int_\infty^\infty e^{2\pi i f t} dt$$

The two expressions are equivalent, though the $\omega$ version is more commonly seen in the context of quantum mechanics. The forms above are commonly called the forward Fourier transforms, though this has no deep meaning. The idea of going 'forward' likely results from our far more pedestrian experience of time as a dimension than frequency; we transform 'forward' to a new domain, and then inverse transform to return to 'normal'. This reasoning holds for the other common Fourier transforms as well; consider that position space is far more familiar to us than momentum space!

Another tool that will be useful is the Dirac delta function, the result of transforming forward, then back, and expecting your original function to return (a reasonable hope):

$$ \int_\infty^\infty e^{i\omega(t - t')} d\omega = 2\pi \delta(t-t')\\ $$

This function has the property of being zero everywhere except where t=t'. It has the following property when integrated with a function:

$$ \int y(t') \delta(t-t') dt' = y(t) $$

It thus has the inverse units of your integration variable.

Discrete FT

Consider discretizing the above definition of the Fourier transform:

$$ \text{For }t = m \Delta t \text{, and} f_n = n/(N\Delta t): \\ y_m = \sum_{n=0}^{N-1} Y_n e^{-2\pi imn/N} \qquad Y_n =\frac{1}{N} \sum_{m=0} y_m e^{2 \pi i m n / N} $$

We can recreate the discrete version of a Dirac delta function with these definitions; the Kronecker delta function:

$$ \sum_{n=0}^{N-1}e^{2 \pi i n(m - m')/N} = N \delta_{m,m'} $$

Note that $\Delta t$ has disappeared from our discrete FT. There are subtleties involving the amount of information contained in the DFT between real and complex cases.

The DFT as implented below loops through the time and frequency indices, m and n, to compute matrix elements $$\bra{n} M \ket{m}$$, where $M$ is the DFT Matrix. Since this matrix has $N\times N$ elements, the computation is of $\mathcal{O}(N^2)$, which will quickly become cumbersome. The Fast Fourier Transform is an algorithm much more efficient at computing the FT, which will be discussed in the next section.

Nyquist Frequency

The Nyquist critical frequency is defined as $$ f_c = \frac{1}{2\Delta t} $$ The simplest case to consider is a sine wave: to properly sample a sine wave, you must take at least 2 points per cycle. This means that you need to sample at a minimum of twice the maximum frequency of your data to preserve information below that frequency. In a real device, such as an experimental setup limited by the response of particular electrical cables, there will probably be natural frequency cutoffs. You need to sample at twice the max frequency you want. For example, human hearing cuts off at around 20kHz. To be able to properly sample audio recordings for future listening (on headphones, for instance), you need to sample at twice this frequency. It's no accident that audio CDs sample at 44.1kHz!

If the datatype is not set to complex, you will be discarding the phase of your DFT; $$ e^{i\theta} = \cos{\theta} + i \sin{\theta} $$ With only half of the information, you will introduce phase shifts into the Fourier transform. This can be easily tested by changing the functions above to read dtype=float.

FFT

The Fast Fourier Transform (FFT) is a ubiquitous algorithm for computing the Fourier transform of a dataset in $\mathcal{O}(N \log N) $ time. Compare the difference in the number of operations needed for the two cases.

$$ N = 8 = 2^3 \\ Y_n = \sum_{m=0}^{N-1} M_{n,m}Y_m \\ M_{nm} = e^{2 \pi i n m / N} \\ ex: Y_0 = \sum_{m=0}^{7} e^{2\pi i (0) m / N}y_m = y_0 + y_1 +...\\ \text{Split into even (e) & odd (o) m:} \\ Y_n = Y^e_n + w^n Y_n^o \qquad w = e^{2 \pi i / N} $$
$$ = \sum_{m'=0}^3 Y_{2m'} w^{2m'n} + w^n \sum_{m'=0}^3 Y_{2m'+1} w^{2m'n} \\ Y_n = Y^{e}_{n} + w^n y_n^o $$

For each level, we can recursively split the evens and odds again and get twice as many terms:

$$ Y^e_n = Y^{ee}_n + w^{2n} Y^{eo}_n \\ Y^o_n = Y^{oe}_n + w^{2n} Y^{oo}_n \\ $$

split until each $Y^{eoeoeoeoeoooeooeo...eoeooooeooooeeeeeoeoee}_n$, etc, term is a single element. For the n=8 example, this means three iterations of dividing the input array.

$$ Y^{ee}_n = Y^{eee}_{(0)} + w^{4n} Y^{eeo}_{(4)} \\ Y^{eo}_n = Y^{eoe}_{(2)} + w^{4n} Y^{eoo}_{(6)} \\ Y^{ee}_n = Y^{ooe}_{(3)} + w^{4n} Y^{ooo}_{(7)} \\ Y^{ee}_n = Y^{oee}_{(1)} + w^{4n} Y^{oeo}_{(5)} \\ $$

Recall what these Y terms are. Each is a broken down component of the initial array $Y_n$, $n \in [0,8).$ A string like eeo, eoe, etc, represents a reversed binary representation of the index to $Y_n$. $e\rightarrow 0, o \rightarrow 1$. If you flip the text string, you recover the binary representation of the index. An example is essential here.

$$ Y^{eoo}: eoo \rightarrow 011; \text{reversed, } 110 \rightarrow Y_6 \\ Y^{eoo} = Y_6 $$

So we can bit-reverse the index of the original array to find the new location of that element after the FFT algorithm begins working.

$$ Y_m = \sum_{n=0}^{N-1} y_n e^{-2\pi i m n / N} \\ $$ $$\order{N \log_2{N}}$$

Bit reversal

The point of bit reversing is a bookkeeping shortcut. Rather than splitting down the levels, sorting even/odd elements until they come out as a new array

$$[y_0, y_4, y_2, y_6, y_1, y_5, y_3, y_7]$$

There is a simple pattern here that is not obvious at first glance. We can recognize that, if we convert the indices to binary, with the number of bits corresponding to $ln_2(len(y))$, we get

$$[y_{000}, y_{100}, y_{010}, y_{110}, y_{001}, y_{101}, y_{011}, y_{111}]$$

And, if we reverse the order of each of these bits, the final result is

$$[y_{000}, y_{001}, y_{010}, y_{011}, y_{100}, y_{101}, y_{110}, y_{111}]$$

a slick trick that can save us time by instantly swapping our input array elements to the place that sorting would have put them already.

Recursive FFT

Below is a sketch of the recursive fft algorithm. Note that this does not utilize the bit reversal technique; deriving this method led to the recognition of bit reversal as a bookkeeping trick, which is implemented in Iterative FFT below.

def _recursive_fft(data):
    """recursively compute the fft of the input data."""
    N = int(len(data))
    if N == 1:
        return data # The FT of a single element is itelf
    w_N = np.exp(-2j * np.pi / N)
    w = 1.0 + 0.0j

    data_even = data[0::2] # starts from 0 and moves two elements until N-2
    data_odd = data[1::2] # starts from 1 and moves two elements until N-1

    y_e = recursive_fft(data_even)
    y_o = recursive_fft(data_odd)

    y = np.zeros([N], dtype=complex)
    for k in range(0, N//2):
        t = w * y_o[k] # butterfly operation
        y[k] = y_e[k] + t
        y[k+N//2] = y_e[k] - t
        w = w * w_N
        del t
    return y

Iterative FFT

The above recursive FFT is messy (I'm not sure where the extra frequency in the fft is coming from), and padding zeros is difficult without significantly more bookkeeping than I've done above. We will move from the bottom up in the tree describing the FFT, following the bit-reversal scheme.

def iterative_fft(cls, data):
    data = data.astype(complex)
    assert (np.log2(len(data)) % 1) == 0
    A = cls.reverse_array(data)
    N = int(len(A))
    S = int(np.log2(N))
    for s in range(1, S+1):
        m = int(2**s) # length of split array at each level
        w_m = np.exp(-2j * np.pi / m)
        for k in range(0, N, m): # by m
            w = 1.0 + 0.0j
            m_2 = (m // 2)
            for j in range(0, m_2):
                t = w * A[k + j + m_2]
                u = A[k + j]
                A[k + j] = u + t
                A[k + j + m_2] = u - t
                w = w * w_m
        return A

Bit reversal

def bit_reverse(num, bits):
    num = int(num)
    return int('{:0{s}b}'.format(num, s=bits)[::-1], 2)