Java Programming Notes # 1510
 Preface
 Preview
 General
Background Information  Practical Examples
 Run the Program
 Summary
 What’s Next?
 References
 Complete Program
Listing
Preface
This is the first installment of a multipart lesson on digital recursive
filtering.
A Recursive Filtering Workbench
The complete collection of installments presents and explains the code for an interactive
Recursive Filtering Workbench (See Figure 1) that can be used to design, experiment with, and
evaluate the behavior of digital recursive filters.
(The digital Recursive Filtering Workbench will be referred to
hereafter simply as the workbench.)
By the end of the last installment, you should have learned how to write a
Java program to create such a workbench.
Hopefully, you will also have gained some understanding of what recursive
filters are, how they behave, and how they fit into the larger overall
technology of Digital Signal Processing (DSP).
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 listings and figures 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.
I also recommend that you pay particular attention to the lessons listed in
the References section of this document.
Preview
Complex poles and zeros
The workbench makes it possible for the user to design a recursive filter by
specifying the locations of
sixteen poles and sixteen
zeros in the complex zplane and then to evaluate the behavior of the recursive filter
for that set of poles and zeros. The user can relocate the poles and zeros
and reevaluate the behavior of the corresponding recursive filter as many times
as may be beneficial without a requirement to restart the program.
An image of the zplane
The GUI in the top center portion of the screen in Figure 1 is an image of the
complex zplane showing the unit circle along with the locations of all the poles and
zeros.
(Because of the reduction in the size of the image, the poles and
zeros are only barely visible in Figure 1. A
fullsize version is shown in Figure 2.)
The user can interactively change the location of any pair of complex poles or zeros
by selecting a specific pair of poles or zeros and clicking the new location
with the mouse in the zplane GUI. This provides a quick and easy way to position
the poles and zeros in the zplane.
Numeric text input
The GUI that is partially showing in the upper right of Figure 1 contains
four text fields and a radio
button for each pair of complex poles and each pair of complex zeros.
(A fullsize version of this GUI is shown in Figure 5.)
This
GUI has several purposes. For example, the user can specify the locations
of pairs of poles or zeros very accurately by entering the real
and imaginary coordinate values corresponding to the locations into the text
fields in this GUI. The user can also view the length and the angle
(relative to the real axis) of an imaginary vector that connects each pole
and each zero to the origin.
This GUI also provides some other features that are used to control the
behavior of the workbench program. These features will be explained later.
The two GUIs are connected
When this GUI is used to relocate a pair of poles or zeros, the graphical
image of the zplane in Figure 2 is automatically updated to
reflect the new location. Similarly, when the mouse is used with the
graphical image of the zplane to relocate a pair of poles or zeros, the numeric
information in the GUI in Figure 5 is automatically updated to
reflect the new location.
Two ways to specify the location of poles and zeros
Thus, these two GUIs provide two different ways that the user can specify the
locations of poles and zeros. Clicking the location with the mouse is very
quick and easy but not particularly accurate. Entering the numeric coordinate values
into the text fields takes a little more effort, but is very accurate.
The two approaches can be used in combination to create a rough design with
the mouse and then to polish the design by entering accurate pole and zero
locations into the text fields.
The behavior of the recursive filter
Once the locations of the poles and zeros that define a recursive filter have
been established using either or both of the two GUIs discussed above, the
leftmost GUI in Figure 1 can be used to display the
following information about the recursive filter. (See
Figure 6 for a fullsize version of this GUI.)
 The impulse response in the time domain.
 The amplitude response in the frequency domain.
 The phase response in the frequency domain.
Two computational approaches
Two versions of the amplitude response and two versions of the phase response
are plotted in the GUI shown in Figure 6. One version is based on the
Fourier Transform of the impulse response. The other version is based on a
vector analysis of the locations of the poles and zeros in the zplane.
Adjusting the plotting parameters
The GUI shown in Figure 6 provides for the input of
seven different plotting parameters (vertical scale, horizontal scale, tic
mark locations, etc.). The user can modify the plotting parameters and replot
the graphs as many times as may be needed to get a good visual representation of
the behavior of the recursive filter.
General Background Information
I will not attempt to teach you the theory behind recursive filtering. Rather, I will assume that you already understand that theory, or
that you are studying a good theoretical resource on recursive filtering
concurrently with the study of this lesson.
(A very good resource on the
theory of digital filtering in general, and recursive digital filtering in
particular, is the Introduction to Digital Filters with Audio Applications by
Julius O. Smith III .)
As I explain the material in this lesson, I will make technical
statements of fact without providing a theoretical justification for those
statements. The theoretical justification for those statements can be
found in Smith’s
material as well as other resources on the web.
What is a recursive filter?
Here is some of what
Wikipedia
has to say about recursive filters:
Recursive filters are
signal processing filters which reuse one or more output(s) of the
filter as inputs. This feedback results in an unending
impulse response characterized by either
exponentially growing, decaying, or sinusoidal signal output components.
A simple recursive filter
Although it is possible to design very complex recursive filter structures,
this lesson will deal with simple recursive filters. As I define it,
a simple recursive filter contains a single feedforward section, which is essentially
equivalent to the convolution filters described in
earlier
lessons. It also contains a single feedback loop, which performs a
convolution operation on the previous output samples before combining that
result with the new input data.
The recursive filter used in this workbench implements a difference equation
very similar to that provided by
Smith.
Most digital filters produce one output sample for every input sample.
The computational cost of a digital filtering operation can usually be measured
in terms of the number of multiplyadd operations (MADs) required for
each output sample.
It is not possible to design a recursive filter to achieve every possible
filtering behavior that may be needed. Given all the possible needs for
digital filters, it is more likely that you can
design a nonrecursive convolution filter to satisfy a given need than that you can
design a recursive filter to satisfy that same need. However, if a recursive
filter can be designed to satisfy a given need, the implementation of the
recursive filter is likely to have a lower computational cost than would be the
case with a nonrecursive convolution filter designed to satisfy the same need.
Greater arithmetic accuracy is required
However, in the computing world as well as the world outside the computer, there is no such thing as a free
lunch. While a recursive filter designed to satisfy a given
filtering need may have a lower computational cost than a nonrecursive convolution
filter design to satisfy the same need, the recursive filter will probably
require greater arithmetic accuracy than the convolution filter. This is
because a recursive filter feeds part of its output back into the input.
As a result, small arithmetic errors can build up to produce large errors.
Phase distortion can be a problem
All filters produce some amount of phase distortion. This means that
some frequency components in the output are delayed more or less than other
frequency components relative to the relationship between those components in
the input. This results in a modification of the output waveform relative
to the input waveform. This may, or may not be a problem, depending on the
application.
A time delay
A nonrecursive convolution filter can be designed in such a way that the phase
distortion manifests itself simply as a time delay between the input and the
output with the waveform being otherwise unmodified by the phase response.
This is the result of designing the filter in such a way that the phase response
is a linear function of frequency. Sometimes this is a very acceptable
form of phase distortion, particularly if all the signals being processed are
subjected to the same amount of time delay.
To my knowledge, it is not possible to design a general recursive filter where the
phase response is always a linear function of frequency. (It is possible, or
almost so for some unique cases of recursive filters.) This may be another
disadvantage of recursive filters relative to nonrecursive convolution filters
if nonlinear phase distortion is a problem.
Advantages and disadvantages of recursive filters
Thus, the primary advantage of a recursive filter relative to a nonrecursive convolution
filter is that it is likely to have a lower computational cost. The
potential disadvantages are the requirement for greater arithmetic accuracy and
a nonlinear phase response.
Design of a nonrecursive convolution filter
If you are not constrained as to the number of filter coefficients required, nonrecursive convolution filters
are relatively easy to design. Although there are
several design approaches available, one of the most straightforward approaches
involves developing the timedomain impulse response of the filter by performing
an inverse Fourier Transform on the desired complex frequency response.
However, if you are constrained as to the number of filter coefficients
required, this often involves making a compromise between convolution filter
length and expectations in the frequency domain. The first attempt may not
produce an acceptable frequency response and an acceptable convolution filter
length. If not, the process can be repeated with
modifications to the frequencydomain expectations until an impulse response having an acceptable
length that matches an acceptable frequency response is achieved. This
process is discussed on one of the
pages of
London’s Global University.
Perform an inverse Fourier Transform on complex
frequency response
With this approach, you specify the desired amplitude and phase response of
the filter as a set of complex samples in the frequency domain. Then you
perform an inverse Fourier Transform on that set of complex samples.
(In this case, it is probably better to use an ordinary inverse
Fourier Transform algorithm instead of an inverse
FFT algorithm due to the
improved flexibility of the nonFFT algorithm relative to the frequency
sampling interval, the timedomain sampling interval, and the length of the
timedomain impulse response.)
If the complex samples are correctly specified, the inverse Fourier Transform
will produce a purelyreal impulse response in the time domain.
Truncate or window the impulse response
The next step is to truncate (or possibly
window) the resulting impulse response to the desired
number of convolution filter coefficients. Then perform a forward Fourier Transform on the truncated (or windowed) impulse response to see how well the frequency
response of the truncated (or windowed) convolution filter matches the desired frequency response.
If the match is good, your design task is finished. Use the truncated
(or windowed) impulse response as your convolution filter.
If the match is not good …
If the match is not good, you must modify your expectations and desires in
order to
make them more realistic. At this point, you can either
 Extend the length of the truncated (or windowed) impulse response until you get a good
