[Research] Fast Fourier Transform，FFT

1. Introduction

Fast Fourier transform (FFT) are various algorithms that compute discrete Fourier transform (DFT) of a sequence, or its inverse (iDFT) that shares the same principle as DFT. Fourier analysis converts a signal from its original domain (often time or space) to a representation in the frequency domain and vice versa. FFT is widely used for many applications in engineering, science, and mathematics. More specifically, in the area of digital signal processing, Fourier transform plays a quite important role.

Actually, Fourier transform has different forms, continuous Fourier transform as well as discrete Fourier transform. When processing the signal, the computer can only work with a discrete signal. That is why I mainly focus on Discrete Fourier transform here. As for the content of this post, I am working with the inference of the DFT taking the polynomial computation as an example first. Then I will give the formula definition of DFT followed by explaining the link between the polynomial computation and DFT. From this, you can have a good understanding of the mathematical logic of the Fourier transform.

2. Background

1) Polynomial

Assume A(x) is a n-1th order polynomial with a coefficient vector (a0, a1, a2,..., an-1). This polynomial is determined uniquely by this coefficient vector and is expressed in the form of:

Formula 1

Let's change our perspective towards this polynomial. Assume the coefficient vector is unknown, in order to figure out these coefficients, we need to establish n independent equations related to these coefficients to get these coefficients via solving equations. If (x0, y0), (x1, y1), (x2, y2),..., (xn-1, yn-1) are n points (xi is different from each other) that make yi=A(xi) established, these points are able to be employed to establish the parameter equations. Once those parameter equations solved, this polynomial is uniquely determined.

Formula 2

Those are two ways to uniquely display a polynomial, one called coefficient representation and the other called point-value representation.

2) complex number

An imaginary number is a number that can be written as a real number, b, multiplied by the imaginary unit i which is defined by its property i^2=-1. And a complex number is formed by a real number, a, added to an imaginary number bi, a+bi. a and b are called real part and imaginary part of the complex number, respectively.
About the operation rules:

• addition: (a+bi)+(c+di)=(a+c)+(b+d)i
• subtraction: (a+bi)-(c+di)=(a-c)+(b-d)i
• multiplication: (a+bi)(c+di)=(ac-bd)+(bc+ad)i
• division: (a+bi)/(c+di)=(ac+bd)/(c2+d2) +((bc-ad)/(c2+d2))i

The complex number can be expressed by a complex plane. Two coordinate systems we often use to display the complex number, the vertical coordinate system and the polar coordinate system. a complex number on the vertical coordinate system can be displayed by z=x+yi, expressed by z=re^iθ on the polar coordinate system.

Figure 1. complex numbers on the vertical coordinate system

Figure 2. a complex number on the polar coordinate system

Particularly, on the polar coordinate system, z is called Magnitude/ Absolute Value/ Modulus. θ is called Argument (arg for short) and it is the angle between the positive direction of the real part axis and line direction consisted of the origin and point z. r is called modular length. we have:

Formula 3

2.1) Root of unity

Roots of unity are used in many branches of mathematics, and are especially important in number theory, the theory of group characters, and the discrete Fourier transform. The most intuitive explanation of the root of unity is: an nth root of unity, where n is a positive integer (i.e. n = 1, 2, 3, …), is a number z satisfying the equation: Z^n=1. Since n has the different integer values, so a group of roots of unity includes n items. The root of unity is often expressed use w as shown in figure 3-2.
On a complex plane, a unit circle is divided into n (=16) parts evenly (figure 3.1).

Figure 3-1. n even division of a unit circle

Figure 3-2. group of root of unities

Depending on the Euler formula (formula. 4 (1), i stands for imaginary unit, θ is an argument), if let ωn=e^(-2πi/n), the group of n-order root of unities can be presented by formula 4 (4).

Formula. 4

2.2) Features of the unit root

Formula. 5

Formula. 6

Formula. 7

Formula. 8

3. FFT in Polynomial Computation Problem

The polynomial computation part I am going to show you is a good example that can tell you all the details about FFT. Actually, the polynomial is the carrier of the FFT of any signals. So it is very important to understand the process of FFT in polynomial computation.

1) The function of FFT in the polynomial computation

Before we explore the FFT, I want to introduce a polynomial multiplication problem to help you understand the function of it.
So problem occurs in polynomial multiplication.

