AN APPLICATION OF ADAPTIVE LINE ENHANCERS TO
OPTICAL PARTICLE COUNTER TECHNOLOGY
by
Edward R. Green
B.S., University of Colorado at Boulder, 1979
A thesis submitted to the
Faculty of the Graduate School of the
University of Colorado at Denver
in partial fulfillment
of the requirements for the degree of
Master of Science
Electrical Engineering
1995
This thesis for the Master of Science
degree by
Edward R. Green
has been approved
by
Af/*j/gr
Date
Boris Stilman
Green, Edward R. (M.S., Electrical Engineering)
An Application of Adaptive Line Enhancers to Optical Particle Counter Technology
Thesis directed by Assistant Professor Miloje S. Radenkovic
ABSTRACT
Adaptive filtering techniques are applied to the task of removing narrowband
interference from a desired signal when a reference correlated with the noise is not available.
The specific structure examined is called the adaptive line enhancer. Various implementations
are explored, including finiteimpulseresponse and infiniteinpulseresponse topologies.
Specific examples are drawn from the field of optical particle counting. Results indicate that
adaptive line enhancers are able to effectively remove narrowband interference from the
narrow pulses arising from particles passing through a laser beam.
This abstract accurately represents the content of the candidates thesis. I recommend its
publication.
Signed
Miloje S. Radenkovic
111
CONTENTS
CHAPTER
1. Introduction............................................ 1
2. Optical Particle Counting............................... 3
2.1 Baseline Shifting..................................... 4
2.2 Ambient light interference ........................... 7
2.3 Vibration............................................. 8
3. Optimal Noise Cancelers.......... ...................... 9
3.1 Meansquared solution................................... 11
3.2 Least Squares Solution ................................ 15
4. Leastmeansquare Adaptive Line Enhancers............... 18
4.1 LMS Algorithm......................................... 18
4.2 The Normalized LMS Algorithm.......................... 23
4.3 An example using the normalized LMS algorithm......... 24
4.4 Coefficient leakage .................................. 30
4.5 LMS normalized by recursive power estimation.......... 32
4.5.1 An example of the powernormalized LMS algorithm ... 33
4.6 An accelerated LMS algorithm.......................... 33
4.7 An example which includes noise................... 34
5.0 Recursive Leastsquare Adaptive Noise Cancelers....... 37
5.1 An example using the RLS algorithm.................... 39
5.2 An example that includes noise ....................... 41
IV
5.3 Exponential weighting..................................... 42
5.3.1 An example of the exponentially weighted RLS algorithm 43
5.4 Comparison to the Accelerated LMS Algorithm............. 47
6. Infinite impulse response adaptive noise cancelers.......... 49
6.1 Basic theory............................................ 49
6.2 Advantages and disadvantages of the HR implementation . 52
6.3 An example using the IIR prediction filter in an ALE .... 53
7.0 Conclusions............................................... 57
Appendix
A. Matlab routines used in this thesis......................... 60
References..................................................... 67
v
1.0 Introduction
A problem commonly found in many signal processing systems is that of
removing narrowband periodic interference from a desired signal. Many methods exist
for eliminating narrowband interference, but most require detailed knowledge about the
desired signal and the corrupting noise. In many applications such knowledge is not
available.
One obvious approach to removing periodic noise is to build a notchfilter at the
frequency, or frequencies, of the interference. A simple nonadaptive notchfilter with a
wide enough stopband to accommodate variations in the nature of the interference from
unit to unit could cause severe degradation in the quality of the signal. For proper
operation, the notchfilter must be precisely tuned to the exact frequency of the
interference. These filters must often posses a very high Q, which requires a complex,
expensive implementation. In addition, many more filters than are actually required for
a given implementation may be needed in order that the system be capable of handling
unforseen conditions.
Traditional techniques for adjusting the location of the stopband in the notch filter
require the difficult task of dynamically measuring the statistics of the interference and
modifying the filter.
1
Other methods of dealing with this type of interference, such as nonlinear polarity
sensing filters (baseline restorers) suffer from problems such as sensitivity to noise and the
lack of the ability to adapt to changing environments. They are also of limited value when
the frequency of the interference is near a frequency of interest in the desired signal.
Techniques for building Adaptive noise cancelers to address this problem have
been proposed by many authors. Glover [Glover, 1977] proposed a configuration for
adaptively eliminating sinusoidal interferences when a signal correlated with the
interference is available. This system is similar to the systems described in this thesis,
except that his solution required a separate reference signal. He was able to show that
when the reference consists of a summation of sinusoids the system converges to a notch
filter, with notches at each of the frequencies contained in the reference.
This thesis describes the theory associated with adaptive line enhancers(ALE), and
presents the results of applying ALE techniques to a specific problem, separating random
gaussian pulses from periodic interference. Specific examples are taken from the field of
optical particle detection where elimination of periodic interference is a serious problem.
2
2.0 Optical Particle Counting
Detecting particles suspended in gases and liquids is an important problem in many
industries. These particles can reduce the quality of many products. Optical particle
counters can be used to detect contamination in semiconductor processes where the
contamination can dramatically reduce yields; they can be used to detect deadly micro
organisms such as giardia in treated drinking water and they can be used to monitor the
quality of pharmaceutical products.
Optical Particle counters detect and quantify the size of particles suspended in air
and liquids by collecting the light scattered from the particles as they pass through a
narrow laser beam. Detecting the small pulses of light resulting from a submicron
particle passing through a narrow laserbeam presents many signal processing problems.
These pulses tend to be approximately gaussian in shape because the intensity profile of
the laser beam is approximately gaussian, and because the bandwidth limitations on the
processing electronics tend to widen the pulses. The pulses vary from a few hundred
nanoseconds to hundreds of microseconds in duration.
Detecting the short, low intensity light pulses from submicron particles passing
through a narrow beam requires the elimination of several sources of periodic
interference. One large problem is associated with ambient light. These particle counters
need to have unrestricted paths from the external environment into the optical chambers
of the benches because any bends or obstructions can remove the particles being
monitored. This direct path allows ambient room light to enter the detection volume.
3
Since many sources of room light have a periodic nature, e.g. florescent lights,
there is a large opportunity for low frequency periodic interference. Figure 2.0.0 shows
an oscilloscope trace of this type of noise. Also shown in this figure are several particle
induced pulses.
Tek fciftim Single Seq 20kS/s
Fig 2.0.0 Periodic noise from florescent light
2.1 Baseline Shifting
One particularly difficult problem is associated with baseline shifting." The
detectors used to collect the scattered light from a particle and convert it to an electrical
4
pulse have aD.C. component associated with darkcurrent or with light coming from the
direct laserbeam in the case of extinction instruments, which is most easily removed by
A.C. coupling the signal. For low pulse rates, this technique works well, but when the
pulse rate is high the A.C. coupling removes the D.C. component associated with the
particle pulses, which results in a D.C. shift of the entire waveform in a direction opposite
to the direction of the pulses. Various methods have been used to restore the low
frequency components to the signal, and thereby compensate for the undesired shift.
Figure 2.1.0 shows an example of a single particle pulse where the baseline of the signal
has been shifted below ground by the A.C. coupling in the signal path.
Tek gBiES Single Seq 2MS/S
5
One method is to detect the negative peak of the signal (assuming positive going
pulses) and to subtract this from the original signal to place the baseline of the signal at
ground. This type of baseline restorer has several problems. One is that it is very sensitive
to negative noise spikes, and the other is that the entire waveform shifts when the
background noise changes in amplitude. A block diagram of such a system is shown in
figure 2.1.1.
Fig 2.1.1
Another method is to switch between two highpass filters depending on whether
a pulse has been detected or not. When the system detects that there is no particle event
currently being detected, a filter tailored to remove as much low frequency signal as
possible but still detect a pulse is switched into the system. When a particle is detected,
this filter is replaced with a filter with a high cutoff frequency to allow the particle signal
to be processed. A diagram of this type of baseline restorer is shown in figure 2.1.2.
6
Fig 2.1.2
2.2 Ambient light interference
Optical particle detectors are often affected by interference from ambient light
sources. Because the laser source used in most modem particle detectors has a
wavelength is in the near infrared (780 nm) portion of the spectrum, it can be very
difficult to design an optical system that is immune to ambient light interference. Many
of the coatings traditionally used to prevent reflections from mechanical parts, e.g. black
anodizing, are reflective to infrared signals making it difficult to design effective light
traps. Making the design of light traps more difficult is the fact that particles bigger than
a few microns are easily removed by bends or obstructions in the flow path. These
considerations make it desirable to be able to detect particles in the presence of some level
of ambient light.
When the light is continuous in intensity, e.g. sun light, the A.C. coupling
described above can be effective in removing the interference. Florescent light, which is
7
common in the environments where these instruments are typically used, poses a more
difficult problem. These Florescent lights switch on and off at twice the line frequency,
which results in an A.C. noise source in the system. Unfortunately, the frequency of the
noise, 120 Hz in the U.S., is located at a frequency where there is substantial power in the
gaussian pulses generated by a particle. Traditional notch filters have proven to not be
a very satisfactory solution because a filter with a broad enough stop band to eliminate
the noise would seriously degrade the quality of the pulses.
2.3 Vibration
Because it is very difficult to eliminate all sources of reflection in the optical
portion of a particle counter, any vibration will cause the light to be scattered in a
different way causes an A.C. signal to enter the detector. This signal is impossible to
predict since the source of the vibration is unknown, but the frequency associated with
such a signal is low, typically less than a few kilohertz. Traditional solutions to this
problem are to painstakingly remove all sources of scattered light in the optical system,
and to try to mechanically isolate the instrument from external sources of vibration.
These solutions are not very desirable because they involve a lot of expensive labor, and
are only partially effective.
8
3.0 Optimal Noise Cancelers
A noise canceler is a signal processing system that accepts one or more inputs
containing a desired signal contaminated with noise, and one or more reference inputs that
are strongly correlated with the noise, but not the signal, and attempts to produce a noise
free version of the signal. In the simplest case there is one input containing the noise
corrupted signal, and one input correlated with the noise, a block diagram of the noise
canceler is shown in figure 3.0.0. It should be noted that this system produces estimates
of both the signal and the noise which makes it possible to use the system when the
reference is correlated with the signal rather than the noise.
Fig 3.0.0
In figure 3.0.0, the x(7)=s(X)+h(0input, where s(t) is a desired signal, and n(t) is
corrupting noise. Input n(t) is a signal which is correlated with n(t), but has little or no
9
correlation with s(t). The noise canceller attempts to extract s(t) by producing, an
estimate of n(t) from n'(t) and subtracting this estimate from x(t). If the statistics of n(t)
and n'(t) are known in advance, a noise canceler can be designed using the results of
Optimal Filtering Theory. However, if the statistics are not well known, or are time
varying, subtracting one noisy signal from another can reduce the signaltonoise ratio at
the output rather than enhance it.
The central component to a noise canceler is the prediction filter. This filter
produces an estimate of the noise component of x(t) which is then used to extract s(t).
This prediction filter can be designed using the results of optimal filter theory.
A filter is considered optimal when it minimizes some predetermined cost function
or performance index. If the cost function is selected correctly, the result is a filter with
good performance, but not necessarily the best or optimal performance in the usual
meanings of those terms. Most of the work with digital optimal filters uses a finite
impulse response (FIR) topology for the filter because such filters are inherently stable
since they contain only forward signal paths. Assuming a FIR structure, the cost function
can then be specified. There are many possible choices for this function, some of which
are [Orfanidis, 1988, p. 158]:
1. Meansquare value of the estimation error
2. Summed squared error
3. Expectation of the absolute value of the estimation error.
10
Of these possible choices for the cost function, which is usually denoted J, the first
two are the most attractive because they lead to tractable mathematics as well as filters
that tend to perform well. These two methods differ in the way that they approach the
problem. The summedsquareerror method leads to a filter that is optimal for some
specific realization of n(t) and n'(t), while when the performance function used is the
meansquarederror the filter adapts to the optimal value for the ensemble average of the
possible realizations for n(t) and n'(t).
3.1 Meansquared solution
The meansquare value of the estimation error is an attractive choice for the cost
function because it leads to manageable mathematics. Specifically it results in a cost
function which has a unique minimum, and which is a simple second order function of the
filter coefficients [Haykin, .1986, p. 163],
The FIR prediction filter in figure 3.0.0 produces an output n(t)which is
calculated from the input n'(t) using the following convolution equation, assuming that
all quantities are real:
L1
n(t) = ^WfpXtk) (3.1.0)
k=0
where wk is the kth filter coefficient, n'(tk) is the delayed filter input presented to the kth
element in the filter and L is the length of the FIR filter. The estimation error is
e(t)=n(f)n(t), and the performance index to be minimized is defined as,
11
Jr=Â£[e(?)2] = Â£[((?) n(t)){n(t) ~n(t))]. This performance function is a quadratic in
L dimensions which, appropriate signals, n'(t) and n(t), has a welldefined minimum. The
optimal values for the filter weights can be determined by solving for the minimum.
Define the gradient operator as:
V,=
dw,
Â£=0,1,2,...
(3.1.1)
then the filter coefficients that minimize the cost function can be found by applying the
gradient operator to the cost function and setting the result to zero,
V^=0 (3.1.2)
After simplification the result is
W = 2 ^n'{tk)eo(t)
= 0
(3.1.3)
where e0 is the error that results from the filter operating at its optimum condition.
This shows that a necessary and sufficient condition for the cost function J to attain its
minimum value is that the estimation error is orthogonal to each input sample that enters
into the calculation of the output at time t.
A corollary to this principle can be seen by calculating Â£[ 7(f)e(r)],
12
L1
E[n(i)e(t)]=E Y wkn Xtk)e(t)
ko
L1 1,1
=YwiAnXttywhYwk=0 (314)
k=0
L1
r
r
Equation (3.1.4) states the output of an optimal filter is uncorrelated with the error.
The above principle of orthogonality states the necessary and sufficient conditions
for optimality, which can be rewritten as:
Ll
E n'{tk)
n{t)~Y w,n Xt~i)
i=0
=0, fc=0,l,2,...,Ll
(3.1.5)
which, after distributing the expectation operator becomes,
Ll
Y/wi^pXt~k)riXt~o} = E^Xtk^t)] (3.i.6)
7=0
The expectation on the left side of equation (3.1.6) is the autocorrelation
function of the input signal for a lag of I k, which is denoted r(Ik). The expectation on
the right side of (3.1.6) is the crosscorrelation between the input signal, n'(tk), and the
desired signal, n(t), for a lag of k. This crosscorrelation is usually denotedp(k). So the
necessary and sufficient conditions for optimality can be written:
13
(3.1.7)
L1
Y,wir{ik) = P(~k\ k0,l,2,...JLl
which are often called the WienerHopf equations. These equations show that when the
autocorrelation function r() and the crosscorrelation function p() are known, and are
time invariant, the filter coefficients for an optimal noise canceler are defined by a system
of linear equations. The solution to equations (3.1.7) is unique when the autocorrelation
matrix is fullrank. Unfortunately, these functions are rarely known a priori, and are
usually changing with time. Since the auto and crosscorrelation functions used are the
ensemble average for the input sequences, the optimal solution for the filter coefficients
found using this performance criterion is optimal for the ensemble average of n'(t) and
n(t).
The WienerHopf equations can be written in matrix notation as:
(3.1.8)
where r.. Zsjft Xki)n /(&/)] = r(ji) = r{i j) is the ijth entry in and
pi = Efyi(k)n Xki)\ is the ith entry in Pms. When /?msis invertible, then the optimum
value for the vector of filter coefficients is given by:
(3.1.9)
14
3.2 Least Squares Solution
While the above meansquareerror method produces a method that is easy to
analyze theoretically, in application the ensemble auto and crosscorrelation functions are
usually not known. They must be approximated by examining a sequence of samples
taken from a single realization.
This fact leads to a second common selection for the performance index, the
sumsquared error measure. Using this index, the solution of the optimal filtering problem
leads to a leastsquares solution. If the number of taps in the FIR structure is L, and the
input sequences are available for 7=0..JV1, then the performance index can be written
as:
J = E [K0)(0)2
t=L1
(3.2.0)
Expanding equation (3.2.0) and using the fact that n(t) = WWit), and assuming all
quantities are real, yields:
N1 AM AM
J = X) K(0I2 2E W^ityn'if) + Y, WtN/(t)Nll{t)W(3.2.1)
t=L\ tL\ t=L1
If the following identifications are made,
15
^E HOI2
t=Ll
Nl
ps,= E
(3.2.2)
t=L1
Nl
E iv'(0JV"(0
/=L1
Equation (3.2.1) becomes
J = D 2WPt + WRW
ss ss ss ss
(3.2.3)
The value of Wthat minimizes the performance index given in (3.2.3) can be found using
matrix calculus by satisfying the following two necessary and sufficient conditions
[Treichler, 1987, p.34 ]:
1. The gradient of J with respect to W is zero.
2. Hyy, the Hessian matrix of J with respect to W is positive definite.
It can be shown that the value of Wthat satisfies the first requirement is:
And the second requirement is satisfied if Rss is positive definite. Also, if Rss is positive
definite, it is invertible.
The form of (3.2.4) is obviously very similar to (3.1.9), The difference being the
interpretations of P and R For the meansquare case the auto and crosscorrelation
(3.2.4)
16
matrices are ensemble averages over the possible realizations of n'(t) and n(t). In the
leastsquare (or sumsquare) case, the values for the correlation matrices are the time
averaged values for one specific realization of n'(t) and n(t).
17
4.0 Leastmeansquare Adaptive Line Enhancers
The optimal noise cancelers described in the previous chapter can be used as a
basis for the development of adaptive systems. The performance indices used to design
the optimal filters provide an ideal mechanism for adapting the filter parameters so as to
produce a filter with good performance. The performance indices are second order
surfaces in Zspace, where Z is the order of the FIR filter. These functions have a well
defined minimum which can be found iteratively. Some of the most common adaptive
systems use gradient search algorithms to find this minimum.
Several authors, e.g. [Zeidler, 1978] and [Glover, 1977], have proposed an
adaptive system based on the Leastmeansquare (LMS) algorithm for detecting and
enhancing periodic signals in the presence of aperiodic noise. These systems produce two
outputs, one which is the predicted narrowband signal and one which is the predicted
broadband signal. Usually the output of interest is the narrowband signal. This thesis
describes the use of these systems to extract the broadband, or aperiodic signal.
4.1 LMS Algorithm
The Least Mean Square (LMS) algorithm as described by Widrow, et al.
[Widrow, 1975] is a method which seeks the minimum of the mean squared error
performance index by iteratively estimating the gradient of the performance index, and
moving a small distance toward the minimum. The advantage of the LMS algorithm is
that it requires relatively little computation at each step, and can be shown to converge
to a solution near the optimal solution when certain requirements are satisfied.
18
(
The key observation that makes the LMS algorithm possible is that the gradient
of J, the performance index, can be estimated from just the error at the output of the
prediction filter, and the input to the filter. The instantaneous value of the performance
index is e 2(t). Taking the gradient with respect to W, the instantaneous impulseresponse
vector for the FIR filter gives,
G(i) = Vje2] = 2e(0V(f)
= 2e(t)vJn(t)WnXt)] (4.1.0)
= MtyvJjr'N'] = 2
The above instantaneous gradient function can be time averaged to yield the gradient of
the sumsquared performance index, or it can be ensemble averaged to get the mean
squared index. The above gradient can be used to move toward the optimal value of W
by the following recursion equations,
W(t + 1) = W(t) + \ie(t)N\t) (4.1.1)
where p is a constant that controls the rate of convergence of the algorithm. If p is
chosen too large, the algorithm will become unstable; if it is chosen too small the
algorithm will take a very long time to converge. An upper bound on p can be calculated
by substituting the analytic expression for the gradient, 2P + 2RW(t), for eft) in
equation (4.1.1) giving,
19
W(t + !) = (/ \iR)W(t) + (IP
(4.1.2)
To study the convergence properties of equation (4.1.2) is useful to apply a
similarity transformation to convert the equation into L uncoupled scalar equations. If Q
is defined as the L by L matrix whose columns are the eigenvalues of R, then the
following identifications can be made,
W = QW'(t\
R = QAQ, (4.1.3)
pf =Qtp
where A is a diagonal matrix containing the eigenvalues of R. Then equation (4.1.2) can
be written as,
W'(t + 1) = (I \iA)W'(t) + 1LP (4.1.4)
which are L scalar equations of the form,
wi(t + 1) = (1  + VPi (4.1.5)
[I A.  < 1. So a condition
(4.1.6)
These difference equations converge to a solution if 11
for convergence of equations (4.1.5) is,
0 < i <
max
20
The time constant associated with each of the uncoupled modes is given by,
1
T.
(4.1.7)
If we combine the limitations on p given by equation (4.1.6) with the expression for the
time constant given by equation (4.1.7), and defining a normalized step size given by
T
(4.1.8)
The convergence of the LMS algorithm is controlled by the spread between the
minimum and maximum eigenvalues of R. This is proportional to the condition number
of the matrix under the 2norm, which shows that when the autocorrelation matrix is ill
conditioned, the filter will take a long time to converge.
Since the original filter coefficients, Wj are linear combinations of the uncoupled
modes, they will converge only as fast as the slowest of the uncoupled modes. Therefore,
the time .constant for convergence of the filter coefficients is bounded from below by
It is possible to examine the behavior of the prediction filter when the auto
correlation matrix R is singular. If R is singular, it can be shown that at least one of the
eigenvalues of R is zero. Then the recursion equation for a transformed mode associated
1
21
with a zero eigenvalue becomes,
wft + 1) = (1 0) wft) + 0 = wfi + 1) (4.1.9)
because the crosscorrelation term associated with this mode is also zero. This shows that
the modes associated with the zero eigenvalues are undriven and undamped. Therefore,
the filter coefficient vector W will never converge.
When the matrix R is singular there are infinitely many values of W which satisfy
the normal equations. This can be shown by the following calculation. If V is an
eigenvector associated with a zero eigenvalue then RU = 0. Then substitute into the
normal equation W+yL/where Wis a solution to the normal equation. This gives:
R(Wq + y U) = RW* + yRU = P + 0 = P (4.1.10)
This shows that any linear combination of the eigenvectors associated with zero
eigenvalues can be added to a solution of the normal equation and the normal equation
will still be satisfied. Another implication of equation (4.1.10) is that the failure of the
undriven modes to converge does not prevent finding a solution that minimizes the
performance criterion.
The rate of convergence when R is singular is given by an equation similar to
equation (4.1.8) where the minimum eigenvalue is replaced by the minimum nonzero
eigenvalue. Specifically the time constant is given by[Treichler, 1987, p.81],
22
X
(4.U1)
4.2 The Normalized LMS Algorithm
One problem with the LMS algorithm is that the adaptation constant can become
too large for the inputs to the filter. When this happens, the filter becomes unstable.
Equation (4.1.6) shows the range of values allowed for the adaption constant. The
normalized LMS algorithm replaces the adaptation constant with a value that is
normalized with an estimate of the largest eigenvalue of R. The new recursion equation
for the filter coefficients is,
(4.2.0)
where a is constrained to be between 0 and 2, and N/tN/ is an estimate of the largest
eigenvalue of R. Usually this equation is modified to insure that when is zero that
a division by zero is avoided. The modified equation is
W(t + 1) = W{t) +
a
e(f)N'(f)
y +N,tN/
(4.2.1)
where y is a small constant.
23
4.3 An example using the normalized LMS algorithm
An example of the normalized LMS algorithm that approximates the type of
problem encountered in particle detection is examined in this section. The narrowband
contamination is generated by combining two sinusoids of different frequencies.
The narrowband signal used is shown in figure 4.3.0.
Narrow Band Signal
2 0
1 5
1 0
5
0
 5
 1 0
 1 5
 2 0
0 50 100 150 200 250 300
Fig 4.3.0
The particle pulse is simulated by a narrow cosinesquared function centered at the
sample which is 3/4 of the maximum index in the input waveform (the length of the
waveform is variable since the rate of adaption is changed from one simulation to the
next, and the filter needs to have converged before the pulse is injected). An example of
the input waveform with the particle pulse is shown in figure 4.3.1.
24
Narrow band signal and pulse
The normalized LMS algorithm has three parameters which can be adjusted to
achieve the desired performance; the length of the filter, the adaptation constant a, and
the amount that the combined signal is delayed before being presented to the filter. Of
these three parameters one of the most important is the length of the prediction filter.
Figure 4.3.2 shows the effect of changing the length. In fig 4.3.2 the plots are of a small
segment of the broadband signal. In each case the plot shows the performance after the
filter has converged. From top to bottom are plots from filters with 20,40,80 and 160
elements in length.
25
In fig 4.3.2, the top plot clearly does not have acceptable performance, whereas the filter
length of the bottom plot is 160 which extracts the pulse with good fidelity. A plot of the
filter weights for the case with 160 tap points is shown in figure 4.3.3.
I________I________I_______I________I________I
6000 6200 6400 6600 6800 7000
Fig 4.3.2
26
Fig 4.3.3
The example signals used in this example produce an autocorrelation matrix, R,
which is poorly conditioned and has many zero eigenvalues. This has a negative effect
A
on the time of convergence for the adaptive system. A plot of the factor, ax, resulting
^min
from the input shown in figure 4.3.0 is shown in figure 4.3.4.
Figure 4.3.4 shows that the filter length can have a very significant effect on the
convergence time of the filter. Specifically there are filter lengths that have a much longer
convergence time than do other lengths. Moreover, the relationship between convergence
time and filter length is not simple.
27
Fig 4.3.4
Since there are many eigenvalues of the R matrix that are zero in this example,
there are many modes that are undriven. Since the weights that are associated with these
modes will not change, it is important to select an appropriate initial value for these
coefficients. In fig 4.3.2, the coefficients were all initialized to zero. To demonstrate the
importance of the initial values in this example, the simulation was rerun with the
coefficients initialized to one. The result of this simulation is shown in figure 4.3.5.
28
2Q 6000 6200 6400 6600____6Â§QQ___ZQ00
6000 6200 6400 6600 6800 7000
Fig 4.3.5
These results show that the coefficients of the prediction filter which dont converge
result in severe degradation in the output signal after the passage of the pulse. The
coefficients for the case where there are 160 taps in the FIR prediction filter are plotted
in fig 4.3.6.
29
Fig 4.3.6
4.4 Coefficient leakage
A problem with the LMS and normalizedLMS algorithms is that the input signal
may not excite all of the modes on the prediction filter, which results in an R matrix with
zero eigenvalues. This results in a filter that does not converge to a unique solution for
W. Several methods have been proposed to deal with this problem such as mixing a small
amount of random noise into the input signal. A better solution is the coefficient leakage
method proposed by Zahm [Zahm, 1973], This method modifies the algorithm by
inserting a term that causes the filter coefficients to slowly leak away to zero in the
absence of an update term from the error. Specifically the update equation for the filter
coefficients becomes,
Wit + 1) = (lpy)fF(O + (4.4.0)
or in the case of the normalized LMS algorithm,
30
a
(4.4.1)
e+N/TN'
W(t + 1) = (1 y\i)W(t) + \ie(t)W(t)
This modifies the R matrix to be R' = yl + R. Which adds the small factor y to each
eigenvalue. This results in a system that converges to a unique solution even when the
input does not excite all of the modes of the filter. The price of this convergence is that
the solution is not optimal. Specifically the coefficient vector converges to (R + yI)~lP
rather than R _1P. A plot of the results of running this modified algorithm under
conditions similar to those in the bottom plot in figure 4.3.5 is shown in figure 4.4.0.
Broadband alpha = 0.01, Len=160
6000 6100 6200 6300
Fig 4.4.0
This plot is very similar to the plot in figure 4.3.5. The addition of coefficient leakage
does not improve the recovery of the filter after the passage of the pulse. A plot of the
filter coefficients for this case is shown in figure 4.4.1.
31
1 .8
Fig 4.4.1
4.5 LMS normalized by recursive power estimation
In the normalized LMS algorithm, the adaption constant was normalized with an
upper bound on the largest eigenvalue of R. The computation of this normalization factor
requires N multiplications and N additions. The value of the normalization factor is given
byfN'^N'k which is an estimate of the power in the input vector at time k. It is possible
to estimate the value of the input power with a recursive equation that is more
computationally efficient. If the estimate of the input power at time k is denoted, n(k)
then a recursive estimate of the input power can be calculated as,
7t(*+l) = (1P)7T(*) + p n/2(k) (4.5.0)
where P is a constant used to control the averaging of the power estimate. With this
modification, the update equation for the filter weights is,
W(k+1) = W(k) + ae(k)N'(k) (4.5.1)
y + 7t(&)
32
if P is chosen as, P = ^ then the powernormalized LMS algorithm will have
performance similar to the normalized LMS algorithm [Treichler, p. 84].
4.5.1 An example of the powernormalized LMS algorithm
When the powernormalized algorithm is applied to the same type of signals used
in the previous examples, the plot shown in figure 4.5.1.0 results.
This plot shows that the powernormalized algorithm produces results very similar to the
results from the LMS algorithm normalized with the upper bound of the largest
Broadband mu =0.0001, Len=160
6000 6100 6200 6300
Fig 4.5.1.0
eigenvalue of R.
4.6 An accelerated LMS algorithm
The LMS algorithms described to this point all move a small distance in the
direction of steepest decent along the gradient of the performance index. While the
normalized algorithms will eventually converge near the optimal, they do not attempt to
33
move directly toward the optimal. The following algorithm describes a modification that
attempts to move directly toward the minimum. While this algorithm is of interest
theoretically, as will be shown in a future section, it is not generally of much practical
significance. This is because the modification requires knowledge about the statistics of
the inputs which are generally not available. If the update equation for the filter weight
vector is modified to include the additional matrix C in the following manor,
W(k+\) = W(k)+\ve(k)CN'{k) (4.6.0)
and C can be chosen to approximate the inverse of R, then the algorithm will have
improved convergence when the eigenvalue disparity of R is large [Treichler, p. 84],
Unfortunately, these adaptive algorithms are of most use when R is not known.
4.7 An example which includes noise
In the previous examples, the input did not include random noise. In any practical
situation, the input signal will contain some amount of random noise from various sources
such as Johnson noise from resistive elements and shotnoise from active elements. To
examine the behavior of the LMS algorithm in the presence of noise, the previous test
signal which includes a periodic signal and a single cosinesquared pulse is modified by
adding gaussian random noise. The results of applying the normalized LMS algorithm to
this noisy signal where the variance of the noise is 1 is shown in figure 4.7.0.
34
Figure 4.7.0 demonstrates that the ALE is able to detect a pulse even in the
presence of a moderate level of noise. Figure 4.7.1 shows the output of the same system
when the variance of the noise is increased to 4.
20
0
20
JL
ill
L
500 1000 1500 2000 2500 3000
500 1000 1500 2000 2500 3000
2200
2300
2400
2500
Fig 4.7.0 topinput, middleoutput, bottomexpanded output
35
20
500 1000 1500 2000 2500 3000
2200
2300
2400
2500
Fig 4.7.1
36
5.0 Recursive Leastsquare Adaptive Noise Cancelers
The gradient search algorithms which were described in chapter 4 form a class of
systems which find great use in application. Another class of algorithms is based on the
leastsquares techniques. In contrast to the LMS algorithms that attempt to move the
state of the system toward the optimal condition at each iteration, the recursive least
square techniques seek to produce an optimal solution, based on the past states of the
system, at each iteration. The RLS method accomplishes this solution by continually
estimating the value of R, and solving the WienerHopf equations. The direct solution
to this problem is very computationally expensive. Fortunately, there exist methods to
minimize the calculations involved in solving the WienerHopf equations.
The goal of producing a solution at each step that is optimal leads to a different
interpretation of the solution. As was pointed out in chapter 3, the LMS type algorithms
lead to a solution that is optimal for the ensemble of possible inputs having the statistical
properties observed in the input functions. The RLS methods produce a solution that is
optimal for the specific averaged statistical properties of the inputs presented to the filter.
The optimal solution for a linear prediction filter, as given by the WienerHopf
equations is:
W = R1 P (5.0.0)
The problem with equation (5.0.0) is that the values of R and P are usually not known.
One obvious technique is to estimate the values of R and P at time k by the
37
following equations:
Rt+1 = Rt + n'n'!
pt+1 = pt + am;
(5.0.1)
and then use the values to compute Wx = Rt+\ Pt+1. Unfortunately, this technique
requires the very computationally expensive action of inverting R at each step. The
solution to this problem is the matrix inversion lemma, which is also known as the
ABCD lemma or Woodburys identity. This important lemma is given in equation
(5.0.2).
(A + BCD) 1 = A1 A lB(DA 1B + C1)'DA_1 (5.0.2)
Making the following identifications:
A = Rt
B = JV/
C = I
D = N?
(5.0.3)
an iterative equation for the inverse of R can be written as,
i + n'/r^n'
(5.0.4)
Several interesting, and simplifying identifications can be made in equation (5.0.4).
38
Define the filtered information vector as:
R<
n;
(5.0.5)
and the a priori output as:
y0(t) N!T Wt
(5.0.6)
and the normalized input power as:
qr = NtTRtiN't = NJzt (5.0.7)
With the above identifications, the update equation for the optimal filter impulse response
vector can be written [Treichler, 1987, p. 95],
TI/0 m y0(t)] zt 
wt =  (5.0.8)
\ + q,
Or in the case of an adaptive line enhancer where d(t) is zero, equation (5.0.8) becomes,
Wt
o
[*(0) y0(t)] zt
1 + q,
(5.0.9)
5.1 An example using the RLS Algorithm
39
Applying the RLS algorithm to data similar to the data that produced the plots in
figure 4.3.2 yields the results shown in figure 5.1.0. From these results it is clear that the
RLS technique produces good results when the filter length is sufficiently long. One
noticeable difference between the results from the LMS and RLS algorithms is the nature
of the initial convergence to a solution. The LMS algorithm has a period when the filter
weights are exponentially converging to their steadystate values. The RLS method
mu =1e005, Len= 20
0
200 400 600
Figure 5.1.0
800 1000
40
produces a solution that is optimal much more rapidly.
5.2 An example that includes noise
In this example, the RLS algorithm is applied to a signal that includes noise in a
mannor similar to the case presented in section 4.7. The magnitude of the noise used in
this example is a2=l. The results of this simulation, shown in Fig 5.2.0, are very
interesting. Not only has the RLS filter successfully separated the pulse from the noise,
but the disturbance of the system after the passage of the pulse is much reduced as
compared to the RLS example without noise. Another notable feature of the results of
the simulation is that the convergence of the algorithm is substantially slower than it was
in the case with no noise.
input mu =1e005, Len= 160
Fig 5.2.0
41
In figure 5.2.1, the simulation was run with the magnitude of the noise doubled to o2= 4.
input mu =1e005, Len= 160
Fig 5.2.1
5.3 Exponential weighting
A problem associated with the RLS algorithm is that as more and more data are
available to the system, the influence of each additional sample is reduced. In a system
where the characteristics of the interference are stationary this Vanishing gain
phenomenon has no negative effects on the operation of the filter. However, when the
characteristics of the interference are changing with time, the filter will be unable to
respond to these changes. On solution to this problem is the exponential weighting
modification. This technique modifies the definitions of R and P to cause the
contributions of old data to be progressively reduced. The modified equations are,
42
(5.3.0)
/=V1
P = E (5.3.1)
Z=AM
where p is a factor included to reduce the effects of old data.
The only changes required to the algorithm given in equations 5.0.5 through 5.0.9
is that the following equations should be substituted for their equivalents in the above
referenced equations,
mzu.
p+q
(5.3.2)
Ki = w;+\
d(k)y0(k)
p+q
(5.3.3)
5.3.1 An example of the exponentially weighted RLS algorithm
The following plot shows the results of running the exponentiallyweighted
algorithm with p=l, which is equivalent to the normal RLS technique. The input consists
of four sinusoids with nearly the same frequency, each of which occupies a sequential
43
section of 500 samples. In addition, the input includes added gaussian random noise with
unity variance. The last one thousand samples consist of only the noise. It is obvious that
the prediction filter has increasing difficulty adapting to the input with each successive
segment. Also the variance of the noise in the final segment is increased from the variance
in the input, or even the first segment. This increase in the variance is caused by the
addition of uncorrelated noise at each sample because the filter is unable to adapt by
driving the filter coefficients to zero.
input rho =1, Len=50
Fig 5.3.1.0
44
A plot of the R matrix after the completion of the algorithm is shown in figure
5.3.1.1. The plot shows that there are diagonals in this matrix that are significantly non
zero. The results for random noise should be a matrix that is nonzero only on the main
diagonal because random noise is uncorrelated for lags greater than zero.
In figure 5.3.1.2 the same data are applied to the weighted RLS algorithm, but
with p=0.99. This value of p will result in a timeconstant for the forgetting factor of
approximately 100 samples.
45
input rho =0.99, Len=50
2011''
2011'11
0 500 1000 1500 2000 2500 3000
output
20......
20 u^1
0 500 1000 1500 2000 2500 3000
Fig 5.3.1.2
Clearly these results show that the forgetting factor has allowed the filter to
adapt to the changing characteristics of the data. Each of the five segments has about the
same ability to eliminate the periodic noise. In figure 5.3.1.3, a plot of the R matrix at the
end of the algorithm is shown.
46
Fig 5.3.1.3
This plot shows that the R matrix has successfully incorporated the new data and shows
the characteristics that would be expected for random noise, i.e., that the R matrix has
substantially nonzero values only on the main diagonal.
In addition to the ability to adapt to changing input statistics, the weighted
algorithm keeps the R and P matrices from becoming progressively larger as more and
more data are included. This has computational benefits when the algorithm is
implemented in hardware with limited dynamic range.
5.4 Comparison to the Accelerated LMS Algorithm
In section 4.6 an accelerated LMS algorithm was presented. It was observed that
if an estimate of the autocorrelation function was available, the convergence properties
of the LMS algorithm could be improved. If equation (4.6.0) which is repeated here for
convenience,
47
W(k+1) = W{k)+\ne(k)CN'(k)
(5.4.0)
is compared to the following update equation for the RLS algorithm, which can be
generated from equations (5.0.5) through (5.0.9),
K, = K * (5A 0
l+q, V >
it can be seen that the RLS algorithm is identical to the LMS algorithm if these
identifications are made,
1
V = 
l+(lt (5.4.2)
C = V
This shows that the RLS algorithm can be interpreted as a method for calculating at each
step the optimal values for p and C [Treichler, 1987, p. 95],
48
6.0 Infinite impulse adaptive noise cancelers
In chapter 5, the theory of finiteimpulseresponse filters as applied to adaptive
line enchantment was presented. Most of the work done with ALEs has been done with
FIR implementations because they are inherently stable, since all of the poles associated
with the prediction filter are located at the z=0 in the ztransform representation.
In contrast, HR filters have both zeros and poles that are adjusted as the adaption
algorithm proceeds. The problem with adjusting the pole locations is that the algorithm
may cause the pole locations to move outside the unit circle, which produces an unstable
filter.
The theory of HR predictors is very complex; the purpose of this section is to
provide an example of one alternative technique to the FIR techniques discussed in
chapter 5. Only one simple algorithm for implementing an HR system is presented.
6.1 Basic theory
The presentation in this section closely follows the material presented in
[Treichler, 1987, pp. 117121], The form of the prediction filter, using an autoregressive,
movingaverage (ARMA) model is:
l m
m = *,(*>>'(*./) (6.1.0)
/=1 j=0
Clearly, this equation has an infiniteimpulseresponse arising from the first summation
in equation (6.1.0).
49
Like the FIR systems, the adaption process relies on minimizing a performance
index. One possible candidate for the performance index is,
= E [<**) *k)P (6.1.1)
k=l
In order to minimize equation (6.1.1) it is necessary to evaluate the partial derivatives of
J with respect to each of the a and b coefficients in equation (6.1.0). Evaluating the
partial derivative with respect to a, yields:
dn'(k)
daft)
 n
'(fco + XX(*)
s=1
dn '(ks)
dafk)
(6.1.2)
This equation shows that the partial derivative of y with respect to a, at time k, depends
on all of the previous values of n/. By assuming that the values of the coefficients are
changing sufficiently slowly, equation (6.1.2) can be written as,
oaft) U 3a/ks)
(6.13)
and in a similar way, and again using the small step size assumption, the partial derivative
with respect to b is,
50
dn\k)
dbfi)
i
* n(k~j) + J2as(k)
S=1
dn '{ks)
dbj(ks)
Then defining a vector called a regressor as,
X(k)=\n '(kI)n '(k1) n{km)n(l$
There is an autoregressively filtered version of X defined as,
tr(Â£) = X(k) + Y as(k)^(ks)
S = 1
Then the algorithm for adapting the filter coefficients is,
W(k+1) = W(k) + A\\?(k)e(k)
in which the components of equation (6.1.7) are defined as,
W{k)^afky.ax{k) bm(k)...b0(k)]'
e(k) = d(k)X\k)W(k) = x(k) X\k)W{k)
A = dia,pm p0]
(6.1.4)
(6.1.5)
(6.1.6)
(6.1.7)
(6.1.8)
(6.1.10)
(6.1.9)
51
and then the output of the predictor is given by,
n(k+l) = X(k+l)W(k+l) (6.1.11)
This algorithm is computationally difficult since it requires l past values of ir be
kept for the autoregressive equation given in (6.1.6). One alternative is to use filtered
versions of the input and output sequences presented to the prediction filter. These
filtered signals are,
i
nF(k) = n'{k)+Y^as(k)rip{ks) (6.1.12)
J=1
nF( = n(k)+J2 affin^ks) (6.1.13)
S1
and then define a new ir as,
^#) = njÂ£kI)njr(kl) nF{km)nF(k^
(6.1.14)
This approximation to ir is not identical to the one given in equation (6.1.6) but given the
assumption that the coefficients of the filter are changing slowly, either can be used.
6.2 Advantages and disadvantages of the 1IR implementation
The advantage of the HR implementation over the FIR method is that the HR
52
method can produce good results with fewer coefficients in the prediction filter. The FIR
method must have enough filter coefficients to reproduce the desired impulse response
of the desired prediction filter. Obviously the HR method can produce a filter with an
impulse response much longer than the length of the coefficient vectors.
One major disadvantage of the DR method is that it can produce an unstable filter.
The required trajectory for achieving an optimal result may take one or more of the
system poles outside of the unit circle. This requires additional computational effort to
detect such excursions, and to correct the pole placement to a stable configuration. This
requires evaluating the position of the poles at each iteration, which can be a
computationally expensive operation.
Although the DR implementation can, in some circumstances, produce a simpler
filter, this not always true. Specifically, the following section will demonstrate that in the
specific case being considered here the HR prediction filter is not simpler than the FIR.
6.3 An example using an HR prediction filter in an adaptive line enhancer
In the following example, the same data that have been used in previous chapters
are applied to an ALE using an HR prediction filter. In figure 6.3.0 the result when 10
poles and 50 zeros are used in the prediction filter is presented. These results show that
the ALE using this type of filter do converge to a solution, but that they are not
substantially better than the equivalent FIR structure.
53
1=10 m =50 mu = 5 e 0 0 6
Figure 6.3.1 shows the results of running an HR implementation which has one pole and
50 zeros.
54
20
I =1 m =50 mu =5e006
The simulation shown in figure 6.3.2 reduces the number of feedforward terms to 25.
The results clearly show that 25 feedforward terms are unable to effectively predict the
narrowband interference. This can be explained by the fact that the input to the
prediction filter has a delay of 16 samples, and therefore there are only nine samples that
are strongly correlated with the desired signal. In this case where the reference is a
delayed version of the input, the increased complexity of the IIR implementation is not
justified.
55
20
I =1 m =25 mu =5e006
Fig 6.3.1
56
7.0 Conclusions
This thesis has presented several adaptive methods to remove narrowband
interference from a desired broadband signal. It has been shown that a delayed version
of the combined signal and interference can be used as the reference input to an adaptive
linear prediction filter, and that the output of the prediction filter can be used to eliminate
the narrowband interference. This type of system is usually called an adaptive line
enhancer. Traditionally, the ALE has been used to extract narrowband signals from
narrowband interference. It has been shown in this thesis that the ALE is equally
valuable when the narrowband signal is the signal of interest. Several different
implementations for the prediction filter used in the ALE have been examined.
The Leastmeansquare (LMS) method in its various forms have been shown to
provide good results in this application. Several modifications to the basic LMS
algorithm have been examined with encouraging results. Computer simulations have been
presented that show that the basic LMS, the normalized LMS and the power estimate
normalized algorithms produce essentially identical results.
The large number of feedforward terms are required to accommodate the delay
in the reference results in an autocorrelation matrix that is singular or poorly conditioned
in many cases. It was shown that while zero eigenvalues for R do not prevent the system
from correctly adapting to the signals, it is very important to select correct initial
conditions for the filter coefficients so that the uncoupled modes in the dynamic equations
do not adversely affect the performance of the filter. Specifically, it appears important to
57
initialize the filter coefficient vector to zero. When the coefficients were initialized to
one, the uncoupled modes that did not converge resulted in filter coefficients that caused
a serious degradation in the output of the system as a particle pulse propagated though
the prediction filter.
The methods usually employed, e.g. coefficient leakage and noise injection, to
force the system to converge to a unique solution when R is singular or poorly
conditioned do not seem to improve the performance of the system when a pulse is
present. These techniques result in a unique solution which is important when the values
of the coefficients after convergence are the desired output of the system. Fortunately
this application does not require a unique solution for W, so when the initial conditions
on the coefficients are correctly chosen the system exhibits acceptable performance.
The RLS technique has been shown to produce results that are very similar in
performance to the LMS techniques. In the special case where the statistics of the signals
are timeinvariant, the RLS algorithm was shown in section 5.4 to be identical to the
accelerated LMS if the parameters of the LMS algorithm are chosen correctly.
The RLS method has the advantage that it will converge much more quickly to
an acceptable solution when the statistics of the signals are stationary. When the statistics
are not stationary, the normal RLS solution has the undesirable inability to adapt to
changes. The exponentially weighted RLS algorithm corrects for this deficiency with little
additional computation.
The IIR solution presented, which is only one of several possibilities, has many
58
complications and does not seem to provide any advantages. It requires as many feed
forward terms as do the FIR techniques in order for it to compensate for the delay in the
reference signal, as was shown by [Etter and Stems, 1978], It is also a much more
complicated algorithm to implement, especially when the test for poles outside the unit
circle in the Zplane are included. In situations where the delay could be reduced the IIR
techniques may prove to be more useful because acceptable performance might be
achieved with fewer coefficients. Given the low cost of DSP processors that are
specifically designed to implement FIR algorithms, the complications associated with HR
implementations make the FIR techniques much more desirable for the narrowly defined
problem of suppressing the periodic interference commonly found in optical particle
detection.
In the past, advances in the sensitivity of optical particle counters have been made
primarily through innovation in laser technology and optical system design. While
advances in these areas will continue, the quantum leaps they have produced in the past
are unlikely to continue. The future advances in these types of instruments will likely
come from improvements in signal processing technology. As the cost of powerful
computational engines for techniques such as those presented in this thesis continue to
decline, they will certainly play a greater role is optical particle detectors, as well as many
other applications.
59
A.l Matlab routines used to generate the examples used in this thesis
All of the calculations in this thesis were accomplished using Matlab version
4.1 for Windows. The routine for the normalized LMS algorithm is accomplished by
two Matlab files. The first is the user interface that queries for needed parameters,
initializes global variables, and plots the results. The second implements the adaption
algorithm.
filt.m:
clear; % clear all variables.
% the following command is needed to insure that matlab will generate
% meta files that are compatible with word perfect
system_dependent( 14,'on');
% initialize the input data string
% creates the following variables:
% len = length of the input data vector
% N= The number of filter coefficients
% noise = standard deviation of the additive noise
% iv=The initial value for the filter coefficients
% x = actual input data
len = input('What is data length?');
N = input('What is filter length?');
mu = input('What is mu?');
noise = input('What is the noise magnitude? ');
iv = input(' What is the initial W value? ');
% set the seed for the random number generator so the data is reproducible
randn('seed',0);
x = noise randn(l,len);
t = 1 den;
% add two sinusoids to generate the narrowband interference
for I = 1 den
x(l,i) = x(l,i) + 10*sin( 1/50 2 pi);
%x(l, i) = x(l,i) + 1/100;
x(l,i) = x(l,i) + 10*sin( I / 100 2 pi);
end;
% generate a pulse 10 high and 9 wide
for 1= 4:4
60
x(l,3*len/4 +1) = x(l,3*len/4 +1) + 10*cos(I/4 pi / 2)A2;
end;
% initialize the filter parameters
% creates the following variables:
% W = FIR filter coefficients
% d = delay
d 16;
gamma = 0.00001;
% run the filter on the data
% creates the following variables:
% NB = vector of the narrowband signal (zero for d+1 locations)
% BB = vector of the broadband signal (zero for d+1 locations)
W = iv ones(l,N);
nlms; % the LMS algorithm
% plot the results
subplot(3,l,l), plot(x,'w'), axis([0 len 20 20]);
subplot(3,l,2), plot(BB,'w'), axis([0 len 20 20]);
subplot(3,l,3), plot(BB>'), axis(rien*3/450 len*3/4+2*N 20 20]);
nlms.m:
% script to run the normalized LMS algorithm
NB = zeros(l,len);
BB = zeros(l,len);
grad = zeros(l,len);
for k = N+d+l:len
X = x(l,kNd+l:kd);
NB(k) =W X';
BB(k) = x(k) NB(k);
g = BB(k) mu .* X / (gamma + X X');
W = W + g;
end;
Other versions of the LMS algorithms are accomplished by replacing the call to nlms
in the filt.m script file with calls to the following scripts,
% lms.m script to run the unnormalized LMS algorithm
% script to run the normalized LMS algorithm
NB = zeros(l,len);
BB = zeros(l,len);
for k = N+d+ l:len
61
X = x(l,kNd+l :kd);
NB(k) =W X1;
BB(k) = x(k) NB(k);
g = BB(k) mu .* X;
W = W + g;
end;________________________________________________________________
% leakage.m
% script to run the normalized LMS algorithm with leakage
NB zeros(l,len);
BB = zeros(l,len);
grad = zeros(l,Ien);
for k = N+d+1 den
X = x(l,kNd+l :kd);
NB(k) =W X';
BB(k) = x(k) NB(k);
mu = alpha / (epsilon + X X');
g = BB(k) mu .* X;
W = (lmu*gamma) .* W + g;
end;____________________;____________________________________________
% power, m
% script to run the normalized LMS algorithm with power normalization
NB = zeros(len,l);
BB = zeros(len,l);
P = 0;
for k = N+d+1 :len
X = x(kNd+l:kd,l);
NB(k) =W X;
BB(k) = x(k) NB(k);
alpha = mu /(gamma + p);
p = p (1beta) + beta XX(k) XX(k);
g = BB(k) alpha X;
W = W + g;
end;________________________________________________________________
% rls.m
% script to run the recursive least square algorithm
NB = zeros(len,l);
BB = zeros(len,l);
Rinv = 10 eye(N);
for k = N+d+1 den
% Form X(k)
62
X = x(kNd+l :kd, 1);
%
%Compute the apriori output
%
yO X' W;
%
%Which is also the Narrowband signal of interest
%
NB(k) = yO;
%
%Compute the apriori error
%
eO = x(k) yO;
%
%Which is also the Broadband signal
%
BB(k) = eO;
%
%Compute the filtered information vector
%
Z Rinv X;
%
%Compute the normalized error power
%
q = X' Z;
%
%Compute the gain constant
%
v = l/(l+q);
%
%Compute the normalized filter information vector
%
Ztilde = v Z;
%
%Update the weight vector
%
W = W + eO Ztilde;
%
%Update the inverse correlation matrix
Rinv = Rinv Ztilde Z';
end;
63
% rls.m
% script to run the recursive least square algorithm with exponential weighting
NB = zeros(len,l);
BB = zeros(len,l);
for k = N+d+1 :len
%
% Form X(k)
X = x(kNd+l :kd, 1);
%
%Compute the apriori output
yO = X1 W;
%
%Which is also the Narrowband signal of interest
NB(k) = yO;
%
%Compute the apriori error
eO = x(k) yO;
%
%Which is also the Broadband signal
BB(k) = eO;
%
%Compute the filtered information vector
Z = Rinv X;
%
%Compute the normalized error power
q = X' Z;
%
%Compute the gain constant
v = l/(rho + q);
%
%Update the weight vector
W = W + eO*v*Z;
%
%Update the inverse correlation matrix
Rinv = (1/rho) (Rinv v Z Z');
end;
The following two files were used to generate the simulation outputs for the
HR algorithm.
64
clear;
% the following command is needed to insure that matlab will generate
% meta files that are compatible with the dimwitted word perfect
sy stem_dep endent( 14,'on');
% initialize the input data string
% creates the following variables:
% len = length of the input data vector
% x = actual input data
len = input('What is data length?');
al = input('How many a terms ?');
bl = inputCHow many b terms ?');
mu = input('What is mu?');
lambda = mu eye(al+bl+l);
noise = input('What is the noise magnitude? ');
% set the seed for the random number generator so the data is reproducible
randn('seed',0);
x = noise randn(len, 1);
t = 1 :len;
% add two sinusoids and a ramp to the random noise
for i = 1 :len
x(i,l) = x(i,l) + 10*sin( i / 50 2 pi);
x(i,l) = x(i,l) + 10*sin( i /100 2 pi);
end;
% generate a pulse 10 high and 9 wide
for i= 4:4
x(3*len/4 + i,l) = x(3*len/4 + i,l) + 10*cos(i/4 pi / 2)A2;
end;
% initialize the filter parameters
% creates the following variables:
% mu = adaptation constant
% W = FIR filter coefficients
% d = delay
% N = number of filter coefficients
d= 16;
W = zeros(al+bl+l,l);
% run the filter on the data
% creates the following variables:
% NB = vector of the narrowband signal (zero for d+1 locations)
% BB = vector of the broadband signal (zero for d+1 locations)
iir;
% plot the results
65
%subplot(2,1,1),
plot(BB,'w'); axis([0 len 20 20]);
title([ '1 num2str(al)'m num2str(bl) mu num2str(mu) ]);
%subplot(2,l,2), plot(BB,'w'), axis([len*3/450 len*3/4+2*N 20 20]);
%title([ 'Expanded Broadband1]); ____________________________________
% iir.m
% script to run the normalized IIR algorithm
NB = zeros(len,l);
BB = zeros(len, 1);
npf = zeros(len,l);
nhf = zeros(len,l);
for k = al+bl+d+1 :lenl
X_k = [NB(kal:kl,l); x(kdbl:kd,l)];
NB(k) = X_k' W;
BB(k) = x(k) NB(k);
npf(k) = x(kd) + sum(W(l:al,l) .* x(kdal:kdl,l));
nhf(k) = NB(k) + sum(W(l:al,l) NB(kal:kl,l));
psi = [ nhf(kal:kl,l); npf(kbl:k,l)];
W = W + lambda psi BB(k);
end;__________________________________________________________________
66
References
Etter, Delores M. and Steams, Samuel D., An Adaptive Noise Canceller with
Applications to Seismic Noise, Conference Record, 12th Annual Asilomar
Conferences on Circuits, Systems, and Computers, 1978.
Glover, John R., Adaptive Noise Canceling Applied to Sinusoidal Interferences,
IEEE Trans. On Acoustics, Speech, and Signal Processing, vol ASSP25, no. 6,
December 1977.
Haykin, Simon, Adaptive Filter Theory, Second Edition, New Jersey: Prentice Hall,
1986.
Orfanidis S. J., Optimum Signal Processing: An Introduction, Second edition, New
York: McGrawHill, 1988.
Trichler, John R. et al., Theory and Design of Adaptive Filters, New York: John
Wiley & Sons, 1987.
Widrow, Bernard, Glover, John R. Jr. et al., Adaptive Noise Cancelling: Principles
and Applications, Proceedings of the IEEE, vol. 63, No. 12, Dec 1975.
Zahm, C. L., Application of Adaptive Arrays to Suppress Strong Jammers in the
Presence of Weak Signals, IEEE Trans. On Aerospace and Electronic Systems, vol.
AES0, pp. 260271, March 1973.
Zeidler, James R., Adaptive Enhancement of Multiple Sinusoids in Uncorrelated
Noise, IEEE Transactions of Acoustics, Speech, and Signal Processing, vol. ASSP
26, No. 3, June 1978
67
