A Gentle Introduction to Discrete Fourier Transform

Cover image
πŸ› A Project Cauchy Op-ed

Named after French mathematician Augustin-Louis Cauchy, Project Cauchy column is where I invite some of the HFI Programming club members to provide neat proofs or explanations about some number theory puzzles. I highly suggest that you read these articles with a pencil and paper so you can sketch things out and scribble solutions to exercises as you come across them. This week, Travor will explain how fundamentals of DFT arises from the heat equation, and, at the end of the day, there will also be my minimal implementation of DFT.

Fourier transform is a commonly used techniques in many different fields. In mathematics, people use Fourier transform to solve differential equations, and in signal processing it is used to study the seasonality of time series. To prevent misuse, we had better understand theoretical background of Fourier transform before diving into applications. Prepare some papers and pencils, our adventure will begin from the heat equation.

The Heat Equation

In the field of partial differential equations (PDEs), heat equation is defined as follows

βˆ‚uβˆ‚t=kβˆ‡2u{\partial u\over\partial t}=k\nabla^2u

where βˆ‡2=βˆ‡β‹…βˆ‡\nabla^2=\nabla\cdot\nabla is the Laplacian operator. To prevent too much digression, we only want to focus on its one-dimensional form:

βˆ‚βˆ‚tu(x,t)=kβˆ‚2βˆ‚x2u(x,t)(1){\partial\over\partial t}u(x,t)=k{\partial^2\over\partial x^2}u(x,t)\tag1

with the initial condition and the boundary condition being u(x,0)=f(x)u(x,0)=f(x) and u(0,t)=u(L,t)=0u(0,t)=u(L,t)=0

We will not go into too much detail about this equation since this is not the topic of this blog.

Solving One-dimensional Heat Equation

Since the conditions and (1) are linear and homogeneous, we are allowed to perform separation of variables. That is, we set

u(x,t)=X(x)T(t)u(x,t)=X(x)T(t)

Plugging this new definition into (1), we get