Figure 4

If we use the way shown in figure 4, the time complexity is O(n^2). As the polynomial has two representation fashions displayed in section 2. We can use the point-value pairs to represent the polynomial as shown in Figure 5. Randomly sample 2-n pairs of points for A(x) and B(x), respectively. Then we can get 2-n pairs of points for C(x). So C(x) can be determined uniquely. And the next step is we just transfer the point-value representation to the coefficient representation.

Figure 5

The effect of randomly sampling is usually low, so another question is to consider the sampling strategy of those 2-n pairs of points. Then someone proposed we could use the roots of unity that transform the real plane to the complex plane. What is the feature of roots of unity?
Before coming to the next sub-part, I first show you the flow chart of the FFT used in polynomial multiplication as shown in Figure 6. The process of DFT is also a process to generate the point-value pairs to transform the coefficient representation to point-value representation.  Figure 5 shows the process of how to realize from A(x), B(x) to C(x). Finally, a iDTF is used to transform the point-value representation to the coefficient representation.

Figure 6

2) Generate the point-value pairs using the root of unity

In order to show you how the root of unity is used in generating point-value pairs, I arbitrarily take a polynomial (coefficient a is known) as an example. We can divide it into two sub-polynomials, the even part as well as the odd part according to the power of x. The lost parts can be replaced by zero. So we have figure 7.

Figure 7

In section 2.2, I already show you the feature of the root of unity. We replace the x with the root of unity as shown in figure 8. For this example, {(w0n, A(w0n)), (w1n, A(w1n)), ..., (wn-1n, A(wn-1n))} is the point-value set we want. [We also call (A(w0n), A(w1n), ..., A(wn-1n)) the DFT of (a0, a1, a2,..., an-1)].

Figure 8

Actually, we can generate different point-value pairs for another polynomial, e.g. B(x) using the same roots of unity. [we also have (B(w0n), B(w1n), ..., B(wn-1n)) for the DFT of (b0, b1, b2,..., bn-1)]. Back to our question, to compute C(x)=A(x)B(x), we utilize these two point-value vectors to generate the point-value set for C(x), {(w0n, A(w0n)B(w0n)), (w1n, A(w1n)B(w1n)), ..., (wn-1n, A(wn-1n)B(wn-1n))}.
This step we generate the point-value sets for A(x), B(x) and C(x), then we need to transform the point-value representation of C(x) to the coefficient representation.

3) Point-value representation to coefficient representation

C(x)=A(x)B(x) and {(w0n, A(w0n)B(w0n)), (w1n, A(w1n)B(w1n)), ..., (wn-1n, A(wn-1n)B(wn-1n))} (expressed by (y0, y1, y2, ..., yn-1)) is the point-value set of C(x). Then we have the formula. 9 to get ck (c0, c1, c2,..., cn-1).

Formula. 9

Now we need to prove whether the ck is the coefficient of C(x).

The idea is that we first assume (a0, a1, a2,..., an-1) the coefficient of C(x) to generate the point-value pairs, and then we use those point-value pairs to recover the coefficient by using Formula 9 to show the relationship between (a0, a1, a2,..., an-1) and ck.

Figure 9

Through figure 9, we obtain the relationship between (a0, a1, a2,..., an-1) and ck, which shows (c0, c1, c2,..., cn-1)=n*(a0, a1, a2,..., an-1). Finally, we successfully transfer the point-value representation to the coefficient representation.

4) Summary

Let's link these three steps together. In order to increase the efficiency of the polynomial computation, we replace the coefficient representation with the point-value representation. Then being able to avoid a low-efficiency randomly sampling, we introduce the group of roots of unity as the independent variables to generate the point-value set. At last, we make a transformation between the point-value representation and the coefficient representation. The reason why the roots of unity are used here is mainly two: increasing the sampling efficiency and reducing the computation load in computing A(w) and recovering the coefficient because of the feature of roots of unity (section 2.2).    Actually, we call the process of using (a0, a1, a2,..., an-1) to calculate (A(w0n), A(w1n), ..., A(wn-1n)) as DFT. And we call the process of using (A(w0n), A(w1n), ..., A(wn-1n)) to recover (a0, a1, a2,..., an-1) as iDFT. The difference between them is using the root of unity to conduct the DFT and using the reciprocal of the root of unity to conduct iDFT (the result should be divided by n). 