match, or  Modify the desired frequency response and start the design process over
in an attempt to get a truncated (or windowed) impulse response of a reasonable length
that matches an acceptable frequency response, or  A combination of the two.
Typically, for the second case, the required modifications will involve
reducing the slope of changes in the amplitude response. Abrupt changes
in the required amplitude response typically require long impulse responses to achieve.
Designing recursive filters
The
design of
recursive filters is considerably more complex than the design of
nonrecursive convolution filters, even for the common cases of
lowpass, bandpass, or highpass filters. (The design of recursive
filters for other more exotic frequencyresponse characteristics is even more
complex.)
There are mathematical techniques
available for
designing recursive filters (many of which derive largely from the design of analog
electronic filters). However, this workbench doesn’t use mathematical
design techniques. Rather, it provides a visual
interactive design approach instead of a mathematical approach.
Not intended for
recursive filter design
Although this workbench provides the capability for designing recursive
filters (even those having exotic frequency response characteristics),
the workbench is intended to serve as a learning tool rather than a design tool.
My hope is that through experimentation with the workbench and examination of
the code, you will gain a
better understanding of the behavior of recursive filters, and an understanding
of how to program them. However, if you also
want to use the workbench as a design tool, be my guest.
(If you do decide to use the workbench for serious recursive filter
design, you might want to expand the size of the zplane GUI shown in
Figure 2. The radius of
the unit circle is currently 200 pixels.
Therefore, you cannot possibly use the
mouse to specify the location of a pole or a zero to an accuracy greater
than about one part in 200. You might want to consider changing the radius
to 500 pixels (or more) and displaying only the top half of the zplane.
This will improve the accuracy with which you can use the mouse to specify
the location of a pole or a zero to about one part in 500.)
A ratio of polynomials in the zplane
The transfer function H(z) for a digital recursive filter can be
thought of as a ratio of two polynomials:
H(z) = B(z)/A(z)
The zplane
The data length
The text field in the top right in Figure 5 shows that this polezero
configuration will be analyzed using an impulse response having a length of 1024
samples, and that the Fourier Transform of the impulse response will be computed
at 1024 points between a frequency of zero and the sampling frequency. We
will see the results of that analysis later.
Numeric pole and zero coordinates
The grid of radio buttons and text fields bordered by green in Figure 5
contains information about eight pairs of complexconjugate zeros that are
created and analyzed by the workbench. The grid of radio button and text
fields bordered by yellow contain information about eight pairs of complex
conjugate poles that are created and analyzed by the workbench.
(A recursive filter can also be described by individual poles and
zeros on the real axis. However, to reduce the complexity of the
program, this workbench
always deals with pole and zeros in complex conjugate pairs. If a pair
is placed on the real axis, as in Figure 2, the two
items in the pair are stacked on top of one
another. You may want to modify the program to relax this requirement.)
A complex conjugate pair
Each row of four text fields and a radio button in Figure 5 corresponds to a complex conjugate
pair of poles or zeros. As you can see, one text field contains the real
coordinate value for the pole or zero. A second text field contains the
imaginary coordinate value. The third text field contains the angle in
degrees relative to the real axis described by a vector drawn from the origin to
one of the poles or zeros in the pair. The fourth text field contains the
length of that vector.
Disabled text fields
The text fields containing angle and length information are disabled insofar
as user input is concerned. Thus, they serve to display the angle and
length values computed on the basis of the contents of the real and imaginary
fields.
The radio buttons
Only one of the radio buttons shown in Figure 5 can be selected at any point
in time. When
one of the radio buttons is selected and the mouse is clicked in the zplane
shown in Figure 2, the selected pole or zero is moved to the location of the
mouse click in the zplane.
The coordinate, angle, and length information in the
corresponding text fields in Figure 5 are automatically updated to reflect the new location for
the pole or zero.
Conversely, if a new numeric value is entered into one
of the text fields in Figure 5 to cause the location of the pole or zero to
change, the image of the corresponding pole or zero is
automatically updated to show the new location in the zplane in Figure 2.
The behavior of the recursive filter
There are five graphs in Figure 6. The top graph shows the first 512 samples of the impulse response
in the time domain of a recursive filter having the polezero configuration
shown in Figure 2.
Figure 6 
The amplitude response
The second and third graphs in Figure 6 show the amplitude response of that
recursive filter (computed two different ways) between a frequency of
zero and the Nyquist folding frequency. The amplitude response is plotted
in decibels.
The decibel scale
Quite a bit of scaling and normalization takes place in the workbench with
regard to the amplitude response graphs. For example, the amplitude
response graphs are normalized so that the peak in the amplitude response always
plots as a value of 100 units.
When the amplitude response graphs are calibrated in such a way as to take
the scaling and normalization into account, the vertical scale comes out to be
approximately 0.2 decibels per plotting unit. Thus, the amplitude response
graphs shown in Figure 6 cover a range of 270 units or 54
decibels.
Each factor of two change in the amplitude response represents a change of
approximately six decibels. Therefore, a range of 54 decibels represents a
dynamic range in the amplitude response of approximately a factor of 512.
The phase response
The fourth and fifth graphs in Figure 6 show the phase response of the
recursive filter (also computed two different ways) over the same
frequency range as the amplituderesponse graphs. The phase response is plotted in degrees
The effect of the pair of poles
The peak in the amplitude at a frequency of zero in Figure 6 is caused by the
pair of poles that reside on the real axis near the intersection of the real
axis and the unit circle in Figure 2.
The effect of the pair of zeros
The notch in the amplitude response and the corresponding rapid
change in the phase response in Figure 6 are caused by the pair of
complex conjugate zeros that
are very close to the right side of the unit circle in Figure 2.
The top half of the unit circle in Figure 2 maps into the frequency axes in
the four bottom graphs in Figure 6.
The rightmost intersection of the
unit circle and the real axis corresponds to a frequency of zero in Figure 6.
The leftmost intersection of the unit circle and the real axis corresponds to
the Nyquist folding frequency (onehalf the sampling frequency) at the
extreme right end of Figure 6.
The topmost intersection of the unit
circle and the imaginary axis in Figure 2 corresponds to onefourth the sampling
frequency (the center of the horizontal frequency axes in Figure 6).
The bottom half
The bottom half of the unit circle in Figure 2 maps into the frequency range
from the Nyquist folding frequency to the sampling frequency. This
corresponds to the frequency range from the center of the horizontal frequency
axes to the rightmost end of the frequency axes in Figure 7.
Figure 7 
(The horizontal scale for Figure 7 was modified relative to Figure 6 so as to
show the entire amplitude and phase response from zero to the sampling
frequency. The Nyquist folding frequency is at the
right end of the frequency axes in Figure 6, and is at
the center of the frequency
axes in Figure 7.)
Represents the complete unit circle
Thus, the frequency axes in Figure 7 represent all of the points on the entire unit
circle in Figure 2.
Note that the amplitude response is symmetric
about the Nyquist folding frequency at the center
of Figure 7, and the phase response is
asymmetric about the Nyquist folding frequency.
Therefore, no new information is imparted by observing the plotted results to the right of the folding frequency relative to the results to the left of the
folding frequency. As a result, I will often ignore the results to the
right of the folding frequency in this lesson.
A point represents a frequency
Now let’s turn our attention back to Figure 6. A specific point somewhere on the top half of the unit circle in Figure 2
represents a specific frequency somewhere along a frequency axis in Figure 6.
The value of the amplitude response at that frequency depends entirely on the
distances from that point on the unit circle to the locations of all the poles
and all the zeros in the zplane.
A nearby zero
For example, a zero on or near the point on the unit circle will cause a low value in the amplitude
response at that frequency unless offset by poles located nearby. Thus,
the zero near the unit circle in Figure 2 causes a notch in the amplitude response
shown at that frequency in Figure 6.
A nearby pole
A pole near the point on the unit circle will cause a high value in the amplitude response
at that frequency unless offset
by zeros located nearby. Thus, the pair of poles located on the real axis
near the intersection of the real axis and the unit circle in Figure 2 causes
the peak at zero frequency in the amplitude response shown in Figure 6.
An unstable recursive filter
Before I forget it, let me mention that the existence of any pole outside the
unit circle in the zplane will cause the recursive filter to be
unstable.
This will cause the output values produced by the recursive filter to increase
without bounds until the arithmetic processor overflows in some manner.
Computing the amplitude response
The amplitude response of the recursive filter at a
given frequency (a point on the unit circle) can be determined by
computing the product of the distances from that point on the unit circle to all
of the zeros and dividing that value by the product of the distances from the
same point on the unit circle to all of the poles.
The amplitude response at all frequencies between zero and the
Nyquist
folding frequency can be determined by performing this computation at all points
on the top half of the unit circle.
(Because there are an infinite number of points on the top half of the
unit circle, we must choose the set of points at which to perform the
computation.)
The amplitude response shown in the third graph in Figure 6 was computed in
this manner. (A little later, I will explain the method that was used to
compute the amplitude response shown in the second graph in Figure 6.)
Poles and zeros at the origin
When poles and zeros are located at the origin, the distance to those
poles and zeros from any point on the unit circle is 1.0. Therefore,
multiplying or dividing by that distance is of no consequence, and poles and
zeros located at the origin have no effect on the amplitude response of the
recursive filter.
Computing the phase angle
The phase angle of the recursive filter at a particular frequency
(represented by a point on the unit circle) can be determined by subtracting
the sum of the angles from the point on the unit circle to all of the poles from
the sum of the angles from the same point on the unit circle to all of the
zeros.
The phase response in the bottom graph of Figure 6 was computed using this
method.
Getting the recursive filter coefficients
The primary objective of a recursive filter design effort is not to identify the pole
and zero locations necessary to produce the desired recursive filtering
behavior. Rather, it is to obtain the feedforward
and feedback coefficients to use in the application of the recursive filtering algorithm
to produce the desired behavior.
So far, we have only discussed the locations of the poles and zeros for the
denominator and numerator polynomials. The next thing that we need to do
is to use those pole and zero locations to determine the values of the
coefficients for the numerator and denominator polynomials.
This is the easy part
The zeros represent the roots of
the numerator polynomial and the poles represent the roots of the denominator
polynomial. All that is necessary to obtain the polynomial coefficients is to multiply the roots.
For example, the real and
imaginary values for the first zero in Figure 5 would require us to multiply the following roots:
(z0.95+0.31j)*(z0.950.31j)
The multiplication of these roots produces a secondorder polynomial
with the FeedForward (Numerator) Coefficients shown in Figure 3.
Similarly, the denominator coefficients are obtained by using the pole values
to construct the roots of the denominator polynomial and to multiply those
roots to produce the coefficients for the denominator polynomial.
Must implement the recursive filtering algorithm
Having obtained the feedforward and feedback coefficients, the next step is
to install those coefficients in a recursive filtering algorithm and to use it
to filter the sampled data of interest. Once we do that, we need a way to confirm
that the coefficients are being used correctly and that the recursive filtering
algorithm is behaving as it should. One way to obtain this confirmation is to:
 Apply the recursive filter to a time series consisting of a single
nonzero value followed by sample values of zero thereafter. (The
output of such a filtering operation is commonly referred to as the
impulse
response of the filter.)  Perform a Fourier Transform on the
impulse response.  Compare the amplitude response and the phase response produced by the
Fourier Transform of the
impulse
response to that produced using the poles
and zeros as described above.
Do they match?
If the amplitude and phase responses produced using the two different approaches match,
this is a very good indication that the filter coefficients are being correctly
applied in the recursive filtering process. If they don’t match, the
recursive filtering process should be subject to question.
The results for this example
The impulse
response is shown in the first graph in Figure 6. The
amplitude response shown in the second graph in Figure 6 was produced by
performing a Fourier Transform on the
impulse
response.
The phase response shown in the fourth graph in Figure 6 was produced by
performing the same Fourier Transform on the
impulse
response.
The amplitude and phases responses produced by applying the Fourier Transform
to the impulse
response match those produced using the pole and zero locations
in the third and fifth graphs in Figure 6. This is a strong indication
that the recursive filtering process in the workbench is implemented correctly.
Practical Examples
Now that you understand how to control the locations of the poles and zeros,
and understand how those locations impact the behavior of the recursive filter,
it’s time to take a look at a few practical examples.
The complete source code for the workbench is provided for your benefit in
Listing 1 near the end of the lesson. If you wish to do so, you can copy, compile, and execute that source
code and follow along for the practical examples that follow.
Summary of the behavior of the workbench program
Before getting into the details of the examples, I will summarize the
behavior of the workbench program.
This program provides a visual, interactive recursive filtering workbench.
The purpose of the program is to make it easy to experiment with the behavior of
recursive filters and to visualize the results of those experiments.
Complex poles and zeros
The program implements a recursive filter having eight pairs of complex
conjugate poles and eight pairs of complex conjugate zeros. The locations
of the pairs of poles and zeros in the zplane are controlled by the user.
Although the pairs of poles and zeros can be colocated on the real axis, the
program does not support the placement of individual poles and zeros on the real
axis.
The user can reduce the number of poles and zeros that impact the behavior of
the recursive filter by moving excess poles and zeros to the origin in the
zplane. This renders them ineffective in the behavior of the recursive
filter.
The impulse response
The program provides three interactive displays on the screen. The
first display (see Figure 6) contains five graphs
in a vertical stack. The first graph in this display shows the
impulse
response of the recursive filter in the time domain.
The amplitude response
The second and third
graphs show the amplitude response of the recursive filter in the frequency
domain computed using two different computational approaches. The two
different computational approaches are provided for comparison purposes.
The first computational approach for computing the amplitude response is to
perform a Fourier Transform on the
impulse
response using an
FFT algorithm.
The quality of the estimate of the amplitude response using this approach is
dependent on the extent to which the entire impulse response is captured in the
set of samples used to perform the
FFT. If the impulse response is
truncated, the estimate will be degraded.
The second approach for computing the amplitude response involves computing
the product of the vector lengths from each point on the unit circle to each of
the poles and each of the zeros in the zplane. This approach provides an
idealized estimate of the amplitude response of the recursive filter, unaffected
by impulseresponse length considerations. This approach produces the same
results that should be produced by performing the
FFT on a set of
impulseresponse samples of sufficient length to guarantee that the values in
the impulse response have damped down to zero (the impulse response is
totally captured in the set of samples on which the
FFT is performed).
The phase response
The fourth and fifth graphs in Figure 6 show the phase response of
the recursive filter computed using similar approaches to those described above
for the amplitude response. (The second approach uses the sum of vector
angles instead of the product of vector lengths.) Once again, the two
approaches are provided for comparison purposes.
Computational data length
By default, the program computes and captures the impulse response for a
length of 1024 samples and performs the Fourier Transform on that length.
However, the length of the captured impulse response and the corresponding
FFT
can be changed by the user to any length between 2 samples and 16384 samples,
provided that the length is an even power of two. (If the length
specified by the user is not an even power of two, it is automatically changed
to an even power of two by the program.)
Interactive plotting parameters
The display shown in Figure 6 is interactive in the sense that there are seven different
plotting parameters (such as vertical scale, horizontal scale, location of
tic marks, etc.) that can be adjusted by the user to produce plots that are
visually useful. The user can modify any of the parameters and then click
the Graph button to cause the graphs to be replotted using the new
parameters.
Poles and zeros in the zplane
The second display (see Figure 2) shows the
locations of all of the poles and all the zeros in the zplane. The user can use
the mouse to change the location of any pair of complex conjugate poles or zeros
by first selecting a specific pair of poles or zeros and then clicking the new
location in the zplane.
This interactive capability makes it possible for the user to modify the
design of the recursive filter in a completely graphic manner by positioning the
poles and zeros in the zplane with the mouse.
Having relocated one or more pairs of poles or zeros in the zplane, the user
can then click the Graph button in Figure 6 to cause
the impulse response, the amplitude response, and the phase
response of the new recursive filter with the modified pole and zero locations
to be computed and displayed.
A general user control panel
The third display (see Figure 5) is a control
panel that uses text fields, ordinary buttons, and radio buttons to allow the
user to perform the following tasks:
 Specify a new length for the impulse response as an even power of two.
(Once again, if the user fails to specify an even power of two, the value
provided by the user is automatically converted to an even power of two by the program.)  Cause all of the poles to be moved to the origin in the zplane.
 Cause all of the zeros to be moved to the origin in the zplane.
 Select a pair of complex conjugate poles or zeros to be relocated using
the mouse in the display of the zplane.  View the angle (relative to the horizontal axis in the zplane) described by an imaginary vector that connects each pole
and each zero to the origin.  View the length of an imaginary vector that connects each pole and each zero
to the origin in the zplane.  Enter (into a text field) a new real or imaginary value
specifying the location of a pair of complex conjugate poles or zeros in the
zplane.
Entering a new location value
When a new real or imaginary value is entered into a text field, the angle
and the length are automatically updated to reflect the new location of the pole
or zero and the display of the pole or zero in the zplane is also updated to
show the new location.
Clicking a new location in the zplane
Conversely, when the mouse is used to relocate a pole or zero in the zplane, the corresponding real and imaginary values in the text fields and the
corresponding angle and length in Figure 5 are automatically updated to reflect the new
location for the pole or zero.
Entry of invalid data
Note that mainly for convenience (but for technical reasons as well),
the angle and length in the text fields and the location of the pole or zero in
the Zplane display are updated as the user enters each character into the text
field.
If the user enters text into the text field that cannot be
converted into a numeric value of type double, (such as the pair of
characters "." for example) the contents of the text field are
automatically converted to a single "0" character. A warning message is
displayed on the commandline screen when this happens.
There is no
longterm harm when this occurs, but the user may need to start over to enter
the new value. Thus, the user should exercise some care regarding the
order in which the characters in the text field are modified when entering new
real and imaginary values.
Coefficients displayed on the commandline screen
Each time the impulse response and the spectral data are plotted, the
seventeen feedforward filter coefficients and the sixteen feedback coefficients
used by the recursive filter to produce the output being plotted are displayed
on the commandline screen.
Usage
This program requires access to the following source code or class files,
which were published in earlier tutorials (see the
References section of this document):
 GraphIntfc01
 ForwardRealToComplexFFT01
 Graph03
