In this article, I show the mechanics of the discrete Fourier transform (DFT) and circular convolution (within the context of time series analysis). The goal of this post is to present these two methods in a practical yet non-superficial way. We will only deal with real, finite length time series $\left\{a_t\right\}$. The primary reference for this article is “Wavelet Methods for Time Series Analysis” by Percival and Walden.

Discrete Fourier Transform (DFT)

The discrete Fourier transform takes a time series from the time domain to the frequency domain. The transform allows us to see the frequency content of a time series.

Suppose we are given a time series

Its discrete Fourier transform is the complex sequence $\left\{A_k\right\}$ given by:

Note that $A_k$ is associated with frequency $f_k \equiv k / N$.

import numpy as np
from matplotlib import pyplot as plt

def naive_finite_fourier_transform(a):
    # a can be a list or 1D numpy array
    N = len(a)
    A = []

    for k in range(N):

        A_k = 0

        for t in range(N):
            A_k += a[t] * np.exp((-1j * 2 * np.pi * t * k)/N)
        
        A.append(A_k)

    A = np.asarray(A)

    f = [k/N for k in range(N)]
    f = np.asarray(f)

    return A, f

We can reconstruct $\left\{a_t\right\}$ from its DFT $\left\{A_k\right\}$ using the inverse DFT:

def naive_inverse_finite_fourier_transform(A):
    # a can be a list or 1D numpy array
    N = len(A)
    a = []

    for t in range(N):

        a_t = 0

        for k in range(N):
            a_t += A[k] * np.exp((1j * 2 * np.pi * t * k)/N)
        
        a.append((1/N)*a_t)

    a = np.asarray(a)

    t = list(range(N))
    t = np.asarray(t)

    return np.real(a), t

Example

a = np.random.random(100)

A, f = naive_finite_fourier_transform(a)

plt.plot(f[1:], np.abs(A[1:])**2)
plt.axvline(x=0.5, color="red", linestyle = "--")
plt.show()

a1, _ = naive_inverse_finite_fourier_transform(A)

np.allclose(a, a1)

True

Circular Convolution/Filtering

Circular convolution (also called circular filtering) is a simple idea, but is cumbersome to read when described with math. For this reason, I will first show how circular convolution works using pictures before presenting any mathematics.

Suppose we have two time series, $\{a_t\}$ and $\{b_t\}$, both of length $N=10$. I won’t specify the value of the time series here (I’ll leave them to the reader’s imagination). I will only show the time indices of each time series in the pictures below. The indices of the first time series $\{a_t\}$ will be shown in red as follows:

and the indices of the second series $\{b_t\}$ will be shown in black to distinguish it from $\{a_t\}$:

The first step is to reverse the order of $\{a_t\}$, which makes its indices look as follows:

Next, we concatenate $\{b_t\}$ front to back and perform a sequence of sliding dot products between $\{a_t\}$ and the newly lengthened version of $\{b_t\}$ as follows:

The above picture shows how the first element of the result of a circular convolution between $\{a_t\}$ and $\{b_t\}$ is constructed:

The next element can be constructed as follows:

and here’s element $N-2$ of the resultant circular convolution time series:

In general, the result of a convolution between two time series (both of length $N$) is a resultant time series also of length $N$.

Mathematically, the $t^{th}$ element of the resultant time series from a circular convolution between two time series $\{a_t\}$ and $\{b_t\}$ (both of length $N$) is:

where $t = 0, 1, \dots, N-1$.

Whenever the subscript on $b$ is negative (which occurs whenever $u$ is greater than $t$), we define $b_{-1}=b_{N-1} ; b_{-2}=b_{N-2}$ and so on.

More compactly (at the expense of readability), we can also write the $t^{th}$ element of the resultant time series of a circular convolution between two time series $\{a_t\}$ and $\{b_t\}$ as:

where $\bmod$ is the modulo operator.

Here’s a coded implementation:

def naive_circular_convolution(a, b):
    assert len(a) == len(b)

    N = len(a)
    c = []
    for t in range(N):
        c_t = 0
        for u in range(N):
            c_t += a[u]*b[(t-u) % N]
        c.append(c_t)
    return np.asarray(c)

Circular Convolution Using DFT

A key fact is that a circular convolution can be computed using the DFT and inverse DFT. Here’s how:

1. Compute the DFT of $\{a_t\}$, denoted by $\{A_k\}$
2. Compute the DFT of $\{b_t\}$, denoted by $\{B_k\}$
3. Take the inverse DFT of $\{A_{k} B_{k}\}$

def naive_circular_convolution_dft(a, b):

    A, _ = naive_finite_fourier_transform(a) # STEP 1
    B, _ = naive_finite_fourier_transform(b) # STEP 2

    # STEP 3
    C = A*B
    c, _ = naive_inverse_finite_fourier_transform(C)
    
    return c

Example

a = np.random.random(100)
b = np.random.random(100)

c = naive_circular_convolution(a, b)
c1 = naive_circular_convolution_dft(a, b)
np.allclose(c, c1)

True

Final Words

An interesting point to mention is that circular convolution can be performed using a filter $\{a_t\}$ that is infinite in length (while the time series $\{b_t\}$ is finite in length). This is made possible by a process known as “filter periodization”. I will make a post in the future about this interesting topic.

Thanks for reading!