4. DFT Model Analysis

1) DFT and iDFT Formula

From section 3, we have an overview of the process of DFT and iDFT. Here, a formula definition is given in Formula 10. Assume f(x) is the discrete sequence of a signal input in the time domain. F(u) is the result of DFT in the frequency domain. The DFT and iDFT can be expressed in the form:

Formula. 10

where:

• f(x) is the discrete sequence of a signal input in the time domain, e.g. the accelerometer signal as the input.
• N is the size of f(x). Usually, N is the integer power of 2, e.g. 32, 64, 128...
• x=0, 1, 2, ..., N-1. It has two meanings: 1) x-th input in the discrete sequence of f(x); 2) the x-th power of the independent variable in a polynomial; 
• F(u) is the result of DFT in the frequency domain. It is a complex number, (a+bi)n or (a+bj)n
• u=0, 1, 2, ..., N-1. it stands for the u-order root of unity W^u.

From Formula 10, we can see that the Discrete Fourier Transform has a close relation with polynomial computation. DFT and iDFT are all concerning polynomial computation.

2) Frequency

Frequency is the number of occurrences of a repeating event per unit of time. But for a DFT, what is the meaning of the frequency to a DFT.
In section 2.1, I introduced the roots of unity. Depending on the definition, n-th root of unity is e^θnj. If let Δθ=2π/N that is the unit angle of the unit circle divided into N parts evenly, θn equals n*Δθ. For each θn, it stands the number of sectors being scanned by the rotation axis from 0 to θn. Therefore, given a fixed time, the θn scans n sectors and θm scans m sectors. The frequency is a fraction of the number of scanned sectors to whole sectors, e.g. fn=n/N.
Then, making a connection between the n-th input in the time domain and the e^θnj in the frequency domain, for each input in the time domain, their frequencies are determined by the number of their scanning sectors, respectively. That is why Fourier transform can generate a graph subject to the frequency when processing a time-domain signal. This process is similar to figure 9-1.

Figure 9-1.

Obviously, the frequency of n-th input is in the form of: where, n=[1, 2, ...., N], fs is the sampling frequence of the signal.

Formula. 11

3) Frequency-Magnitude graph

As the result of DFT is expressed in the form of the complex number, so the amplitude of z (=a+bi), the magnitude/modular length is |z|=sqrt(a^2+b^2), the argument (also called phrase) is sita=atan(y/x).

Figure 11, DFT-magnitude and frequency

Obviously, each magnitude has its corresponding sub-frequency, we can use a specific frequency to find its magnitude. And we can also use a magnitude to find the frequency.

5. Source Codes

import numpy as np
import matplotlib.pyplot as plt
x=np.linspace(0,2*np.pi,64)
#print(x)
wave=np.sin(x)
#print(wave)

transformed=np.fft.fft(wave) #DFT
plt.plot(np.abs(transformed)) 

plt.plot(np.fft.ifft(transformed))
plt.show()

Figure 11, DFT-magnitude and sub-frequency

Figure 12, the result of DFT

Figure 13. iDFT

5. Conclusion

This post is mainly working on the understanding and inference of FFT (DFT). Two most of important things are 1) understanding the implementation of DFT from the example of polynomial computation; 2) understanding how DFT generates the frequency-amplitude graph from a time-domain input (this is also the essence of Fourier Analysis).
Actually, Fourier transform is closely related to signal process, including signal generating,  propagation, FIR filter design and so on. besides, we also need an understanding of Convolution, which I will show you in my next post.

Reference

https://oi.men.ci/fft-notes/
https://zhuanlan.zhihu.com/p/31584464
https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5796454/
https://www.zhihu.com/question/20460630
https://www.zhihu.com/question/25023139
http://www.cnblogs.com/v-July-v/archive/2011/02/20/1983676.html
http://blog.miskcoo.com/2015/04/polynomial-multiplication-and-fast-fourier-transform
http://people.chu.edu.tw/~cclu/ComplexVariableAnalysis/Summary%20309.pdf
Each signal can be decomposed into numerous sin waves: https://www.youtube.com/watch?v=r18Gi8lSkfM
http://download.ni.com/evaluation/pxi/Understanding%20FFTs%20and%20Windowing.pdf