Discrete Fourier Transform

The integral form of the Fourier transform equations is only useful if you have a mathematical function that characterizes the signal to be analyzed . A more common situation is to have a series of measurements of a given signal taken at a discrete time or frequency increments . The DFT is a way to apply a Fourier transform to a discrete set of data.

Let's say the amplitude of a signal is measured N times at equal time increments over a total sampling time of T . You would then have N measurements of x ( t ), [ x , x ₁ , x ₂ , ..., x _{N “ 1} ]. Any one discrete x value would be given by the following equation.

Equation 22.8

Eq. (22.8) assumes that the sampling started at time t = 0. We can discretize the frequency domain in a similar manner. The overall sampling frequency is F = N / T . We will take N discrete values of the frequency amplitude, [ X , X ₁ , X ₂ , ..., X _{N “ 1} ]; where any individual frequency amplitude is given by the expression in Eq. (22.9).

Equation 22.9

Eqs. (22.6), (22.8), and (22.9) can be used to derive a forward transform for a discrete value of X _k along the sampling range, Eq. (22.10).

Equation 22.10

Using the same process, an equation can also be derived for the inverse transform.

Equation 22.11

The D f term to the left of the summation in Eq. (22.11) has been replaced with D f = 1/( N D T ). The standard form of the DFT equations can be obtained by using the following relations.

Equation 22.12

Using the variables from Eq. (22.12), the DFT equations take their commonly expressed form shown in Eq. (22.13) and Eq. (22.14).

Equation 22.13

Equation 22.14

Because Eq. (22.13) is a summation, the magnitude of the X _k terms will increase with increasing sample number. If you are taking millions of samples, this can make the frequency response output difficult to plot. To avoid this sit uation, the 1/ N term is sometimes incorporated into the forward transform equation by using the expression X ( f _k ) = N D tX _k .

As it turns out, the Java implementation of the DFT equations is quite simple. We will write a method named discreteFT() that will perform a forward or inverse DFT analysis on an array of data. The method will be public and static and will be declared in a class named Fourier . The Fourier class will be placed in the TechJava.MathLib package.

Before we can implement the discreteFT() method, we need to convert Eq. (22.13) and Eq. (22.14) back into trigonometric form. We will do this conversion using the following relations.

Equation 22.15

The other important thing to consider is that both the initial and transformed data can have real and imaginary components .

Equation 22.16

When Eq. (22.15) and Eq. (22.16) are incorporated into Eq. (22.13) and Eq. (22.14), we come up with the following expressions for a forward and inverse DFT. We will use the 1/ N term in the forward transform.

Equation 22.17

Equation 22.18

If you look at Eq. (22.17) and Eq. (22.18), you will see that the equations are very similar. The only differences are the 1/ N term and the leading edge coefficient on the sin() terms. This similarity means that we can implement the forward and inverse Fourier transforms with a single method. You should also note that, generally speaking, the transformed values will have both real and imaginary components whether the initial data is real, imaginary, or both.

The discreteFT() method is really very simple. It takes three arguments. The first is the array of data to be transformed. This array will contain real-imaginary data pairs for each sample point. The second argument is the number of sample points. The final argument is a boolean that is true if this is a forward transform calculation (from time domain to frequency domain) or false if it is a reverse transform (from frequency to time).

The method initializes a temporary variable omega that represents w = 2 p / N . If this calculation is for an inverse transform, a negative sign is placed in front of the variable. This will cause the Math.sin() method evaluations later on to flip their signs. The discreteFT() method then goes through the summations from Eq. (22.17) and Eq. (22.18). If this is a forward transform, the resulting X _k values are scaled by 1/ N .

