Spectrum Analysis Using Java
Java Programming, Notes # 1485
Preface
Spectral analysis
A previous lesson entitled Fun with Java, How and Why Spectral Analysis Works explained some of the fundamentals regarding spectral analysis.
The lesson entitled Spectrum Analysis using Java, Sampling Frequency, Folding Frequency, and the FFT Algorithm presented and explained several Java programs for doing spectral analysis, including both DFT programs and FFT programs. That lesson illustrated the fundamental aspects of spectral analysis that center around the sampling frequency and the Nyquist folding frequency.
The lesson entitled Spectrum Analysis using Java, Frequency Resolution versus Data Length used similar Java programs to explain spectral frequency resolution.
The lesson entitled Spectrum Analysis using Java, Complex Spectrum and Phase Angle explained issues involving the complex spectrum, the phase angle, and shifts in the time domain.
This lesson will illustrate and explain forward and inverse Fourier transforms using both DFT and FFT algorithms. I will also illustrate and explain the implementation of frequency filtering by modifying the complex spectrum in the frequency domain and then transforming the modified complex spectra back into the time domain.
Viewing tip
You may find it useful to open another copy of this lesson in a separate browser window. That will make it easier for you to scroll back and forth among the different figures and listings while you are reading about them.
Supplementary material
I recommend that you also study the other lessons in my extensive collection of online Java tutorials. You will find those lessons published at Gamelan.com. However, as of the date of this writing, Gamelan doesn't maintain a consolidated index of my Java tutorial lessons, and sometimes they are difficult to locate there. You will find a consolidated index at www.DickBaldwin.com.
Preview
In this lesson, I will present and explain the following new programs:
- Dsp035 - Illustrates the reversible nature of the Fourier transform. This program transforms a real time series into a complex spectrum, and then reproduces the real time series by performing an inverse Fourier transform on the complex spectrum. This is accomplished using a DFT algorithm.
- InverseComplexToReal01 - Class that implements an inverse DFT algorithm for transforming a complex spectrum into a real time series.
- Dsp036 - Replicates the behavior of the program named Dsp035 but uses an FFT algorithm instead of a DFT algorithm.
- InverseComplexToRealFFT01 - Class that implements an inverse FFT algorithm for transforming a complex spectrum into a real time series.
- Dsp037 - Illustrates filtering in the frequency domain. Uses an FFT algorithm to transform a time-domain impulse into the frequency domain. Modifies the complex spectrum, eliminating energy within a specific band of frequencies. Uses an inverse FFT algorithm to produce the filtered version of the impulse in the time domain.
In addition, I will use the following programs that I explained in the lesson entitled Spectrum Analysis using Java, Sampling Frequency, Folding Frequency, and the FFT Algorithm.
- ForwardRealToComplex01 - Class that implements a forward DFT algorithm for transforming a real time series into a complex spectrum.
- ForwardRealToComplexFFT01 - Class that implements a forward FFT algorithm for transforming a real time series into a complex spectrum.
- Graph03 - Used to display various types of data. (The concepts were explained in an earlier lesson.)
- Graph06 - Also used to display various types of data in a somewhat different format. (The concepts were also explained in an earlier lesson.)
Discussion and Sample Code
Description of the program named Dsp035
The program named Dsp035 illustrates forward and inverse Fourier transforms using DFT algorithms.
The program performs spectral analysis on a time series consisting of pulses and a sinusoid. Then it passes the resulting real and complex parts of the spectrum to an inverse Fourier transform program. This program performs an inverse Fourier transform on the complex spectral data to reconstruct the original time series.
This program can be run with either Graph03 or Graph06 in order to plot the results. Enter the following at the command-line prompt to run the program with Graph03:
java Graph03 Dsp035
The program was tested using J2SE 1.4.2 under WinXP.
The order of the plotted results
When the data is plotted (see Figure 1) using the programs Graph03 or Graph06, the plots appear in the following order from top to bottom:
- The input time series
- The real spectrum of the input time series
- The imaginary spectrum of the input time series
- The amplitude spectrum of the input time series
- The output time series produced by the inverse Fourier transform
The format of the plots
There were 256 values plotted horizontally in each section. I plotted the values on a grid that is 270 units wide to make it easier to view the plots on the rightmost end. This leaves some blank space on the rightmost end to contain the numbers, preventing the numbers from being mixed in with the plotted values. The last actual data value coincides with the rightmost tick mark on each plot.
The forward Fourier transform
A static method named transform belonging to the class named ForwardRealToComplex01 was used to perform the forward Fourier transform.
(I explained this class and the transform method in the earlier lesson entitled Spectrum Analysis using Java, Sampling Frequency, Folding Frequency, and the FFT Algorithm.)
The method named transform does not implement an FFT algorithm. Rather, it implements a DFT algorithm, which is more general than, but much slower than an FFT algorithm.
(See the program named Dsp036 later in the lesson for the use of an FFT algorithm.)
The inverse Fourier transform
A static method named inverseTransform belonging to the class named InvereComplexToReal01 was used to perform the inverse Fourier transform. I will explain this method later in this lesson.
Let's see some results
Before getting into the technical details of the program, let's take a look at the results shown in Figure 1.
The top plot in Figure 1 shows the input time series used in this experiment.
Figure 1 Forward and inverse transform of a time series using DFT algorithm |
Length is a power of two
The time series is 256 samples long. Although the DFT algorithm can accommodate time series of arbitrary lengths, I set the length of this time series to a power of two so that I can compare the results with results produced by an FFT algorithm later in the lesson.
(Recall that most FFT algorithms are restricted to input data lengths that are a power of two.)
The input time series
As you can see, the input time series consists of three concatenated pulses separated by blank spaces. The pulse on the leftmost end consists simply of some values that I entered into the time series to create a pulse with an interesting shape.
The middle pulse is a truncated sinusoid.
The rightmost pulse is a truncated square wave.
The objective
The objective of the experiment is to confirm that it is possible to transform this time series into the frequency domain using a forward Fourier transform, and then to recreate the time series by using an inverse Fourier transform to transform the complex spectrum back into the time domain.
The real part of the spectrum is symmetrical
The real part of the complex spectrum is shown in the second plot from the top in Figure 1. It will become important later to note that the real part of the spectrum is symmetrical about the folding frequency near the center of the plot (at the eighth tick mark).
Without attempting to explain why, I will simply tell you that the real part of the Fourier transform of a complex series whose imaginary part is all zeros is always symmetrical about the folding frequency.
The imaginary part of the spectrum is asymmetrical
The imaginary part of the complex spectrum is shown in the third plot from the top. Again, it will become important later to note that the imaginary part of the spectrum is asymmetrical about the folding frequency.
Once again, without attempting to explain why, the imaginary part of the Fourier transform of a complex series whose imaginary part is all zeros is always asymmetrical about the folding frequency.
The converse is also true
It is also true that the values of the imaginary part of the Fourier transform of a complex spectrum whose real part is symmetrical about the folding frequency and whose imaginary part is asymmetrical about the folding frequency will be all be zero. I will take advantage of these facts later to simplify the computing and plotting process.
The amplitude spectrum
The amplitude spectrum is shown in the fourth plot down from the top. Recall from previous lessons that the amplitude values are always positive, consisting of the square root of the sum of the squares of the real and imaginary parts.
The output time series
The output time series, produced by performing an inverse Fourier transform on the complex spectrum is shown in the bottom plot in Figure 1. Compare the bottom plot to the top plot. As you can see, they are the same, demonstrating the reversible nature of the Fourier transform.
Let's see some code
I will discuss this program in fragments. A complete listing of the program is provided in Listing 14 near the end of the lesson.
The beginning of the class for Dsp035, including the declaration of some variables and the creation of some array objects is shown in Listing 1. This code is straightforward.
class Dsp035 implements GraphIntfc01{ final double pi = Math.PI; int len = 256; double[] timeDataIn = new double[len]; double[] realSpect = new double[len]; double[] imagSpect = new double[len]; double[] angle = new double[len];//unused double[] magnitude = new double[len]; double[] timeDataOut = new double[len]; int zero = 0; Listing 1 |
The constructor
The constructor begins in Listing 2. The code in Listing 2 creates the input time series data shown in the top plot of Figure 1.
public Dsp035(){//constructor //Create the raw data pulses timeDataIn[0] = 0; timeDataIn[1] = 50; //... //code deleted for brevity //... timeDataIn[254] = -80; timeDataIn[255] = -80; //Create raw data sinusoid for(int x = len/3;x < 3*len/4;x++){ timeDataIn[x] = 80.0 * Math.sin( 2*pi*(x)*1.0/20.0); }//end for loop Listing 2 |
Note that I deleted much of the code from Listing 2 for brevity. You can view the missing code in Listing 14 near the end of the lesson.
Compute the complex spectrum
The code in Listing 3 invokes the static transform method of the ForwardRealToComplex01 class to compute and save the complex spectrum of the time series shown in the top plot of Figure 1.
ForwardRealToComplex01.transform(timeDataIn, realSpect, imagSpect, angle, magnitude, zero, 0.0, 1.0); Listing 3 |
The method parameters
I explained the transform method in the earlier lesson entitled Spectrum Analysis using Java, Sampling Frequency, Folding Frequency, and the FFT Algorithm. The three middle plots in Figure 1 are plots of the data returned in the arrays referred to by realSpect, imagSpect, and magnitude by the transform method.
The angle results returned by the transform program are not used in this lesson.
One of the parameters (zero) establishes that the first sample in the time series array referred to by timeDataIn represents the zero time origin.
The parameters also specify that the complex spectrum is to be computed at a set of equally spaced frequencies ranging from zero (0.0) to the sampling frequency (1.0).
Perform the inverse Fourier transform
The code in Listing 4 invokes the static inverseTransform method of the InverseComplexToReal01 class to perform an inverse Fourier transform on the complex spectral data, producing the output time series shown in the bottom plot in Figure 1.
InverseComplexToReal01.inverseTransform( realSpect, imagSpect, timeDataOut); }//end constructor Listing 4 |
I will explain the inverseTransform method later.
An object of the class Dsp035
Listing 4 also signals the end of the constructor. Once the constructor has completed executing, an object of the Dsp035 class exists. The array objects belonging to the object have been populated with the original time series, the complex spectrum of the original time series, and the output time series produced by performing an inverse Fourier transform on that complex spectrum. This data is ready for plotting.
The interface methods
All of the remaining code in Dsp035 consists of the six methods necessary to satisfy the interface named GraphIntfc01. Those methods are required to provide data to the plotting program, as explained in earlier lessons in this series.
If you have studied the previous lessons in this series, you probably don't want to hear any more about those methods, so I won't discuss them further. You can view the six interface methods in Listing 14 near the end of the lessons.
The inverseTransform method of the InverseComplexToReal01 class
The static method named inverseTransform performs a complex-to-real inverse discrete Fourier transform returning a real result only. In other words, the method transforms a complex input to a real output.
There are more efficient ways to write this method taking known symmetry and asymmetry conditions into account. However, I wrote the method the way that I did because I wanted it to mimic the behavior of an FFT algorithm. Therefore, the complex input must extend from zero to the sampling frequency.
The method does not implement an FFT algorithm. Rather, the inverseTransform method implements a straight forward sampled-data version of the continuous inverse Fourier transform that is defined using integral calculus.
The parameters
The parameters to the inverseTransform are:
- double[] realIn - incoming real data
- double[] imagIn - incoming imag data
- double[] realOut - outgoing real data
The method considers the data length to be realIn.length, and considers the computational time increment to be 1.0/realIn.length.
Assumptions
The method returns a number of points equal to the data length. It assumes that the real input consists of positive frequency points for a symmetric real frequency function. That is, the real input is assumed to be symmetric about the folding frequency. The method does not test this assumption.
The method assumes that imaginary input consists of positive frequency points for an asymmetric imaginary frequency function. That is, the imaginary input is assumed to be asymmetric about the folding frequency. Once again, the method does not test this assumption.
A real output
A symmetric real part and an
asymmetric imaginary part guarantee that the
imaginary output will be all zero values. Having made that assumption, the program makes no attempt
to compute an imaginary output. If the assumptions described above are not valid, the
results won't be valid.
The program was tested using J2SE v1.4.2 under WinXP.
The inverseTransform method
The beginning of the class and the beginning of the static inverseTransform method is shown in Listing 5.
public class InverseComplexToReal01{ public static void inverseTransform( double[] realIn, double[] imagIn, double[] realOut){ int dataLen = realIn.length; double delT = 1.0/realIn.length; double startTime = 0.0; Listing 5 |
Listing 5 declares and initializes some variables that will be used later.
The inverse transform computation
Listing 6 contains a pair of nested for loops that perform the actual inverse transform computation.
//Outer loop interates on time domain // values. for(int i=0; i < dataLen;i++){ double time = startTime + i*delT; double real = 0; //Inner loop iterates on frequency // domain values. for(int j=0; j < dataLen; j++){ real += realIn[j]* Math.cos(2*Math.PI*time*j) + imagIn[j]* Math.sin(2*Math.PI*time*j); }//end inner loop realOut[i] = real; }//end outer loop }//end inverseTransform }//end class InverseComplexToReal01 Listing 6 |
If you have been studying the earlier lessons in this series, you should be able to understand the code in Listing 6 without further explanation. Pay particular attention to the comments that describe the two for loops.
The program named Dsp036
The program named Dsp036 replicates the behavior of the program named Dsp035, except that it uses an FFT algorithm to perform the inverse Fourier transform instead of using a DFT algorithm as in Dsp035.
The output from Dsp036
The output produced by running the program named Dsp036 and plotting the output using the program named Graph03 is shown in Figure 2.
Figure 2 Forward and inverse transform of a time series using FFT algorithm |
Compare Figure 2 with Figure 1. The two should be identical. The program named Dsp036 was designed to use an FFT algorithm for the inverse Fourier transform and to replicate the behavior of the program named Dsp035, which uses a DFT algorithm for the inverse Fourier transform. In addition, the same plotting parameters were used for both figures.
The code
I'm only going to show you one short code fragment from the program named Dsp036. Listing 7 shows the code that invokes the methods to perform the forward and inverse Fourier transforms using the FFT algorithm. A complete listing of the program named Dsp036 is shown in Listing 16 near the end of the lesson.
//Compute FFT of the time data and save it in // the output arrays. ForwardRealToComplexFFT01.transform( timeDataIn, realSpect, imagSpect, angle, magnitude); //Compute inverse FFT of the spectral data InverseComplexToRealFFT01. inverseTransform( realSpect, imagSpect, timeOut); Listing 7 |
The forward Fourier transform
The transform method used to perform the forward Fourier transform in Listing 7 was discussed in an earlier lesson, so I won't discuss it further here.
The inverse Fourier transform
The static inverseTransform method of the InverseComplexToRealFFT01 class was used to perform the inverse Fourier transform in Listing 7. You can view this method in Listing 17 near the end of the lesson.
I'm not going to discuss this method in detail either, because it is very similar to the method named InverseComplexToReal01 discussed earlier in conjunction with Listing 4 and the listings following that one.
A couple of things to note
There are a couple of things, however, that I do want to point out.
The transform method and the inverseTransform method each invoke a method named complexToComplex to actually perform the Fourier transform. This method implements a classical FFT algorithm accepting complex input data and producing complex output data. The restriction of real-to-complex and complex-to-real is imposed by the methods named transform and inverseTransform.
(The method named complexToComplex is also suitable for use if you have a need to perform complex-to-complex Fourier transforms.)
The signature of the complexToComplex method
The signature for the complexToComplex method is shown in Figure 3.
public static void complexToComplex( int sign, int len, double real[], double imag[]){ Figure 3 |
The complexToComplex method can be used to perform either a forward or an inverse transform. The value of the first parameter determines whether the method performs a forward or an inverse Fourier transform.
The first parameter of the complexToComplex method
A value of +1 for the first parameter causes the complexToComplex method to perform a forward Fourier transform.
A value of -1 for the first parameter causes the complexToComplex method to perform an inverse Fourier transform.
The forward transform
Although I didn't include the code in this lesson, (because it was shown in an earlier lesson), the transform method in Figure 7 passes a value of +1 to the complexToComplex method to cause it to perform a forward Fourier transform.
The inverse transform
Similarly, the inverseTransform method shown in Listing 17 passes a value of -1 to the complexToComplex method to cause it to perform an inverse Fourier transform.
FFT and DFT produce equivalent results
As evidenced in Figure 1 and Figure 2, the program named Dsp035, which uses a DFT algorithm, produces the same results as the program named Dsp036, which uses an FFT algorithm. However, if you were to put a timer on each of the programs, you would find that Dsp036 runs faster due to the improved speed of the FFT algorithm over the DFT algorithm.
Using a Fourier transform to perform frequency filtering
The program named Dsp037 illustrates frequency filtering accomplished by modifying the complex spectrum in the frequency domain and then performing an inverse Fourier transform on the modified frequency-domain data. The results are shown in Figure 4.
Figure 4 Illustrates filtering in the frequency domain |
Operation of the program
The program begins by using an FFT algorithm to perform a forward Fourier transform on a single impulse in the time domain.
(A DFT algorithm could have been used equally as well, but it would have been slower.)
The impulse is shown as the input time series in the topmost plot in Figure 4.
(Although I didn't show the complex spectrum of the impulse, we know that the magnitude of the spectrum of an impulse is constant across all frequencies. In other words, the magnitude spectrum of an impulse is a flat line from zero to the sampling frequency and above.)
Modify the complex spectrum
Then the program eliminates all energy between one-sixth and five-sixths of the sampling frequency by setting the real and imaginary parts of the FFT output to zero.
The second, third, and fourth plots in Figure 4 show the real part, imaginary part, and amplitude respectively of the modified complex spectrum.
(The two boxes in the fourth plot in Figure 4 show what's left of the spectral energy after the energy in the middle of the band has been eliminated.)
The folding frequency in these three plots is near the center of the plot at the eighth tick mark.
The plotting format
The input data length was 256 samples. All but one of the input data values was set to zero resulting in a single impulse in the input time series near the second tick mark in Figure 4.
(The real and complex parts of the frequency spectrum were computed at 256 frequencies between zero and the sampling frequency.)
There were 256 values plotted horizontally in each separate plot. Once again, to make it easier to view the plots on the rightmost end, I plotted the values on a grid that is 270 units wide. This leaves some blank space on the rightmost end to contain the numbers, thus preventing the numbers from being mixed in with the plotted values. The last actual data value coincides with the rightmost tick mark on each plot.
Perform an inverse Fourier transform
After modifying the complex spectrum as described above, the program performs an inverse Fourier transform on the modified complex spectrum to produce the filtered impulse.
The filtered impulse
The filtered impulse is shown as the bottom plot in Figure 4. As you can see, the pulse is smeared out in time relative to the input pulse in the top plot. This is the typical result of reducing the bandwidth of a pulse.
(This particular modification of the complex spectrum resulted in a filtered pulse that has the waveform of a SIN(X)/X function. A different modification of the complex spectrum would have resulted in a filtered pulse with a different waveform.
This example also illustrates one of the miracles of digital signal processing. Energy appears in the output before it occurs in the input. Obviously that is not possible in the real world of analog systems, but many things are possible in the digital world that are not possible in the real world.)
The code for Dsp037
Listing 8 shows the beginning of the class definition for the program named Dsp037.
class Dsp037 implements GraphIntfc01{ final double pi = Math.PI; int len = 256; double[] timeDataIn = new double[len]; double[] realSpect = new double[len]; double[] imagSpect = new double[len]; double[] angle = new double[len];//unused double[] magnitude = new double[len]; double[] timeOut = new double[len]; Listing 8 |
Listing 8 simply declares and initializes some variables that will be used later.
The constructor
The constructor begins in Listing 9.
public Dsp037(){//constructor timeDataIn[32] = 90; Listing 9 |
Listing 9 creates the raw pulse data shown in the topmost plot in Figure 4.
When the array object referred to by timeDataIn is created, the values of all array elements are set to zero by default. Listing 9 modifies one of the elements to have a value of 90. This results in a single impulse at an index of 32.
Compute the Fourier transform
Still in the constructor, the code in Listing 10 uses an FFT algorithm in the method named transform (discussed earlier) to compute the Fourier transform of the impulse.
//Compute FFT of the time data and save it in // the output arrays. ForwardRealToComplexFFT01.transform( timeDataIn, realSpect, imagSpect, angle, magnitude); Listing 10 |
The results of the Fourier transform are stored in the array objects referred to by realSpect, imagSpect, and magnitude.
(The phase angle is also computed but is of no interest in this example.)
Apply the filter to the frequency data
Listing 11 applies the filter by setting sample values in a portion of the real and imaginary parts of the complex spectrum to zero.
for(int cnt = len/6;cnt < 5*len/6;cnt++){ realSpect[cnt] = 0.0; imagSpect[cnt] = 0.0; }//end for loop Listing 11 |
This code eliminates all energy between one-sixth and five-sixths of the sampling frequency. The modified data for the real and imaginary parts of the complex spectrum are shown in the second and third plots in Figure 4.
Recompute the magnitude
Listing 12 recomputes the magnitude values for the modified real and imaginary values of the complex spectrum.
//Recompute the magnitude based on the // modified real and imaginary spectra. for(int cnt = 0;cnt < len;cnt++){ magnitude[cnt] = (Math.sqrt( realSpect[cnt]*realSpect[cnt] + imagSpect[cnt]*imagSpect[cnt])/len); }//end for loop Listing 12 |
The modified data for the amplitude of the complex spectrum are shown in the fourth plot in Figure 4.
Compute the inverse Fourier transform
Listing 13 uses the inverseTransform method to compute the inverse Fourier transform of the modified complex spectrum stored in realSpect and imagSpect. The results of the inverse transform are stored in timeOut.
InverseComplexToRealFFT01.inverseTransform( realSpect, imagSpect, timeOut); }//end constructor Listing 13 |
The results of the inverse transform are shown in the bottom plot in Figure 4.
Listing 13 also signals the end of the constructor.
Display the results
Once the constructor returns, all of the data that is to be plotted has been stored in the various array objects. The remaining code in the program consists of the definition of the six methods required by the interface named GraphIntfc01. These methods are required to make it possible to use the program named Graph03 to plot the results as shown in Figure 4.
I have discussed these methods on numerous previous occasions, and won't repeat that discussion here.
One more example, Dsp038
Figure 5 illustrates one more example of performing frequency filtering by modifying the complex spectrum and then performing an inverse transform on the modified spectrum.
While discussing the program named Dsp037, I told you that performing a different modification on the complex spectrum would result in a different waveform for the filtered impulse. The program named Dsp038 applies a different modification to the complex spectrum, but is otherwise the same as Dsp037.
(Because of the similarity of the two programs, I won't discuss the code in Dsp038. You can view that code in Listing 19 near the end of the lesson.)
Figure 5 shows the output produced by the program named Dsp038.
Figure 5 Illustrates filtering in the frequency domain |
Compare Figure 5 with Figure 4
The basic plotting format of Figure 5 is the same as Figure 4.
The first difference to note between the two figures is that I moved the impulse in the input time series in the topmost plot sixteen samples further to the right in Dsp038.
(This has no impact on the final result, which you can verify by modifying the program to move the impulse to a different position and then compiling and running the modified program.)
Compare the bandwidth of the pass band
The second difference to note is shown in the modified amplitude spectrum in the fourth plot in the two figures. The bandwidth of the pass band is significantly narrower in Figure 5 than in Figure 4. Also, the pass band in Figure 4 extends all the way down to zero frequency, while Figure 5 eliminates all energy below a frequency of three thirty-seconds of the sampling frequency.
Waveforms of filtered impulse
Finally, note the waveforms of the two filtered impulses. The overall amplitude of the filtered impulse in Figure 5 is less than in Figure 4, simply because it contains less total energy. In addition, the filtered impulse in Figure 5 is broader than the filtered impulse in Figure 4. This is because it has a narrower bandwidth.
(Pulses that are narrow in terms of time duration require a wider bandwidth than pulses that have a longer time duration. The time duration of the pulse tends to be inversely related to the bandwidth of the pulse.)
Run the Programs
I encourage you to copy, compile, and run the programs provided in this lesson. Experiment with them, making changes and observing the results of your changes.
Create more complex experiments. For example, use more complex input time series when experimenting with frequency filtering. Apply different modifications to the complex spectrum when experimenting with frequency filtering.
Most of all enjoy yourself and learn something in the process.
Summary
This lesson illustrates and explains forward and inverse Fourier transforms using both DFT and FFT algorithms.
The lesson also illustrates and explains the implementation of frequency filtering by modifying the complex spectrum in the frequency domain and then transforming the modified complex spectrum back into the time domain.
Complete Program Listings
Complete listings of the programs discussed in this lesson follow.
import java.util.*; class Dsp035 implements GraphIntfc01{ final double pi = Math.PI; int len = 256; double[] timeDataIn = new double[len]; double[] realSpect = new double[len]; double[] imagSpect = new double[len]; double[] angle = new double[len];//unused double[] magnitude = new double[len]; double[] timeDataOut = new double[len]; int zero = 0; public Dsp035(){//constructor //Create the raw data pulses timeDataIn[0] = 0; timeDataIn[1] = 50; timeDataIn[2] = 75; timeDataIn[3] = 80; timeDataIn[4] = 75; timeDataIn[5] = 50; timeDataIn[6] = 25; timeDataIn[7] = 0; timeDataIn[8] = -25; timeDataIn[9] = -50; timeDataIn[10] = -75; timeDataIn[11] = -80; timeDataIn[12] = -60; timeDataIn[13] = -40; timeDataIn[14] = -26; timeDataIn[15] = -17; timeDataIn[16] = -11; timeDataIn[17] = -8; timeDataIn[18] = -5; timeDataIn[19] = -3; timeDataIn[20] = -2; timeDataIn[21] = -1; timeDataIn[240] = 80; timeDataIn[241] = 80; timeDataIn[242] = 80; timeDataIn[243] = 80; timeDataIn[244] = -80; timeDataIn[245] = -80; timeDataIn[246] = -80; timeDataIn[247] = -80; timeDataIn[248] = 80; timeDataIn[249] = 80; timeDataIn[250] = 80; timeDataIn[251] = 80; timeDataIn[252] = -80; timeDataIn[253] = -80; timeDataIn[254] = -80; timeDataIn[255] = -80; //Create raw data sinusoid for(int x = len/3;x < 3*len/4;x++){ timeDataIn[x] = 80.0 * Math.sin( 2*pi*(x)*1.0/20.0); }//end for loop //Compute DFT of the time data and save it in // the output arrays. ForwardRealToComplex01.transform(timeDataIn, realSpect, imagSpect, angle, magnitude, zero, 0.0, 1.0); //Compute inverse DFT of spectral data and // save output time data in output array InverseComplexToReal01.inverseTransform( realSpect, imagSpect, timeDataOut); }//end constructor //-------------------------------------------// //The following six methods are required by the // interface named GraphIntfc01. public int getNmbr(){ //Return number of curves to plot. Must not // exceed 5. return 5; }//end getNmbr //-------------------------------------------// public double f1(double x){ int index = (int)Math.round(x); if(index < 0 || index > timeDataIn.length-1){ return 0; }else{ return timeDataIn[index]; }//end else }//end function //-------------------------------------------// public double f2(double x){ int index = (int)Math.round(x); if(index < 0 || index > realSpect.length-1){ return 0; }else{ //scale for convenient viewing return 5*realSpect[index]; }//end else }//end function //-------------------------------------------// public double f3(double x){ int index = (int)Math.round(x); if(index < 0 || index > imagSpect.length-1){ return 0; }else{ //scale for convenient viewing return 5*imagSpect[index]; }//end else }//end function //-------------------------------------------// public double f4(double x){ int index = (int)Math.round(x); if(index < 0 || index > magnitude.length-1){ return 0; }else{ //scale for convenient viewing return 5*magnitude[index]; }//end else }//end function //-------------------------------------------// public double f5(double x){ int index = (int)Math.round(x); if(index < 0 || index > timeDataOut.length-1){ return 0; }else{ return timeDataOut[index]; }//end else }//end function }//end sample class Dsp035 Listing 14 |
/*File InverseComplexToReal01.java Copyright 2004, R.G.Baldwin Rev 5/24/04 Although there are more efficient ways to write this program, it was written the way it was to mimic the behavior of an FFT algorithm. Therefore, the complex input must extend from zero to the sampling frequency. The static method named inverseTransform performs a complex to real inverse discrete Fourier transform returning a real result only. In other words, the method transforms a complex input to a real output. Does not implement the FFT algorithm. Implements a straight-forward sampled-data version of the continuous inverse Fourier transform defined using integral calculus. The parameters are: double[] realIn - incoming real data double[] imagIn - incoming imag data double[] realOut - outgoing real data Considers the data length to be realIn.length Computational time increment is 1.0/realIn.length Returns a number of points equal to the data length. Assumes real input consists of positive frequency points for a symmetric real frequency function. That is, the real input is assumed to be symmetric about the folding frequency. Does not test this assumption. Assumes imaginary input consists of positive frequency points for an asymmetric imaginary frequency function. That is, the imaginary input is assumed to be asymmetric about the folding frequency. Does not test this assumption. The assumption of a symmetric real part and an asymmetric imaginary part guarantees that the imaginary output would be all zero if it were to be computed. Thus the program makes no attempt to compute an imaginary output. Tested using J2SE v1.4.2 under WinXP. ************************************************/ public class InverseComplexToReal01{ public static void inverseTransform( double[] realIn, double[] imagIn, double[] realOut){ int dataLen = realIn.length; double delT = 1.0/realIn.length; double startTime = 0.0; //Outer loop interates on time domain // values. for(int i=0; i < dataLen;i++){ double time = startTime + i*delT; double real = 0; //Inner loop iterates on frequency // domain values. for(int j=0; j < dataLen; j++){ real += realIn[j]* Math.cos(2*Math.PI*time*j) + imagIn[j]* Math.sin(2*Math.PI*time*j); }//end inner loop realOut[i] = real; }//end outer loop }//end inverseTransform }//end class InverseComplexToReal01 Listing 15 |
/* File Dsp036.java Copyright 2004, R.G.Baldwin Revised 5/24/04 Illustrates forward and inverse Fourier transforms using FFT algorithms. Performs spectral analysis on a time series consisting of pulses and a sinusoid. Passes resulting real and complex parts to inverse Fourier transform program to reconstruct the original time series. Run with Graph03. Tested using J2SE 1.4.2 under WinXP. ************************************************/ import java.util.*; class Dsp036 implements GraphIntfc01{ final double pi = Math.PI; int len = 256; double[] timeDataIn = new double[len]; double[] realSpect = new double[len]; double[] imagSpect = new double[len]; double[] angle = new double[len];//unused double[] magnitude = new double[len]; double[] timeOut = new double[len]; public Dsp036(){//constructor //Create the raw data pulses timeDataIn[0] = 0; timeDataIn[1] = 50; timeDataIn[2] = 75; timeDataIn[3] = 80; timeDataIn[4] = 75; timeDataIn[5] = 50; timeDataIn[6] = 25; timeDataIn[7] = 0; timeDataIn[8] = -25; timeDataIn[9] = -50; timeDataIn[10] = -75; timeDataIn[11] = -80; timeDataIn[12] = -60; timeDataIn[13] = -40; timeDataIn[14] = -26; timeDataIn[15] = -17; timeDataIn[16] = -11; timeDataIn[17] = -8; timeDataIn[18] = -5; timeDataIn[19] = -3; timeDataIn[20] = -2; timeDataIn[21] = -1; timeDataIn[240] = 80; timeDataIn[241] = 80; timeDataIn[242] = 80; timeDataIn[243] = 80; timeDataIn[244] = -80; timeDataIn[245] = -80; timeDataIn[246] = -80; timeDataIn[247] = -80; timeDataIn[248] = 80; timeDataIn[249] = 80; timeDataIn[250] = 80; timeDataIn[251] = 80; timeDataIn[252] = -80; timeDataIn[253] = -80; timeDataIn[254] = -80; timeDataIn[255] = -80; //Create raw data sinusoid for(int x = len/3;x < 3*len/4;x++){ timeDataIn[x] = 80.0 * Math.sin( 2*pi*(x)*1.0/20.0); }//end for loop //Compute FFT of the time data and save it in // the output arrays. ForwardRealToComplexFFT01.transform( timeDataIn, realSpect, imagSpect, angle, magnitude); //Compute inverse FFT of spectral data InverseComplexToRealFFT01. inverseTransform( realSpect, imagSpect, timeOut); }//end constructor //-------------------------------------------// //The following six methods are required by the // interface named GraphIntfc01. public int getNmbr(){ //Return number of curves to plot. Must not // exceed 5. return 5; }//end getNmbr //-------------------------------------------// public double f1(double x){ int index = (int)Math.round(x); if(index < 0 || index > timeDataIn.length-1){ return 0; }else{ return timeDataIn[index]; }//end else }//end function //-------------------------------------------// public double f2(double x){ int index = (int)Math.round(x); if(index < 0 || index > realSpect.length-1){ return 0; }else{ //scale for convenient viewing return 5*realSpect[index]/len; }//end else }//end function //-------------------------------------------// public double f3(double x){ int index = (int)Math.round(x); if(index < 0 || index > imagSpect.length-1){ return 0; }else{ //scale for convenient viewing return 5*imagSpect[index]/len; }//end else }//end function //-------------------------------------------// public double f4(double x){ int index = (int)Math.round(x); if(index < 0 || index > magnitude.length-1){ return 0; }else{ //scale for convenient viewing return 5*magnitude[index]; }//end else }//end function //-------------------------------------------// public double f5(double x){ int index = (int)Math.round(x); if(index < 0 || index > timeOut.length-1){ return 0; }else{ //scale for convenient viewing return timeOut[index]/len; }//end else }//end function }//end sample class Dsp036 Listing 16 |
/*File InverseComplexToRealFFT01.java Copyright 2004, R.G.Baldwin Rev 5/24/04 The static method named inverseTransform performs a complex to real Fourier transform using a complex-to-complex FFT algorithm. A specific parameter is passed to the FFT algorithm that causes this to be an inverse Fourier transform. See InverseComplexToReal01 for a version that does not use an FFT algorithm but uses a DFT algorithm instead. Incoming parameters are: double[] realIn - incoming real data double[] imagIn - incoming imaginary data double[] realOut - outgoing real data Requires spectral input data extending from zero to the sampling frequency. Assumes real input consists of positive frequency points for a symmetric real frequency function. That is, the real input is assumed to be symmetric about the folding frequency. Does not test this assumption. Assumes imaginary input consists of positive frequency points for an asymmetric imaginary frequency function. That is, the imaginary input is assumed to be asymmetric about the folding frequency. Does not test this assumption. The assumption of a symmetric real part and an asymmetric imaginary part guarantees that the imaginary output is all zeros. Thus, the program does not return an imaginary output. Does not test the assumption that the imaginary is all zeros. CAUTION: THE INCOMING DATA LENGTH MUST BE A POWER OF TWO. OTHERWISE, THIS PROGRAM WILL FAIL TO RUN PROPERLY. Returns a number of points equal to the incoming data length. Those points are uniformly distributed beginning at zero. ************************************************/ public class InverseComplexToRealFFT01{ public static void inverseTransform( double[] realIn, double[] imagIn, double[] realOut){ double pi = Math.PI;//for convenience int dataLen = realIn.length; double[] imagOut = new double[dataLen]; //The complexToComplex FFT method does an // in-place transform causing the output // complex data to be stored in the arrays // containing the input complex data. // Therefore, it is necessary to copy the // input data into the output arrays before // passing them to the FFT algorithm. System.arraycopy(realIn,0,realOut,0,dataLen); System.arraycopy(imagIn,0,imagOut,0,dataLen); //Perform the spectral analysis. The results // are stored in realOut and imagOut. Note // that the -1 value for the first // parameter causes the transform to be an // inverse transform. A +1 value would cause // it to be a forward transform. complexToComplex(-1,dataLen,realOut,imagOut); }//end inverseTransform method //-------------------------------------------// //This method computes a complex-to-complex // FFT. The value of sign must be 1 for a // forward FFT and -1 for an inverse FFT. public static void complexToComplex( int sign, int len, double real[], double imag[]){ double scale = 1.0; //Reorder the input data into reverse binary // order. int i,j; for (i=j=0; i < len; ++i) { if (j>=i) { double tempr = real[j]*scale; double tempi = imag[j]*scale; real[j] = real[i]*scale; imag[j] = imag[i]*scale; real[i] = tempr; imag[i] = tempi; }//end if int m = len/2; while (m>=1 && j>=m) { j -= m; m /= 2; }//end while loop j += m; }//end for loop //Input data has been reordered. int stage = 0; int maxSpectraForStage,stepSize; //Loop once for each stage in the spectral // recombination process. for(maxSpectraForStage = 1, stepSize = 2*maxSpectraForStage; maxSpectraForStage < len; maxSpectraForStage = stepSize, stepSize = 2*maxSpectraForStage){ double deltaAngle = sign*Math.PI/maxSpectraForStage; //Loop once for each individual spectra for (int spectraCnt = 0; spectraCnt < maxSpectraForStage; ++spectraCnt){ double angle = spectraCnt*deltaAngle; double realCorrection = Math.cos(angle); double imagCorrection = Math.sin(angle); int right = 0; for (int left = spectraCnt; left < len;left += stepSize){ right = left + maxSpectraForStage; double tempReal = realCorrection*real[right] - imagCorrection*imag[right]; double tempImag = realCorrection*imag[right] + imagCorrection*real[right]; real[right] = real[left]-tempReal; imag[right] = imag[left]-tempImag; real[left] += tempReal; imag[left] += tempImag; }//end for loop }//end for loop for individual spectra maxSpectraForStage = stepSize; }//end for loop for stages }//end complexToComplex method }//end class InverseComplexToRealFFT01 Listing 17 |
/* File Dsp037.java Copyright 2004, R.G.Baldwin Revised 5/24/04 Illustrates filtering in the frequency domain. Performs FFT on an impulse. Eliminates all energy between one-sixth and five-sixths of the sampling frequency by modifying the real and imaginary parts of the FFT output. Then performs inverse FFT to produce the filtered impulse. Run with Graph03. Tested using J2SE 1.4.2 under WinXP. ************************************************/ import java.util.*; class Dsp037 implements GraphIntfc01{ final double pi = Math.PI; int len = 256; double[] timeDataIn = new double[len]; double[] realSpect = new double[len]; double[] imagSpect = new double[len]; double[] angle = new double[len];//unused double[] magnitude = new double[len]; double[] timeOut = new double[len]; public Dsp037(){//constructor //Create the raw data pulse timeDataIn[32] = 90; //Compute FFT of the time data and save it in // the output arrays. ForwardRealToComplexFFT01.transform( timeDataIn, realSpect, imagSpect, angle, magnitude); //Apply the frequency filter eliminating all // energy between one-sixth and five-sixths // of the sampling frequency by modifying the // real and imaginary parts of the spectrum. for(int cnt = len/6;cnt < 5*len/6;cnt++){ realSpect[cnt] = 0.0; imagSpect[cnt] = 0.0; }//end for loop //Recompute the magnitude based on the // modified real and imaginary spectra. for(int cnt = 0;cnt < len;cnt++){ magnitude[cnt] = (Math.sqrt( realSpect[cnt]*realSpect[cnt] + imagSpect[cnt]*imagSpect[cnt])/len); }//end for loop //Compute inverse FFT of modified spectral // data. InverseComplexToRealFFT01.inverseTransform( realSpect, imagSpect, timeOut); }//end constructor //-------------------------------------------// //The following six methods are required by the // interface named GraphIntfc01. public int getNmbr(){ //Return number of curves to plot. Must not // exceed 5. return 5; }//end getNmbr //-------------------------------------------// public double f1(double x){ int index = (int)Math.round(x); if(index < 0 || index > timeDataIn.length-1){ return 0; }else{ return timeDataIn[index]; }//end else }//end function //-------------------------------------------// public double f2(double x){ int index = (int)Math.round(x); if(index < 0 || index > realSpect.length-1){ return 0; }else{ return realSpect[index]; }//end else }//end function //-------------------------------------------// public double f3(double x){ int index = (int)Math.round(x); if(index < 0 || index > imagSpect.length-1){ return 0; }else{ return imagSpect[index]; }//end else }//end function //-------------------------------------------// public double f4(double x){ int index = (int)Math.round(x); if(index < 0 || index > magnitude.length-1){ return 0; }else{ //scale for convenient viewing return len*magnitude[index]; }//end else }//end function //-------------------------------------------// public double f5(double x){ int index = (int)Math.round(x); if(index < 0 || index > timeOut.length-1){ return 0; }else{ //scale for convenient viewing return 3.0*timeOut[index]/len; }//end else }//end function }//end sample class Dsp037 Listing 18 |
/* File Dsp038.java Copyright 2004, R.G.Baldwin Revised 5/24/04 Illustrates filtering in the frequency domain. Performs FFT on an impulse. Modifies the complex spectrum. Then performs inverse FFT to produce the filtered impulse. Run with Graph03. Tested using J2SE 1.4.2 under WinXP. ************************************************/ import java.util.*; class Dsp038 implements GraphIntfc01{ final double pi = Math.PI; int len = 256; double[] timeDataIn = new double[len]; double[] realSpect = new double[len]; double[] imagSpect = new double[len]; double[] angle = new double[len];//unused double[] magnitude = new double[len]; double[] timeOut = new double[len]; public Dsp038(){//constructor //Create the raw data pulse timeDataIn[64] = 90; //Compute FFT of the time data and save it in // the output arrays. ForwardRealToComplexFFT01.transform( timeDataIn, realSpect, imagSpect, angle, magnitude); //Apply the frequency filter. for(int cnt = 0;cnt <= len/2;cnt++){ if(cnt < 3*len/32){ realSpect[cnt] = 0; imagSpect[cnt] = 0; }//end if if(cnt > 5*len/32){ realSpect[cnt] = 0; imagSpect[cnt] = 0; }//end if //Fold complex spectral data if(cnt > 0){ realSpect[len - cnt] = realSpect[cnt]; }//end if if(cnt > 0){ imagSpect[len - cnt] = -imagSpect[cnt]; }//end if }//end for loop //Recompute the magnitude based on the // modified real and imaginary spectra. for(int cnt = 0;cnt < len;cnt++){ magnitude[cnt] = (Math.sqrt( realSpect[cnt]*realSpect[cnt] + imagSpect[cnt]*imagSpect[cnt])/len); }//end for loop //Compute inverse FFT of modified spectral // data. InverseComplexToRealFFT01.inverseTransform( realSpect, imagSpect, timeOut); }//end constructor //-------------------------------------------// //The following six methods are required by the // interface named GraphIntfc01. public int getNmbr(){ //Return number of curves to plot. Must not // exceed 5. return 5; }//end getNmbr //-------------------------------------------// public double f1(double x){ int index = (int)Math.round(x); if(index < 0 || index > timeDataIn.length-1){ return 0; }else{ return timeDataIn[index]; }//end else }//end function //-------------------------------------------// public double f2(double x){ int index = (int)Math.round(x); if(index < 0 || index > realSpect.length-1){ return 0; }else{ return realSpect[index]; }//end else }//end function //-------------------------------------------// public double f3(double x){ int index = (int)Math.round(x); if(index < 0 || index > imagSpect.length-1){ return 0; }else{ return imagSpect[index]; }//end else }//end function //-------------------------------------------// public double f4(double x){ int index = (int)Math.round(x); if(index < 0 || index > magnitude.length-1){ return 0; }else{ //scale for convenient viewing return len*magnitude[index]; }//end else }//end function //-------------------------------------------// public double f5(double x){ int index = (int)Math.round(x); if(index < 0 || index > timeOut.length-1){ return 0; }else{ //scale for convenient viewing return 3.0*timeOut[index]/len; }//end else }//end function }//end sample class Dsp038 Listing 19 |
Copyright 2004, Richard G. Baldwin. Reproduction in whole or in part in any form or medium without express written permission from Richard Baldwin is prohibited.
About the author
Richard Baldwin is a college professor (at Austin Community College in Austin, TX) and private consultant whose primary focus is a combination of Java, C#, and XML. In addition to the many platform and/or language independent benefits of Java and C# applications, he believes that a combination of Java, C#, and XML will become the primary driving force in the delivery of structured information on the Web.Richard has participated in numerous consulting projects, and he frequently provides onsite training at the high-tech companies located in and around Austin, Texas. He is the author of Baldwin's Programming Tutorials, which has gained a worldwide following among experienced and aspiring programmers. He has also published articles in JavaPro magazine.
In addition to his programming expertise, Richard has many years of practical experience in Digital Signal Processing (DSP). His first job after he earned his Bachelor's degree was doing DSP in the Seismic Research Department of Texas Instruments. (TI is still a world leader in DSP.) In the following years, he applied his programming and DSP expertise to other interesting areas including sonar and underwater acoustics.
Richard holds an MSEE degree from Southern Methodist University and has many years of experience in the application of computer technology to real-world problems.
-end-
This article was originally published on November 16, 2004