A Gentle Introduction to Discrete Fourier Transform
Posted December 25. 2020.
3023 words.
14 min read.
π 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
βtβuβ=kβ2u
where β2=ββ β is the Laplacian operator. To prevent too much digression, we only want to focus on its one-dimensional form:
βtββu(x,t)=kβx2β2βu(x,t)(1)
with the initial condition and the boundary condition being u(x,0)=f(x) and u(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)
Plugging this new definition into (1), we get
TTβ²β(t)=kXXβ²β²β(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:
k1βTTβ²β(t)=XXβ²β²β(x)=βΞ»2
which allows us to separate this PDE into two ordinary differential equations (ODEs):
{Tβ²(t)Xβ²β²(x)β=βΞ»2kT(t)=βΞ»2X(x)β
and using techniques learned from calculus, we obtain a special solution for u(x,t):
u(x,t)=AeβΞ»2ktsin(Ξ»x+Ο)
Now, if we were to plug in the boundary conditions of this problem, we get
unβ(x,t)=AeβΞ»n2βktsin(Ξ»nβx)
with Ξ»nβ being nΟ/L. Finally, by the linearity of (1), we get its general solution:
which, applied to (2), gives a formula to determine Anβ:
Anβ=L2ββ«0Lβf(x)sin(Ξ»nβ)dx(3)
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] can be represented as superposition of sinusoids:
with Ξ»nβ=2Ο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:
Oftentimes, Fourier expansion of certain functions help us evaluate values of certain series. For instance, let's set T0β=0, T=2Ο, and f(x)=x2; we have
Equation (5) still looks a bit complicated, but Euler's formula for trigonometric functions can help us simplify:
cosΞΈ=2eiΞΈ+eβiΞΈβsinΞΈ=2ieiΞΈβeβiΞΈβ
As a result, we obtain the following series
f(x)=nβZββcnβe2Οinx/T(6)
and by similar means in previous sections, we acquire the following formula for cnβ:
cnβ=T1ββ«T0βT0β+Tβf(x)eβ2Οinx/Tdx
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 (ββ,β)? In our case, set T0β=βT/2, so that if we take the limit Tββ the Fourier series can express f(x) entirely.
wherein f^β(ΞΎ) is the Fourier transform of f(x) and f(x) is called the inverse Fourier transform of f^β(ΞΎ). Moreover, because Fourier transform is regarded as continuous analog of Fourier series, the variable ΞΎ is usually called the frequency, and f^β(ΞΎ) gives the amplitude of the sinusoid with frequency ΞΎ.
For these reasons, in the field of signal processing, f(t) is often called the time domain, and f^β(ΞΎ) 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:
As a result, if we were to denote Y as the Fourier transform of y, then the Fourier transform on (11) results in
βΟ2YβiYβ²=0
Using elementary techniques, we obtain the solution as follows:
Y(Ο)=Y0βeiΟ3/3
which gives us the spectrum of the square-integrable special solution to (11). Subsequently, inverse Fourier transform on both side gives
y=2ΟY0βββ«ββββeiΟxeΟ3/3dΟ
Because the complex exponential makes the solution look too intimidating, we may consider splitting the interval of integration to get a more friendly version:
and if we specify Y0β=1, we get the Airy's function of the first kind:
Ai(x)=Ο1ββ«0ββcos(Οx+3Ο3β)dΟ
As (11) suggests, Airy's equation is a second order ODE, meaning there exists another branch of solutions that are linearly independent to Y0β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/T, we get
fTβ(t)=f(t)n=βββββΞ΄(tβnT)
The summation of delta functions are often called the Dirac comb.
Now, using the property of Dirac delta function that
β«ββββΞ΄(tβa)g(t)dt=g(a)
We deduce the formula for discrete-time Fourier transform (DTFT).
F(Ο)=n=βββββf(nT)eβiΟ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 N:
f(nT)={xnβ0β0β€nβ€Nβ1otherwiseβ
Then using this definition, we obtain the DTFT for discrete and finite signals:
F(Ο)=n=0βNβ1βxnβeβiΟnT(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 NT, as a result its discretized frequency can be written as
Ο=NT2Οkβ
which gives us
F(NT2Οkβ)=n=0βNβ1βxnβeβ2Οikn/N
By the periodicity of complex exponential, we have that kβ{0,1,β¦,Nβ1}. To recover the original signal, we apply the (10) with the new discretized version of the frequency domain:
Now, if we were to define the Fourier coefficientsXkβ as the amplitude of discrete frequencies F(2Ο/NT), then we get the following equations for discrete Fourier transform (DFT) and its inversion:
Eventually, we discover that the Fourier transform of the greatest common divisor function gives us the Euler's totient function. In addition, because Ο(n) is real, we can take the real component on the expression of X1β, yielding
Ο(n)=r=1βnβgcd(n,r)cos(n2Οrβ)
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
voiddft(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<typenameIter>voiddft(Iter X, Iter last){constauto 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
defdft(X):
N =len(X)
temp =[0]* N
for i inrange(N):for k inrange(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 0 to Nβ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!
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?