The discreteFT() method source code is shown here.

 package TechJava.MathLib; public class Fourier {   //  Discrete Fourier transform   public static double[] discreteFT(double[] data,                            int N, boolean forward) {     double X[] = new double[2*N];     double omega;     int k, ki, kr, n;     //  If this is a inverse transform, reverse the     //  sign of the angle so the sin() terms will     //  change sign.     if (forward) {       omega = 2.0*Math.PI/N;     } else {       omega = -2.0*Math.PI/N;     }     //  Perform the discrete Fourier transform.     //  The real and imaginary data are stored in the     //  x[] and X[] vectors as pairs, one real and     //  one imaginary value for each sample point.     for(k=0; k<N; ++k) {       kr = 2*k;       ki = 2*k + 1;       X[kr] = 0.0;       X[ki] = 0.0;       for(n=0; n<N; ++n) {         X[kr] += data[2*n]*Math.cos(omega*n*k) +                  data[2*n+1]*Math.sin(omega*n*k);         X[ki] += -data[2*n]*Math.sin(omega*n*k) +                   data[2*n+1]*Math.cos(omega*n*k);       }     }     //  If this is a forward transform, multiply     //  in the 1/N terms     if ( forward ) {       for(k=0; k<N; ++k) {         X[2*k] /= N;         X[2*k + 1] /= N;       }     }     //  Return the transformed data.     return X;   } } 

Let's apply the discreteFT() method to compute an FFT on a 2 Hz cosine wave. The amplitude as a function of time for this wave is given by the equation below.

Equation 22.19

We will take 64 data samples over a 2-second sample period. The discreteFT() method will first compute a FFT to obtain the frequency spectrum for a 2 Hz cosine wave. The discreteFT() method will then be used to perform an inverse transform that will reconstruct a 2 Hz cosine wave from its frequency spectrum. The program that will do all this is named TestDFT.java . Its code listing follows . When the transformed data is written to the console, the sample indexes are converted to units of time or frequency.

 import TechJava.MathLib.*; public class TestDFT {   public static void main(String args[]) {     int N = 64;     double T = 2.0;     double tn, fk;     double data[] = new double[2*N];     //  A 2 Hz cosine function that is sampled     //  for two seconds. First load the amplitude data     for(int i=0; i<N; ++i) {       data[2*i] = Math.cos(4.0*Math.PI*i*T/N);       data[2*i+1] = 0.0;     }     //  Compute the DFT     double X[] = Fourier.discreteFT(data, N, true);     //  Print out the frequency spectrum     for(int k=0; k<N; ++k) {       fk = k/T;       System.out.println("f["+k+"] = "+fk+"Xr["+k+             "] = "+X[2*k]+ "  Xi["+k+"] = "+X[2*k + 1]);     }     //  Reconstruct a 2 Hz cosine wave from its     //  frequency spectrum. First load the frequency     //  data.     for(int i=0; i<N; ++i) {       data[2*i] = 0.0;       data[2*i+1] = 0.0;       if (i == 4  i == N-4 ) {         data[2*i] = 0.5;       }     }     //  Compute the DFT     double x[] = Fourier.discreteFT(data, N, false);     //  Print out the amplitude vs time data.     System.out.println();     for(int n=0; n<N; ++n) {       tn = n*T/N;       System.out.println("t["+n+"] = "+tn+"xr["+n+               "] = "+x[2*n]+"  xi["+n+"] ="+x[2*n + 1]);     }     } } 

The results of the FFT are shown in Figure 22.1. The frequency amplitude data was converted to decibels using the Math2.log10() method we developed in Chapter 17. There are two peaks in the data. The first is at 2 Hz, the frequency of the cosine signal. The second is at 30 Hz. This may seem a bit confusing but is the result of starting our summations at n = 0. We took samples at a rate of 32 samples a second for this analysis. The DFT actually computes the spectrum from “16 to 16 Hz, but we have shifted it over to start at 0 so the negative part of the spectrum has been mapped onto the frequency range 16 to 32 Hz. The peak at 30 Hz on the figure really corresponds to a peak at “2 Hz. Having peaks at ±2 Hz is the correct solution. The imaginary part of the frequency part of the spectrum as computed by the discreteFT() method was uniformly zero as it should be since a cosine wave is a real, even function.

Figure 22.1. Frequency response of 2 Hz cosine wave

The results of the inverse Fourier transform are shown in Figure 22.2. In this case we reconstructed a 2 Hz cosine wave from a frequency spectrum using the discreteFT() method. The DFT method did an excellent job of reconstructing the cosine wave. The computed points lie right on the exact solution curve.

Figure 22.2. Cosine wave reconstruction