Tuesday, December 21, 2010

Introduction to Adaptive Filtering

The Kalman filter is just one of many adaptive filtering (or estimation) algorithms. Despite its elegant derivation and often excellent performance, the Kalman filter has two drawbacks:
  1. The derivation and hence performance of the Kalman filter depends on the accuracy of the a priori assumptions. The performance can be less than impressive if the assumptions are erroneous.
  2. The Kalman filter is fairly computationally demanding, requiring O(P2) operations per sample. This can limit the utility of Kalman filters in high rate real time applications.
As a popular alternative to the Kalman filter, we will investigate the so-called least-mean-square (LMS) adaptive filtering algorithm.
The principle advantages of LMS are
  1. No prior assumptions are made regarding the signal to be estimated.
  2. Computationally, LMS is very efficient, requiring O(P) per sample.
The price we pay with LMS instead of a Kalman filter is that the rate of convergence and adaptation to sudden changes is slower for LMS than for the Kalman filter (with correct prior assumptions).

Adaptive Filtering Applications

Channel/System Identification

Figure 1
Channel/System Identification
Channel/System Identification (csid.png)

Noise Cancellation

Suppression of maternal ECG component in fetal ECG (Figure 2).
Figure 2: Cancelling maternal heartbeat in fetal electrocardiography (ECG): position of leads.
Figure 2 (.png)
Figure 3
Figure 3 (ecg.png)
y is an estimate of the maternal ECG signal present in the abdominal signal (Figure 4).
Figure 4: Results of fetal ECG experiment (bandwidth, 3-35Hz; sampling rate, 256Hz): (a)reference input (chest lead); (b)primary input (abdominal lead); (c)noise-canceller output.
Figure 4 (.png)

Channel Equalization

Figure 5
Channel Equalization
Channel Equalization (ce.png)

Adaptive Controller

Figure 6: Here, the reference signal is the desired output. The adaptive controller adjusts the controller gains (filter weights) to keep them appropriate to the system as it changes over time.
Adaptive Controller
Adaptive Controller (ac.png)

Iterative Minimization

Most adaptive filtering alogrithms (LMS included) are modifications of standard iterative procedures for solving minimization problems in a real-time or on-line fashion. Therefore, before deriving the LMS algorithm we will look at iterative methods of minimizing error criteria such as MSE.
Conider the following set-up: xk:observation yk:signal to be estimated

Linear estimator

yk=w1xk+w2x<apply><minus></minus>k1</apply>++wpx<apply><plus></plus><apply><minus></minus>kp</apply>1</apply>(1)
Figure 7
Figure 7 (fir.png)
Impulse response of the filter: ,0,0,w1,w2,wp,0,0,

Vector notation

yk=xkTw(2)
Where xk=(
xk
x<apply><minus></minus>k1</apply>
x<apply><plus></plus><apply><minus></minus>kp</apply>1</apply>
)
 and w=(
w1
w2
wp
)

Error signal

ek=ykyk
=ykxkTw
(3)

Assumptions

(xk,yk) are jointly stationary with zero-mean.

MSE

E[ek2]=E[(ykxkTw)2]
=E[yk2]2wTE[xkyk]+wTE[xkxkT]w
=Ryy2wTRxy+wTRxxw
(4)
Where Ryy is the variance of yk2Rxx is the covariance matrix of xk, and Rxy=E[xkyk] is the cross-covariance between xk and yk

NOTE: 

The MSE is quadratic in W which implies the MSE surface is "bowl" shaped with a unique minimum point (Figure 8).
Figure 8
Figure 8 (mse.png)

Optimum Filter

Minimize MSE:
w
 (E[ek2])
=-2Rxy+2Rxxw=0
wopt=Rxx-1Rxy
(5)
Notice that we can re-write Equation 5 as
E[xkxkTw]=E[xkyk](6)
or
E[xk(ykxkTw)]=E[xkek]
=0
(7)
Which shows that the error signal is orthogonal to the input xk (by the orthogonality principle of minimum MSE estimator).

Steepest Descent

