Fourier Transform - Part 1

There are plenty of literature explaining Fourier Transform. However, at least most of it is fairly terse and, as I feel, are difficult to understand. I will try to give a very simple and easy to follow introduction. I will explain the idea in three different ways, giving three different ways to look at the concept. The first way is through the lens of polynomial interpolation. The second is the way of actually visualizing how the waves interfere to produce the function, and the third is seeing functions as vectors. Let us start with polynomial interpolation This would be a series of blogs, the first starting this one. We will cover the interpolation view for the discrete Fourier transform in this article.

Polynomial Interpolation

Let us say we have an unknown polynomial \(f(x)\), that we know would evaluate to \(y_i\) at the points \(x_i\) for \(i\) between \(1\) and some integer \(n\). For convience of not having to reiterate fringe cases, we will assume that \(n\) is a pretty large number than, say bigger than \(10\). We want to find this polynomial \(f(x)\) so that \(y_i=f(x_i)\) for all such \(i\). When can we do it? If we expand the polynomial \(f(x)\) as \(f(x)=\sum_{j=0}^{n-1}a_j x^j\) as an \(n-1\) degree polynomial, we can get \(n\) linear equations like \(\sum_{j=0}^{n-1} a_jx_i^j = y_i\). Notice that in the linear equations, \(x_i\) and \(y_i\) are known and the coefficients \(a_j\) are the unknows that we have to solve for. There will be \(n\) such linear equations in the variables \(a_j\), which we can then solve to find the coefficients \(a_j\). However, whether the set of linear equations are solvable depends on whether the equations are linearly independent. While we can prove that that would be the case in here, there is a simpler way to find the solution.

Suppose we can somehow find some polynomials \(f_i(x)\) for each \(i\) such that it evaluates to \(1\) at \(x_i\) and evaluates to \(0\) everywhere else. That is to say \(f_i(x_i)=1\) and \(f_i(x_j)=0\) for all \(j \neq i\) where \(1 \le j \le n\). Then we can say \(f(x)=\sum_{i=1}^n y_i f_i(x)\). Why? because for any evaluation point \(x_i\), only \(f_i(x_i)\) evaluates to \(1\) and for all the other \(f_j(x)\) for \(j\neq i\) evalutes to \(0\). So, \(f(x)\) as the sum of \(y_if_i(x)\) must evaluate to \(y_i\).

How do we get such a polynomial \(f_i(x)\)? Well, first thing is that it has to evaluate to \(0\) at every point \(j \neq i\) for \(1\le j \le n\). So, all those \(x_j\) are the roots of \(f_i(x)\). The simplest way to have such a polynomial \(g_i(x)\) is \(g_i(x) = \prod_{j\neq i} (x-x_j)\). This polynomial is only zero at those roots by construction. The only thing remaining is making the value evaluate to \(1\) at \(x_i\). This is very simple, we simply divide by whatever value it takes at \(x_i\). Notice that it cannot be zero at \(x_i\) because we already know all the zeros. So, we simply define \(f_i(x) = g_i(x)/g_i(x_i)\). Or explanding, we get the following – \[ f(x) = \sum_{i=1}^n y_i f_i(x) = \sum_{i=1}^n y_i \frac{g_i(x)}{g_i(x_i)} = \sum_{i=1}^n y_i \frac{\prod_{j\neq i} (x-x_j)}{\prod_{j\neq i} (x_i-x_j)}\] A few things to notice. Each \(f_i\) is a degree \(n-1\) polynomial. Since \(f(x)\) is a linear sum of all the \(f_i(x)\), it can have a maximum degree of \(n-1\), as we would like. Also note that we do not need any other complicated conditions to prove for this to work out, we already covered why the denominator cannot be zero and that’s all we need. This is called Lagrange interpolation method.

Discrete Fourier Transform

Some further simplication happens when we fix the evaluation points at the roots of \(1\). Let me explain. We are going to assume that \(f(x)\) is a polynomial with complex coefficients and complex \(x\). So we are trying to interpolate \(f(x)\) evaluating at complex \(x_i\) to give a complex \(y_i\). They can be real numbers as well since real numbers are also complex numbers, i.e. the set of real number is a subset of the set of complex numbers. But the point is they can be all complex. Now, we want \(x_i\) to be one of the \(n^{th}\) roots of \(1\). The number \(1\) has \(n\) such roots, one of those is \(1\) itself. We can write them as \(e^{2 \pi \hat{i}/n }\) where \(\hat{i}\) is the square root of \(-1\) [so we do not confuse it with our index \(i\)] . If we plot then on the complex plain, the roots of \(1\) make vectors separated by an equal angle of \(2\pi/n\). For our convinience, we can assign \(\omega=e^{2\pi\hat{i}/n}\), and then the roots of \(1\) are \(\omega^0, \omega^1,\cdots, \omega^{n-1}\).

Now we want to simplify the Langrange interpolation as a polynomial. This is a lot of algebraic manipulation, but we can simplify our lives by using a few symmetric expressions. We first define \(\sigma_k\) to be the sum of all combinations of \(k\) of those roots. So there are \(n\) roots, we can choose \(k\) of them in \(^nC_k\) ways. So this sum has \(^nC_k\) terms. For example \(\sigma_1=\sum_{i=0}^{n-1} \omega^i\), \(\sigma_2= \sum_{i=0}^{n-2}\sum_{j=i+1}^{n-1} \omega^i \omega^j\), \(\sigma_3= \sum_{i=0}^{n-3}\sum_{j=i+1}^{n-2}\sum_{k=j+1}^{n-1} \omega^i \omega^j \omega^k\) and so on. Note that all \(\omega^i\) are roots of \(1\), so they are roots of the polynomial \(x^n-1\). Therefore, \(x^n-1 = \prod_{i=0}^{n-1} (x-\omega^i)\). Comparing the coeffficients on both sides, we can see that \(\sigma_i\) is the coefficient of \(x^{n-i}\) in the polynomial \(x^n-1\) multiplied by \((-1)^i\). That means \(\sigma_i=0\) for all \(1 \le i \le n-1\) and \(\sigma_n = (-1)^{n+1}\).