Tβ€²T(t)=kXβ€²β€²X(x){T'\over T}(t)=k{X''\over X}(x)

By the properties of this equation, we observe that any partial derivatives taken on this equation will give zero, so both side must equal to a constant. In order for this constant to be physically meaningful, we set it negative:

1kTβ€²T(t)=Xβ€²β€²X(x)=βˆ’Ξ»2\frac1k{T'\over T}(t)={X''\over X}(x)=-\lambda^2

which allows us to separate this PDE into two ordinary differential equations (ODEs):

{Tβ€²(t)=βˆ’Ξ»2kT(t)Xβ€²β€²(x)=βˆ’Ξ»2X(x)\begin{cases} T'(t)&=-\lambda^2kT(t) \\ X''(x)&=-\lambda^2X(x) \end{cases}

and using techniques learned from calculus, we obtain a special solution for u(x,t)u(x,t):

u(x,t)=Aeβˆ’Ξ»2ktsin⁑(Ξ»x+Ο†)u(x,t)=Ae^{-\lambda^2kt}\sin\left(\lambda x+\varphi\right)

Now, if we were to plug in the boundary conditions of this problem, we get

un(x,t)=Aeβˆ’Ξ»n2ktsin⁑(Ξ»nx)u_n(x,t)=Ae^{-\lambda_n^2kt}\sin(\lambda_nx)

with Ξ»n\lambda_n being nΟ€/Ln\pi/L. Finally, by the linearity of (1), we get its general solution:

u(x,t)=βˆ‘n=1∞un(x,t)=βˆ‘n=1∞Aneβˆ’Ξ»n2ktsin⁑(Ξ»nx)u(x,t)=\sum_{n=1}^\infty u_n(x,t)=\sum_{n=1}^\infty A_ne^{-\lambda_n^2kt}\sin(\lambda_nx)

We are almost there. Plugging in the initial condition gives

f(x)=βˆ‘n=1∞Ansin⁑(Ξ»nx)(2)f(x)=\sum_{n=1}^\infty A_n\sin(\lambda_nx)\tag2

All left is to determine the coefficients. To begin with, we consider the following integral

Im,n=∫0Lsin⁑(λmx)sin⁑(λnx)dxI_{m,n}=\int_0^L\sin(\lambda_mx)\sin(\lambda_nx)\mathrm dx

For m=nm=n, we have

Im,m=∫0Lsin⁑2(Ξ»mx)dx=∫0L1βˆ’cos⁑(2Ξ»mx)2dx=L2\begin{aligned} I_{m,m} &=\int_0^L\sin^2(\lambda_mx)\mathrm dx \\ &=\int_0^L{1-\cos(2\lambda_mx)\over2}\mathrm dx \\ &=\frac L2 \end{aligned}

but when m≠nm\ne n, we can use the fact that

sin⁑αsin⁑β=cos⁑(Ξ±βˆ’Ξ²)βˆ’cos⁑(Ξ±+Ξ²)2\sin\alpha\sin\beta={\cos(\alpha-\beta)-\cos(\alpha+\beta)\over2}

to get

Im,n=12∫0L[cos⁑(Ξ»mβˆ’nx)βˆ’cos⁑(Ξ»m+nx)]dx=0I_{m,n}=\frac12\int_0^L[\cos(\lambda_{m-n}x)-\cos(\lambda_{m+n}x)]\mathrm dx=0

Combining these results, we get

∫0Lsin⁑(λmx)sin⁑(λnx)dx=L2δmn\int_0^L\sin(\lambda_mx)\sin(\lambda_nx)\mathrm dx=\frac L2\delta_{mn}

which, applied to (2), gives a formula to determine AnA_n:

An=2L∫0Lf(x)sin⁑(λn)dx(3)A_n=\frac2L\int_0^Lf(x)\sin(\lambda_n)\mathrm dx\tag3

Heat Equation and Fourier Series

As shown in (2), the solution to the heat equation requires us to determine the coefficients of a trigonometric series, and this type of series is known as the Fourier series to acknowledge Joseph Fourier for his contribution in the study of heat equation. Particularly, Fourier hypothesizes all functions defined on an arbitrary [T0,T0+T][T_0,T_0+T] can be represented as superposition of sinusoids:

f(x)=A0+βˆ‘n=1∞Ansin⁑(Ξ»nx+Ο†n)(4)f(x)=A_0+\sum_{n=1}^\infty A_n\sin(\lambda_nx+\varphi_n)\tag4

with Ξ»n=2Ο€n/T\lambda_n=2\pi n/T. However, (4) is not a good version of Fourier series to work with, so we may consider applying angle-sum identity of sine function:

f(x)=A0+βˆ‘n=1∞An[sin⁑φncos⁑(Ξ»nx)+cos⁑φnsin⁑(Ξ»nx)]=A0+βˆ‘n=1∞[ancos⁑(Ξ»nx)+bnsin⁑(Ξ»nx)]\begin{aligned} f(x) &=A_0+\sum_{n=1}^\infty A_n[\sin\varphi_n\cos(\lambda_nx)+\cos\varphi_n\sin(\lambda_nx)] \\ &=A_0+\sum_{n=1}^\infty[a_n\cos(\lambda_nx)+b_n\sin(\lambda_nx)] \end{aligned}

Hence, we are able to develop similar means in (3) to get

A0=1T∫T0T0+Tf(x)dxan=2T∫T0T0+Tf(x)cos⁑(λnx)dxbn=2T∫T0T0+Tf(x)sin⁑(λnx)dx\begin{aligned} A_0&=\frac1T\int_{T_0}^{T_0+T}f(x)\mathrm dx \\ a_n&=\frac2T\int_{T_0}^{T_0+T}f(x)\cos(\lambda_nx)\mathrm dx \\ b_n&=\frac2T\int_{T_0}^{T_0+T}f(x)\sin(\lambda_nx)\mathrm dx \end{aligned}

Because the formula for A0A_0 and ana_n are too similar, we can often see some texts used the following fashion to define Fourier series:

f(x)=a02+βˆ‘n=1∞[ancos⁑(2Ο€nxT)+bnsin⁑(2Ο€nxT)](5)f(x)={a_0\over2}+\sum_{n=1}^\infty\left[a_n\cos\left(2\pi nx\over T\right)+b_n\sin\left(2\pi nx\over T\right)\right]\tag5

Application of Fourier Series - the Basel Problem

Oftentimes, Fourier expansion of certain functions help us evaluate values of certain series. For instance, let's set T0=0T_0=0, T=2Ο€T=2\pi, and f(x)=x2f(x)=x^2; we have

a02=12Ο€βˆ«02Ο€x2dx=4Ο€23{a_0\over2}={1\over2\pi}\int_0^{2\pi} x^2\mathrm dx={4\pi^2\over3}

and for nβ‰₯1n\ge1:

an=1Ο€βˆ«02Ο€x2cos⁑(nx)dx=1Ο€[1nx2sin⁑(nx)∣02Ο€βˆ’1n∫02Ο€2xsin⁑(nx)dx]=1Ο€[1n22xcos⁑(nx)∣02Ο€βˆ’2n2∫02Ο€cos⁑(nx)dx]=4n2\begin{aligned} a_n &=\frac1\pi\int_0^{2\pi} x^2\cos(nx)\mathrm dx \\ &=\frac1\pi\left[\frac1nx^2\sin(nx)|_0^{2\pi}-\frac1n\int_0^{2\pi}2x\sin(nx)\mathrm dx\right] \\ &=\frac1\pi\left[{1\over n^2}2x\cos(nx)|_0^{2\pi}-{2\over n^2}\int_0^{2\pi}\cos(nx)\mathrm dx\right] \\ &={4\over n^2} \end{aligned}

which gives

x2=8Ο€33+βˆ‘n=1∞[4n2cos⁑(nx)+bnsin⁑(nx)]x^2={8\pi^3\over3}+\sum_{n=1}^\infty\left[{4\over n^2}\cos(nx)+b_n\sin(nx)\right]

set x=Ο€x=\pi, we have

Ο€2=4Ο€23+4βˆ‘n=1∞(βˆ’1)nn2\pi^2={4\pi^2\over3}+4\sum_{n=1}^\infty{(-1)^n\over n^2}

resulting in

βˆ‘n=1∞(βˆ’1)nn2=βˆ’Ο€212\sum_{n=1}^\infty{(-1)^n\over n^2}=-{\pi^2\over12}

which is an equivalent form of the Basel problem.

Fourier Series in a Complex Sense

Equation (5) still looks a bit complicated, but Euler's formula for trigonometric functions can help us simplify:

cos⁑θ=eiΞΈ+eβˆ’iΞΈ2sin⁑θ=eiΞΈβˆ’eβˆ’iΞΈ2i\cos\theta={e^{i\theta}+e^{-i\theta}\over2} \\ \sin\theta={e^{i\theta}-e^{-i\theta}\over2i}

As a result, we obtain the following series

f(x)=βˆ‘n∈Zcne2Ο€inx/T(6)f(x)=\sum_{n\in\mathbb Z}c_ne^{2\pi inx/T}\tag6

and by similar means in previous sections, we acquire the following formula for cnc_n:

cn=1T∫T0T0+Tf(x)eβˆ’2Ο€inx/Tdxc_n=\frac1T\int_{T_0}^{T_0+T}f(x)e^{-2\pi inx/T}\mathrm dx

Fourier Transform

Since sinusoids are periodic, Fourier series virtually serve to produce a trigonometric representations of periodic functions. Nonetheless, in mathematics most functions we study are not periodic, indicating that a stronger tool is needed.

If a function is aperiodic, then why not consider its period is the entire real line (βˆ’βˆž,∞)(-\infty,\infty)? In our case, set T0=βˆ’T/2T_0=-T/2, so that if we take the limit Tβ†’βˆžT\to\infty the Fourier series can express f(x)f(x) entirely.

In addition, define

f^(ΞΎ)=Tcn=βˆ«βˆ’T/2T/2f(x)eβˆ’2Ο€ixΞΎdx(7)\hat f\left(\xi\right)=Tc_n=\int_{-T/2}^{T/2}f(x)e^{-2\pi ix\xi}\mathrm dx\tag7

so that

f(x)=lim⁑Mβ†’βˆž1Tβˆ‘βˆ’MT/2≀n≀MT/2f^(nT)e2Ο€ix(n/T)f(x)=\lim_{M\to\infty}\frac1T\sum_{-MT/2\le n\le MT/2}\hat f\left(\frac nT\right)e^{2\pi ix(n/T)}

Now, let's consider dividing the summation region into subintervals:

βˆ’M≀⋯<ΞΎ0<ΞΎ1<⋯≀M-M\le\cdots<\xi_0<\xi_1<\cdots\le M

and require that the length of each subintervals to be 1/T1/T, then we have

f(x)=lim⁑Mβ†’βˆžβˆ‘nf^(ΞΎn)e2Ο€ixΞΎnΔξn(8)f(x)=\lim_{M\to\infty}\sum_n\hat f(\xi_n)e^{2\pi ix\xi_n}\Delta\xi_n\tag8

and if we take Tβ†’βˆžT\to\infty then (8) and (7) becomes

f(x)=βˆ«βˆ’βˆžβˆžf^(ΞΎ)e2Ο€ixΞΎdΞΎf(x)=\int_{-\infty}^\infty\hat f(\xi)e^{2\pi ix\xi}\mathrm d\xi f^(ΞΎ)=βˆ«βˆ’βˆžβˆžf(x)eβˆ’2Ο€ixΞΎdx\hat f(\xi)=\int_{-\infty}^\infty f(x)e^{-2\pi ix\xi}\mathrm dx

wherein f^(ΞΎ)\hat f(\xi) is the Fourier transform of f(x)f(x) and f(x)f(x) is called the inverse Fourier transform of f^(ΞΎ)\hat f(\xi). Moreover, because Fourier transform is regarded as continuous analog of Fourier series, the variable ΞΎ\xi is usually called the frequency, and f^(ΞΎ)\hat f(\xi) gives the amplitude of the sinusoid with frequency ΞΎ\xi.

For these reasons, in the field of signal processing, f(t)f(t) is often called the time domain, and f^(ΞΎ)\hat f(\xi) is often referred as the frequency domain so that Fourier transform became a tool to connect them.

Combining the two identities above gives us the Fourier inversion theorem:

f(t)=βˆ«βˆ’βˆžβˆže2Ο€iΞΎtβˆ«βˆ’βˆžβˆžf(x)eβˆ’2Ο€iΞΎxdxdΞΎf(t)=\int_{-\infty}^\infty e^{2\pi i\xi t}\int_{-\infty}^\infty f(x)e^{-2\pi i\xi x}\mathrm dx\mathrm d\xi

For simplicity in the exponential terms, we define the angular frequency Ο‰=2πξ\omega=2\pi\xi, resulting in an alternative version of Fourier inversion formula:

f(t)=12Ο€βˆ«βˆ’βˆžβˆžeiΟ‰tβˆ«βˆ’βˆžβˆžeβˆ’iΟ‰xf(x)dxdΟ‰f(t)={1\over2\pi}\int_{-\infty}^\infty e^{i\omega t}\int_{-\infty}^\infty e^{-i\omega x}f(x)\mathrm dx\mathrm d\omega

This implies a pair of Fourier transform equations using angular frequencies:

F(Ο‰)=βˆ«βˆ’βˆžβˆžeβˆ’iΟ‰tf(t)dt(9)F(\omega)=\int_{-\infty}^\infty e^{-i\omega t}f(t)\mathrm dt\tag9 f(t)=12Ο€βˆ«βˆ’βˆžβˆžeiΟ‰tF(Ο‰)dΟ‰(10)f(t)={1\over2\pi}\int_{-\infty}^\infty e^{i\omega t}F(\omega)\mathrm d\omega\tag{10}

In order for F(Ο‰)F(\omega) to exist, we ensure that f(t)f(t) is square-integrable (i.e. this integral ∫∣f∣2\int|f|^2 must converge).

Application of Fourier Transform - Airy's Equation

Oftentimes in quantum mechanics, Schrodinger equations in particular problems were simplified into Airy's equation:

yβ€²β€²βˆ’xy=0(11)y''-xy=0\tag{11}

If we differentiate (9) with respect to Ο‰\omega, we get

Fβ€²(Ο‰)=βˆ’iβˆ«βˆ’βˆžβˆžeβˆ’iΟ‰ttf(t)dtF'(\omega)=-i\int_{-\infty}^\infty e^{-i\omega t}tf(t)\mathrm dt

Moreover, if we differentiate (10) with respect to tt twice, we get

fβ€²β€²(t)=12Ο€βˆ«βˆ’βˆžβˆžeiΟ‰t[βˆ’Ο‰2F(Ο‰)]dΟ‰f''(t)={1\over2\pi}\int_{-\infty}^\infty e^{i\omega t}[-\omega^2F(\omega)]\mathrm d\omega

As a result, if we were to denote YY as the Fourier transform of yy, then the Fourier transform on (11) results in

βˆ’Ο‰2Yβˆ’iYβ€²=0-\omega^2Y-iY'=0

Using elementary techniques, we obtain the solution as follows:

Y(ω)=Y0eiω3/3Y(\omega)=Y_0e^{i\omega^3/3}

which gives us the spectrum of the square-integrable special solution to (11). Subsequently, inverse Fourier transform on both side gives

y=Y02Ο€βˆ«βˆ’βˆžβˆžeiΟ‰xeΟ‰3/3dΟ‰y={Y_0\over2\pi}\int_{-\infty}^\infty e^{i\omega x}e^{\omega^3/3}\mathrm d\omega

Because the complex exponential makes the solution look too intimidating, we may consider splitting the interval of integration to get a more friendly version:

y=Y02Ο€[βˆ«βˆ’βˆž0ei(Ο‰x+Ο‰3/3)dΟ‰βŸΟ‰β†¦βˆ’Ο‰+∫0∞ei(Ο‰x+Ο‰3/3)dΟ‰]=Y0Ο€βˆ«0∞ei(Ο‰x+Ο‰3/3)βˆ’eβˆ’i(Ο‰x+Ο‰3/3)2dΟ‰=Y0Ο€βˆ«0∞cos⁑(Ο‰x+Ο‰33)dΟ‰\begin{aligned} y &={Y_0\over2\pi}\left[\underbrace{\int_{-\infty}^0 e^{i(\omega x+\omega^3/3)}\mathrm d\omega}_{\omega\mapsto-\omega}+\int_0^\infty e^{i(\omega x+\omega^3/3)}\mathrm d\omega\right] \\ &={Y_0\over\pi}\int_0^\infty{e^{i(\omega x+\omega^3/3)}-e^{-i(\omega x+\omega^3/3)}\over2}\mathrm d\omega \\ &={Y_0\over\pi}\int_0^\infty\cos\left(\omega x+{\omega^3\over3}\right)\mathrm d\omega \end{aligned}

and if we specify Y0=1Y_0=1, we get the Airy's function of the first kind:

Ai⁑(x)=1Ο€βˆ«0∞cos⁑(Ο‰x+Ο‰33)dΟ‰\operatorname{Ai}(x)=\frac1\pi\int_0^\infty\cos\left(\omega x+{\omega^3\over3}\right)\mathrm d\omega

As (11) suggests, Airy's equation is a second order ODE, meaning there exists another branch of solutions that are linearly independent to Y0Ai⁑(x)Y_0\operatorname{Ai}(x). Consequently, Fourier transform only gives the square-integrable branch of the general solution.

Discretizing Fourier Transform

Although integrals look beautiful, it is not easy for computers to evaluate integrals, especially when it is integrating over the entire real line. As a result, we may consider discretizing the problem for computer use.

Sampling the Time Domain - Discrete-time Fourier Transform

In reality, we do not capture signals as a continuous flow but instead a discrete sequence, so discretizing the time domain will help us apply Fourier transform to discrete signals. Particularly, if we were to discretize the time-domain with sampling period 1/T1/T, we get

fT(t)=f(t)βˆ‘n=βˆ’βˆžβˆžΞ΄(tβˆ’nT)f_T(t)=f(t)\sum_{n=-\infty}^\infty\delta(t-nT)

The summation of delta functions are often called the Dirac comb.

Consequently, its Fourier transform became

F(Ο‰)=βˆ«βˆ’βˆžβˆžfT(t)eβˆ’iΟ‰tdt=βˆ«βˆ’βˆžβˆžf(t)βˆ‘n=βˆ’βˆžβˆžΞ΄(tβˆ’nT)eβˆ’iΟ‰tdt=βˆ‘n=βˆ’βˆžβˆžβˆ«βˆ’βˆžβˆžΞ΄(tβˆ’nT)f(t)eβˆ’iΟ‰tdt\begin{aligned} F(\omega) &=\int_{-\infty}^\infty f_T(t)e^{-i\omega t}\mathrm dt\\ &=\int_{-\infty}^\infty f(t)\sum_{n=-\infty}^\infty\delta(t-nT)e^{-i\omega t}\mathrm dt \\ &=\sum_{n=-\infty}^\infty\int_{-\infty}^\infty\delta(t-nT)f(t)e^{-i\omega t}\mathrm dt \\ \end{aligned}

Now, using the property of Dirac delta function that

βˆ«βˆ’βˆžβˆžΞ΄(tβˆ’a)g(t)dt=g(a)\int_{-\infty}^\infty\delta(t-a)g(t)\mathrm dt=g(a)

We deduce the formula for discrete-time Fourier transform (DTFT).

F(Ο‰)=βˆ‘n=βˆ’βˆžβˆžf(nT)eβˆ’iΟ‰nTF(\omega)=\sum_{n=-\infty}^\infty f(nT)e^{-i\omega nT}

wherein the inversion formula is the same as (10).

Sampling the Frequency Domain - Discrete Fourier Transform

Thanks to delta function, we turn integrals into summation. However, the inversion formula for DTFT is still in terms of integrals, so it is needed for us to manipulate this expression to feet real-life use. Because in reality we process finitely sized signal, the discrete signals are only finite sequences. As a result, let's consider DTFT on an discrete signal with length NN:

f(nT)={xn0≀n≀Nβˆ’10otherwisef(nT)= \begin{cases} x_n & 0\le n\le N-1 \\ 0 & \text{otherwise} \end{cases}

Then using this definition, we obtain the DTFT for discrete and finite signals:

F(Ο‰)=βˆ‘n=0Nβˆ’1xneβˆ’iΟ‰nT(12)F(\omega)=\sum_{n=0}^{N-1}x_ne^{-i\omega nT}\tag{12}

Since discrete signal only contains discrete periods, we only need to consider discrete frequencies. According to the definition of the sequence, its period is NTNT, as a result its discretized frequency can be written as

Ο‰=2Ο€kNT\omega={2\pi k\over NT}

which gives us

F(2Ο€kNT)=βˆ‘n=0Nβˆ’1xneβˆ’2Ο€ikn/NF\left(2\pi k\over NT\right)=\sum_{n=0}^{N-1}x_ne^{-2\pi ikn/N}

By the periodicity of complex exponential, we have that k∈{0,1,…,Nβˆ’1}k\in\{0,1,\dots,N-1\}. To recover the original signal, we apply the (10) with the new discretized version of the frequency domain:

xn=f(nT)=12Ο€βˆ«βˆ’βˆžβˆžeiΟ‰nTβˆ‘k=0Nβˆ’1F(Ο‰)Ξ΄(NΟ‰2Ο€βˆ’kT)dΟ‰=1Nβˆ‘k=0Nβˆ’1e2Ο€ikn/NF(2Ο€kNT)\begin{aligned} x_n &=f(nT) \\ &={1\over2\pi}\int_{-\infty}^\infty e^{i\omega nT}\sum_{k=0}^{N-1}F(\omega)\delta\left({N\omega\over2\pi}-{k\over T}\right)\mathrm d\omega \\ &=\frac1N\sum_{k=0}^{N-1}e^{2\pi ikn/N}F\left(2\pi k\over NT\right) \end{aligned}

Now, if we were to define the Fourier coefficients XkX_k as the amplitude of discrete frequencies F(2Ο€/NT)F(2\pi /NT), then we get the following equations for discrete Fourier transform (DFT) and its inversion:

Xk=βˆ‘n=0Nβˆ’1xneβˆ’2Ο€ink/N(13)X_k=\sum_{n=0}^{N-1}x_ne^{-2\pi ink/N}\tag{13} xn=1Nβˆ‘k=0Nβˆ’1Xke2Ο€ikn/N(14)x_n=\frac1N\sum_{k=0}^{N-1}X_ke^{2\pi ikn/N}\tag{14}

It can be verified that (13) and (14) conforms to the definition of DFT in numpy and scipy.

Application of DFT - Spectrum of the Greatest Common Divisor

It is well-known that DFT is frequently used in signal processing, and now let's try analyzing the spectrum of the greatest common divisor function:

Xk=βˆ‘m=0nβˆ’1gcd⁑(n,m)eβˆ’2Ο€imk/N=βˆ‘r=1ngcd⁑(n,r)eβˆ’2Ο€irk/N\begin{aligned} X_k &=\sum_{m=0}^{n-1}\gcd(n,m)e^{-2\pi imk/N} \\ &=\sum_{r=1}^n\gcd(n,r)e^{-2\pi irk/N} \end{aligned}

Since gcd⁑(n,r)\gcd(n,r) divides nn, we can consider summing the summands using the value of gcd⁑(n,r)\gcd(n,r):

Xk=βˆ‘d∣ndβˆ‘r=jdgcd⁑(j,n/d)=1eβˆ’2Ο€ijdk/N=βˆ‘d∣ndβˆ‘1≀j≀n/dgcd⁑(j,n/d)=1eβˆ’2Ο€ijk/(n/d)β‰œβˆ‘d∣ndβ‹…cn/d(k)\begin{aligned} X_k &=\sum_{d|n}d\sum_{r=jd\atop\gcd(j,n/d)=1}e^{-2\pi ijdk/N} \\ &=\sum_{d|n}d\sum_{1\le j\le n/d\atop\gcd(j,n/d)=1}e^{-2\pi ijk/(n/d)} \\ &\triangleq\sum_{d|n}d\cdot c_{n/d}(k) \end{aligned}

Now, let's turn our focus to cm(k)c_m(k):

cm(k)=βˆ‘1≀j≀mgcd⁑(j,m)=1eβˆ’2Ο€ijk/m=βˆ‘j=1meβˆ’2Ο€ijk/m⌊1gcd⁑(j,m)βŒ‹\begin{aligned} c_m(k) &=\sum_{1\le j\le m\atop\gcd(j,m)=1}e^{-2\pi ijk/m} \\ &=\sum_{j=1}^me^{-2\pi ijk/m}\left\lfloor1\over\gcd(j,m)\right\rfloor \end{aligned}

Now, apply Mobius inversion formula to the floor function, so we get

cm(k)=βˆ‘j=1meβˆ’2Ο€ijk/mβˆ‘d∣j,d∣mΞΌ(d)=βˆ‘d∣mΞΌ(d)βˆ‘j=qdeβˆ’2Ο€ijk/m=βˆ‘d∣mΞΌ(d)βˆ‘q=1m/deβˆ’2Ο€iqk/(m/d)=βˆ‘d∣mΞΌ(md)βˆ‘q=1deβˆ’2Ο€iqk/d\begin{aligned} c_m(k) &=\sum_{j=1}^me^{-2\pi ijk/m}\sum_{d|j,d|m}\mu(d) \\ &=\sum_{d|m}\mu(d)\sum_{j=qd}e^{-2\pi ijk/m} \\ &=\sum_{d|m}\mu(d)\sum_{q=1}^{m/d}e^{-2\pi iqk/(m/d)} \\ &=\sum_{d|m}\mu\left(\frac md\right)\sum_{q=1}^de^{-2\pi iqk/d} \end{aligned}

If d∣kd|k, the inner summation becomes dd, but when d∀kd\nmid k we have

βˆ‘q=1deβˆ’2Ο€iqk/d=eβˆ’2Ο€iq/dβ‹…eβˆ’2Ο€iqk/dβ‹…dβˆ’1eβˆ’2Ο€iq/dβˆ’1=0\sum_{q=1}^de^{-2\pi iqk/d}=e^{-2\pi iq/d}\cdot{e^{-2\pi iqk/d\cdot d}-1\over e^{-2\pi iq/d}-1}=0

As a result, we obtain a nice expression for cm(k)c_m(k):

cm(k)=βˆ‘d∣m,d∣kΞΌ(md)dc_m(k)=\sum_{d|m,d|k}\mu\left(\frac md\right)d

In particular, if k=1k=1 this expression becomes the mobius function. Accordingly, plugging k=1k=1 in to the XkX_k expression gives

X1=βˆ‘d∣nΞΌ(nd)d=cn(0)=βˆ‘1≀j≀ngcd⁑(j,n)=11=Ο†(n)X_1=\sum_{d|n}\mu\left(\frac nd\right)d=c_n(0)=\sum_{1\le j\le n\atop\gcd(j,n)=1}1=\varphi(n)

Eventually, we discover that the Fourier transform of the greatest common divisor function gives us the Euler's totient function. In addition, because Ο†(n)\varphi(n) is real, we can take the real component on the expression of X1X_1, yielding

Ο†(n)=βˆ‘r=1ngcd⁑(n,r)cos⁑(2Ο€rn)\varphi(n)=\sum_{r=1}^n\gcd(n,r)\cos\left(2\pi r\over n\right)

Summary

In this article, we began our exploration on a physics concepts. Then after we derived the continuous form of Fourier transform, we apply it to solving differential equation. In addition, we also discover application of discrete Fourier transform in number theory. Hence, we can see that Fourier transform has applications in different areas of mathematics because in the derivation we brought together different areas of mathematics: mathematical physics, mathematical analysis, and number theory!

Implementation

If you have no prior knowledge in MVC and find it difficult to understand the maths behind DFT. You can always refer to Jason's implementation of DFT below.

C

void dft(double complex *X, const size_t N) {
    double complex tmp[N];
    for (size_t i = 0; i < N; ++i) {
        tmp[i] = 0;
        for (size_t j = 0; j < N; ++j) {
            tmp[i] += X[j] * cexp(-2.0 * M_PI * I * j * i / N);
        }
    }

    memcpy(X, tmp, N * sizeof(*X));
}

C++

template <typename Iter>
void dft(Iter X, Iter last) {
  const auto N = last - X;
  std::vector<complex> tmp(N);
  for (auto i = 0; i < N; ++i) {
    for (auto j = 0; j < N; ++j) {
      tmp[i] += X[j] * exp(complex(0, -2.0 * M_PI * i * j / N));
    }
  }
  std::copy(std::begin(tmp), std::end(tmp), X);
}

Python

def dft(X):
    N = len(X)
    temp = [0] * N
    for i in range(N):
        for k in range(N):
            temp[i] += X[k] * exp(-2.0j * pi * i * k / N)
    return temp

In this function, we define n to be a set of integers from 00 to Nβˆ’1Nβˆ’1 and arrange them in a column. We then set k to be the same thing, but in a row. This means that when we multiply them together, we get a matrix, but not just any matrix! This matrix is the heart to the transformation itself!

M = [1.0+0.0im  1.0+0.0im           1.0+0.0im          1.0+0.0im;
     1.0+0.0im  6.12323e-17-1.0im  -1.0-1.22465e-16im -1.83697e-16+1.0im;
     1.0+0.0im -1.0-1.22465e-16im   1.0+2.44929e-16im -1.0-3.67394e-16im;
     1.0+0.0im -1.83697e-16+1.0im  -1.0-3.67394e-16im  5.51091e-16-1.0im]

It was amazing to us when we saw the transform for what it truly was: an actual transformation matrix! That said, the Discrete Fourier Transform is slow -- primarily because matrix multiplication is slow, and as mentioned before, slow code is not particularly useful.

So what was the trick that everyone used to go from a Discrete Fourier Transform to a Fast Fourier Transform?

Recursion! More about that next week.

β—€ MU-GAN RevisitedAnalytic Continuation of the Riemann Zeta Function β–Ά