Although we can easily determine wopt by solving the system of equations
Rxxw=Rxy(8)
Let's look at an iterative procedure for solving this problem. This will set the stage for our adaptive filtering algorithm.
We want to minimize the MSE. The idea is simple. Starting at some initial weight vector w0, iteratively adjust the values to decrease the MSE (Figure 9).
Figure 9
In One-Dimension
In One-Dimension (iterative.png)
We want to move w0 towards the optimal vector wopt. In order to move in the correct direction, we must move downhill or in the direction opposite to the gradient of the MSE surface at the point w0. Thus, a natural and simple adjustment takes the form
w1=w0
1
2
 
μ
w
 (E[ek2])
|w=w0
(9)
Where μ is the step size and tells us how far to move in the negative gradient direction (Figure 10).
Figure 10
Figure 10 (adjust.png)
Generalizing this idea to an iterative strategy, we get
wk=wk1
1
2
 
μ
w
 (E[ek2])
|w=wk1
(10)
and we can repeatedly update ww0,w1,,wk. Hopefully each subsequent wk is closer to wopt. Does the procedure converge? Can we adapt it to an on-line, real-time, dynamic situation in which the signals may not be stationary?

Content actions

GIVE FEEDBACK:

System Identification for the Rectilinear Plant

Summary: The objective of this lab is to review of the behavior of second-order systems. Students will use gain a better understanding of the importance system identification. Students will also develop a hands-on understanding of the concept of hardware gain and why it will play a crucial role in controller design. System Identification will be implemented in LabVIEW using the System Identification Toolkit.





System Identification for the Rectilinear Plant

Objectives

  • Review of the behavior of second-order systems
  • Understand the importance of system identification
  • Understand the concept of hardware gain and why it will play a crucial role in controller design

Pre-Lab

  1. Derive the equations of motion for the 1DOF spring-mass system with friction damping shown below and find the transfer function. Simulate the unforced, natural response in LabVIEW to an initial displacement of 3.0cm (look for the Initial Response function in the Control Design Toolkit). Determine the natural frequency from the plot. Use the following parameters: M=3.0kg k=350N/m c=4.0Nsec/m
    Figure 1: 1DOF Spring-Mass System with Friction Damping
    Figure 1 (spring-mass.jpg)

