http://www.developer.com/
Adaptive Filtering in Java, Getting StartedSeptember 20, 2005 Java Programming, Notes # 2350
PrefaceDSP and adaptive filtering With the decrease in cost and the increase in speed of digital devices, Digital Signal Processing (DSP) is showing up in everything from cell phones to hearing aids to rock concerts. Many applications of DSP are static. That is, the characteristics of the digital processor don't change with time or circumstances. However, a particularly interesting branch of DSP is adaptive filtering. This is a situation where the characteristics of the digital processor change with time, circumstances, or both. First in a series This is the first lesson in a series designed to teach you about adaptive filtering in Java. This lesson will introduce you to the topic by showing you how to write a Java program to solve a relatively simple time-adaptive filtering problem for which the correct solution is well known in advance. This will make it possible to check the adaptive solution against the known correct solution. An adaptive whitening filter The next lesson will show you how to write an adaptive whitening filter program in Java, which is conceptually more difficult than the filter that I will explain in this lesson. The next lesson will also show you how to use the whitening filter to extract wide band signal from a channel in which the signal is corrupted by one or more components of narrow band noise. More general adaptive filtering considerations Following that, the lessons in the series will become somewhat more general. I plan to publish lessons that explain and provide examples for the four common scenarios in which adaptive filtering is used:
Somewhere along the way I will probably also publish a lesson that explains and illustrates the difference between least mean square (LMS) and recursive least squares (RLS) adaptive algorithms. 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. PreviewA time-delay filter with a flat amplitude response The program that I will present and explain in this lesson illustrates an aspect of adaptive filtering for which the correct solution is already well known. The program adaptively designs a time-delay filter with a flat amplitude response and a linear phase response in the frequency domain. A straightforward but useful scenario for learning Although this is a relatively straightforward scenario, it is also a useful scenario for learning purposes. To begin with, the program illustrates the use of a least mean square (LMS) adaptive algorithm in a relatively simple setting, making it easy to understand what the algorithm is doing. In addition, the program teaches you about the use of digital delay lines, a topic with which you may not yet be familiar. Beyond that, it is easy to confirm that the adaptive solution matches the known correct solution. User experimentation is encouraged When running this program, the user provides several parameters that have an impact on the adaptive process. This allows the user to experiment with the adaptive process comparing results for different input parameters. Two channels of input data Two sampled time series, chanA and chanB, are presented to the adaptive processing system. Each time series consists of the same wide band signal plus white noise that is uncorrelated between the two channels. A time shift between the two channels On the basis of user input, the signal in chanB can be delayed or advanced by up to six samples relative to the signal in chanA. In other words, the time base for chanB can be shifted in either direction causing chanB to lead chanA in time, or causing chanB to lag chanA in time.
Filtering chanA A nine-point convolution operator is developed adaptively and applied to chanA. The purpose of the adaptive process is to cause the filtered output to be in time registration with chanB.
How do you measure success? When the adaptive process converges successfully, the time series produced by applying the convolution operator to chanA matches the signal on chanB. It is already well known that the correct solution to this problem is a finite impulse response (FIR) convolution filter in which one coefficient has a value of 1 and all the other coefficients have a value of 0. The location of the coefficient having the value of 1 is such as to cause the result of filtering chanA to be either advanced or delayed in time by a number of samples that causes it to be in time registration with chanB. It is also already well known that the Fourier transform of the filter described above will have a flat amplitude response and a linear phase response. If the adaptive process converges to this result, it is successful. User inputs The user provides the following information as command line parameters:
The first example was run using default parameters. This produced the following output on the command line screen. (Note that a line break was manually entered into the first line to force it to fit into this narrow publication format.) Usage: java Adapt01 timeShift feedbackGain noiseLevel numberIterations Negative timeShift is delay Using -4 sample shift by default Using 0.001 feedbackGain by default noiseLevel is a decimal fraction Using 0.0 by default numberIterations is an int Using 100 by default Graphic output Figure 1 shows the first of three graphic outputs that are produced by this program. (The other two graphic outputs are shown later in Figure 2.) The following four time series are plotted in color in Figure 1 showing the convergence (or lack thereof) of the adaptive algorithm:
Figure 1 Traces wrap around and down When the top four traces reach the right end of the plotting area in Figure 1, they wrap around and down resulting in four new traces further down the page. Thus, the bottom four traces in Figure 1 are the continuation of the right end of the top four traces. Was the adaptive process successful? If the adaptive process is successful for this problem, the green (error) trace should go to zero, and the red (filter output) trace should match the blue (target) trace. As you can see, the adaptive process converges to this solution about half way across the top four traces in Figure 1. Thus, the adaptive process was successful. I will have more to say about Figure 1 later when I explain the code. Impulse response plots The second of the three graphic outputs produced by this program is shown in the left column of Figure 2. The impulse response of the adaptive convolution filter at the beginning and at the end of every tenth iteration is shown in the left column of Figure 2. Thus, the changes in the shape of the impulse response can be viewed from the beginning to the end of the adaptive process. The progressive stages of the impulse response are shown from top to bottom in the left column in Figure 2. Figure 2
Initial impulse response at the top The filter is initialized with a single coefficient value of 1 at the center and 0 for all of the other eight coefficient values. The correct solution is a single coefficient value of 1 at a location in the convolution filter that matches the time shift between chanA and the target, chanB. For the case shown in Figure 2, the target was delayed by four samples relative to chanA. The peak shifts to the left As you can see from Figure 2, the impulse response begins at the top with a value of 1 in the center and values of 0 elsewhere. As the adaptive process progresses down the page, the peak value in the impulse response progresses from the center of the impulse response to the correct location at the left end. In other words, the impulse response modifies itself such that after about 70 iterations (the eighth impulse response plot), it has a value of 1 at the left end and zeros elsewhere. The frequency response The third graphic output produced by this program is shown in the right column of Figure 2. The frequency response of the convolution filter at the beginning and at the end of every tenth iteration is computed and displayed in the right column of Figure 2. Both the amplitude response and the phase response of the filter are displayed, with the amplitude response being plotted above the phase response in the right column of Figure 2. The frequency-domain plots in the right column extend from zero frequency on the left to the Nyquist folding frequency, (which is one-half the sampling frequency), on the right. Frequency response on the right corresponds to impulse response on the left Each pair of plots in the right column immediately to the right of a single impulse response in the left column consists of the amplitude and phase responses of the corresponding impulse response. The amplitude response is plotted in black and the phase response is plotted in red below the amplitude response. Frequency response starts out flat At the beginning (at the top), the impulse response consists of a single impulse in the center of the convolution filter.
At that point in the adaptive process, both the amplitude response and the phase response of the convolution filter are flat across the entire frequency spectrum. Deviations appear by the twentieth iteration By the third impulse response going down the page (twenty adaptive iterations), the impulse response is no longer represented by a single impulse, and some deviation from flatness is apparent in the corresponding amplitude response. By the fourth impulse response (thirty adaptive iterations), the shape of the impulse response has changed considerably and quite a lot of activity is apparent in the corresponding amplitude response and phase response. An adaptive solution in 90 iterations By the eighth impulse response (70 iterations), the impulse response has become a single time-shifted impulse, the amplitude response has returned to being flat across the frequency spectrum, and the phase response has taken on a saw tooth character. This is the correct solution as described above. A flat amplitude response The fact that the impulse response is once again a single impulse means that the amplitude response of the convolution filter is flat across the entire frequency spectrum. The new location (relative to the zero time origin) of the single impulse in the convolution filter causes the output of the filter to be shifted in time relative to its input. As evidenced in Figure 1, this causes the convolution filter output (red) to be in time registration with the target (blue). At that point, the error (green) has converged to zero. A linear phase response A convolution operator that produces a simple time shift is represented in the frequency domain by a flat amplitude response and a linear relationship between phase and frequency. In other words, the phase is a straight line that goes through the zero frequency origin. The slope of the line indicates the direction of the time shift. The magnitude of the slope indicates the amount of the time shift. Because the phase response is plotted in the range from -180 degrees to +180 degrees, and wraps around whenever it exceeds either of those limits, a linear phase shift has a saw tooth appearance when plotted as shown near the bottom of Figure 2. Testing This program was tested using J2SE 5.0 and WinXP. J2SE 5.0 or later is required. Discussion and Sample CodeThe program named Adapt01 I will discuss and explain this program in fragments. A complete listing of the program is provided in Listing 32 near the end of the lesson. The beginning of the class and the beginning of the main method is shown in Listing 1.
The code in the main method in Listing 1 deals with the command line parameters. If the user enters command line parameters, those parameters are used in the adaptive process. If the user doesn't enter command line parameters, a set of default parameters are used in the adaptive process. The code in Listing 1 is completely straightforward and shouldn't require explanation. Perform the adaptive process The code in Listing 2 instantiates an object of the Adapt01 class and invokes the process method to cause the adaptive process to be performed. As described above, the parameters passed to the process method are either provided by the user as command line parameters, or provided by the program by default if the user doesn't enter command line parameters.
Listing 2 also signals the end of the main method. The parameters for Figure 1 and Figure 2 The case that produced the output shown in Figure 1 and Figure 2 was run without entering command line parameters. Hence, the adaptive process was executed using the following default parameters, which were displayed by the code in the main method: Using -4 sample shift by default Using 0.001 feedbackGain by default noiseLevel is a decimal fraction Using 0.0 by default numberIterations is an int Using 100 by default The process method of the Adapt01 class The process method of the Adapt01 class begins in Listing 3. The code in Listing 3 begins by declaring and populating a nine-element array of type double containing the initial convolution filter. This is the filter that is shown by the impulse response at the top of the left column in Figure 2. The contents of this array will be adaptively modified as the program executes.
Create three tapped delay line objects Then the code in Listing 3 declares three array objects, which will be used as tapped delay lines for the rawData, chanA, and chanB. A delay line is similar to a queue. For example, think of a short section of hose that is full of colored marbles. When you insert a new marble in one end, a marble gets pushed out and discarded from the other end. Now think of cutting a series of small holes along the hose through which you can see the color of the marble next to each hole. You can think of these holes as taps from which you can extract the color of the marbles as they progress past the holes. Each of the delay lines in this program is implemented using an array object. Data enters the delay line at the topmost element and moves to the next lower element once during each iteration. The data value that is in element 0 during a particular iteration is discarded and replaced by the value from element 1 during the next iteration. The ability for the code to access any individual element provides taps by which the contents at any stage in the delay line can be retrieved by the code. Instantiate a plotting object for time series data The code in listing 4 instantiates an object of the class named PlotALot05, which provides the ability to plot the time series data in the format shown in Figure 1.
The class named PlotALot05 is a simple extension of the class named PlotALot04, which I explained in the lesson entitled Plotting Large Quantities of Data using Java. I will refer you to that lesson for a general explanation of the class.
Output on the command-line screen The parameters passed to the constructor for the class caused the constructor to display the following information about the plotting object on the command line screen. Title: Time Series Frame width: 398 Frame height: 250 Page width: 390 Page height: 223 Trace spacing: 25 Sample spacing: 5 Traces per page: 8 Samples per page: 156 Instantiate a plotting object for frequency response data The code in Listing 5 instantiates a plotting object for two channels of frequency response data. One channel is used to plot the amplitude response in db and the other channel is used to plot the phase on a scale that extends from -180 degrees to +180 degrees.
The class named PlotALot03 is one of the classes that I explained in the earlier lesson entitled Plotting Large Quantities of Data using Java. (I will refer you back to that lesson for a copy of the source code for this class.) The parameters that describe this plotting object are given below.
Title: Freq Frame width: 264 Frame height: 487 Page width: 256 Page height: 460 Trace spacing: 20 Sample spacing: 2 Traces per page: 22 Samples per page: 1408 Instantiate a plotting object for the impulse response data The code in Listing 6 instantiates a plotting object to display the impulse response of the convolution filter at intervals during the adaptive process. I explained the class named PlotALot01 in the earlier lesson entitled Plotting Large Quantities of Data using Java. (Once again, I will refer you back to that lesson for a copy of the source code for this class.)
The actual plotting object parameters The actual parameters that describe this plotting object are: Title: Filter Frame width: 112 Frame height: 487 Page width: 104 Page height: 460 Trace spacing: 40 Sample spacing: 4 Traces per page: 11 Samples per page: 286 The Frame width is different Note that the actual Frame width of 112 is different from the specified Frame width of 44 in Listing 6. This is because the minimum allowable width for an AWT Frame object is 112 pixels under WinXP regardless of the specified width.
Because the actual Frame width doesn't match the specified Frame width, the code in Listing 6 won't cause the plotting process to synchronize properly and to plot a single impulse response on each axis for filter lengths less than 25 coefficients. To compensate for this problem, the code that feeds the filter data to the plotting object later in the program extends the length of the filter to cause it to synchronize and to plot one impulse response on each axis.
Display frequency response of the filter Listing 7 invokes the method named displayFreqResponse to compute and display the frequency response of the initial convolution filter at 128 points between zero frequency and the Nyquist folding frequency.
At this point, I am going to set the discussion of the process method aside momentarily while I explain the method named displayFreqResponse. I will return to the discussion of the process method shortly. The displayFreqResponse method The displayFreqResponse method, which begins in Listing 8, receives a reference to a double array containing a convolution filter along with a reference to a plotting object capable of plotting two channels of data. The method also receives a value specifying the number of frequencies at which a discrete Fourier transform (DFT) is to be performed on the filter, along with the sample number that represents the zero time location in the filter. The method uses this information to perform a DFT on the filter from zero to the Nyquist folding frequency. It feeds the resulting amplitude spectrum and the phase spectrum to the plotting object for plotting later.
Performing the discrete Fourier transform The process method uses a static method named transform belonging to a class named ForwardRealToComplex01 to perform the Fourier transform. I explained that class and method in detail in the earlier lesson entitled Spectrum Analysis using Java, Sampling Frequency, Folding Frequency, and the FFT Algorithm. Description of the transform method The static method named transform performs a real to complex Fourier transform. The method does not implement the FFT algorithm. Rather, it implements a straightforward sampled data version of the continuous Fourier transform defined using integral calculus. The return values The method returns the following:
The transform method parameters The method parameters for the transform method are:
Frequency increment, magnitude spectrum, and number of returned values The computational frequency increment is the difference between the high and low limits divided by the length of the magnitude array. The magnitude (amplitude) is computed as the square root of the sum of the squares of the real and imaginary parts. This value is divided by the incoming data length, which is given by data.length. The method returns a number of points in the frequency domain equal to the incoming data length regardless of the high and low frequency limits. Prepare the required array objects Listing 9 creates the five array objects required as input parameters by the transform method. Then Listing 9 uses the arraycopy method of the System class to copy the incoming filter data into the array object referred to by timeDataIn.
Note that the length of the timeDataIn array is about fourteen times greater than the length of the actual convolution filter. This results in a small frequency increment value producing smooth representations of the frequency response curves in Figure 2. Compute the frequency response of the convolution filter Listing 10 invokes the static transform method of the ForwardRealToComplex01 class to compute and return the frequency response of the convolution operator from zero to the Nyquist folding frequency.
You can compare the parameter values in Listing 10 with the description of the parameters given earlier. We will discard all of the transform results other than the amplitude response returned in the array referred to by magnitude, and the phase response returned in the array referred to by angle. Prepare to convert amplitude response to log base 10 That's all there is to getting the amplitude and phase response of the convolution filter. The remaining code in the displayFreqResponse method is concerned with plotting the data in the format shown in the right column of Figure 2. The amplitude response will be converted to decibels before plotting. This requires converting the amplitude data to log base 10. I will use the log10 method of the Math class to perform the conversion. Some possible amplitude response values are not compatible with the log10 method. I invite you to examine Sun's documentation for that method to see what I mean. The code in Listing 11 converts those values to compatible values if they exist in the amplitude response.
Convert amplitude response to log base 10 The code in Listing 12 invokes the log10 method in a for loop to convert each amplitude response value to log base 10.
Because this is an in-place conversion, all future references to the array referred to by magnitude will be referring to data values that have been converted to log base 10. This data can be thought of scaled decibels. Normalize the data for plotting One of the difficulties that is always encountered in plotting large quantities of data is the problem of scaling the data so that it fits in the range allocated to the plotting space. The code in Listing 13 begins by finding the absolute peak value of the log base 10 amplitude response data. Then it scales that data by a factor of 50, divided by the peak value, and adds a constant value of 50 producing the results shown in Figure 2.
Feed the data to the plotting object The code in Listing 14 feeds the normalized decibel data and the phase response data to the plotting object where it will be plotted later. The phase data ranges from -180 to +180 degrees. This data is divided by a factor of 20 to make it compatible with the plotting format being used.
Listing 14 signals the end of the displayFreqResponse method. At this point, I will resume my discussion of the process method of the Adapt01 class. Display impulse response of the initial convolution filter Back in the process method, Listing 15 feeds the values of the initial convolution filter (see Listing 3) to one of the plotting objects (see Listing 6) for plotting as shown at the top of the left column in Figure 2.
Extend the convolution filter Recall from the lesson entitled Plotting Large Quantities of Data using Java that this plotting object knows nothing about the end of one convolution filter and the beginning of the next. Rather, the plotting object is designed to simply accept large quantities of data values to be plotted and to cause the plot to automatically wrap from the current axis down to the next axis when the right side of the current axis is encountered. Adjust width of the plotting surface to match length of impulse response Normally, to cause each impulse response to appear on a separate axis, I would adjust the width of the plotting surface (Page width) to cause it to match the length of the impulse response. However, an AWT Frame object has a minimum allowable width of 112 pixels under WinXP, regardless of the width that you specify in the setSize method. Therefore, it was not possible to cause the width of the plotting surface to match the length of the impulse response in this case. Extend the length of the impulse response Therefore, the code in Listing 16 (working in conjunction with the code in Listing 6) extends the length of each impulse response to cause it to match the width of the plotting surface.
A constant value of 2.5 I extended the impulse response with a constant value of 2.5 to cause the extended portion to be visually separable from the actual impulse response values. The extended portion is shown by the broad flat areas to the right of the impulse response plots in the left column of Figure 2. Declare and initialize adaptive variables Listing 17 declares and initializes some variables that are used later in the adaptive process.
Execute the adaptive process Listing 18 shows the beginning of a for loop that creates the data and performs the adaptive process on that data. The body of this loop executes once for each adaptive iteration specified by the user as numberIterations (see Listing 1).
Create the wide bandwidth signal data Listing 18 uses a random number generator to get a random value to serve as a wide bandwidth signal sample. The random number generator produces random values having a flat distribution from 0 to 1.0. Listing 18 subtracts 0.5 from the random value causing the resulting distribution to range from -0.5 to +0.5. Insert data into the rawData delay line Listing 18 passes the random value to the method named flowLine. Listing 18 also passes a reference to an array object named rawData to the flowLine method.
At this point, I will put the discussion of the process method aside for momentarily and explain the method named flowLine. I will return to the explanation of the process method shortly. The flowLine method The flowLine method, which is shown in its entirety in Listing 19, is a simple utility method that causes an array object to behave as a tapped delay line.
The flowLine method receives a reference to an array and a reference to a value of type double. It discards the value at index 0 of the array, moves all the other values down by one element toward index 0, and inserts the new value at the top of the array. That's all that there is to the flowLine method, so I will resume my discussion of the process method at this point. Getting raw data for chanA The purpose of the rawData delay line is to make it possible to get wide band data representing chanA and chanB where chanB can be shifted in time by up to plus or minus six samples relative to chanA. The array used for the rawData delay line has a length of thirteen elements (see Listing 3). In this case, the value at element index 6 is considered to represent the zero time origin. Therefore, the value at index 0 represents a time delay of six samples and the value at index 12 represents a time advance of six samples. Feed signal plus noise into chanA Listing 20 extracts the middle sample (at zero time) from the rawData delay line, adds some random noise to it, and inserts it into a nine-element delay line (see Listing 3) containing the data for chanA.
The amount of additive random noise is controlled by the value of noiseLevel, which is provided by the user as an input parameter. Feed signal plus noise into chanB Listing 21 extracts data with a timeShift from the rawData delay line, adds some random noise, and inserts it into a delay line (see Listing 3) containing the data for chanB.
Once again, the level of the additive noise is controlled by the value of noiseLevel described earlier. The time shift is controlled by the value of timeShift, which is provided by the user as a command line parameter.
Get the input (chanA) data for plotting Listing 22 gets and saves the value at the middle of the chanA delay line for plotting later. These are the values that are plotted in black in Figure 1.
Time is relative When processing sampled data in a computer, (assuming that you are willing to suffer an overall time delay), time is relative. You can consider any sample in a time series to represent zero time. At this point, we will shift our thinking and consider the center taps for the chanA and chanB delay lines to represent zero time. Thus, each of those delay lines contains the sample value at the current time, four samples of historical data, and four samples of future data. Why should we think about it this way? The primary motivation for thinking about it this way is to cause the output from a convolution filter with a single impulse at the center to represent a time series with no time shift in either direction. Not possible for a real-time system Obviously, in a real time system, the delay line could not possibly contain future data. Considering the center tap of the chanA delay line to represent zero time actually inserts an overall time delay of four samples into the entire process. This is evident in Figure 5 for which you see four samples with a value of zero at the beginning of the black trace. In this case, noise is inserted into the chanA delay line during the first iteration, but it doesn't emerge from the center tap on the delay line until the fifth iteration. Figure 1, on the other hand, shows ten iterations before anything emerges from the center tap on the chanA delay line. This is because Figure 1 shows a noise free case where the only input to the chanA delay line is obtained from the center tap on the rawData delay line. Thus, the rawData delay line causes a six-sample delay in addition to the four-sample delay caused by the chanA delay line. The bottom line The bottom line is that the value that is extracted and saved for plotting in Listing 22 represents our concept of the value at zero time from the input chanA (but in actuality, an overall time delay has been incurred). Apply the convolution filter to chanA Normally, if I were going to apply a fixed convolution filter to a time series, I would write a method that receives the filter and the time series, performs the convolution, and returns the filtered time series. For example, that is what I did in Listing 13 in the earlier lesson entitled Convolution and Frequency Filtering in Java. However, the situation here is different. The convolution filter is not fixed. Rather, I need to interrupt the convolution operation at the end of each iteration and modify the filter coefficients. Listing 23 invokes the dotProduct method to apply the filter coefficients to the current data in the chanA delay line returning a single value of type double as a result.
I will explain the dotProduct method at this point and return to the explanation of the process method shortly. The dotProduct method This method receives two arrays and treats the first n elements in each array as a pair of vectors. It computes and returns the vector dot product of the two vectors. If the length of one array is greater than the length of the other array, it considers the number of dimensions of the vectors to be equal to the length of the smaller array.
As you can see, the vector dot product is simply the sum of the products of the values that describe each of the two vectors. Get the adaptive target Returning now to the explanation of the process method, the code in Listing 25 gets and saves the middle sample from the chanB delay line to be used as the adaptive target.
In other words, the adaptive process will attempt to cause the filtered version of chanA to match the value in the middle of the chanB delay line.
Compute the error value As you will see shortly, the adaptive process makes an adjustment to each of the filter coefficients based on the value of an error.
As you can see in Listing 26, the error value is the difference between the output produced by applying the filter coefficients to the contents of the chanA delay line (see Listing 23) and the value of the target from Listing 25.
Update the filter coefficients We've now arrived at the heart of the matter for an LMS adaptive algorithm. The code in Listing 27 uses a for loop to modify each of the filter coefficient values by computing a correction value for each coefficient and subtracting that correction value from the coefficient.
The correction value The correction value for each coefficient value consists of the product of the following:
As you can see from Figure 1 and Listing 27, as the value of the error approaches zero, the magnitude of the correction values also approaches zero and the coefficient values converge to fixed values. A theoretical justification I'm not going to attempt to justify this adaptive algorithm theoretically. There are hundreds of articles on the web that provide such justification. If you are interested in a justification, I recommend that you use Google to search them out and read them. For example, you might search for the keywords LMS Adaptive Algorithm or for the keywords Steepest Descent. An important concept This is an important concept because I will be using this adaptive algorithm in several future lessons. As mentioned earlier, I wanted to present the algorithm here in a relatively simple case so that when we get to the more difficult cases, you will already have an understanding of the adaptive algorithm and can concentrate on the bigger picture. The end of the adaptive code This is the end of the adaptive process. All of the remaining code in the process method is used to display information about the adaptive process. Feed the time series plotting object Recall that we are still in a for loop processing one input sample during each iteration of the loop. The code in Listing 28 feeds the values from four time series to the plotting object that will eventually be used to plot the time series in Figure 1.
The order of the input parameters to the feedData in Listing 28 corresponds to the black, red, blue, and green traces in Figure 1. Compute and display frequency response data Listing 29 shows the beginning of an if statement that is used to plot the impulse response and the frequency response of the convolution operator at the end of every tenth iteration.
The code in Listing 29 invokes the displayFreqResponse method discussed earlier to produce the amplitude and phase response plots shown in the right column of Figure 2. Extend and display impulse response data The code in Listing 30 deals with the need to extend the impulse response data in order to plot the impulse response at the end of every tenth iteration as shown in the left column of Figure 2.
Extending the impulse response for plotting I discussed the need for extending impulse response for plotting earlier and won't discuss it further here. Suffice it to say that the code in Listing 30 feeds the extended impulse response data to the plotting object so that it can be plotted later. Listing 30 also signals the end of the if statement that is used to cause its body to be executed every tenth iteration. Listing 30 also signals the end of the for loop whose body is executed once for every adaptive iteration. Plot the results Up to this point in time, a lot of data has been fed into the three plotting objects responsible for plotting the results shown in Figure 1 and Figure 2. However, nothing has actually been plotted yet. By invoking the plotData method on each of the three plotting objects, the code in Listing 31 causes them to plot their data at the screen coordinates passed as parameters to the plotData method.
The end of the process method Listing 31 also signals the end of the process method and the end of the Adapt01 class. Some more sample results Before ending this lesson, I want to show you the results for a few more cases. I will begin with a case having no uncorrelated noise on chanA and chanB but with chanB having a time shift of +3 samples relative to chanA. I will run this case using the following command on the command line: java Adapt01 3 0.001 0.0 100 This produces the following output on the command line: timeShift: 3 feedbackGain: 0.0010 noiseLevel: 0.0 numberIterations: 100 Compare this output with the command line output for the first example, particularly with respect to the value of timeShift. The graphic output The time series output for this case is shown in Figure 3.
Figure 3 If you compare Figure 3 with Figure 1, you will see that the target (blue) data in Figure 3 leads the input (black) data in Figure 3, which is just the reverse of Figure 1 where the target data lags the input data. In both cases, the error is decreased to near zero about midway across the top four traces, and the filter output (red) data converges to become a good match for the target (blue) data at that point. The impulse and frequency response The impulse response of the convolution filter is shown in the left column of Figure 4. The frequency response is shown in the right column of Figure 4. Figure 4 Compare Figure 4 with Figure 2 If you compare Figure 4 with Figure 2, you will see that the adaptive process causes the peak in the impulse response to move three samples to the right in Figure 4, as compared to four samples to the left in Figure 2. This is the difference between a leading target for this example and a lagging target for the previous example. In addition, the slope of the phase response in positive in Figure 4 and negative in Figure 2. A positive slope indicates a leading target for this example whereas a negative slope indicates a lagging target for the previous example. Finally, the magnitude of the slope of the phase response in Figure 4 is less than the magnitude of the slope in Figure 2. This is an indication that the relative time shift of three samples in this example is less than the relative time shift of four samples in the previous example. One final example For this case, white random noise was added to both chanA and chanB at a level equal to twenty-five percent of the signal level on each channel. The noise was uncorrelated between chanA and chanB. Some changes to the input parameters were required to achieve reasonably good results. The following command was entered on the command line for this case: java Adapt01 3 0.0005 0.25 210 Compare with the previous example Compare this with the command line input for the previous example and you will see that:
These parameters produced the following output on the command line: timeShift: 3 feedbackGain: 5.0E-4 noiseLevel: 0.25 numberIterations: 210 Note that the time shift of the target relative to chanA is still +3 samples as in the previous example. The time series output The time series output for this case is shown in Figure 5.
Figure 5 Error never goes to zero One thing worth noting is that the error (green) never goes to zero in Figure 5 regardless of the quality of the adaptive solution. The target contains white random noise. There is nothing that the convolution filter being applied to chanA can do to eliminate that noise and it passes straight through from the target to the error in the subtraction process shown in Listing 26. There is also white random noise on chanA. If the convolution filter maintains a flat amplitude response as intended, this noise will be passed through to the output also and will contribute to the noise on the error signal in Figure 5.
Despite the uncorrelated noise on both channels, the filtered output (red) is a reasonably good replica of the target (blue) by the beginning of the traces in the bottom page of Figure 5. The impulse response The impulse response at the beginning and at the end of every tenth adaptive iteration is shown in the left column in the two pages in Figure 6. As you can see, this impulse response is converted from a single impulse in the center of the convolution filter to a single impulse (plus a few very small non-zero values) three samples to the right of center after about 140 iterations. As you know, this is the correct solution, and it was achieved despite the presence of random noise on the two channels.
Figure 6 What this means is that after about 140 iterations, the output from the convolution filter has been shifted so as to register with the target. Most of the error that we see in Figure 5 after this point in time is simply the result of the random white noise being passed through from chanA and the target to the error. The frequency response After about 140 iterations, the amplitude response shown in the right column of Figure 6 is very flat, and the phase response has a positive slope and a linear shape indicating that the convolution filter has converged to the simple time-shift filter that we know to be the correct solution to the problem. Run the ProgramsI encourage you to copy, compile and run the program provided in Listing 32 below. You will need some other classes in addition to the program in Listing 32. I have provided the source code for the class named PlotALot05 in Listing 33. You will need to go to the previous lesson entitled Plotting Large Quantities of Data using Java to get the source code for the other required PlotALot classes. In addition, you will need to go to the lesson entitled Spectrum Analysis using Java, Sampling Frequency, Folding Frequency, and the FFT Algorithm to get the source code for the class named ForwardRealToComplex01. Have fun and learn Above all, have fun and use this program to learn as much as you can about the basics of adaptive filtering. SummaryIn this lesson, I explained adaptive filtering using an LMS adaptive algorithm in a relatively simple scenario. What's Next?The next lesson will tackle a considerably more complicated scenario for adaptive filtering. In the next lesson, I will teach you how to write a whitening filter program for the extraction of wide band signals corrupted by narrow band noise. Following that, the lessons in the series will become somewhat more general. I plan to publish lessons that explain and provide examples for the four common scenarios in which adaptive filtering is used:
Somewhere along the way I may publish a lesson that explains and illustrates the difference between least mean square (LMS) and recursive least squares (RLS) adaptive algorithms. Complete Program ListingsComplete listings of two of the programs discussed in this lesson are provided in Listing 32 and Listing 33 below.
Copyright 2005, 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 authorRichard 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. |