Having compiled all the source code, enter the following command at the command prompt to run this program:
java Graph03 Dsp046
The program was tested using J2SE 5.0 under WinXP. J2SE 5.0 or later is required due to the
use of static imports and the printf method.
The behavior at startup
The workbench establishes a set of default pole and zero locations on
startup. Among other reasons, this is done to confirm that the workbench
program is operating properly.
When you first start the program and make a few adjustments to
the plotting parameters, the three GUIs that appear on your screen should match
those shown in Figure 8, Figure 9, and Figure 10.
(To adjust the plotting parameters, modify the values in the text
fields at the bottom and click the Graph button Figure 10.)
Sixteen poles and sixteen zeros
The poles and zeros in Figure 8 are located at a distance of 0.995 from the
origin. This is very close to, but not on the unit circle.
Located at incremental angles
The amplitude and phase response
The remaining four graphs in Figure 10 show the amplitude response and the phase response
computed using both methods discussed earlier and displayed from zero frequency
to the sampling frequency.
If you examine the amplitude response carefully, you will see sixteen peaks
corresponding to the locations of the sixteen poles shown in Figure 8. You
will also see sixteen notches corresponding to the sixteen zeros shown in Figure 8.
Similarly, the big changes in the phase response occur at the frequency locations of
the poles and the zeros.
Insufficient data length
Note also that the xMax plotting parameter in Figure 12 was reduced
from 1024 to 128 so that only 128 samples are displayed horizontally in the
graphs.
A truncated impulse response
As I mentioned earlier, about 700 samples are required for the impulse response to damp down to near zero. As you can see in the top graph of
Figure 12, that portion of the impulse response captured and analyzed for this
example was only a small portion of the total impulse response.
Figure 12 
A poor match to the actual amplitude and phase
response
The Fourier Transform of the truncated impulse response provided a poor
estimate of the behavior of the recursive filter.
The estimate of the
amplitude response produced by the Fourier Transform shown in the second
graph in Figure 12 is a poor match for the amplitude response shown in the
third graph.
Similarly, the estimate of the phase response produced by the Fourier Transform shown in the fourth graph is a poor match for the phase
response shown in the fifth graph.
Not a problem with the recursive filter
This doesn’t mean that there was something wrong with the recursive filter.
It didn’t change. Rather, the change was made in one of the computational
approaches used
to estimate the amplitude and phase response of the recursive filter. The amount of data captured
and subjected to Fourier analysis was simply insufficient to provide good
results.
Don’t truncate the impulse response
The bottom line is that when using Fourier analysis to confirm proper
operation of a recursive filter, you must be careful to capture the entire
impulse response instead of truncating it and performing Fourier analysis on a
truncated version of the impulse response.
Moving poles and zeros to the origin
The coordinate values
The real and imaginary coordinate values for this pole are shown in
Figure
14.
Figure 14 
A very long impulse response
(After finding the approximate locations of the poles and zeros using the
mouse, I used a calculator to polish the real and imaginary coordinate
values that specify the locations.)
How does this work?
The zero on the unit circle in Figure 18 produces the notch at the desired
frequency. Unfortunately, this zero also tends to drag down the response
on both sides of the notch.
The pair of poles in Figure 18 helps to bring the response on each side of
the notch back up to achieve a narrower notch and a flatter overall response.
The zero between the poles and the origin also tends to cause the overall
response to be flatter outside the vicinity of the notch. This is because
at points on the unit circle some distance from the notch, there are the same
number of poles and zeros all at approximately the same distance from the point
on the unit circle. As a result, the product of the distances to the poles
is approximately the same as the product of the distances to the zeros.
When the product of zero distances is divided by the product of pole distances,
the quotient is approximately unity.
Locations of the poles and zeros
The actual real and imaginary coordinate values for the poles and the zeros
in this filter are shown in Figure 19. I used the angle and length information shown
in Figure 19 in conjunction with the mouse to locate the poles and the zeros.
The angle and length information helped me to create a pattern of poles and zeros that were
nearly symmetrical with respect to a line at an angle of 60 degrees relative to
the real axis. As mentioned earlier, I then went back with a calculator
and polished the coordinate values producing the values shown in
Figure 19.
Figure 19 
The filter coefficients
The feedforward and feedback filter coefficients that resulted
from this placement of poles and zeros are shown in Figure 20. This
recursive filter has a total of nine nonzero coefficients. Therefore,
this filter could be realized at a computational cost of nine
MADs per output
sample.
FeedForward (Numerator) Coefficients 0 1.000 1 1.739 2 2.285 3 1.285 4 0.546 Feedback (Denominator) Coefficients 1 1.582 2 1.874 3 0.995 4 0.395 Figure 20 
A narrower notch filter
Now assume that you need the notch to be even narrower than that shown in
Figure 17 and that you are willing to sacrifice some flatness in the overall
amplitude response to achieve the narrower notch as shown in
Figure 21.
Figure 21 
New locations for the poles
This can be accomplished by leaving the zeros alone and pushing the pair of
complex conjugate poles in Figure 18 closer to the unit circle and closer to an
imaginary line that connects the zero on the unit circle to the origin.
This is shown in the zplane in Figure 22.
Figure 22 
The poles in Figure 22 were pushed so close to the imaginary line mentioned
above that they are actually stacked on top of one another on that imaginary
line.
Coordinates of the poles and zeros
The numeric coordinates for the poles and zeros in Figure 22 are shown in
Figure 23.
Figure 23 
The filter coefficients
Finally, the feedforward and feedback coefficient values for this filter are
shown in Figure 24.
FeedForward (Numerator) Coefficients 0 1.000 1 1.739 2 2.285 3 1.285 4 0.546 Feedback (Denominator) Coefficients 1 1.840 2 2.539 3 1.557 4 0.716 Figure 24 
As before, this filter has a total of nine nonzero coefficients.
Therefore, the filter could be realized at a computational cost of nine
MADs per output sample.
If you compare Figure 21 with Figure 17, you will see that the notch for this
filter is considerably narrower than the notch for the filter shown in Figure 17.
When the poles move to the origin …
The next thing that I will show you is what happens when all of the poles
are moved to the origin. To do this demonstration, I will restart the
workbench program causing the poles and zeros to be placed in their default
locations as shown in Figure 8. Then I will click the button labeled
Move Poles to Origin to cause all of the poles to be moved to the origin.
This leaves the zeros positioned as shown graphically in
Figure 25.
Figure 25 
The numeric pole and zero locations
The numeric values for the pole and zero locations are shown in
Figure 26.
Figure 26 
The filter coefficients
As you learned earlier and can confirm from Figure 27, when all of the poles
are moved to the origin, all of the feedback coefficients go to zero. This
means that none of the output is fed back into the input.
FeedForward (Numerator) Coefficients 0 1.000 1 1.960 2 2.851 3 3.646 4 4.324 5 4.864 6 5.251 7 5.476 8 5.532 9 5.421 10 5.147 11 4.720 12 4.154 13 3.468 14 2.684 15 1.827 16 0.923 Feedback (Denominator) Coefficients 1 0.000 2 0.000 3 0.000 4 0.000 5 0.000 6 0.000 7 0.000 8 0.000 9 0.000 10 0.000 11 0.000 12 0.000 13 0.000 14 0.000 15 0.000 16 0.000 Figure 27 
A nonrecursive convolution filter
The result of eliminating the feedback in the recursive filter algorithm is
to cause that algorithm to behave like a nonrecursive convolution filter.
The impulse response of a nonrecursive convolution is the same length as the
number of coefficients in the convolution filter. Once the impulse moves through all
of the filter coefficients, the output goes to zero where it remains from that
time forward.
The impulse response of the filter
The top graph in Figure 28 shows that is what happens in this case.
(Simply ignore the other graphs in Figure 28 at this point.)
Figure 28 
As you can see from Figure 27, there are 17 coefficients in the feedforward
loop. The horizontal scale in Figure 28 was expanded to show the
details of the filter output. If you count the number of nonzero values
in the output shown in the top graph of Figure 28, you will see that it matches
the number of nonzero values in Figure 27. If you estimate the amplitudes
of the nonzero values in Figure 28, you will see that they match the values of
the feedforward coefficients in Figure 27.
The amplitude and phase response
If you count the number of notches in the amplitude response in Figure 29,
you will see that the number matches the number of zeros shown in the zplane in
Figure 25. In addition, the locations of the notches in Figure 29 match
the locations of the zeros near the unit circle in the zplane shown in Figure 25.
Add a pole
Now I will show you what happens to the impulse response when a single pair of
complex conjugate poles is moved away from the origin, causing some of the
output to be fed back into the input. To do this, I will leave the zeros
in the same locations as in Figure 25 and revert back to the same horizontal
plotting scale factor that I used in Figure 28.
I moved one pair of complex conjugate poles from the origin to a location
near the unit circle at about twenty degrees off the real axis as shown in
Figure 30.
Figure 30 
The numeric location of the pole
The actual location of the pole in numeric terms is shown for the first pole
in the bottom chart in Figure 31.
Figure 31 
The feedback coefficients
Because the locations of the zeros were not modified, the values of the
feedforward coefficients did not change. However, there are now two
nonzero feedback coefficients as shown in Figure 32.
Feedback (Denominator) Coefficients 1 1.790 2 0.910 Figure 32 
The new impulse response
One of the results of moving this pole away from the origin in the zplane is
shown by the top graph in Figure 33.
Figure 33 
As you can see, the filter output no longer goes to zero after seventeen
samples. Rather, the value that is fed back into the input as a result of
the nonzero feedback coefficients continues to circulate in the filter causing
the output to resemble a damped lowfrequency sinusoid.
(Although I didn’t
show it here, by compressing the horizontal scale and making the vertical scale
more sensitive for the first graph shown in Figure 33, it can be shown that at
least 750 output samples are required to cause the
impulse response for this
filter to go to near zero.)
The amplitude and phase response
As you can see by comparing the amplitude response graphs, moving the pair of complex conjugate poles from the origin to
near the unit circle produces a peak at a low frequency in the amplitude
response that isn’t there in Figure 29.
A pole very close to the unit circle
I told you earlier that the recursive filter is unstable if any of the poles
are outside the unit circle. What happens if a pole is on the unit circle?
The answer is that the recursive filter turns into an oscillator. I will
illustrate this with several graphs.
Let’s begin with the polezero configuration shown in
Figure 35.
Figure 35 
As you can see, all of the zeros in Figure 35 are located at the origin. A single
pair of complex conjugate poles is placed very near the unit circle at an angle of 45
degrees relative to the real axis. The distance from the pole to the origin is 0.9998489885977782,
so you can see that the pole is very close to the unit circle.
(The real and imaginary coordinate values for this case were 0.707.)
The beginning of the impulse response
The top graph in Figure 36 shows the first 100 samples of the impulse response of the filter. This will give you some idea of the general
waveform of the impulse response.
Figure 36 
As you can see, the top graph is Figure 36 looks a lot
like a sinusoid.
The complete impulse response
The top graph in Figure 37 shows the first 16384 samples
of the impulse
response.
Figure 37 
This graph shows that although it takes a long time to do so, the impulse response does finally damp down to zero,
at least that appears to be the case for the vertical scale used for the graph.
Move the pole outside the unit circle
Now let’s move the pole to a location that is only slightly outside the unit
circle. We will leave it at an angle of 45 degrees relative to the real
axis and cause its distance from the origin to be 1.0000611206321341.
(The real and imaginary coordinate values for this case were 0.70715.)
The top graph in Figure 38 shows the first 16384 samples
of the impulse
response.
Figure 38 
An unstable recursive filter
As you can see, this impulse response continues to grow over time.
Eventually the values in the impulse response will grow so large the double
precision numeric format in my computer won’t be able to handle them and the
output will become unpredictable. This is what I mean when I refer to an
unstable recursive
filter.
Put the pole on the unit circle
Now I am going to set the real and imaginary coordinate values for the pole
to 0.70710678118. This will keep the pole at an angle of 45 degrees
relative to the real axis and cause the distance from the origin to the pole to
be 0.9999999999907404 (displayed as a value of type double).
This pole is as close to being on the unit circle as I am able to position it.
The top graph in Figure 39 shows the first 16384 samples of the impulse
response.
Figure 39 
Is the filter output growing or shrinking?
I can’t say with absolute certainty that the amplitude of the filter output is neither growing nor
shrinking. However, if the amplitude is changing over time, the magnitude of the
change is very small.
A single frequency oscillator
This is what I meant when I said that placing a pole
on the unit circle will cause the recursive filter to become an
oscillator.
Ideally, the output for this filter will continue producing the sinusoidal
waveform shown in Figure 36 indefinitely without growing or shrinking.
(If you wanted an oscillator with a more complex waveform, you could
put more than one pole on the unit circle.)
The amplitude response
The amplitude response for this filter as shown in Figure 39 is also quite
interesting. As you can see, it consists of a very narrow spike.
No matter how much I increase the vertical sensitivity of the plot, I can’t see
any grass growing around the bottom of the spike.
(However, this probably means that the estimate of the amplitude
response is going through zero at the computation points, and if I were to
compute at half the frequency interval without changing the length of the
sample, I would probably see some grass.)
In any event, this strongly suggests that if this filter receives any input
at all, it will "ring" at a constant frequency for a very long time,
which is a characteristic of a very narrow band filter. If it rings
without growing or shrinking forever, it is an oscillator.
A realworld analogy
The previous three examples are analogous to a child playing on a swing at the playground. As
long as the child continues pumping energy into the process, the magnitude of
the excursions of the swing will continue to grow. This is analogous to
having a pole outside the unit circle.
(Of course there are physical factors at work on the playground that
will limit the size of the excursions at some point, no matter how much
energy the child attempts to pump into the process.)
As every child knows, if they stop pumping, the excursions will eventually
damp down to zero, although that may take quite a while. This is analogous
to having a pole inside the unit circle and hitting the filter with an impulse.
If the child pumps just the right amount of energy into the process during
each cycle, the excursions will become constant and neither grow nor shrink.
This is analogous to having a pole exactly on the unit circle.
The pendulum of a clock
This is
also analogous to the pendulum in a clock where the spring mechanism pumps a
little bit of energy into the process during each cycle of the pendulum, keeping
it swinging at a constant rate until the energy in the spring is dissipated.
Conscientious clock owners rarely let that happen so that the pendulums in many
clocks oscillate at the same rate for years on end.
A bandpass filter
Now I will demonstrate how to use the workbench to design a
bandpass filter.
(Again, however, I want to emphasize that the workbench is intended
primarily as a learning tool and not as a design tool. The purpose of
the workbench is to help you develop an understanding of recursive
filtering so that you will understand what you are doing when you use other
design techniques.)
Will develop in two stages
I will develop this filter in two stages. In the first stage, I will
position a set of sixteen poles in the pass band to illustrate how the poles contribute to the
behavior of the bandpass filter. In the second stage, I will add sixteen
zeros to show how zeros can be used to improve the outofband reject
characteristics of the filter as well as to improve the flatness in the pass
band.
The filter specifications
Assume that the specifications call for a filter that will pass energy at all
frequencies between 18.3percent and 26percent of the folding frequency and
reject energy at all other frequencies. When viewed in the zplane, this
represents the angular region of the unit circle extending from 33 degrees to 47
degrees. Therefore, we know at the outset that we will primarily be
dealing with that angular range of the unit circle.
The poles and the zeros
Given the above specification, when I started to develop this filter, I knew generally where the poles and
zeros should be to produce the desired behavior. However, I didn’t know exactly where
they should be or how they should be arranged.
An interactive graphic design process
I developed the initial polezero configuration for the filter
by clicking in the zplane to locate the poles and zeros and then observing the
results of that placement in the frequency domain. Once I had determined
approximately where the poles and
zeros should be using the graphic approach, I resorted to the use of a calculator to polish those locations
producing the kind of symmetry that I wanted.
(As mentioned earlier, the ability to design a filter by clicking in the zplane could be
improved by enlarging the image of the zplane. This would make it
possible to use the mouse to locate the poles and zeros more accurately.)
The locations of the poles
The design rationale for poles
As mentioned above, this was an iterative click and evaluate process, but it
definitely wasn’t a process of randomly positioning the poles until something
good was produced. Generally the rationale
for this pole configuration included two important considerations:
 Each edge of the pass band should be well defined. Therefore, I
