Consider the following integral:
(2.3) |
(2.4) |
(2.5) |
The discrete Fourier transform matrix shown in equation (2.6) is a complex Vandermonde matrix. Such matrices appear in polynomial interpolation. It is easy to write a simple Fortran program that builds a discrete Fourier transform matrix. But instead of doing it with Fortran, this time we'll use our interactive Fortran, Octave, to build a discrete Fourier transform matrix.
Octave, version 2.0.13 (sparc-sun-solaris2.6). Copyright (C) 1996, 1997, 1998 John W. Eaton. This is free software with ABSOLUTELY NO WARRANTY. For details, type `warranty'. octave:1> w = exp(-2 * pi * i / 4) w = -i octave:2> F = ones(4,4) F = 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 octave:3> (0:3) ans = 0 1 2 3 octave:4> F(2,:) = w .^ (0:3) F = 1 1 1 1 1 -i -1 i 1 1 1 1 1 1 1 1 octave:5> F(3,:) = w .^ (2 * (0:3)) F = 1 1 1 1 1 -i -1 i 1 -1 1 -1 1 1 1 1 octave:6> F(4,:) = w .^ (3 * (0:3)) F = 1 1 1 1 1 -i -1 i 1 -1 1 -1 1 i -1 -i octave:7>
Another way to accomplish the same would be as follows:
octave:8> F = ones(4,4) F = 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 octave:9> F(2,:) = w .^ (0:3) F = 1 1 1 1 1 -i -1 i 1 1 1 1 1 1 1 1 octave:10> F(3,:) = F(2,:) .* F(2,:) F = 1 1 1 1 1 -i -1 i 1 -1 1 -1 1 1 1 1 octave:11> F(4,:) = F(2,:) .* F(3,:) F = 1 1 1 1 1 -i -1 i 1 -1 1 -1 1 i -1 -i octave:13>This is a faster way because multiplication is cheaper than exponentiation.
Once we have built the discrete Fourier transform matrix we can attempt to find a numerical Fourier transform of, say, a sine function:
octave:39> x = [ 0, 2*pi/4, 2*2*pi/4, 3*2*pi/4 ] x = 0.00000 1.57080 3.14159 4.71239 octave:40> x = sin(x) x = 0.00000 1.00000 0.00000 -1.00000 octave:41> 2 * pi / 4 * F * x' ans = 0.00000 + 0.00000i 0.00000 - 3.14159i 0.00000 + 0.00000i -0.00000 + 3.14159i octave:42>Our numerical result, neglecting the last harmonic, is that
alg' integ(sin(x)^2, x, 0, 2 * pi) pibut remember to do it in the radians mode.
What about that last harmonic? Repeating the same computation for n = 8 returns:
octave:43> w = exp(-2 * pi * i / 8) w = 0.70711 - 0.70711i octave:44> F = ones(8,8); octave:45> for k = 2:8 > F(k,:) = w .^ ((k - 1) * (0:7)); > end octave:46> x = 2 * pi / 8 * (0:7); octave:47> x = sin(x); octave:48> 2 * pi / 8 * F * x' ans = -0.00000 + 0.00000i 0.00000 - 3.14159i 0.00000 + 0.00000i 0.00000 + 0.00000i 0.00000 + 0.00000i 0.00000 + 0.00000i 0.00000 + 0.00000i -0.00000 + 3.14159i octave:49>So we are still getting the same result, , and the only difference is that the last harmonic has been moved further away.
Let us try a signal with higher harmonics:
octave:45> x = 2 * pi / 8 * (0:7); octave:46> y = sin(x) + 0.5 * sin (2 * x); octave:47> 2 * pi / 8 * F * y' ans = 0.00000 + 0.00000i 0.00000 - 3.14159i 0.00000 - 1.57080i -0.00000 + 0.00000i 0.00000 + 0.00000i 0.00000 + 0.00000i -0.00000 + 1.57080i -0.00000 + 3.14159i octave:48> y = cos(x) + sin(2 * x) - cos(3 * x); octave:49> 2 * pi / 8 * F * y' ans = 0.00000 + 0.00000i 3.14159 + 0.00000i 0.00000 - 3.14159i -3.14159 - 0.00000i -0.00000 + 0.00000i -3.14159 - 0.00000i -0.00000 + 3.14159i 3.14159 + 0.00000i octave:50>What is happening is that harmonics are inserted from both ends. The reason for this is that the discrete Fourier transform returns amplitudes that correspond to both positive and negative wave numbers. Amplitudes for positive wave numbers are returned in slots labeled with . Amplitudes for negative wave numbers are returned in reverse order in slots and the Nyquist freqency amplitude
(2.7) |
(2.8) |
In our case , k_{c} = 4, and the Nyquist frequency Amplitude is -0.0 + 0.0 i.
How did those negative frequencies creep into our sum:
F_{n - k} = F_{-k} | (2.9) |
What is going to happen if your signal contains harmonics that are higher than the Nyquist frequency k_{c}? It turns out that all of the power spectral density that lies outside of the range [-k_{c}, +k_{c}] will be sneakily moved into that range, i.e., any amplitude that corresponds to a frequency outside that range is going to be aliased into that range.
There is no way to overcome aliasing other than to know the bandwidth limit of a signal in advance and then sample sufficiently densely, so that the Nyquist frequency is higher than the bandwidth limit.