Lab Procedure

  1. Power on the PXI, and while it's booting up configure the plant as explained in the next step.
  2. Clamp the second and third mass carriages, and attach the medium stiffness spring between the first and second carriage as shown in Fig. 2. Secure four 0.5kgbrass weights to the first and second carriages.
    Figure 2: Rectilinear Plant System ID Configurations
    Figure 2 (plantconfigurations.jpg)
  3. Start LabVIEW on your host PC and target the RT system.
  4. Open the Lab 2 - Rectilinear System Identification.vi file, enter a loop time of 0.005 seconds, and then run the VI. Also, turn on the amplifier now.
  5. Reset the encoders if necessary and enter the parameters to execute a zero amplitude step with a dwell time of 4000ms. Prepare to manually displace the first mass carriage. When you are ready, execute the step command, and manually displace the first carriage approximately 3cm, then release it. Observe the natural response of the carriage, and wait for the system to finish executing the trajectory and acquiring the data.
  6. Recall that for underdamped systems, the damped natural frequency is related to the undamped natural frequency through the following relationship:
    ωd=ωn\1ζ2(1)
    and it follows that
    ωn=
    ωd
    \1ζ2
    (2)
    Therefore, for light damping (small values of ζ)
    ωnωd(3)
    For the purposes of plant identification, you may assume that the damping due to friction in the system is sufficiently small to satisfy the above condition. Plot the encoder 1 position data and calculate the natural frequency of this system in both Hz and rad/sec . Call this value ωnm11 (the subscript m11 denotes mass #1 trial #1). Use only oscillations with amplitude greater than 1000 counts in your calculation. This is because smaller amplitudes start to become dominated by nonlinear friction effects.
  7. Remove the four brass weights and repeat steps 5 and 6 to obtain ωnm12 for the unloaded carriage. Decrease the execution time if necessary.
  8. Measure the initial cycle amplitude X0 and the last cycle amplitude Xn for the n cycles measured in Step 6. Using relationships associated with the logarithmic decrement find the damping ratio (call it ζm12) as
    ζ
    \1ζ2
     =
    1
    2pin
     ln
    X0
    Xn
     ζ
    1
    2pin
     ln
    X0
    Xn
    (4)
  9. Configure the plant as shown in Fig. 2b where the first mass carriage is now clamped and the second is free. Repeat steps 5 through 8 for the second mass carriage to obtain ωnm21 , ωnm22 and ζm22 . How does this damping ratio compare with that for the first mass? Be sure to save this plotted data as it will be used in the next experiment.
  10. Connect the mass carriage extension bracket and dashpot to the second mass as shown in Fig. 2c, and load four brass weights onto the second mass carriage. Open the damping (air flow) adjustment knob 2.0 turns from the fully closed position. Perform the necessary steps to obtain ζd where the "d" subscript denotes the dashpot. For this calculation use amplitudes  500 counts.
  11. Each brass weight has a mass of 0.5±0.01kg. (You may weigh the pieces if a more precise value is desired). Call the mass of the four weights combined mw. Recall that the dynamics for a second-order system in terms of the Laplace variable s is given by s2+2ζωns+ωn2 Compare the above with the equations of motion you derived in terms of mass, damping coefficient, and spring constant ( mc, and k). Use your knowledge of second-order systems along with your recorded values ofωnm11 and ωnm12 to solve for the unloaded carriage mass mc1, and spring constant kRepeat these calculations for the trials involving the second mass carriage and spring, mc2 and k. Since the same spring was used in the experiments, the values for k that result from the above calculations should be very close. You may use the average of the two for all future calculations and experiments. Call this value kmed to denote the medium stiffness spring. Now with your recorded values for ζcalculate c1, the damping coefficient for the first mass carriage. Repeat this calculation to find the damping coefficient for the second mass carriage, c2Finally, calculate the damping coefficient of the dashpot, cd.
  12. Remove the carriage extension bracket and dashpot from the second mass carriage, replace the medium stiffness spring with the high stiffness spring, and perform the necessary steps to obtain the resulting natural frequency ωm23. Repeat this frequency measurement using the least stiff spring to obtain ωm24.Now you can compute khigh and klowAll dynamic parameters have been identified! Values for mc1 and mc2 for any configuration of masses may be found by adding the calculated mass contribution of the weights to that of the unloaded carriages. Be sure that you have filled in all of the system parameters in the table provided before proceeding.

    NOTE: 

    The following is necessary to establish the hardware gain for control modeling purposes.
  13. Remove the spring connecting the first and second carriages and secure four 0.5kg brass weights onto the first carriage. Using the stop blocks, secure the second mass carriage clear from the first, and remove any stop blocks that are constraining the first carriage. Verify that the masses are secure and that the carriage slides freely. Position the first mass approximately 3cm to the left (negative x1 position) of its center of travel.
  14. Perform a 2.0V bidirectional step with a dwell time of 75ms (recall that on the PXI, ±215 counts corresponds to a voltage of ±10V).
  15. Plot the velocity data and observe the four segments of the profile with nominal shapes of: linear increase (constant acceleration), nearly constant velocity (zero acceleration), linear decrease (constant deceleration), and constant. Obtain the acceleration from the linear, positive-sloped segment of the velocity profile. Repeat this for the negative-sloped segment. You may call these two values x1eacc and x1edec, respectively. Be sure that you account for the scaling of velocity data that is shown on the front panel graph indicator, and remember to include units in your calculation.Calculate the average of the magnitudes of the positive and negative accelerations (call this value x1e) to be used in obtaining khwbelow.
  16. Recall that the hardware gain for this system is given by khw=kckaktkmpkepkewhere kc, the DAC gain, = 10V/32,768 DAC counts ka, the servo amp gain, = to be determined Amps/V kt, the servo motor torque constant, = to be determined Nm/Amp kmp, the motor pinion pitch radius inverse, = 26.25m1 kep, the encoder pinion pitch radius inverse, = 89m1 ke, the encoder gain, = 16,000 pulses/2 πradians Using what you have learned about hardware gain, calculate a value for the product kakt. Use this value to compute the hardware gain khw.
  17. Save any files or plots of interest. Stop and exit the VI, power off the amplifier and PXI, and return all plant materials to the instructor.

Post-Lab

  1. Complete the table below making sure to include units.
    Figure 3: Post Lab Table
    Figure 3 (postlab.jpg)
  2. How do your values of mc1 and c1 compare with mc2 and c2. Can you explain why they may be different?
  3. What were your values for ζm12 and ζm22 Are these values sufficiently small for the approximations made in steps 6 and 8 to be valid?
  4. Explain the importance of accurately identifying the plant parameters and hardware gain in terms of how they will affect controller design.