placed one pole at each edge of the pass band close to the unit circle
(four poles total when you count those in the bottom and the top halves of
the unit circle.).  As you traverse the pass band along the unit circle, the product of the distance to all of the
poles should remain relatively constant. This caused me to come up
with the pole configuration shown in Figure 40.
The amplitude response for poles only
With all the zeros located at the origin, this pole configuration resulted in
the amplitude response shown in Figure 41.
Figure 41 
Now add the zeros
Zeros near the edge of the pass band
As you can see, I added two zeros on the unit circle very close to each of the
four poles that define the edges of the pass band. These eight zeros are
close to, but outside the pass band. There were two purposes for
adding them:
 To pull down the horns at the edges of the pass band in Figure 41.
 To cause the sides of the pass band to be more vertical than is the case
in Figure 41.
The remaining eight zeros
Beyond that, I added eight zeros at 15degree intervals on each side of the pass band. This
was done to hold the response down in regions far removed from the
pass band.
Was this a successful design?
You can judge for yourself how successful this was. The amplitude
response of the resulting recursive filter is shown in the second and third
graphs of Figure 43.
Figure 43 
The amplitude response
The pass band is reasonably flat. By working with the plotting
parameters in Figure 43, I determined that the pass band is
flat to within plus or minus about 1.8 decibels.
As you can see in Figure 43, the pass band has nice vertical sides.
Again, by working with the plotting parameters, I determined that the first
lobe outside the pass band is down by about 24 decibels. Other than the
first lobe, the maximum response outside the pass band is down by about 36
decibels.
The impulse response
As you can see from the first graph in Figure 43, the impulse response for
this filter has a very long tail. Approximately 2000 output samples are
required for the impulse response to damp down to near zero. This suggests
that a nonrecursive convolution filter capable of providing this same bandpass
filter would be quite long, requiring thousands of
MADs
per output sample. On the other hand, this recursive filter requires only 33
MADs per output
sample.
The numeric coordinates
The numeric coordinates of the poles and zeros shown in Figure 42 are listed
in Figure 44.
Figure 44 
The filter coefficients
The recursive filter coefficients are shown in Figure 45.
FeedForward (Numerator) Coefficients 0 1.000 1 11.176 2 61.512 3 221.102 4 580.629 5 1180.833 6 1923.497 7 2559.773 8 2812.494 9 2559.773 10 1923.497 11 1180.833 12 580.629 13 221.102 14 61.512 15 11.176 16 1.000 Feedback (Denominator) Coefficients 1 11.140 2 60.899 3 215.617 4 551.413 5 1077.473 6 1661.374 7 2059.823 8 2074.043 9 1701.302 10 1133.309 11 606.977 12 256.488 13 82.798 14 19.303 15 2.914 16 0.216 Figure 45 
Run the Program
I encourage you to copy the code from Listing 1 into your text
editor, compile it, and execute it. Experiment with it, making
changes, and observing the results of your changes.
Remember that in addition to the code from Listing 1,
this workbench program requires access to the following source code or class
files, which were published in earlier tutorials (see the
References section of this document):
 GraphIntfc01
 ForwardRealToComplexFFT01
 Graph03
(You can also find the source code for these classes by searching for
the class names along with the keywords baldwin java on
Google.)
Having compiled all the source code, enter the following command at the
command prompt to run this program:
java Graph03 Dsp046
Summary
This document, which is the first part of a multipart lesson on recursive
filtering, has provided an overview of an interactive Recursive Filtering
Workbench that can be used to design, experiment with, and evaluate the
behavior of digital recursive filters.
Numerous practical examples of recursive filtering were presented and
explained. Hopefully the study of this lesson has helped you to gain a
better understanding of the behavior of digital recursive filters.
What’s Next?
An understanding of this workbench program requires an understanding of the following
Java classes, which are new
to this program (plus some anonymous classes not listed here):
 Dsp046 – Driver class for the workbench.
 InputGUI – Provides user input capability (mainly text)
for the workbench as shown in Figure 9.  InputGUI$MyTextListener – Processes text input to the workbench, updating the angle output and updating the
zplane
display.  InputGUI$ZPlane – Provides the zplane display along with graphic
mouse input capability for the workbench as shown in
Figure 8.
In addition, an understanding of this program requires an understanding of
the following Java classes, which were presented and explained in earlier lessons:
 ForwardRealToComplexFFT01 – Used to perform an
FFT on
the impulse response.  Graph03 – Plotting program as shown in Figure 10.
 GraphIntfc01 – Required to use Graph03 for plotting.
 GUI – Required to use Graph03 for plotting.
 GUI$MyCanvas – Required to use Graph03 for plotting.
This installment of this multipart lesson has provided an overview of the
workbench.
Although it is possible that I may change the schedule as I write and publish
the future installments of this lesson, here are my plans at this point in time.
The second installment will present and explain the class named Dsp046,
and will also explain its relationship to the classes in the second list above.
The third installment will present and explain the class named InputGUI
along with the inner class named InputGUI$MyTextListener.
The fourth installment will present and explain the class named
InputGUI$MyTextListener.
References
An understanding of the material in the following previously published
lessons will be very helpful to you in understanding the material in this
lesson:

1468 Plotting Engineering and Scientific Data using Java  100 Periodic
Motion and Sinusoids  104 Sampled
Time Series  108
Averaging Time Series  1478
Fun with Java, How and Why Spectral Analysis Works  1482
Spectrum Analysis using Java, Sampling Frequency, Folding Frequency, and the
FFT Algorithm  1483
Spectrum Analysis using Java, Frequency Resolution versus Data Length  1484
Spectrum Analysis using Java, Complex Spectrum and Phase Angle  1485
Spectrum Analysis using Java, Forward and Inverse Transforms, Filtering in
the Frequency Domain  1486
Fun with Java, Understanding the Fast Fourier Transform (FFT) Algorithm
 1487
