Discrete Fourier Transforms
Fourier transforms have many important applications in all branches
of pure and applied mathematics. Given an arbitrary function f(t),
the basic idea is to decompose that function (over some finite
interval) optimally into a sum of pure sine waves, and to find the
frequencies, amplitudes, and phases of those waves.
Recall that the quantity A*exp(theta i) can be regarded as a
vector on the complex plane, with length A and directed at an angle
theta up from the positive real axis. The real and imaginary parts
of this quantity are A*cos(theta) and A*sin(theta) respectively.
Multiplying any complex number by A*exp(theta i) has the effect of
scaling the magnitude of that number by A and "rotating" it in the
complex plane through the angle theta.
Suppose we're given the sequence of values taken on by the function
f(t) for t=0,1,2,..,n-1. (We will refer to the independent variable
"t" as the time, although it not literally be time.) It is most
convenient to fit these points to a sum of exponentials
sqrt(n) f(t) = c_0 + c_1 exp(-1qti) +...+ c_(n-1) exp(-(n-1)qti)
where q = (2PI/n). (We've applied the factor sqrt(n) because it
makes the inverse transform work out nicely, as shown below.) We
have n data points and n unknown coefficients, so we can solve
for the coefficients to give an exact fit of our data points. Of
course, this will give us a complex function of t, but the values
will be purely real at t=0,1,2,..,n-1, so the imaginary parts
vanish at these points, which implies that if we simply omit
the imaginary parts altogether we are left with a real function
f(t) that passes through each of the original data points and
consists of a sum of n cosine terms, each with its own amplitude,
frequency, and phase.
The coefficients c_j will generally be complex, which means they
can be written in the form c_j = C_j exp(phi_j i) for some angle
phi_j. Thus the entire jth term in the fitted expression for f(t)
is C_j exp((jqt + phi_j)i), so phi_j represents the phase angle of
this term. Omitting the imaginary part, the term becomes simply
C_j cos(jqt + phi_j).
The sequence of coefficients c_j represents the Fourier transform
of the time sequence f(t). The jth term corresponds to the
frequency w = jq = j(2PI/n) rad/sec (assuming the units of t are
seconds). Thus our frequencies come in integer multiples of 2PI/n,
so the more often we sample, the better resolution of frequencies
we will achieve.
The straight-forward way of computing the coefficients c_j is by
simply solving the system of simultaneous equations for the n data
points. To illustrate, if we had just four data points, f(0), f(1),
f(2), and f(3), we would solve this system of equations
_ _ _ _ _ _
| | | | | |
| 1 1 1 1 | | c_0 | | f(0) |
| | | | | |
1 | 1 exp(-qi) exp(-2qi) exp(-3qi) | | c_1 | | f(1) |
------- | | | | = | |
sqrt(n) | 1 exp(-2qi) exp(-4qi) exp(-6qi) | | c_2 | | f(2) |
| | | | | |
| 1 exp(-3qi) exp(-6qi) exp(-9qi) | | c_3 | | f(3) |
|_ _| |_ _| |_ _|
The inverse of the square matrix is given by simply negating each of
the exponential arguments, and then multiplying the entire array by
1/n, as can be verified from the fact that the element in the jth row
of the kth column of that matrix times its inverse is the sum
exp(0(j-k)qi) + exp(1(j-k)qi) + exp(2(j-k)qi) + exp(3(j-k)qi)
which equals n when j=k, and otherwise equals 0. Therefore, the
inverse of the above matrix equation is
_ _ _ _ _ _
| | | | | |
| 1 1 1 1 | | f(0) | | c_0 |
| | | | | |
1 | 1 exp(qi) exp(2qi) exp(3qi) | | f(1) | | c_1 |
------- | | | | = | |
sqrt(n) | 1 exp(2qi) exp(4qi) exp(6qi) | | f(2) | | c_2 |
| | | | | |
| 1 exp(3qi) exp(6qi) exp(9qi) | | f(3) | | c_3 |
|_ _| |_ _| |_ _|
The sign in the exponentials just determines the direction of rotation
of these periodic terms, so we see that the sequences f(t) and c_j
are essentially duals of each other, meaning that each is the Fourier
transform of the other. (It should be clear now why we divided by
sqrt(n) earlier, to make these transforms formally identical, up to
the direction of rotation.) Let's call the square matrix multiplied
by 1/sqrt(n) in this equation E. Then the elements of the E array
can be specified simply as
exp(jk q i)
E(j,k) = ----------- where q = 2PI/n
sqrt(n)
So, given the sequence of values f(t) as a column vector f we can
compute the Fourier transform c(w), which is column vector of n
elements, by the equation
c = Ef
Obviously the inverse transform is simply f = E^-1 c. The values of
c represent the (complex) amplitudes of a set of n cosine functions
with the frequencies
|w| = 0, 2PI/n, 2(2PI/n),..., (n-1)(2PI/n)
radians per [unit of t]. The power spectral density of the signal
f(t) at a particular frequency is simply the squared amplitude
at that frequency. In other words, the power spectrum can be
represented as a column vector with the values
P(k) = Re(c(k))^2 + Im(c(k))^2
Notice that this covers the power of the signal over (1/n)th of the
total frequency range. If we want to represent a smaller or larger
band of frequencies, we would scale the amplitude accordingly.
The phase spectrum is similar to the power spectrum, except that it
represents the phase angle corresponding to a given frequency. Thus
is can be expressed as the column vector
/ Im(c(k) \
phi(k) = arg(c(k)) = atan( -------- ) [+- PI as needed]
\ Re(c(k) /
As an example, suppose we are given the sequence of values
f(t) = 3,5,7,2 for t = 0,1,2,3 sec
Multiplying the column vector of these values by the 4th-order E
matrix gives the Fourier transform
c(w) = 17/2, -2+(3/2)i, 3/2, -2-(3/2)i
for w = 0, PI/2, PI, and 3PI/2 radians/sec. The power spectrum is
P = 72.25, 6.25, 2.25, 6.25
at the frequencies 0, PI/2, PI, and 3PI/2 radians/sec, so the
predominant component is the DC signal. The amplitudes and phases
of the four components are
A = 8.5, 2.5, 1.5, 2.5 phi = 0, 2.498, 0, -2.498 radians
Therefore, we have fit the function f(t) to the continuous function
f(t) = [8.5 + 2.5cos(-qt+2.498) + 1.5cos(-2qt) + 2.5cos(-3qt-2.498)]/2
where q = PI/2. A plot of this function is shown below:
Incidentally, the E matrices for n=2,3,4 are as shown below
_ _
1 | 1 1 |
E_2 = ------- | |
sqrt(2) |_ 1 -1 _|
_ _
1 | 1 1 1 | -1 + sqrt(-3)
E_3 = ------- | | r = ------------
sqrt(3) | 1 r r* | 2
| |
|_ 1 r* r _|
_ _
| 1 1 1 1 |
| |
1 | 1 i -1 -i |
E_4 = --- | |
2 | 1 -1 1 -1 |
| |
|_ 1 -i -1 i _|
This shows that the familiar reciprocal transform
x + y x - y
X = ------- Y = -------
sqrt(2) sqrt(2)
is a primitive example of a Fourier transform. These simple
discrete transforms occur quite often in quantum mechanics and
the theory of relativity.
Another interesting aspect of the continuous representations
given by the Fourier series for a linear sequence of discrete
data points is illustrated by the case f(t) = t for t=0,1,..,19.
This is just a simple ramp function, but the 20th order discrete
transform interprets it as an oscillating function as shown below:
This "smudging" of a linear ramp of data points is quite general,
and applies to any number of data points, provided we perform an
nth order Fourier transform on n points. However, it's worth noting
that we don't need to use an nth order Fourier transform with n
data points. Given a time sequence with many thousands of data
points (say) we can apply multiple linear regression to fit a
Fourier transform of any order less than or equal to the number
of samples. (See the note Multiple Linear Regression for details.)
Continuous Fourier transforms are essentially the same as the
discrete transforms described above, except that we are given
a continuous function f(t) instead of a finite set of discrete
samples, and we replace the summations with integrations.
Return to MathPages Main Menu
Сайт управляется системой
uCoz