Notes for Discrete Fourier Transform
We will use [] to denote superscripts, while () will generally be subscripts.
1. Grid Functions and Periodicity
We will present two equivalent views of what discrete periodic functions are. In different situations one encounters these different views so it is good to have thought about them.
1.A. Finite Grids and Grid Functions:
A grid of size n is a set of n equally spaced points on the line. We can label (number) the points of the grid as 0, 1, 2, ..., n-1 and we will refer to this grid as G[n] (read [n] as a superscript).
A grid function is a function that has values at the points of a grid. In other words a grid function is a set of n numbers, one assigned to each grid point. A grid function will then be a set of n numbers u = {u(0),u(1),...u(n-1)}. We will consider complex functions u where each u(i) is a complex number. Clearly real functions are a special case. The set of all grid functions on G[n] is the n-dimensional vector space C[n] of complex n-component vectors.
1.B. Periodic Grids and Periodic Functions:
A periodic grid of period n is an infinite set of points which repeat consecutively after every n. Thus we might describe it as a set of points labeled as: ...., n-2, n-1, 0, 1, 2, ...., n-2, n-1 ,0, 1, 2, ......
The best example is a clock: after you reach 12 you return to 1,2,..... Other names for a periodic grid are a ring or a torus.
A periodic function U on a periodic grid of period n is an infinite set of numbers, one for each grid point, which has the same value at the "same point". Thus there are only n distinct values: .., U(n-2), U(n-1), U(0), U(1), U(2), ..., U(n-2), U(n-1), U(0), U(1),U(2),...
In other words a periodic function on a grid G[n] of period n is an infinite set of numbers U={U(i)} with the property that:
U(i) = U(i pmod n ), for any i.
Here the pmod function is the modulo (remainder) function defined by:
Let i = cn + r, with c integer and 0<=r<n. Then i pmod n = r.
Note that the % operator in C (i%n) and the Fortran mod() function (mod(i,n)) do not satisfy this definition if i is negative - so be careful when writing code in either language.
Since the values taken by U are just a repeat of the values of the grid function of size n: u={U(0),..U(n-1)}, we usually identify the two and use the name u for U.
Example of a periodic function:
Define U(i) = sin(2*pi*i/n), -inf < i < inf,
If i = cn + r, as above, with c an integer, then sin(2*pi*i/n) = sin(2*pi*c + 2*pi*r/n) = sin(2*pi*r/n) since adding a multiple of 2*pi does not change the sine of an angle. This shows that U(i) really only takes n distinct values: U(0), ..., U(n-1) and otherwise simple repeats these forever.
1.C. Extending Grid Functions to Periodic Functions:
Given any grid function u={u(0),...,u(n-1)} on G[n] we can imagine extending the set to an infinite set of numbers U, indexed from -infinity to +infinity, by defining: U(i) = u(i pmod n).
Then the set of numbers U is a periodic grid function of period n. We will frequently do this, and sometimes not even say we are doing it. Furthermore we'll typically use the same name for either function.
2. Complex Vector Spaces
2.1 Scalars:
We will consider vectors whose elements are taken from either the set of all real numbers (R) or the set of all complex numbers (C). The quantities from this space of elements are called scalars, and we denote the space of scalars by S. We use * to denote complex conjugate. If s = a+ib, then s* = a-ib. Of course if the scalars S = R then s* = s since they are real.
2.2: Vectors:
For vector spaces, read Halmos, "Finite Dimensional Vector Spaces", or any of 1000 other sources. Here is a quick review.
An n-component vector v is a set of n scalars: v = (v(0),v(1),...v(n-1)). When discussing vectors, the number of components n will be fixed.
We define addition of vectors, and multiplication of vectors by scalars: u+v is the vector with components u(i)+v(i) su is the vector with components su(i)
We define the operation of inner product (scalar or dot product) of vectors: <u,v> = u(0)*v(0) + u(1)*v(1) + .... + u(n-1)*v(n-1) where the *'s denote complex conjugate.
We define the norm (length) of a vector as: ||u|| = sqrt(<u,u>).
Based on these concrete examples we define the abstract concept of a vector space by:
Definition: A vector space V over a scalar space S is a set of quantities (vectors) on which operations of addition, multiplication by scalars and inner product (scalar product) are defined and where: u,v in V => u+v in V u in V, s in S => su in S u,v in V => <u,v> in S. <u,av+bw> = a<u,v> + b<u,w>
Example 1: The space of real n-dimensional vectors, often called R[n]. The inner product of two vectors is their dot product: <u,v> = u(0)v(0) + u(1)v(1) + .... u(n-1)v(n-1)
Example 2: The space of complex n-dimensional vectors, often called C[n]. The inner product of two vectors is their conjugate dot product: <u,v> = u(0)* v(0) + u(1)* v(1) + .... u(n-1)* v(n-1)
2.3 Linear Independence:
Definition: A set of m vectors v[i], 1<=i<=m, in a vector space is linearly dependent if there exist constants a(i), not all zero, such that a(1)v[1] + a(2)v[2] + .... +a(m)v[m] = 0 If no such relation exists then the v[i] are called linearly independent.
Corollary: if v[i] are linearly dependent then one v is a linear combination of the others.
Proof: Suppose a[m]!=0 above. Then clearly: v[m] = b(1)v[1] + b(2)v[2] + .. b(m-1)v(m-1), where b(i) = a(i)/a(m).
2.4 Basis Sets and Dimensionality:
Definition: A set of m linearly independent vectors E = {e[i]} in a vector space is called a basis set for V if EVERY vector in V can be written as a linear combination of vectors from E:
v = c(0) e[0] + c(1) e[1] + ... c(m-1) e[m-1] or
v = SUM c(k) e[k] (2.1) 0<=k<m
We can also write this out in component form as:
v(j) = SUM c(k) e[k](j) 0<=k<m
Definition: The dimension of a vector space is defined to be the size of the smallest basis set - i.e. the smallest m for which the above is true for all v in V, for some set E of size m. A vector space is called finite dimensional if it has a basis with a finite number of elements.
Example:
One basis for either R[n] or C[n] is the set of n "unit coordinate" vectors E: e[0] = (1,0,0,....,0) e[1] = (0,1,0,....,0) .... e[n-1] = (0,0,0,....,1)
Alternatively we might say that the e[k] are defined by: e[k](j) = 1 if j=k, and = 0 if j!=k.
Clearly any vector v = (v(0),v(1),...,v(n-1)) can be written as:
v = v(0)e[0] + v(1)e[1] + .... + v(n-1)e[n-1]
showing that E is a basis.
One can show quite easily that there is no smaller basis than E. Therefore both R[n] and C[n] are n-dimensional vector spaces.
The basic fact of finite dimensional vector spaces is that:
Theorem: Every basis for a finite dimensional vector space V has the same number of elements n. Conversely, if there is a basis of n elements, then every set of n linearly independent vectors is a basis for the space. This number n is called the dimension of the vector space.
2.5: Orthonormal Sets and Bases:
Any set of vectors E are called orthonormal, if each is of length 1 and is perpendicular to all the others:
<e[k],e[l]> = 1 if k=l, and = 0 if k!=l.
where <a,b> is the inner product on the vector space.
Lemma: Any orthonormal set is linearly independent
Proof: Suppose not. Then one v[i] can be expressed in terms of the others: v[m] = a(1)v[1] + a(2)v[2] + ... + a(m-1)v[m-1] Now take the scalar product with v[m]: 1 = <v[m],v[m]> = a[1]<v[m},v1]> + ... + a(m-1)<v[m],v[m-1]> = 0 , using the orthonormality, which is a contradiction.
Corollary: Any set of n orthonormal vectors in an n-dimensional vector space is a basis.
Proof: By the lemma, the set are linearly independent. Since there are n of them, they must be a basis, by the theorem.
Comment: This gives a nice way to prove that a set of n vectors is a basis. Just show it is an orthonormal set, and it has to be a basis. Proving directly that ANY vector can be expanded in terms of the set of vectors would be much harder.
If the e[k] are orthonormal, and a vector v is a linear combination of the e[k], then we can get a formula for the coefficients c(k) above by taking the dot product of eqn (2.1) by e[l]:
<e[l],v> = SUM c(k) <e[l],e[k]> = c[l]
0<=k<n
where we used the orthonormality of the e[k]. So the converse of equation (2.1) is: c(k) = <e[k],v> (2.2)
In particular these equations apply if the set E is an orthonormal basis. Equations (2.1) and (2.2) are complementary equations and fully describe the basis expansion. These equations allow us to move back and forth between vectors and the coefficients c[k] of basis expansions.
2.6: Representations of a Vector:
For each basis E, a vector v will have a set of n coefficients c(k) and the vector can be recovered from the coefficients. For this reason we call the coefficients c = {c(k)} a representation of the vector v. In fact the vector itself is a representation of itself - in the basis E of unit coordinate vectors, since in that case c(k) = v(k).
Key point: Usually we distinguish between a vector v and its infinitely many representations in terms of basis expansions. The vector v is the important quantity, and we may switch basis to get a convenient view of it.
Example: It may be useful to view objects (vectors) in 1D or 3D from a rotated viewpoint (rotated basis) rather than from the usual x,y,z view. In the simplest case, when rock climbing (1D) we may define the basic direction as up, whereas for deep sea diving we may define the basic direction as down. From a point on the (rotating) surface of the earth, we choose basis vectors that are parallel and perpendicular to the surface at that point - even though these are different from one point to the next. So what's up to an Australian is down to us.
2.7. Complete Example (Revision):
Consider the vector space C(n) of complex n-component vectors where we define the scalar product of u,v as:
<u,v> = SUM u(k)* v(k) 0<=k<n
where c* denotes the complex conjugate of a number c. The norm of a vector is:
||u|| = sqrt(<u,u>) = sqrt { SUM |u(k)|**2 } 0<=k<n
Introduce the basis functions defined by:
e[k](j) = 1, if j=k, e[k](j) = 0, if j!=k
Thus e[0] = {1,0,0,...}, e[1] = {0,1,0,.... }, .. e[n-1] = {0,0,...,1} . These are the familiar unit vectors along coordinate axes and are clearly orthonormal.
Any vector v = {v(0),v(1),....v(n-1)} can then be written as (see eqn. 2.1):
v = SUM c(k) e[k], 0<=k<n
and the inverse equation (2.2) has the form:
c(k) = <e[k],v> = v[k]
So for this basis the coefficients c(k) are very simple - they are exactly the elements of v, as was obvious!
However we could have used any other basis f[] instead of e[]. We would then have had the same relations (2.1) and (2.2) holding, except that the fluky conclusion that c(k) = v(k) would not be true - that one is true only for the basis of unit coordinate vectors.
We will next introduce the Discrete Fourier Transform and learn that it is nothing more than a particular representation of a vector in terms of a special basis (Fourier basis) for C[n]. Both the vectors components (in the unit coordinate basis, ie as described in this example) and its Fourier Transform components (i.e. those in the Fourier basis to be described below) are different ways of looking at the same thing - the vector.
3. The Fourier Basis
On the same vector space C(n) as above, consider the basis functions e[k] defined by:
e[k](j) = 1/sqrt(n) exp( -i 2pi jk/n), 0 <=j < n.
These are oscillating trigonometric functions and we will analyze their shapes and properties later below in section 4, 5. This particular basis is called the Fourier basis for C[n].
Lemma 1: The e[k] are a set of n orthonormal vectors in C(n).
Corollary: The e[k] are a basis for C[n].
Proof of Lemma 1:
<e[k],e[l]> = SUM e[k](j)* e[l](j) 0<=j<n
= 1/n SUM exp( i 2pi (k-l)j/n) 0<=j<n
= 1/n SUM w(k-l)**j, where w(k-l) = exp(i 2pi (k-l)/n) 0<=j<n
Note that w(k-l) is an nth root of unity: w(k-l)**n = exp(i 2pi (k-l)) = 1.
In fact if k=l, then w(k-l) = 1, and in this case we get immediately:
<e[k],e[k]> = 1/n SUM 1 = 1/n n = 1 0<=j<n
For the remaining cases, k!=l, we can compute the SUM above using the algebraic identity:
Lemma 2: SUM a**j = 1 + a + a**2 + .. + a**(n-1) = (1-a**n)/(1-a) 0<=j<n
Proof: Multiply the right side by 1-a and expand out:
(1-a) (1 + a + a**2 + ... + a**(n-1)) =
1 + a + a**2 + ... + a**(n-1) - a - a**2 - ... - a**(n-1) - a**n = 1 - a**n
from which the Lemma follows on division of both sides by (1-a).
Applying the lemma 2 we conclude that for k!=l:
<e[k],e[l]> = 1/n {1 - w(k-l)**n}/{1-w(k-l)} = 0
since as we have seen, w(k-l)**n = 1!
So the e[k] are an orthonormal basis for C[n].
#####
As a result, any vector v in C[n] can be written as an expansion:
v = SUM c(k) e[k] 0<=k<n
and the c(k) are given by the inverse relation:
c(k) = <e[k],v>
The first equation is called the Discrete Fourier Transform (DFT) and the second is called the Inverse Discrete Fourier Transform (IDFT). We can write these out explicitly, i.e. in terms of the individual coordinates v(j) of v, by inserting the definition of e[k] in each equation:
v(j) = 1/sqrt(n) SUM c(k) exp(-i 2pi jk/n) 0<=k<n and c(k) = 1/sqrt(n) SUM v(j) exp(+i 2pi jk/n) 0<=j<n
Notice the remarkable similarity of the two equations. In fact the IDFT differs from the DFT only in the sign of the exponent. In practice the coefficient c(k) is often denoted v^(k) to remind us that it depends on v. Either the n coefficients {v()} or the n coefficients {c()} may be used equally to describe the vector v. They give views of v from different points of view. Sometimes one view is more useful than the other and vice versa.
4. Shapes of the basis functions:
For k near 0 or n the Fourier basis functions e[k] are very slowly varying, while for k near n/2 they vary fastest. For example we can consider the following cases:
k=0: The function e[0] is: e[0](j) = exp(0) = 1, 0 <= j < n.
So e[0] is the constant function, equal to 1 at every grid point.
k=1: The function e[1] is: e[1](j) = exp(-i 2pi j/n) The real part of this function takes the values: Re e[1](0) = 1; Re e[1](n/4) = Re exp(-i pi/2) = Re cos(pi/2) - i sin(pi/2) = 0 Re e[1](n/2) = Re exp(-i pi) = Re cos(pi) - i sin(pi) = -1 Re e[1](3n/2) = Re exp(-i 3pi/2) = Re cos(3pi/2) - i sin(3pi/2) = 0 Re e[1](n) = Re exp(-i 2pi) = 1
So Re u[1] goes from 1 to 0 to -1 to 0 back to 1 as j goes from 0 to n. Between these four values it takes intermediate cosine values. We can describe e[1] by saying that it is the slowest varying cosine on [0,n] which is periodic (value at n equals value at 0). The imaginary part of e[1] is similarly a slowly varying sine function. Because Re e[1] just performs a single cosine oscillation between 0 and n, we call it a "frequency 1" periodic trigonometric function.
If we examined Re e[2] we would find that it varies twice as fast. It runs through the values 1,0,-1,0,1,0,-1,0,1 as j goes from 0 to n - completing two oscillations. We call this a "frequency 2" periodic trigonometric function..
Similarly e[3] will be found to oscillate 3 times between 0 and n.
We find that e[n/2] oscillates n/2 times from 0 to n. In fact:
e[n/2](j) = exp(-i 2pi j/2) = exp(-i pi j) = (-1)**j
So e[n/2] takes the values 1,-1,1,-1,1,...,-1,1 at 0,1,2,3,4,....,n-1,n. This assumes n is even. This is the fastest varying periodic trigonometric function - it jumps from 1 to -1 between each pair of points of the grid. Nothing could jump about faster. This is the frequency n/2 periodic trigonometric function.
As we look at larger k than n/2, we find that e[k] starts to decrease in frequency. For example e[n/2+1] has the same frequency as e[n/2-1] and so on, up to e[n-1] which is similar to e[1]. In fact
e[n-r](j) = exp(-i 2pi (n-r)j/n) = exp (i 2pi rj/n) = {e[r](j)}* where the * denotes complex conjugate. This shows that e[n-r] has frequency r.
Thus we see that our set of functions e[k] are all a) trigonometric b) periodic on grid [0,n] c) oscillate K times on [0,n] where K=k, k<= n/2,and K=n-k, k >=n/2.
In the sequel we will call k near 0 or n-1 low frequency, while k near n/2 is called high frequency. In practice we might consider: High frequency: n/4 <= k <= 3n/4 Low frequency: everything else.
5. Some Properties of the Fourier Basis:
5.1 Translation of Fourier Basis:
Perhaps the most useful feature of the Fourier basis is that each function has simple behavior when translated. Translation means that you look at vectors from another place on the grid. For example if v is a vector with components (in the unit coord basis) of v(k), then we want to know how v would look if we regarded say point d as the origin. Remember that we are on a periodic grid (think circle) and so the choice of origin is really quite arbitrary.
Suppose we move our viewpoint forward by d grid points. Then what we would call v(k) there would really be v(k+d) relative to the old origin. In general we define the translate of a vector v as v[d] ([d] is a superscript) where v[d] is the vector such that: v[d](j) = v(d+j), 0 <= j < n
Lets translate the Fourier basis functions: e[k] becomes e[k][d] where:
e[k][d](j) = e[k](j+d) = 1/sqrt(n) exp( -i 2pi (j+d)k/n) = exp(-i 2pi dk/n) e[k](j)
Since this is true for all components j, we may rewrite it as a vector equation:
e[k][d] = exp(-i 2pi dk/n) e[k] (5.3)
So the translated e[k] is just e[k] multiplied by a complex coefficient.
The above formula (5.3) is extremely useful. It captures completely the shape aspects of the Fourier basis and can be used to simplify lots of equations that would be a nightmare to write out in detail.
5.2 Fourier Transform of a translated vector:
Consider a vector v and its translation v[d]. How are the Fourier transforms of v[d] and v related? We denote the Fourier transform of v by v^(k) and of v[d] by v[d]^(k).
Theorem: v[d]^(k) = exp(-i 2pi kd/n) v^(k) (5.2)
Proof:
v[d]^(k) = SUM e[k](j)* v[d](j) 0<=j<n
= SUM e[k](j)* v(d+j) 0<=j<n
= SUM e[k](j'-d)* v(j'), j' = j+d d<=j'<n+d
= SUM e[k][-d](j')* v(j') (since v and e[k] are periodic) 0<=j'<n
= exp(-i 2pi kd/n) SUM e[k](j')* v(j')
= exp(-i 2pi kd/n) v^(k)
6. Discrete Fourier Analysis Revisited (ignore if you like)
Here we will simply REDERIVE the DFT and IDFT. This time we wont use vector space terminology. As a result the result is more mysterious but the treatment is more elementary. If you understood section 3 then you can skip this section.
We will show that for any given set of n numbers u(0), .., u(n-1), there are a set of quantities u^(k), k = 0, 1, ..., n-1, such that we may write:
u(j) = 1/sqrt(n) SUM exp(- i 2pi jk/n) u^(k), 0<=j<n. (6.1) 0<=k<n
This formula is called the Discrete Fourier Transform. We will prove this formula is true and we will determine the coefficients u^(k).
The u^(k) are given by a remarkably similar formula:
u^(k) = 1/sqrt(n) SUM exp(+ i 2pi jk/n) u(j), 0<=k<n. (6.2) 0<=j<n
which differs only in the sign of the exponential. The latter formula is referred to as the inverse discrete fourier transform.
6.1. Proof of the Discrete Fourier Transform:
Let us DEFINE constants u^(k) by the formula (6.2) above. We will then PROVE formula (6.1) is true. Let us denote the RHS of Equation (6.1) by RHS1. Then
RHS1 = 1/sqrt(n) SUM exp(- i 2pi jk/n) u^(k), 0<=k<n
= 1/n SUM exp(-i 2pi jk/n) SUM exp(i 2 pi kl/n) u(l) 0<=k<n 0<=l<n
= 1/n SUM { SUM exp (i 2pi (l-j)k/n } u(l) } 0<=l<n 0<=k<n
Define w(p) = exp(i 2 pi p/n), for p any integer.
Note that w(p) is an nth root of unity since w(p)**n = exp (i 2pi p) = 1 since p is an integer.
The inner sum i the { } above then has the form:
{ } = SUM w(l-j)**k 0<=k<n
Now for any a, a simple algebraic identity states that:
1 + a + a**2 + ... a**(n-1) = (1 - a**n)/(1-a)
So the above sum { } has the form:
{ } = (1 - w(l-j)**n)/(1-w(l-j))
However w(l-j)**n = 1 since w is a nth root of unity, and therefore { } = 0. The only exception is the case where l-j=0 in which case w(l-j)=1. In that case each term in the { } SUM is just 1. So in all:
{ } = 0 if j!=l, and = n if j=l.
It follows on substituting into the RHS1 formula that:
RHS1 = u(j)
which is what we set out to prove.
7. Description in 2D and 3D
In 2D we deal with the grid G[n1,n2] of n1*n2 points and the set of functions C[n1,n2] on this grid which is a vector space of dimension n1*n2. We choose a basis e[k1,k2] for this space of the form:
e[k1,k2] = e[k1] e[k2],
where e[k] are the Fourier basis functions introduced in 1D. (To be precise we should label them as e[k1;n1], and e[k2;n2] since they depend on n). We call these basis functions high frequency if either k1 or k2 is high frequency, otherwise we call them low frequency. It is easily checked that this basis is orthonormal. As a result we can expand any function in C[n1,n2] as:
u(j1,j2) = SUM SUM c(k1,k2) e[k1](j1) e[k2](j2) 0<=k1<n1 0<=k2<n2
and the c(k1,k2) are then given by:
c(k1,k2) = <e[k1,k2], u>.
These are the 2D DFT and IDFT respectively. We will often denote c(k1,k2) by u^(k1,k2) to remind us of its dependence on u.
Inserting the expression for e[k] we can also write these as:
u(j1,j2) = 1/sqrt(n1) 1/sqrt(n2) SUM SUM 0<=k1<n1 0<=k2<n2 exp (-i 2pi j1*k1/n1) exp (-i 2pi j2*k2/n2) u^(k1,k2) where
u^(k1,k2) = 1/sqrt(n1) 1/sqrt(n2) SUM SUM 0<=k1<n1 0<=k2<n2 exp (+i 2pi j1*k1/n1) exp (+i 2pi j2*k2/n2) u(j1,j2)
Notice that these explicit representations are far more ugly than the versions using vector notation.
The extension to 3D is analogous.
For 2D there is an obvious analogue of the lemma we proved in section 5.2 for a translated vector. We can now translate in two dimensions. Let u[d1,d2] denote u translated by d1 in x and d2 in y: u[d1,d2](j1,j2) = u(d1+j1,d2+j2)
Lemma: u[d1,d2]^(k1,k2) = exp(-i 2pi k1d1/n1 ) exp(-i 2pi k2d2/n2 ) u^(k1,k2)
Proof of Lemma:
Follows immediately from the 1D version since the 2D transform u^ is just the product of two 1D transforms, one in k1, the other in k2. Each of these picks up a factor of exp(-i 2pi ...), so the whole expression picks up the product.
or, brute force proof (ignore unless you love indices):
u[d1,d2]^(k1,k2) = 1/n SUM e[k1](j1) e[k2](j2) u[d1,d2](j1,j2) 0<=j1<n1 0<=j2<n2
= 1/n SUM exp(i 2pi j1k1/n1) exp(i 2pi j2k2/n2) u(j1+d1,j2+d2) 0<=j1<n1 0<=j2<n2
= exp(-i 2pi d1k1/n1) exp(-i 2pi d2k2/n2) x 1/n SUM exp(i 2pi (j1+d1)k1/n1) exp(i 2pi (j2+d2)k2/n2) u(j1+d1,j2+d2) 0<=j1<n1 0<=j2<n2
= exp(-i 2pi d1k1/n1) exp(-i 2pi d2k2/n2) x 1/n SUM exp(i 2pi i1k1/n1) exp(i 2pi i2k2/n2) u(i1,i2) -d1<=i1<n1-d1 -d2<=i2<n2-d2
= exp(-i 2pi d1k1/n1) exp(-i 2pi d2k2/n2) x u^(k1,k2)
where in the last step we use the periodicity of u(i1,i2) to change the sum limits to 0<=i1<n1 and 0<=i2<n2. This completes the proof of the Lemma.