Convolution and Frequency Filtering in Java  1488
Convolution and Matched Filtering in Java  1492
Plotting Large Quantities of Data using Java  Other
previouslypublished lessons on DSP including adaptive processing and image
processing
In addition, there are numerous good references on DSP available on the web.
For example, good references can be found at the following URLs:
Complete Program Listing
A complete listing of the program discussed in this lesson is shown in
Listing 1 below.
/* File Dsp046.java Copyright 2006, R.G.Baldwin This program provides a visual, interactive recursive filtering workbench. The purpose of the program is to make it easy to experiment with the behavior of recursive filters and to visualize the results of those experiments. The program implements a recursive filter having eight pairs of complex conjugate poles and eight pairs of complex conjugate zeros. The locations of the pairs of poles and zeros in the zplane are controlled by the user. Although the pairs of poles and zeros can be colocated on the real axis, the program does not support the placement of individual poles and zeros on the real axis. The user can reduce the number of poles and zeros used by the recursive filter by moving excess poles and zeros to the origin in the zplane, rendering them ineffective in the behavior of the recursive filter. The program provides three interactive displays on the screen. The first display (in the leftmost position on the screen) contains five graphs. The first graph in this display shows the impulse response of the recursive filter in the time domain. The second and third graphs in the first display show the amplitude response of the recursive filter in the frequency domain computed using two different approaches. The two different computational approaches are provided for comparison purposes. The first computational approach for computing the amplitude response is to perform a Fourier transform on the impulse response using an FFT algorithm. The quality of the estimate of the amplitude response using this approach is dependent on the extent to which the entire impulse response is captured in the set of samples used to perform the FFT. If the impulse response is truncated, the estimate will be degraded. The second approach for computing the amplitude response involves computing the product of the vector lengths from each point on the unit circle to each of the poles and each of the zeros in the zplane. This approach provides an idealized estimate of the amplitude response of the recursive filter, unaffected by impulseresponse considerations. This approach provides the same results that should be produced by performing the FFT on a set of impulseresponse samples of sufficient length to guarantee that the values in the impulse response have been damped down to zero (the impulse response is totally captured in the set of samples on which the FFT is performed). The fourth and fifth graphs in the first display show the phase response of the recursive filter computed using the same (or similar) approaches described above for the amplitude response. (The second approach uses the sum of vector angles instead of the product of vector lengths.) Once again, the two approaches are provided for comparison purposes. By default, the program computes and captures the impulse response for a length of 1024 samples and performs the Fourier transform on that length. However, the length of the captured impulse response and the corresponding FFT can be changed by the user to any length between 2 samples and 16384 samples, provided that the length is an even power of two. (If the length specified by the user is not an even power of two, it is automatically changed to an even power of two by the program.) The first display is interactive in the sense that there are seven different plotting parameters that can be adjusted by the user in order to produce plots that are visually useful in terms of the vertical scale, the horizontal scale, the location of tic marks, etc. The user can modify any of the parameters and then click a Graph button to have the graphs replotted using the new parameters. The second display (which appears in the upper center of the screen) shows the locations of all of the poles and zeros in the zplane. The user can use the mouse to change the location of any pair of complex conjgate poles or zeros by first selecting a specific pair of poles or zeros and then clicking the new location in the zplane. This interactive capability makes it possible for the user to modify the design of the recursive filter in a completely graphic manner by positioning the poles and zeros in the zplane with the mouse. Having relocated one or more pairs of poles or zeros in the zplane, the user can then click the Graph button in the first display described earlier to cause the new impulse response, the new amplitude response, and the new phase response of the new recursive filter with the modified pole and zero locations to be computed and displayed. The third display (that appears in the upperright of the screen) is a control panel that uses text fields, ordinary buttons, and radio buttons to allow the user to perform the following tasks. 1. Specify a new length for the impulse response as an even power of two. (Once again, if the user fails to specify an even power of two, the value provided by the user is converted to an even power of two by the program.) 2. Cause all of the poles to be moved to the origin in the zplane. 3. Cause all of the zeros to be moved to the origin in the zplane. 4. Select a particular pair of complex conjugate poles or zeros to be relocated using the mouse in the display of the zplane. 5. View the angle described by each pole and zero relative to the origin and the horizontal axis in the zplane. 6. View the length of an imaginary vector that connects each pole and zero to the origin. 7. Enter (into a text field) a new real or imaginary value specifying the location of a pair of complex conjugate poles or zeros in the zplane. When a new real or imaginary value is entered, the angle and the length are automatically updated to reflect the new location of the pole or zero and the display of the pole or zero in the zplane is also updated to show the new location. Conversely, when the mouse is used to relocate a pole or zero in the zplane display, the corresponding real and imaginary values in the text fields and the corresponding angle and length are automatically updated to reflect the new location for the pole or zero. Note that mainly for convenience (but for technical reasons as well), the angle and length in the text fields and the location of the pole or zero in the Zplane display are updated as the user enters each character into the text field. If the user enters text into the text field that cannot be converted into a numeric value of type double, (such as the pair of characters "." for example) the contents of the text field are automatically converted to a single "0" character. A warning message is displayed on the commandline screen when this happens. There is no longterm harm when this occurs, but the user may need to start over to enter the new value. Thus, the user should exercise some care regarding the order in which the characters in the text field are modified when entering new real and imaginary values. Each time the impulse response and the spectral data are plotted, the seventeen feedforward filter coefficients and the sixteen feedback coefficients used by the recursive filter to produce the output being plotted are displayed on the commandline screen. Usage: This program requires access to the following source code or class files, which were published in earlier tutorials: GraphIntfc01 ForwardRealToComplexFFT01 Graph03 Enter the following command at the command prompt to run the program: java Graph03 Dsp046. Tested using J2SE 5.0 under WinXP. J2SE 5.0 or later is required due to the use of static imports and printf. **********************************************************/ import static java.lang.Math.*; import java.awt.*; import java.awt.event.*; import javax.swing.*; class Dsp046 implements GraphIntfc01{ //The value stored in the following variable specifies // the number of samples of the impulse response that are // captured. The impulse response serves as the input to // an FFT for the purpose of estimating the amplitude and // phase response of the recursive filter. This data // length must be a power of 2 for the FFT program to // work correctly. If the user enters a value for the // data length that is not a power of two, the value is // automatically converted to a power of two by the // program. int dataLength; double[] impulseResponse; double[] fourierAmplitudeRespnse; double[] fourierPhaseAngle; double[] vectorAmplitudeResponse; double[] vectorPhaseAngle; //This variable contains a reference to a user input GUI // containing buttons, radio buttons, and text fields. InputGUI inputGUI = null; //// Dsp046(){//constructor //If the InputGUI object doesn't already exist, // create it. However, if it already exists, retrieve // the reference to the object from a static variable // belonging to the InputGUI class. This is // necessary because a new object of the Dsp046 class // is instantiated each time the user clicks the Graph // button on the main (Graph03) GUI. However, the // InputGUI object needs to persist across many // clicks of that button because it stores the state of // the poles and zeros designed by the user. When the // InputGUI object is first created, its pole and // zero text fields are initialized with a set of // default pole and zero data values. Its data length // text field is initialized to 1024 samples. if(InputGUI.refToObj == null){ //Instantiate a new InputGUI object. inputGUI = new InputGUI(); //Initialize the length of the array objects that // will contain pole and zero data to a value that is // maintained in the InputGUI object. double[] defaultPoleReal = new double[inputGUI.numberPoles/2]; double[] defaultPoleImag = new double[inputGUI.numberPoles/2]; double[] defaultZeroReal = new double[inputGUI.numberZeros/2]; double[] defaultZeroImag = new double[inputGUI.numberZeros/2]; //Create the default data for eight pairs of complex // conjugate poles spaced at 20degree intervals // around the unit circle. These are the locations // of the poles in the complex zplane. The poles // are barely (0.995) inside the unit circle. for(int cnt = 0;cnt < inputGUI.numberPoles/2;cnt++){ defaultPoleReal[cnt] = 0.995*cos((20+20*cnt)*PI/180.0); defaultPoleImag[cnt] = 0.995*sin((20+20*cnt)*PI/180.0); }//end for loop //Create the default data for eight pairs of complex // conjugate zeros spaced at 20degree intervals // around the unit circle. These are the locations // of the zeros in the complex zplane. The zero // positions are half way between the pole positions. for(int cnt = 0;cnt < inputGUI.numberZeros/2;cnt++){ defaultZeroReal[cnt] = 0.995*cos((10+20*cnt)*PI/180.0); defaultZeroImag[cnt] = 0.995*sin((10+20*cnt)*PI/180.0); }//end for loop //At various points in the program, you may notice // that I have performed separate iterations on // poles and zeros even though the number of poles // is the same as the number of zeros, and I could // have combined them. I did this to make it // possible to modify the number of poles or the // number of zeros later without the requirement for // a major overhaul of the program source code. //Initialize the real and imaginary text fields in // the InputGUI object with the default real and // imaginary pole and zero values. //Start by setting the pole values. for(int cnt = 0;cnt < inputGUI.numberPoles/2;cnt++){ inputGUI.poleReal[cnt].setText( String.valueOf(defaultPoleReal[cnt])); inputGUI.poleImag[cnt].setText( String.valueOf(defaultPoleImag[cnt])); }//end for loop //Now set the zero values. for(int cnt = 0;cnt < inputGUI.numberZeros/2;cnt++){ inputGUI.zeroReal[cnt]. setText(String.valueOf(defaultZeroReal[cnt])); inputGUI.zeroImag[cnt].setText( String.valueOf(defaultZeroImag[cnt])); }//end for loop //Get the default data length from the new object. // There is no requirement to convert it to a power // of two because a power of two is hardcoded into // the program as the default value. dataLength = Integer.parseInt( inputGUI.dataLengthField.getText()); }else{//An InputGUI object already exists. //Retrieve the reference to the existing object that // was saved earlier. inputGUI = InputGUI.refToObj; //Get the current data length from the object. This // value may have been modified by the user and may // not be an even power of two. Convert it to an // even power of two and store the converted value // back into the text field in the existing object. dataLength = Integer.parseInt( inputGUI.dataLengthField.getText()); dataLength = convertToPowerOfTwo(dataLength); inputGUI.dataLengthField.setText("" + dataLength); }//end if //Establish the length of some arrays based on the // current data length impulseResponse = new double[dataLength]; fourierAmplitudeRespnse = new double[dataLength]; fourierPhaseAngle = new double[dataLength]; //Compute and save the impulse response of the filter. // The following code implements a recursive filtering // operation based on the poles and zeros previously // established. However, the conversion from poles // and zeros to feedForward and feedback coefficients // are not routinely involved in the application of a // recursive filter to input data. Therefore, there is // more code here than would be needed for a routine // recursive filtering operation. //Invoke the captureTextFieldData method on the // InputGUI object to cause the current values in the // text fields describing the poles and zeros to be // converted into double values and stored in arrays. inputGUI.captureTextFieldData(); //Create the feedback coefficient array based on the // values stored in the text fields of the InputGUI // object. The process is to first multiply the roots // corresponding to each of eight pairs of complex // conjugate poles. This produces eight secondorder // polynomials. These secondorder polynomials are // multiplied in pairs to produce four fourthorder // polynomials. These four fourthorder polynomials // are multipled in pairs to produce two eighthorder // polynomials. The two eighthorder polynomials are // multiplied to produce one sixteenthorder // polynomial. //The algbraic sign of the real and imag values were // changed to make them match the format of a root // located at a+jb. The format of the root is // (xajb) //Multiply two pairs of complex conjugate roots to // produce two secondorder polynomials. double[] temp1 = conjToSecondOrder( inputGUI.poleRealData[0],inputGUI.poleImagData[0]); double[] temp2 = conjToSecondOrder( inputGUI.poleRealData[1],inputGUI.poleImagData[1]); //Multiply a pair of secondorder polynomials to // produce a fourthorder polynomial. double[] temp3 = secondToFourthOrder( temp1[1],temp1[2],temp2[1],temp2[2]); //Multiply two pairs of complex conjugate roots to // produce two secondorder polynomials. double[] temp4 = conjToSecondOrder( inputGUI.poleRealData[2],inputGUI.poleImagData[2]); double[] temp5 = conjToSecondOrder( inputGUI.poleRealData[3],inputGUI.poleImagData[3]); //Multiply a pair of secondorder polynomials to // produce a fourthorder polynomial. double[] temp6 = secondToFourthOrder( temp4[1],temp4[2],temp5[1],temp5[2]); //Multiply a pair of fourthorder polynomials to // produce an eighthorder polynomial. double[] temp7 = fourthToEighthOrder( temp3[1],temp3[2],temp3[3],temp3[4], temp6[1],temp6[2],temp6[3],temp6[4]); //Multiply two pairs of complex conjugate roots to // produce two secondorder polynomials. double[] temp11 = conjToSecondOrder( inputGUI.poleRealData[4],inputGUI.poleImagData[4]); double[] temp12 = conjToSecondOrder( inputGUI.poleRealData[5],inputGUI.poleImagData[5]); //Multiply a pair of secondorder polynomials to // produce a fourthorder polynomial. double[] temp13 = secondToFourthOrder( temp11[1],temp11[2],temp12[1],temp12[2]); //Multiply two pairs of complex conjugate roots to // produce two secondorder polynomials. double[] temp14 = conjToSecondOrder( inputGUI.poleRealData[6],inputGUI.poleImagData[6]); double[] temp15 = conjToSecondOrder( inputGUI.poleRealData[7],inputGUI.poleImagData[7]); //Multiply a pair of secondorder polynomials to // produce a fourthorder polynomial. double[] temp16 = secondToFourthOrder( temp14[1],temp14[2],temp15[1],temp15[2]); //Multiply a pair of fourthorder polynomials to // produce an eighthorder polynomial. double[] temp17 = fourthToEighthOrder( temp13[1],temp13[2],temp13[3],temp13[4], temp16[1],temp16[2],temp16[3],temp16[4]); //Perform the final polynomial multiplation, // multiplying a pair of eighthorder polynomials to // produce a sixteenthorder polynomial. Place the // coefficients of the sixteenthorder polynomial in // the feedback coefficient array. double[] feedbackCoefficientArray = eighthToSixteenthOrder( temp7[1],temp7[2],temp7[3],temp7[4], temp7[5],temp7[6],temp7[7],temp7[8], temp17[1],temp17[2],temp17[3],temp17[4], temp17[5],temp17[6],temp17[7],temp17[8]); //Determine the length of the delay line required to // perform the feedback arithmetic. int feedbackDelayLineLength = feedbackCoefficientArray.length; //Create the feedForward coefficient array based on // the values stored in the text fields of the // InputGUI object. The process is the same as // described above for the poles. //Multiply two pairs of complex conjugate roots to // produce two secondorder polynomials. double[] temp21 = conjToSecondOrder( inputGUI.zeroRealData[0],inputGUI.zeroImagData[0]); double[] temp22 = conjToSecondOrder( inputGUI.zeroRealData[1],inputGUI.zeroImagData[1]); //Multiply a pair of secondorder polynomials to // produce a fourthorder polynomial. double[] temp23 = secondToFourthOrder( temp21[1],temp21[2],temp22[1],temp22[2]); //Multiply two pairs of complex conjugate roots to // produce two secondorder polynomials. double[] temp24 = conjToSecondOrder( inputGUI.zeroRealData[2],inputGUI.zeroImagData[2]); double[] temp25 = conjToSecondOrder( inputGUI.zeroRealData[3],inputGUI.zeroImagData[3]); //Multiply a pair of secondorder polynomials to // produce a fourthorder polynomial. double[] temp26 = secondToFourthOrder( temp24[1],temp24[2],temp25[1],temp25[2]); //Multiply a pair of fourthorder polynomials to // produce an eighthorder polynomial. double[] temp27 = fourthToEighthOrder( temp23[1],temp23[2],temp23[3],temp23[4], temp26[1],temp26[2],temp26[3],temp26[4]); //Multiply two pairs of complex conjugate roots to // produce two secondorder polynomials. double[] temp31 = conjToSecondOrder( inputGUI.zeroRealData[4],inputGUI.zeroImagData[4]); double[] temp32 = conjToSecondOrder( inputGUI.zeroRealData[5],inputGUI.zeroImagData[5]); //Multiply a pair of secondorder polynomials to // produce a fourthorder polynomial. double[] temp33 = secondToFourthOrder( temp31[1],temp31[2],temp32[1],temp32[2]); //Multiply two pairs of complex conjugate roots to // produce two secondorder polynomials. double[] temp34 = conjToSecondOrder( inputGUI.zeroRealData[6],inputGUI.zeroImagData[6]); double[] temp35 = conjToSecondOrder( inputGUI.zeroRealData[7],inputGUI.zeroImagData[7]); //Multiply a pair of secondorder polynomials to // produce a fourthorder polynomial. double[] temp36 = secondToFourthOrder( temp34[1],temp34[2],temp35[1],temp35[2]); //Multiply a pair of fourthorder polynomials to // produce an eighthorder polynomial. double[] temp37 = fourthToEighthOrder( temp33[1],temp33[2],temp33[3],temp33[4], temp36[1],temp36[2],temp36[3],temp36[4]); //Perform the final polynomial multiplation, // multiplying a pair of eighthorder polynomials to // produce a sixteenthorder polynomial. Place the // coefficients of the sixteenthorder polynomial in // the feedForward coefficient array. double[] feedForwardCoefficientArray = eighthToSixteenthOrder( temp27[1],temp27[2],temp27[3],temp27[4], temp27[5],temp27[6],temp27[7],temp27[8], temp37[1],temp37[2],temp37[3],temp37[4], temp37[5],temp37[6],temp37[7],temp37[8]); //Determine the length of the delay line required to // perform the feedForward arithmetic. int feedForwardDelayLineLength = feedForwardCoefficientArray.length; //Display the feedForward and feedback coefficients System.out.println( "FeedForward (Numerator) Coefficients"); for(int cnt = 0; cnt < feedForwardCoefficientArray.length; cnt++){ System.out.printf ("%2d % 5.3f%n",cnt, feedForwardCoefficientArray[cnt]); }//end for loop System.out.println( "nFeedback (Denominator) Coefficients"); for(int cnt = 1; cnt < feedbackCoefficientArray.length; cnt++){ System.out.printf ("%2d % 5.3f%n",cnt, feedbackCoefficientArray[cnt]); }//end for loop System.out.println(); //Create the data delay lines used for feedForward // and feedback arithmetic. double[] feedForwardDelayLine = new double[feedForwardDelayLineLength]; double[] feedbackDelayLine = new double[feedbackDelayLineLength]; //Initial input and output data values double filterInputSample = 0;//input data double filterOutputSample = 0;//output data //Compute the output values and populate the output // array for further analysis such as FFT analysis. // This is the code that actually applies the // recursive filter to the input data given the // feedForward and feedback coefficients. for(int dataLenCnt = 0;dataLenCnt < dataLength; dataLenCnt++){ //Create the input samples consisting of a single // impulse at time zero and sample values of 0 // thereafter. if(dataLenCnt == 0){ filterInputSample = 100.0; }else{ filterInputSample = 0.0; }//end else //*************************************************// //This is the beginning of one cycle of the actual // recursive filtering process. //Shift the data in the delay lines. Oldest value // has the highest index value. for(int cnt = feedForwardDelayLineLength1; cnt > 0;cnt){ feedForwardDelayLine[cnt] = feedForwardDelayLine[cnt1]; }//end for loop for(int cnt = feedbackDelayLineLength1; cnt > 0;cnt){ feedbackDelayLine[cnt] = feedbackDelayLine[cnt1]; }//end for loop //Insert the input signal into the delay line at // zero index. feedForwardDelayLine[0] = filterInputSample; //Compute sum of products for input signal and // feedForward coefficients from 0 to // feedForwardDelayLineLength1. double xTemp = 0; for(int cnt = 0;cnt < feedForwardDelayLineLength; cnt++){ xTemp += feedForwardCoefficientArray[cnt]* feedForwardDelayLine[cnt]; }//end for loop //Compute sum of products for previous output values // and feedback coefficients from 1 to // feedbackDelayLineLength1. double yTemp = 0; for(int cnt = 1;cnt < feedbackDelayLineLength;cnt++){ yTemp += feedbackCoefficientArray[cnt]* feedbackDelayLine[cnt]; }//end for loop //Compute new output value as the difference. filterOutputSample = xTemp  yTemp; //Save the output value in the array containing the // impulse response. impulseResponse[dataLenCnt] = filterOutputSample; //Insert the output signal into the delay line at // zero index. feedbackDelayLine[0] = filterOutputSample; //This is the end of one cycle of the recursive // filtering process. //*************************************************// }//end for loop //Now compute the Fourier Transform of the impulse // response, placing the magnitude result from the FFT // program into the fourierAmplitudeRespnse array and // the phase angle result in the fourierPhaseAngle // array, each to be plotted later. ForwardRealToComplexFFT01.transform( impulseResponse, new double[dataLength], new double[dataLength], fourierPhaseAngle, fourierAmplitudeRespnse); /*DO NOT DELETE THIS CODE //Note that this normalization code is redundant // because of the normalization that takes place in // the method named convertToDB later. However, if // the conversion to decibels is disabled, the // following code should be enabled. //Scale the fourierAmplitudeRespnse to compensate for // the differences in data length. for(int cnt = 0;cnt < fourierAmplitudeRespnse.length; cnt++){ fourierAmplitudeRespnse[cnt] = fourierAmplitudeRespnse[cnt]*dataLength/16384.0; }//end for loop */ //Compute the amplitude response based on the ratio of // the products of the pole and zero vectors. vectorAmplitudeResponse = getVectorAmplitudeResponse( inputGUI.poleRealData, inputGUI.poleImagData, inputGUI.zeroRealData, inputGUI.zeroImagData); //Compute the phase angle based on sum and difference // of the angles of the pole and zero vectors. vectorPhaseAngle = getVectorPhaseAngle( inputGUI.poleRealData, inputGUI.poleImagData, inputGUI.zeroRealData, inputGUI.zeroImagData); //Convert the fourierAmplitudeRespnse data to decibels. // Disable the following statement to disable the // conversion to decibels. convertToDB(fourierAmplitudeRespnse); //Convert the vectorAmplitudeResponse to decibels. // Disable the following statement to disable the // conversion to decibels. convertToDB(vectorAmplitudeResponse); }//end constructor //// //The purpose of this method is to convert an incoming // array containing amplitude response data to decibels. void convertToDB(double[] magnitude){ //Eliminate or modify all values that are incompatible // with conversion to log base 10 and the log10 method. // Also limit small values to be no less than 0.0001. for(int cnt = 0;cnt < magnitude.length;cnt++){ if((magnitude[cnt] == Double.NaN)  (magnitude[cnt] <= 0.0001)){ magnitude[cnt] = 0.0001; }else if(magnitude[cnt] == Double.POSITIVE_INFINITY){ magnitude[cnt] = 9999999999.0; }//end else if }//end for loop //Find the peak value for use in normalization. double peak = 9999999999.0; for(int cnt = 0;cnt < magnitude.length;cnt++){ if(peak < abs(magnitude[cnt])){ peak = abs(magnitude[cnt]); }//end if }//end for loop //Normalize to the peak value to make the values easier // to plot with regard to scaling. for(int cnt = 0;cnt < magnitude.length;cnt++){ magnitude[cnt] = magnitude[cnt]/peak; }//end for loop //Now convert normalized magnitude data to log base 10. for(int cnt = 0;cnt < magnitude.length;cnt++){ magnitude[cnt] = log10(magnitude[cnt]); }//end for loop }//end convertToDB //// //This method makes certain that the incoming value is a // nonzero positive power of two that is less than or // equal to 16384. If the input is not equal to either a // power of two or one less than a power of two, it is // truncated to the next lower power of two. If it is // either a power of two or one less than a power of two, // the returned value is that power of two. Negative // input values are converted to positive values before // making the conversion. int convertToPowerOfTwo(int dataLength){ //Eliminate negative values dataLength = abs(dataLength); //Make sure the data length is not zero by adding a // value of 1. dataLength++; //Make sure the data length is not greater than 16384. if(dataLength > 16384){ dataLength = 16384; }//end if int cnt = 0; int mask = 0x4000; //Loop and left shift left until the msb of the data // length value matches 0x4000. Count the number of // shifts required to make the match. No shifts are // required for a data length of 16384. while((dataLength & mask) == 0){ cnt++; dataLength = dataLength << 1; }//end while loop //Now shift the mask to the right the same number of // places as were required to make the above match. // Because the mask consists of a single bit, this // guarantees that the resulting value is an even // power of two. dataLength = mask >> cnt; return dataLength; }//end convertToPowerOfTwo //// //The amplitude response of the recursive filter at a // given frequency (a point on the unit circle) can be // estimated by dividing the product of the lengths of // the vectors connecting that point on the unit circle // to all of the zeros by the product of the lengths of // the vectors connecting the same point on the unit // circle to all of the poles. The amplitude response // at all frequencies between zero and the Nyquist // folding frequency can be estimated by performing this // calculation at all points in the top half of the unit // circle. //The bottom half of the unit circle provides the // amplitude response for frequencies between the Nyquist // folding frequency and the sampling frequency. These // values are redundant because they are the same as the // values below the folding frequency. //Despite the redundancy, this method computes the // amplitude response at a set of frequencies between // zero and the sampling frequency where the number of // frequencies in the set is equal to the data length. // This causes the amplitude response to be computed at // a set of frequencies that matches the set of // frequencies for which the amplitude response is // computed using the FFT algorithm, making it easier // to plot the two for comparison purposes. double[] getVectorAmplitudeResponse( double[] poleRealData, double[] poleImagData, double[] zeroRealData, double[] zeroImagData){ double[] amplitudeResponse = new double[dataLength]; double freqAngle = 0; double freqHorizComponent = 0; double freqVertComponent = 0; //Divide the unit circle into dataLength frequencies // and compute the amplitude response at each // frequency. for(int freqCnt = 0;freqCnt < dataLength;freqCnt++){ //Get the angle from the origin to the point on the // unit circle relative to the horizontal axis. freqAngle = freqCnt*2*PI/dataLength; //Get the horizontal and vertical components of the // distance from the origin to the point on the unit // circle. freqHorizComponent = cos(freqAngle); freqVertComponent = sin(freqAngle); //Compute the product of the lengths from the point // on the unit circle to each of the zeros. //Get the distance from the point on the unit circle // to each zero as the square root of the sum of the // squares of the sides of a right triangle formed // by the zero and the point on the unit circle with // the base of the triangle being parallel to the // horiontal axis. //Declare some working variables. double base;//Base of the right triangle double height;//Height of the right triangle double hypo;//Hypotenuse of the right triangle double zeroProduct = 1.0;//Initialize the product. //Loop and process all complex conjugate zeros for(int cnt = 0;cnt < zeroRealData.length;cnt++){ //First compute the product for a zero in the // upper half of the zplane //Get base of triangle base = freqHorizComponent  zeroRealData[cnt]; //Get height of triangle height = freqVertComponent  zeroImagData[cnt]; //Get hypotenuse of triangle hypo = sqrt(base*base + height*height); //Compute the running product. zeroProduct *= hypo; //Continue computing the running product using the // conjugate zero in the lower half of the zplane. // Note the sign change on the imaginary // part. base = freqHorizComponent  zeroRealData[cnt]; height = freqVertComponent + zeroImagData[cnt]; hypo = sqrt(base*base + height*height); zeroProduct *= hypo; }//end for loop  all zeros have been processed //Now compute the product of the lengths to the // poles. double poleProduct = 1.0;//Initialize the product. for(int cnt = 0;cnt < poleRealData.length;cnt++){ //Begin with the pole in the upper half of the // zplane. base = freqHorizComponent  poleRealData[cnt]; height = freqVertComponent  poleImagData[cnt]; hypo = sqrt(base*base + height*height); //Compute the running product. poleProduct *= hypo; //Continue computing the running product using the // conjugate pole in the lower half of the zplane. // Note the sign change on the imaginary part. base = freqHorizComponent  poleRealData[cnt]; height = freqVertComponent + poleImagData[cnt]; hypo = sqrt(base*base + height*height); poleProduct *= hypo;//product of lengths }//end for loop //Divide the zeroProduct by the poleProduct. //Compute and save the amplitudeResponse for this // frequency and then go back to the top of the loop // and compute the amplitudeResponse for the next // frequency. amplitudeResponse[freqCnt] = zeroProduct/poleProduct; }//end for loop on data length  all frequencies done return amplitudeResponse; }//end getVectorAmplitudeResponse //// //The phase angle of the recursive filter at a // particular frequency (represented by a point on the // unit circle) can be determined by subtracting the sum // of the angles from the point on the unit circle to all // of the poles from the sum of the angles from the same // point on the unit circle to all of the zeros. //The phase angle for an equallyspaced set of // frequencies between zero and the sampling frequency is // computed in radians and then converted to degrees in // the range from 180 degrees to +180 degrees for return // to the calling program.. double[] getVectorPhaseAngle(double[] poleRealData, double[] poleImagData, double[] zeroRealData, double[] zeroImagData){ double[] phaseResponse = new double[dataLength]; double freqAngle = 0; double freqHorizComponent = 0; double freqVertComponent = 0; //Divide the unit circle into dataLength frequencies // and compute the phase angle at each frequency. for(int freqCnt = 0;freqCnt < dataLength;freqCnt++){ //Note that the following reference to an angle is // not a reference to the phase angle. Rather, it // is a reference to the angle of a point on the unit // circle relative to the horizontal axis and the // origin of the zplane. freqAngle = freqCnt * 2 * PI/dataLength; //Get the horizontal and vertical components of the // distance from the origin to the point on the unit // circle. freqHorizComponent = cos(freqAngle); freqVertComponent = sin(freqAngle); //Begin by processing all of the complex conjugate // zeros. //Compute the angle from the point on the unit circle // to each of the zeros. Retain as an angle between // zero and 2*PI (360 degrees). //Declare some working variables. double base;//Base of a right triangle double height;//Height of a right triangle double zeroAngle; double zeroAngleSum = 0; //Loop and process all complex conjugate zeros for(int cnt = 0;cnt < zeroRealData.length;cnt++){ //Compute using the zero in upper the half of the // zplane. //Get base of triangle base = (freqHorizComponent  zeroRealData[cnt]); //Get height of triangle height = (freqVertComponent  zeroImagData[cnt]); if(base == 0){//Avoid division by zero. zeroAngle = PI/2.0;//90 degees }else{//Compute the angle. zeroAngle = atan(height/base); }//end else //Adjust for negative coordinates if((base < 0) && (height > 0)){ zeroAngle = PI + zeroAngle; }else if((base < 0) && (height < 0)){ zeroAngle = PI + zeroAngle; }else if((base > 0) && (height < 0)){ zeroAngle = 2*PI + zeroAngle; }//end else //Compute the running sum of the angles. zeroAngleSum += zeroAngle; //Continue computing the running sum of angles // using the conjugate zero in the lower half of // the zplane. Note the sign change on the // imaginary part. base = (freqHorizComponent  zeroRealData[cnt]); height = (freqVertComponent + zeroImagData[cnt]); if(base == 0){ zeroAngle = 3*PI/2.0;//270 degees }else{ zeroAngle = atan(height/base); }//end else //Adjust for negative coordinates if((base < 0) && (height > 0)){ zeroAngle = PI + zeroAngle; }else if((base < 0) && (height < 0)){ zeroAngle = PI + zeroAngle; }else if((base > 0) && (height < 0)){ zeroAngle = 2*PI + zeroAngle; }//end else //Add the angle into the running sum of zero // angles. zeroAngleSum += zeroAngle; }//end for loop  all zeros have been processed //Now compute the sum of the angles from the point // on the unit circle to each of the poles. double poleAngle; double poleAngleSum = 0; //Loop and process all complex conjugate poles for(int cnt = 0;cnt < poleRealData.length;cnt++){ //Compute using pole in the upper half of the // zplane base = (freqHorizComponent  poleRealData[cnt]); height = (freqVertComponent  poleImagData[cnt]); if(base == 0){//Avoid division by zero poleAngle = PI/2.0;//90 degees }else{ poleAngle = atan(height/base); }//end else //Adjust for negative coordinates if((base < 0) && (height > 0)){ poleAngle = PI + poleAngle; }else if((base < 0) && (height < 0)){ poleAngle = PI + poleAngle; }else if((base > 0) && (height < 0)){ poleAngle = 2*PI + poleAngle; }//end else //Compute the running sum of the angles. poleAngleSum += poleAngle; //Continue computing the sum of angles using the // conjugate pole in the lower half of the Zplane. // Note the sign change on the imaginary part. base = (freqHorizComponent  poleRealData[cnt]); height = (freqVertComponent + poleImagData[cnt]); if(base == 0){//Avoid division by 0. poleAngle = 3*PI/2.0;//270 degees }else{ poleAngle = atan(height/base); }//end else //Adjust for negative coordinates if((base < 0) && (height > 0)){ poleAngle = PI + poleAngle; }else if((base < 0) && (height < 0)){ poleAngle = PI + poleAngle; }else if((base > 0) && (height < 0)){ poleAngle = 2*PI + poleAngle; }//end else poleAngleSum += poleAngle; }//end for loop  all poles have been processed //Subtract the sum of the pole angles from the sum of // the zero angles. Convert the angle from radians // to degrees in the process. //Note that the minus sign in the following // expression is required to cause the sign of the // angle computed using this approach to match the // sign of the angle computed by the FFT algorithm. // This indicates that either this computation or the // FFT computation is producing a phase angle having // the wrong sign. double netAngle = (zeroAngleSum  poleAngleSum)*180/PI; //Normalize the angle to the range from 180 degrees // to +180 degrees to make it easier to plot. if(netAngle > 180){ while(netAngle > 180){ netAngle = 360; }//end while }else if(netAngle < 180){ while(netAngle < 180){ netAngle += 360; }//end while }//end else if //Save the phase angle for this frequency and then go // back to the top of the loop and compute the phase // angle for the next frequency. phaseResponse[freqCnt] = netAngle; }//end for loop on data length  all frequencies done return phaseResponse; }//end getVectorPhaseAngle //// //Receives the complex conjugate roots of a secondorder // polynomial in the form (a+jb)(ajb). Multiplies the // roots and returns the coefficients of the secondorder // polynomial as x*x + 2*a*x + (a*a + b*b) in a three // element array of type double. double[] conjToSecondOrder(double a,double b){ double[] result = new double[]{1,2*a,(a*a + b*b)}; return result; }//end conjToSecondOrder //// //Receives the coefficients of a pair of secondorder // polynomials in the form: // x*x + a*x + b // x*x + c*x + d //Multiplies the polynomials and returns the coefficients // of a fourthorder polynomial in a fiveelement array // of type double. double[] secondToFourthOrder(double a,double b, double c,double d){ double[] result = new double[]{1, a + c, b + c*a + d, c*b + d*a, d*b}; return result; }//end secondToFourthOrder //// //Receives the coeficients of a pair of fourth order // polynomials in the form: // x*x*x*x + a*x*x*x + b*x*x + c*x + d // x*x*x*x + e*x*x*x + f*x*x + g*x + h //Multiplies the polynomials and returns the coefficients // of an eighthorder polynomial in a nineelement // array of type double. double[] fourthToEighthOrder(double a,double b, double c,double d, double e,double f, double g,double h){ double[] result = new double[]{ 1, a + e, b + e*a + f, c + e*b + f*a + g, d + e*c + f*b + g*a + h, e*d + f*c + g*b + h*a, f*d + g*c + h*b, g*d + h*c, h*d}; return result; }//end fourthToEighthOrder //// //Receives the coefficients of a pair of eighth order // polynomials in the following form where xn indicates // x to the nth power: // x8 + ax7 + bx6 + cx5 + dx4 + ex3 + fx2 + gx + h // x8 + ix7 + jx6 + kx5 + lx4 + mx3 + nx2 + ox + p //Multiplies the polynomials and returns the coefficients // of a sixteenthorder polynomial in a 17element // array of type double double[] eighthToSixteenthOrder(double a,double b, double c,double d, double e,double f, double g,double h, double i,double j, double k,double l, double m,double n, double o,double p){ double[] result = new double[]{ 1, a+i, b+i*a+j, c+i*b+j*a+k, d+i*c+j*b+k*a+l, e+i*d+j*c+k*b+l*a+m, f+i*e+j*d+k*c+l*b+m*a+n, g+i*f+j*e+k*d+l*c+m*b+n*a+o, h+i*g+j*f+k*e+l*d+m*c+n*b+o*a+p, i*h+j*g+k*f+l*e+m*d+n*c+o*b+p*a, j*h+k*g+l*f+m*e+n*d+o*c+p*b, k*h+l*g+m*f+n*e+o*d+p*c, l*h+m*g+n*f+o*e+p*d, m*h+n*g+o*f+p*e, n*h+o*g+p*f, o*h+p*g, p*h}; return result; }//end eiththToSixteenthOrder //// //The following six methods are declared in the interface // named GraphIntfc01, and are required by the plotting // program named Graph03. //// //This method specifies the number of functions that will // be plotted by the program named Graph03. public int getNmbr(){ //Return number of functions to // process. Must not exceed 5. return 5; }//end getNmbr //// //This method returns the values that will be plotted in // the first graph by the program named Graph03. public double f1(double x){ //Return the impulse response of the filter. if(((int)x >= 0) && ((int)x < impulseResponse.length)){ return impulseResponse[(int)x]; }else{ return 0; }//end else }//end f1 //// //This method returns the values that will be plotted in // the second graph by the program named Graph03. public double f2(double x){ //Return the amplitude response of the recursive filter // obtained by performing a Fourier Transform on the // impulse response and converting the result to // decibels. Recall that adding a constant to a // decibel plot is equivalent to multiplying the // original data by the constant. if(((int)x >= 0) && ((int)x < fourierAmplitudeRespnse.length)){ return 100 + (100.0 * fourierAmplitudeRespnse[(int)x]); }else{ return 0; }//end else }//end f2 //// //This method returns the values that will be plotted in // the third graph by the program named Graph03. public double f3(double x){ //Return the amplitude response of the recursive filter // obtained by dividing the product of the zero vector // lengths by the product of the pole vector lengths // and converting the result to decibels. Recall that // adding a constant to a decibel plot is equivalent to // multiplying the original data by the constant. if(((int)x >= 0) && ((int)x < vectorAmplitudeResponse.length)){ return 100 + (100.0 * vectorAmplitudeResponse[(int)x]); }else{ return 0; }//end else }//end f3 //// //This method returns the values that will be plotted in // the fourth graph by the program named Graph03. public double f4(double x){ //Return the phase response of the recursive filter // obtained by performing a Fourier Transform on the // impulse response. if(((int)x >= 0) && ((int)x < fourierPhaseAngle.length)){ return fourierPhaseAngle[(int)x]; }else{ return 0; }//end else }//end f4 //// //This method returns the values that will be plotted in // the fifth graph by the program named Graph03. public double f5(double x){ //Return the phase response of the recursive filter // obtained by subtracting the sum of the pole vector // angles from the sum of the zero vector angles. if(((int)x >= 0) && ((int)x < vectorPhaseAngle.length)){ return vectorPhaseAngle[(int)x]; }else{ return 0; }//end else }//end f5 }//end class Dsp046 //=======================================================// //An object of this class stores and displays the real and // imaginary parts of sixteen complex poles and sixteen // complex zeros. The sixteen poles and sixteen zeros // form eight conjugate pairs. Thus, there are eight // complex conjugate pairs of poles and eight complex // conjugate pairs of zeros. //The object also computes and displays the angle in // degrees for each pole and each zero relative to the // origin of the zplane. It also computes and displays // the length of an imaginary vector connecting each // pole and zero to the origin. //The real and imaginary part for each pole or zero is // displayed in a TextField object. The user can modify // the values by entering new values into the text fields. // The user can also modify the real and imaginary values // by selecting a radio button associated with a specific // pole or zero and then clicking a new location for that // pole or zero in an auxiliary display that shows the // complex zplane and the unit circle in that plane. When // the user clicks in the zplane, the corresponding values // in the text fields for the selected pole or zero are // automatically updated. When the user changes a real // or imaginary value in a text field, the radio button // for that pole or zero is automatically selected, and // the image of that pole or zero in the zplane is // automatically changed to show the new position of the // pole or zero. //Regardless of which method causes the contents of a text // field to be modified, a TextListener that is registered // on the text field causes the displayed angle and length // to be updated to match the new real and imaginary // values. class InputGUI{ //A reference to an object of this class is stored in the // following static variable. This makes it possible for // the original object that created this object to cease // to exist without this object becoming eligible for // garbage collection. When that original object is // replaced by a new object, the new object can assume // ownership of this object by getting its reference from // the static variable. Thus, ownership of this object // can be passed along from one object to the next. static InputGUI refToObj = null; //The following ButtonGroup object is used to group // radio buttons to cause them to behave in a mututlly // exclusive way. The zero buttons and the pole buttons // are all placed in the same group so that only one zero // or one pole can be selected at any point in time. ButtonGroup buttonGroup = new ButtonGroup(); //A reference to an auxiliary display showing the // zplane is stored in the following instance variable. ZPlane refToZPlane; //This is the default number of poles and zeros // including the conjugates. These values cannot be // changed by the user. However, the effective number of // poles or zeros can be reduced by moving poles and/or // zeros to the origin in the zplane, rendering them // ineffective in the recursive filtering process. // Moving poles and zeros to the origin causes feedback // and feedforward zeros to go to zero. int numberPoles = 16; int numberZeros = 16; //The following variables refer to a pair of buttons used // to make it possible for the user to place all the // poles and zeros at the origin. In effect, these are // reset buttons relative to the pole and zero locations. JButton clearPolesButton = new JButton("Move Poles to Origin"); JButton clearZerosButton = new JButton("Move Zeros to Origin"); //The following text field stores the default data length // at startup. This value can be changed by the user to // investigate the impact of changes to the data length // for a given set of poles and zeros. TextField dataLengthField = new TextField("1024"); //The following array objects get populated with numeric // values from the text fields by the method named // captureTextFieldData. That method should be called // to populate them with the most current text field // data when the most current data is needed. double[] poleRealData = new double[numberPoles/2]; double[] poleImagData = new double[numberPoles/2]; double[] zeroRealData = new double[numberZeros/2]; double[] zeroImagData = new double[numberZeros/2]; //The following arrays are populated with references to // radio buttons, each of which is associated with a // specific pair of complex conjugate poles or zeros. JRadioButton[] poleRadioButtons = new JRadioButton[numberPoles/2]; JRadioButton[] zeroRadioButtons = new JRadioButton[numberZeros/2]; //The following arrays are populated with references to // text fields, each of which is associated with a // specific pair of complex conjugate poles or zeros. // The text fields contain the real values, imaginary // values, and the values of the angle and the length // specified by the real and imaginary values. //The size of the following arrays is only half the // number of poles and zeros because the conjugate is // generated on the fly when it is needed. //I wanted to use JTextField objects, but JTextField // doesn't have an addTextListener method. I needed to // register a TextListener object on each real and // imaginary text field to compute the angle and length // each time the contents of a text field changes, so I // used the AWT TextField class instead. TextField[] poleReal = new TextField[numberPoles/2]; TextField[] poleImag = new TextField[numberZeros/2]; TextField[] poleAngle = new TextField[numberZeros/2]; TextField[] poleLength = new TextField[numberZeros/2]; TextField[] zeroReal = new TextField[numberZeros/2]; TextField[] zeroImag = new TextField[numberZeros/2]; TextField[] zeroAngle = new TextField[numberZeros/2]; TextField[] zeroLength = new TextField[numberZeros/2]; //// InputGUI(){//constructor //Instantiate a new JFrame object and condition its // close button. JFrame guiFrame = new JFrame("Copyright 2006 R.G.Baldwin"); guiFrame.setDefaultCloseOperation( JFrame.EXIT_ON_CLOSE); //The following JPanel contains a JLabel, a TextField, // and two JButton objects. The label simply provides // instructions to the user regarding the entry of a // new data length. The text field is used for entry // of a new data length value by the user at runtime. // The buttons are used to cause the poles and zeros // to be moved to the origin in two groups. Moving // both the poles and the zeros to the origin // effectively converts the recursive filter to an // allpass filter. //This panel is placed in the NORTH location of the // JFrame resulting in the name northControlPanel JPanel northControlPanel = new JPanel(); northControlPanel.setLayout(new GridLayout(0,2)); northControlPanel. add(new JLabel("Data Length as Power of 2")); northControlPanel.add(dataLengthField); northControlPanel.add(clearPolesButton); northControlPanel.add(clearZerosButton); guiFrame.add(northControlPanel,BorderLayout.NORTH); //Register an action listener on the clearPolesButton // to set the contents of the text fields that // represent the locations of the poles to 0. This // also causes the poles to move to the origin in the // display of the zplane. clearPolesButton.addActionListener( new ActionListener(){ public void actionPerformed(ActionEvent e){ for(int cnt = 0;cnt < numberPoles/2;cnt++){ poleReal[cnt].setText("0"); poleImag[cnt].setText("0"); }//end for loop }//end actionPerformed }//end new ActionListener );//end addActionListener( //Register action listener on the clearZerosButton to // set the contents of the text fields that // represent the locations of the zeros to 0 This // also causes the zeros to move to the origin in the // display of the zplane. clearZerosButton.addActionListener( new ActionListener(){ public void actionPerformed(ActionEvent e){ for(int cnt = 0;cnt < numberZeros/2;cnt++){ zeroReal[cnt].setText("0"); zeroImag[cnt].setText("0"); }//end for loop }//end actionPerformed }//end new ActionListener );//end addActionListener( //The following JPanel object contains two other JPanel // objects, one for zeros and one for poles. They // are the same size with one located above the other. // The panel containing zero data is green. The panel // containing pole data is yellow. This panel is // placed in the CENTER of the JFrame object. Hence // the name centerControlPanel. JPanel centerControlPanel = new JPanel(); centerControlPanel.setLayout(new GridLayout(2,1)); //The following JPanel is populated with text fields // and radio buttons that represent the zeros. JPanel zeroPanel = new JPanel(); zeroPanel.setBackground(Color.GREEN); //The following JPanel is populated with text fields // and radio buttons that represent the poles. JPanel polePanel = new JPanel(); polePanel.setBackground(Color.YELLOW); //Add the panels containing textfields and readio // buttons to the larger centerControlPanel. centerControlPanel.add(zeroPanel); centerControlPanel.add(polePanel); //Add the centerControlPanel to the CENTER of the // JFrame object. guiFrame.getContentPane().add( centerControlPanel,BorderLayout.CENTER); //Instantiate a text listener that will be registered // on each of the text fields containing real and // imaginary values, each pair of which specifies the // location of a pole or a zero. MyTextListener textListener = new MyTextListener(); //A great deal is accomplished in each of the following // two for loops. //Begin by populating the arrays described earlier with // radio buttons and text fields. //Then place the radio buttons associated with the // poles and zeros in the same group to make them // behave in a mutually exclusive manner. //Next, place the components on the polePanel and the // zero panel. //Then disable the text fields containing the angle // and the length to prevent the user from entering // data into them. Note that this does not prevent // the program from writing text into the text fields // containing the angle and the length. The text // fields are disabled only insofar as manual input by // the user is concerned. //After that, set the name property for each of the // real and imaginary text fields. These name property // values will be used later by a common TextListener // object to determine which text field fired a // TextEvent. //Finally, register a common TextListener object on the // real and imaginary text fields to cause the angle // and the length to be computed and displayed when the // text value changes for any reason. //Deal with the poles. polePanel.setLayout(new GridLayout(0,5)); //Place a row of column headers polePanel.add(new JLabel("Poles")); polePanel.add(new JLabel("Real")); polePanel.add(new JLabel("Imag")); polePanel.add(new JLabel("Angle (deg)")); polePanel.add(new JLabel("Length")); //Take the actions described above with respect to the // poles. for(int cnt = 0;cnt < numberPoles/2;cnt++){ poleRadioButtons[cnt] = new JRadioButton("" + cnt); poleReal[cnt] = new TextField("0"); poleImag[cnt] = new TextField("0"); poleAngle[cnt] = new TextField("0"); poleLength[cnt] = new TextField("0"); buttonGroup.add(poleRadioButtons[cnt]); polePanel.add(poleRadioButtons[cnt]); polePanel.add(poleReal[cnt]); polePanel.add(poleImag[cnt]); polePanel.add(poleAngle[cnt]); poleAngle[cnt].setEnabled(false); polePanel.add(poleLength[cnt]); poleLength[cnt].setEnabled(false); poleReal[cnt].setName("poleReal" + cnt); poleImag[cnt].setName("poleImag" + cnt); poleReal[cnt].addTextListener(textListener); poleImag[cnt].addTextListener(textListener); }//end for loop //Deal with the zeros. zeroPanel.setLayout(new GridLayout(0,5)); //Place a row of column headers zeroPanel.add(new JLabel("Zeros")); zeroPanel.add(new JLabel("Real")); zeroPanel.add(new JLabel("Imag")); zeroPanel.add(new JLabel("Angle (deg)")); zeroPanel.add(new JLabel("Length")); //Now take the actions described above with respect to // the zeros. for(int cnt = 0;cnt < numberZeros/2;cnt++){ zeroRadioButtons[cnt] = new JRadioButton("" + cnt); zeroReal[cnt] = new TextField("0"); zeroImag[cnt] = new TextField("0"); zeroAngle[cnt] = new TextField("0"); zeroLength[cnt] = new TextField("0"); buttonGroup.add(zeroRadioButtons[cnt]); zeroPanel.add(zeroRadioButtons[cnt]); zeroPanel.add(zeroReal[cnt]); zeroPanel.add(zeroImag[cnt]); zeroPanel.add(zeroAngle[cnt]); zeroAngle[cnt].setEnabled(false); zeroPanel.add(zeroLength[cnt]); zeroLength[cnt].setEnabled(false); zeroReal[cnt].setName("zeroReal" + cnt); zeroImag[cnt].setName("zeroImag" + cnt); zeroReal[cnt].addTextListener(textListener); zeroImag[cnt].addTextListener(textListener); }//end for loop //Now create an auxiliary display of the zplane // showing a unit circle. The user locates poles and // zeros on it by first selecting the radio button that // specifies a particular pole or zero, and then // clicking in the zplane with the mouse. (Note that // locating a pole outside the unit circle should // result in an unstable recursive filter with an // output that continues to grow with time.) refToZPlane = new ZPlane(); //Register an anonymous MouseListerer object on the // zplane. refToZPlane.addMouseListener( new MouseAdapter(){ public void mousePressed(MouseEvent e){ //Get and save the coordinates of the mouse click // relative to an orgin that has been translated // from the upperleft corner to a point near the // center of the frame. //Change the sign on the vertical coordinate to // cause the result to match our expectation of // positive vertical values going up the screen // instead of going down the screen. int realCoor = e.getX()  refToZPlane. translateOffsetHoriz; int imagCoor = (e.getY()  refToZPlane. translateOffsetVert); //The new coordinate values are deposited in the // real and imaginary text fields associated with // the selected radio button. //Examine the radio buttons to identify the pair // of real and imaginary text fields into which // the new coordinate values should be deposited. //Note that one radio button is always selected, // so don't click in the zplane unless you // really do want to modify the coordinate values // in the text fields associated with the // selected radio button. //Note that the integer coordinate values are // converted to fractional coordinate values by // dividing the integer coordinate values by the // radius (in pixels) of the unit circle as // displayed on the zplane. //Examine the pole buttons first. boolean selectedFlag = false; for(int cnt = 0;cnt < numberPoles/2;cnt++){ if(poleRadioButtons[cnt].isSelected()){ poleReal[cnt].setText("" + (realCoor/ (double)(refToZPlane.unitCircleRadius))); poleImag[cnt].setText("" + abs(imagCoor/ (double)(refToZPlane.unitCircleRadius))); //Set the selectedFlag to prevent the zero // radio buttons from being examined. selectedFlag = true; //No other button can be selected. They are // mutually exclusive. break; }//end if }//end for loop if(!selectedFlag){//Skip if selectedFlag is true. //Examine the zero buttons for(int cnt = 0;cnt < numberZeros/2;cnt++){ if(zeroRadioButtons[cnt].isSelected()){ zeroReal[cnt].setText("" + (realCoor/ (double)(refToZPlane.unitCircleRadius))); zeroImag[cnt].setText("" + abs(imagCoor/ (double)(refToZPlane.unitCircleRadius))); break;//No other button can be selected }//end if }//end for loop }//end if on selectedFlag //Cause the display of the zplane to be // repainted showing the new location for the // pole or zero. refToZPlane.repaint(); }//end mousePressed }//end new class );//end addMouseListener //Set the size and location of the InputGUI (JFrame) // object on the screen. Position it in the upper // right corner of a 1024x768 screen. guiFrame.setBounds(1024472,0,472,400); //Cause two displays to become visible. Prevent the // user from resizing them. guiFrame.setResizable(false); guiFrame.setVisible(true); refToZPlane.setResizable(false); refToZPlane.setVisible(true); //Save the reference to this GUI object so that it can // be recovered later after a new instance of the // Dsp046 class is instantiated. InputGUI.refToObj = this; }//end constructor //// //This is a utility method used to capture the latest // text field data, convert it into type double, and // store the numeric values into arrays. //The trycatch handlers are designed to deal with the // possibility that a text field contains a text value // that cannot be converted to a double value when the // method is invoked. In that case, the value is // replaced by 0 and an error message is displayed on // the commandline screen. This is likely to happen, // for example, if the user deletes the contents of a // text field in preparation for entering a new value. // In that case, the TextListener will invoke this // method in an attempt to compute and display new // values for the angle and the length. //One way to enter a new value in the text field is to // highlight the old value before starting to type the // new value. Although this isn't ideal, it is the best // that I could come up with in order to cause the angle // and the length to be automatically computed and // displayed each time a new value is entered into the // text field. void captureTextFieldData(){ //Encapsulate the pole data in an array object. for(int cnt = 0;cnt < numberPoles/2;cnt++){ try{ poleRealData[cnt] = Double.parseDouble(poleReal[cnt].getText()); }catch(NumberFormatException e){ //The text in the text field could not be converted // to type double. poleReal[cnt].setText("0"); System.out.println("Warning: Illegal entry for " + "poleReal[" + cnt + "], " + e.getMessage()); }//end catch try{ poleImagData[cnt] = Double.parseDouble(poleImag[cnt].getText()); }catch(NumberFormatException e){ poleImag[cnt].setText("0"); System.out.println("Warning: Illegal entry for " + "poleImag[" + cnt + "], " + e.getMessage()); }//end catch }//end for loop //Encapsulate the zero data in an array object. for(int cnt = 0;cnt < numberZeros/2;cnt++){ try{ zeroRealData[cnt] = Double.parseDouble(zeroReal[cnt].getText()); }catch(NumberFormatException e){ zeroReal[cnt].setText("0"); System.out.println("Warning: Illegal entry for " + "zeroReal[" + cnt + "], " + e.getMessage()); }//end catch try{ zeroImagData[cnt] = Double.parseDouble(zeroImag[cnt].getText()); }catch(NumberFormatException e){ zeroImag[cnt].setText("0"); System.out.println("Warning: Illegal entry for " + "zeroImag[" + cnt + "], " + e.getMessage()); }//end catch }//end for loop }//end captureTextFieldData method //// //=======================================================// //This is an inner text listener class. A registered // object of the class is notified whenever the text value // for any of the real or imaginary values in the pole and // zero text fields changes for any reason. When the // event handler is notified, it computes and displays the // angle specified by the ratio of the imaginary part to // the real part. All angles are expressed in degrees // between 0 and 359.9 inclusive. //Also, when notified, the event handler computes the // length of an imaginary vector connecting the pole or // zero to the origin as the square root of the sum of the // squares of the real and imaginary parts. //Finally, the event handler also causes the zplane to be // repainted to display the new location for the pole or // zero represented by the text field that fired the event. class MyTextListener implements TextListener{ public void textValueChanged(TextEvent e){ //Invoke the captureTextFieldData method to cause the // latest values in the text fields to be converted to // numeric double values and stored in arrays. captureTextFieldData(); //Identify the text field that fired the event and // respond appropriately. Cause the radio button // associated with the modified text field to become // selected. boolean firingObjFound = false; String name = ((Component)e.getSource()).getName(); for(int cnt = 0;cnt < numberPoles/2;cnt++){ if((name.equals("poleReal" + cnt))  (name.equals("poleImag" + cnt))){ //Compute and set the angle to the pole. poleAngle[cnt].setText(getAngle( poleRealData[cnt],poleImagData[cnt])); //Compute and set the length of an imaginary vector // connecting the pole to the origin. poleLength[cnt].setText("" + sqrt( poleRealData[cnt]*poleRealData[cnt] + poleImagData[cnt]*poleImagData[cnt])); //Select the radio button. poleRadioButtons[cnt].setSelected(true); firingObjFound = true;//Avoid testing zeros break; }//end if }//end for loop if(!firingObjFound){ for(int cnt = 0;cnt < numberZeros/2;cnt++){ if((name.equals("zeroReal" + cnt))  (name.equals("zeroImag" + cnt))){ //Compute and set the angle to the zero. zeroAngle[cnt].setText(getAngle( zeroRealData[cnt],zeroImagData[cnt])); //Compute and set the length of an imaginary // vector connecting the zero to the origin. zeroLength[cnt].setText("" + sqrt( zeroRealData[cnt]*zeroRealData[cnt] + zeroImagData[cnt]*zeroImagData[cnt])); //Select the radio button. zeroRadioButtons[cnt].setSelected(true); break; }//end if }//end for loop }//end if on firingObjFound //Repaint the zplane to show the new pole or zero // location. refToZPlane.repaint(); }//end textValueChanged //// //This method returns the angle in degrees indicated by // the incoming real and imaginary values in the range // from 0 to 359.9 degrees. String getAngle(double realVal,double imagVal){ String result = ""; //Avoid division by 0 if((realVal == 0.0) && (imagVal >= 0.0)){ result = "" + 90; }else if((realVal == 0.0) && (imagVal < 0.0)){ result = "" + 270; }else{ //Compute the angle in radians. double angle = atan(imagVal/realVal); //Adjust for negative coordinates if((realVal < 0) && (imagVal == 0.0)){ angle = PI; }else if((realVal < 0) && (imagVal > 0)){ angle = PI + angle; }else if((realVal < 0) && (imagVal < 0)){ angle = PI + angle; }else if((realVal > 0) && (imagVal < 0)){ angle = 2*PI + angle; }//end else //Convert from radians to degrees angle = angle*180/PI; //Convert the angle from double to String. String temp1 = "" + angle; if(temp1.length() >= 5){ result = temp1.substring(0,5); }else{ result = temp1; }//end else }//end else return result; }//end getAngle //// }//end inner class MyTextListener //=======================================================// //This is an inner class. An object of this class is // an auxiliary display that represents the zplane. class ZPlane extends Frame{ Insets insets; int totalWidth; int totalHeight; //Set the size of the display area in pixels. int workingWidth = 464; int workingHeight = 464; int unitCircleRadius = 200;//Radius in pixels. int translateOffsetHoriz; int translateOffsetVert; ZPlane(){//constructor //Get the size of the borders and the banner. Set the // overall size to accommodate them and still provide // a display area whose size is specified by // workingWidth and workingHeight. //Make the frame visible long enough to get the values // of the insets. setVisible(true); insets = getInsets(); setVisible(false); totalWidth = workingWidth + insets.left + insets.right; totalHeight = workingHeight + insets.top + insets.bottom; setTitle("Copyright 2006, R.G.Baldwin"); //Set the size of the new Frame object so that it will // have a working area that is specified by // workingWidth and workingHeight. Elsewhere in the // program, the resizable property of the Frame is set // to false so that the user cannot modify the size. //Locate the Frame object in the uppercenter of the // screen setBounds(408,0,totalWidth,totalHeight); //Move the origin to the center of the working area. // Note, however, that the direction of positivey is // down the screen. This will be compensated for // elsewhere in the program. translateOffsetHoriz = workingWidth/2 + insets.left; translateOffsetVert = workingHeight/2 + insets.top; //Register a window listener that can be used to // terminate the program by clicking the Xbutton in // the upper right corner of the frame. addWindowListener( new WindowAdapter(){ public void windowClosing(WindowEvent e){ System.exit(0);//terminate the program }//end windowClosing }//end class def );//end addWindowListener }//end constructor //// //This overridden paint method is used to repaint the // ZPlane object when the program requests a repaint on // the object. public void paint(Graphics g){ //Translate the origin to the center of the working // area in the frame. g.translate(translateOffsetHoriz,translateOffsetVert); //Draw a round oval to represent the unit circle in the // zplane. g.drawOval(unitCircleRadius,unitCircleRadius, 2*unitCircleRadius,2*unitCircleRadius); //Draw horizontal and vertical axes at the new origin. g.drawLine(workingWidth/2,0,workingWidth/2,0); g.drawLine(0,workingHeight/2,0,workingHeight/2); //Invoke the captureTextFieldData method to cause the // latest data in the text fields to be converted to // double numeric values and stored in arrays. captureTextFieldData(); //Draw the poles in red using the data retrieved from // the text fields. Note that the unitCircleRadius in // pixels is used to convert the locations of the // poles from double values to screen pixels. g.setColor(Color.RED); for(int cnt = 0;cnt < poleRealData.length;cnt++){ //Draw the conjugate pair of poles as small red // squares centered on the poles. int realInt = (int)(poleRealData[cnt]*unitCircleRadius); int imagInt = (int)(poleImagData[cnt]*unitCircleRadius); //Draw the pair of conjugate poles. g.drawRect(realInt2,imagInt2,4,4); g.drawRect(realInt2,imagInt2,4,4); }//end for loop //Draw the zeros in black using the data retrieved from // the text fields. Note that the unitCircleRadius in // pixels is used to convert the locations of the // zeros from double values to screen pixels. g.setColor(Color.BLACK);//Restore color to black for(int cnt = 0;cnt < zeroRealData.length;cnt++){ //Draw the conjugate pair of zeros as small black // circles centered on the zeros. int realInt = (int)(zeroRealData[cnt]*unitCircleRadius); int imagInt = (int)(zeroImagData[cnt]*unitCircleRadius); //Draw the pair of conjugate zeros. g.drawOval(realInt2,imagInt2,4,4); g.drawOval(realInt2,imagInt2,4,4); }//end for loop }//end overridden paint method }//end inner class ZPlane //=======================================================// }//end class InputGUI 
Copyright 2006, 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 hightech companies located in and
around Austin, Texas. He is the author of Baldwin’s Programming
Tutorials, which have 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 realworld
problems.