Now we consider another type of sum \(\sigma_{k,i}\) which has all the terms the same as \(\sigma_k\) except all the terms that have \(\omega^i\) it them. We can compute the value simply by substrating the sum of all the terms that have \(\omega^i\) in them from \(\sigma_k\). But for all the terms that do have \(\omega^i\), we can take that outside the parenthesis and multiply by the sum of all the terms of \(\sigma_{k-1}\) without any containing \(\omega^k\) (because we already took that out of the parentheses). Therefore, \(\sigma_{k,i} = \sigma_k - \omega^i\sigma_{k-1, i}\), unless \(k=1\), in which case \(\sigma_{1,i} = -\omega^i\) since the sum of all the roots is \(0\).

Now, we can compute the following -

1. \(\sigma_{1,i} = \sigma_1 - \omega^i = -\omega^i\)
2. \(\sigma_{2,i} = \sigma_2 - \omega^i\sigma_{1,i} = -\omega^i\sigma_{1,i}=\omega^{2i}\)
3. \(\sigma_{3,i} = \sigma_3 - \omega^i\sigma_{2,i} = -\omega^i\sigma_{2,i}=-\omega^{3i}\) \(\vdots\)

11. \(\sigma_{k,i} = \sigma_k - \omega^i\sigma_{k-1,i} = -\omega^i\sigma_{k-1,i}=(-1)^k\omega^{ki}\)

up to \(k=n-1\).

Now, let us consider \(g_i(x) = \prod_{j\neq i} (x-\omega^j)\) . A quick consideration gives that the coefficient of \(x^{n-k-1}\) in this polynomial is \((-1)^{k}\sigma_{k,i} = \omega^{ki}\) . And the coefficient of \(x^{n-1}\) is \(1\). Therefore,

\[ g_i(x) =x^{n-1}+\sum_{k=1}^{n-1} \omega^{ki} x^{n-1-k}=\sum_{k=0}^{n-1} \omega^{ki} x^{n-1-k} \] And, we can compute [Note that we use the fact that \(\omega^n=1\)] – \[ g_i(\omega^i) = \sum_{k=0}^{n-1} \omega^{ki} \omega^{i(n-k-1)} = \sum_{k=0}^{n-1} \omega^{i(n-1)} = n\omega^{-i} \] Hence, \[ \begin{aligned} f_i(x) \\ &= g_i(x)/g_i(\omega^i) \\ &= \frac{1}{n\omega^{-i}}\sum_{k=0}^{n-1} \omega^{ki} x^{n-k-1} \\ &= \frac{1}{n\omega^{-i}}\sum_{l=0}^{n-1} \omega^{(n-l-1)i} x^{l} \textit{ [setting $l=n-k-1$]}\\ &=\frac{1}{n}\sum_{l=0}^{n-1} \omega^{(n-l)i} x^{l}\\ &= \frac{1}{n}\sum_{l=0}^{n-1} \omega^{-li} x^{l}\\ &= \frac{1}{n}\sum_{k=0}^{n-1} \omega^{-ki} x^{k} \end{aligned} \]

Therefore, \[ \begin{aligned} f(x) &= \sum_{i=1}^n y_i f_i(x) \\ &= \frac{1}{n}\sum_{i=1}^n y_i\sum_{k=0}^{n-1} \omega^{-ki} x^{k} \\ &= \frac{1}{n}\sum_{k=0}^{n-1}\sum_{i=1}^n y_i \omega^{-ki} x^{k} \\ &=\frac{1}{n}\sum_{k=0}^{n-1}x^{k} \sum_{i=1}^n y_i \omega^{-ki} \end{aligned} \] The coefficient of \(x^k\) in \(f(x)\) is therefore \[ \frac{1}{n} \sum_{i=1}^n y_i \omega^{-ki} \] We can express it as a function of \(k\) as \(Y_k = \frac{1}{n} \sum_{i=1}^n y_i \omega^{-ki}\). We call \(Y_k\) to be the discrete Fourier transform of values \(y_i\). If we want our \(y_i\) back given \(Y_k\), we simply evaluate the polynomial at \(\omega^i\), since the whole point of polynomial interpolation was for it to evaluate to \(y_i\) at \(\omega^i\). That gives - \[ y_i = \sum_{k=0}^{n-1} Y_k \omega^{ik} \] This is the inversion of the polynomial interpolation and is called the inverse discrete Fourier transform of \(Y_k\).

The frequency interpretation

Notice again that \(\omega^{-ik}\) rotates in the complex plane as \(k\) increases. The higher is the value of \(i\) the faster it rotates. In fact, it rotates by an angle of \(2 \pi i/n\) for an increase of \(1\) in the value of \(k\). This suggests that \(i\) is some kind of frequency measure, it is proportional to the speed of circling around the input. In the next article, I will expand on this in more detail with geometrical view allowing one to see exactly what is happening. For now, we take a break.