Adaptive optimal fitlering using the logarithmic number system

MISSING IMAGE

Material Information

Title:
Adaptive optimal fitlering using the logarithmic number system
Added title page title:
Adaptive optimal filtering using the logarithmic number system
Physical Description:
viii, 144 leaves : ill. ; 28 cm.
Language:
English
Creator:
Papadourakis, George Michael, 1959-
Publisher:
s.n.
Publication Date:

Subjects

Subjects / Keywords:
Kalman filtering   ( lcsh )
Stochastic processes   ( lcsh )
Electrical Engineering thesis Ph. D
Dissertations, Academic -- Electrical Engineering -- UF
Genre:
bibliography   ( marcgt )
non-fiction   ( marcgt )

Notes

Thesis:
Thesis (Ph. D.)--University of Florida, 1986.
Bibliography:
Bibliography: leaves 136-143.
General Note:
Typescript.
General Note:
Vita.
Statement of Responsibility:
by George Michael Papadourakis.

Record Information

Source Institution:
University of Florida
Rights Management:
All applicable rights reserved by the source institution and holding location.
Resource Identifier:
aleph - 030573201
oclc - 17647560
System ID:
AA00022722:00001

Table of Contents
    Title Page
        Page i
    Dedication
        Page ii
    Acknowledgement
        Page iii
    Table of Contents
        Page iv
        Page v
        Page vi
    Abstract
        Page vii
        Page viii
    Chapter 1. Introduction
        Page 1
        Page 2
        Page 3
        Page 4
    Chapter 2. Survey of literature
        Page 5
        Page 6
        Page 7
        Page 8
        Page 9
        Page 10
        Page 11
        Page 12
        Page 13
        Page 14
        Page 15
        Page 16
        Page 17
        Page 18
        Page 19
        Page 20
    Chapter 3. Kalman filters
        Page 21
        Page 22
        Page 23
        Page 24
        Page 25
        Page 26
        Page 27
        Page 28
        Page 29
        Page 30
        Page 31
        Page 32
        Page 33
        Page 34
        Page 35
    Chapter 4. Logarithmic number system
        Page 36
        Page 37
        Page 38
        Page 39
        Page 40
        Page 41
        Page 42
        Page 43
        Page 44
    Chapter 5. Error models for LNS arithmetic
        Page 45
        Page 46
        Page 47
        Page 48
        Page 49
        Page 50
        Page 51
        Page 52
        Page 53
        Page 54
        Page 55
        Page 56
    Chapter 6. Logarithmic Kalman filters
        Page 57
        Page 58
        Page 59
        Page 60
        Page 61
        Page 62
        Page 63
        Page 64
        Page 65
        Page 66
        Page 67
        Page 68
        Page 69
        Page 70
        Page 71
    Chapter 7. Adaptive Kalman filters
        Page 72
        Page 73
        Page 74
        Page 75
        Page 76
        Page 77
        Page 78
        Page 79
        Page 80
        Page 81
        Page 82
        Page 83
        Page 84
        Page 85
        Page 86
        Page 87
        Page 88
        Page 89
        Page 90
        Page 91
        Page 92
        Page 93
        Page 94
        Page 95
        Page 96
        Page 97
        Page 98
        Page 99
        Page 100
        Page 101
        Page 102
        Page 103
        Page 104
        Page 105
        Page 106
        Page 107
    Chapter 8. Implementation of Kalman filters using systolic architectures
        Page 108
        Page 109
        Page 110
        Page 111
        Page 112
        Page 113
        Page 114
        Page 115
        Page 116
        Page 117
        Page 118
        Page 119
        Page 120
        Page 121
        Page 122
        Page 123
        Page 124
        Page 125
        Page 126
        Page 127
        Page 128
        Page 129
        Page 130
        Page 131
        Page 132
    Chapter 9. Summary and conclusions
        Page 133
        Page 134
        Page 135
    Bibliography
        Page 136
        Page 137
        Page 138
        Page 139
        Page 140
        Page 141
        Page 142
        Page 143
    Biographical sketch
        Page 144
        Page 145
        Page 146
Full Text













ADAPTIVE OPTIMAL FITLERING USING
THE LOGARITHMIC NUMBER SYSTEM










By

GEORGE MICHAEL PAPADOURAKIS


A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL
OF THE UNIVERSITY OF FLORIDA IN
PARTIAL FULFILLMENT OF THE REQUIREMENTS
FOR THE DEGREE OF DOCTOR OF PHILOSOPHY




UNIVERSITY OF FLORIDA


1986




























DEDICATED TO

EAEU9EPIA














ACKNOWLEDGMENTS


I would like to thank my advisor Dr. Fred J. Taylor

for his valuable guidance and assistance during my

graduate studies. All these years he has provided me

the much needed support and encouragement to complete my

doctoral work.

I would also like to thank Dr. Donald G. Childers,

Dr. Antonio A. Arroyo, Dr. George Logothetis, and

Dr. Jask R. Smith for serving on my supervisory

committee. In addition, I would like to thank

Dr. Jose C. Principe for his useful comments.

Finally, I would like to thank my fellow graduate

students and especially A. Skavantzos and A. Stouraitis

for their help.














TABLE OF CONTENTS


PAGE
ACKNOWLEDGMENTS .................................. ....... iii

ABSTRACT...... ............................................ vii

CHAPTER

ONE INTRODUCTION................ ......................... 1

TWO SURVEY OF LITERATURE........................o....... 5

Kalman Filters............ ...................... ... 5
Correlated Process and Measurement Noise.......... 5
Wiener Filter................... .................. .6
Random Sampling Times..................... ........ 6
Effect of Inaccurate Mathematical Models and
Statistics...... ................................ 6
Performance (Sensitivity) Analysis............... .7
Square Root Filtering................. ........... .8
Suboptimal Filtering ................................ 8
Compensation of Linear Model Inadequacies.........9
Linear Smoothing Problem......................... 10
Observers.............. ....... .................... 10
Adaptive Estimation.............................. 11
Applications............ .......................... 11
Logarithmic Number System.......................... 14
Logarithmic Conversion........................... 15
Addition and Subtraction Algorithms.............. 16
Multiplication and Division Algorithms..........17
Table Reduction Techniques.,..................... 17
Extended Precision LNS........................... 18
Digital Filtering Using LNS....................... 18
FFT Implementation with LNS ......................19
Other Applications with LNS ..................... 19

THREE KALMAN FILTERS .................................... 21

Statement of Problem............................... 21
Least Square Estimation............................ 22












Estimation of Parameters Using
Weighted Least-Squares........................... 25
Recursive Filters .................................. 27
Discrete Kalman Filters ............................28
Kalman Filters with Deterministic Inputs...........35

FOUR LOGARITHMIC NUMBER SYSTEM .......................... 36

LNS Representation ................................. 37
LNS Arithmetic Operations .......................... 38
Multiplication................................... .38
Addition and Subtraction ......................... 39
A Hybrid Floating Point Logarithmic
Arithmetic Unit .................................. 40
Floating Point Format ............................ 41
(FU) Square Unit................................. .41

FIVE ERROR MODELS FOR LNS ARITHMETIC .................... 45

Input Quantization ................................. 45
Coefficient Quantization........................... 54
Quantization in Arithmetic Operations.............. 54

SIX LOGARITHMIC KALMAN FILTERS ......................... 57

LNS Kalman Filter .................................. 57
Theoretical Error Analysis in LNS ..................58
Mean Error Analysis .............................. 62
Variance Error Analysis ..........................63
Simulation Studies .................................64

SEVEN ADAPTIVE KALMAN FILTERS ............................ 72

The Effect of Erroneous Models on the
Kalman filter Response .......................... 72
The Effect of Input and Output Noise Covariances
on the Kalman Gains ..............................75
Adaptive Kalman Filtering ..........................84
Algorithm Implementation............................86
Kalman Filter .................................... 86
Autocorrelation.................................. 90
Search Algorithm ..................................90
Algorithm Integration ..............................94
Mehra's Adaptive Algorithm .........................94
Carew's and Balanger's Adaptive Algorithm..........97
Experimental Results ............................... 99












EIGHT IMPLEMENTATION OF KALMAN FILTERS USING
SYSTOLIC ARCHITECTURES .......................... 108

Systolic Architectures ............................108
Orthogonal Systolic Arrays ........................115
Kung's Systolic Array ...........................116
Liu's and Young's Systolic array ................118
Pipelined Systolic Array ........................121
Kalman Filtering Using Systolic Arrays ............123

NINE SUMMARY AND CONCLUSIONS ...........................133

BIBLIOGRAPHY............................................... 136

BIOGRAPHICAL SKETCH ...................................... .144















Abstract of Dissertation Presented to the Graduate School
of the University of Florida in Partial Fulfillment of the
Requirements for the Degree of Doctor of Philosophy





ADAPTIVE OPTIMAL FILTERING USING
THE LOGARITHMIC NUMBER SYSTEM


By

GEORGE MICHAEL PAPADOURAKIS


August 1986



Chairman: Dr. Fred J. Taylor
Major Department: Electrical Engineering


The Kalman filter has been one of the most widely

applied techniques in the area of modern control, signal

processing and communication applications. The

implementation of a Kalman filter using the logarithmic

number system (LNS) is considered. The choice of LNS is

highly attractive in this digital filtering setting since it

offers high precision, wide dynamic range, attractive

architecture, and ease of arithmetic operations. A

theoretical error analysis is presented concerning the mean

and variance of the actual estimations due to a


vii













finite-precision logarithmic implementation. Simulation

studies are utilized to show how closely the analytical

predictions agree with actual results. A comparison between

the finite-precision conventional Kalman filter and the

proposed scheme is performed in terms of speed and accuracy.

In addition, an adaptive Kalman filter is considered

for which the input and output noise covariances are

unknown. A new procedure is presented for the

identification of the unknown covariances using the

autocorrelation information contained in the innovation

error sequence and a modified Fibonacci search. Real-time

performance is facilitated through the use of the presented

efficient high speed digital autocorrelation algorithms and

search procedures.

Finally, the implementation of the Kalman filter using

systolic array architectures is investigated. Systolic

architectures have received much attention recently as a

means of achieving high data rates. Systolic designs have

been used to perform matrix operations making them suitable

for Kalman filtering processing. Using an orthogonal array

processor, a new algorithm is presented which increases the

throughput of matrix operations. This new algorithm will be

compared with other systolic architectures on the

implementation of Kalman filters.


viii

















CHAPTER ONE
INTRODUCTION





The theory of optimal filtering has application to a

broad range of meaningful problems. Within the family of

known optimal algorithms, the filtering technique known as

the Kalman filter has been the most visible. It has been

shown to be well suited for digital computer implementation

provided high throughput is not a prerequisite. This class

of filter has penetrated the areas of control,

communication, guidance-navigation, and data processing.

To successfully design a Kalman filter, and its

counterparts, the engineer/scientist must address both

modeling and computation problems. The Kalman filter is

intrinsically a model reference filter, and as a result, can

suffer from degraded performance if the system model is

improperly configured. In order to desensitize the filter

to model dependent errors, an adaptive filter policy should

be considered. This leads to the second problem. It has

been assumed that the Kalman filter would lend itself well

to digital computation due to its recursive structure. This











has not always been the case. The computational burdens

imposed by real-time, high speed, high-precision Kalman

filters can exceed a computers throughput capabilities. In

addition, there is always a controversy on whether floating

point or fixed point arithmetic should be used. Furthermore

there is generally the overlooked problem of burst error

management of filters operating in hostile environments

which must be attended to. The additional computational

demands of an adaptive filter often exceed the performance

capability of conventional computer architectures. This is

especially true in those cases where small wordlength,

limited capability processors are found (viz. the

microprocessor). Rather than developing theory over a real

number field, we are proposing to represent numbers over

finite fields and compressed numbering systems.

Recent developments in the area of logarithmic number

systems (LNS) have proven it to offer some significant

advantages over other number weighted systems like the fixed

or floating point. The main advantages of LNS are the

extended dynamic range on which they operate and the high

speed of operations accompanied by a remarkably regular data

flow. Another remarkable feature of the LNS architecture is

its ability to be efficiently transferred into VLSI.

Furthermore LNS arithmetic units are ideally suited to

systolic array data processing, thus increasing throughput.











The organization of the dissertation is as follows:

Chapter Two contains a comprehensive survey of the

literature on Kalman filters and logarithmic number system.

In Chapter Three, the derivation of the discrete-time

(sampled data) optimal estimator, known as the Kalman

filter, will be presented. Chapter Four contains a

description of the logarithmic number system and its

properties. In addition a logarithmic number processor,

(FU)2, will be developed. Chapter Five describes the

quantization-error models used in LNS arithmetic, namely,

input, coefficient and arithmetic operations. In Chapter

Six, the implementation of a Kalman filter using the

logarithmic number system (LNS) will be considered. A

theoretical error analysis will be presented concerning the

mean and variance of the actual estimation due to a

finite-precision logarithmic implementation. Simulation

studies will be utilized to show how closely the analytical

predictions agree with actual results. In addition

comparison between the finite-precision conventional Kalman

filter and the proposed scheme will be performed in terms of

speed and accuracy. In Chapter Seven, an adaptive Kalman

filter will be considered for which the input and output

noise covariances are unknown. A new procedure will be

presented for the identification of the unknown covariances

using the autocorrelation information contained in the

innovation error sequence and a modified Fibonacci search.








4


In Chapter Eight, implementation of a Kalman filter using a

new systolic design will be investigated and compared with

other systolic architectures. Finally in Chapter Nine, by

way of summarizing the proceeding, the significance of the

described research will be brought out and some directions

for future research will be indicated.

















CHAPTER TWO
SURVEY OF LITERATURE






Kalman Filters

In 1960 Kalman [1,2] developed the linear estimation

theory hereafter referred to as the Kalman filter. Others

[3,4] developed the theory even further and gave more

insights to the linear estimation problem. This estimation

theory can be formulated in terms of either a time-discrete

or time continuous model. Additional subjects on Kalman

filters that have been treated in the voluminous literature

will be presented.



Correlated Process and Measurement Noise

The white noise sequences have been assumed to be

uncorrelated in the formulation of the Kalman filter. This

restriction is not necessary. However, the resulting

equations are considerably more complicated [5-7].











Wiener Filter

Underlying Wiener filter design [8] is the so called

Wiener-Horf equation, its solution through spectral

factorization and the practical problems of synthesizing the

theoretically optimal filter from its impulse response. In

the case of time-invariant systems models and stationary

noises, the Wiener filter is shown to be equivalent to the

steady state Kalman filter appropriate to the given problem

[8-10].


Random Sampling Times

The measurement data have been assumed to occur at

prespecified times. It can happen that the time to which a

measurement is related is random. Also, the transition

matrix and/or the observation matrix might contain

parameters that are random. These cases have been

considered by several investigators including Tou [Ii],

Rauch [12], Gunckel [13] and others.



Effect of Inaccurate Mathematical Models and Statistics

The validity of the estimates provided by this

estimation technique require an accurate mathematical

description of the dynamical and observational processes.

The omission of terms that actually enter these processes

can result in extremely poor estimates. If the statistics

(i.e., R, Q) are not accurate descriptions of the











second-order moments of the noise processes [14], the

covariance matrix P might not be a realistic measure of the

accuracy of the estimate.




Performance (Sensitivity) Analysis

The estimates of the state are random variables. To

determine the validity of the estimates and the significance

of P, it is necessary to perform a Monte Carlo simulation of

the physical system under consideration [15].

Unfortunately, this is a costly and time consuming process.

More efficient means of generating the statistical

information is the mean analysis [16] and the covariance

analysis (171. One particular use of a covariance

performance (sensitivity) analysis is a systematic approach

to the filter tuning [18-24] process. The basic objective

of filter tuning is to achieve the best possible estimation

performance from a filter of specified structural form,

i.e., totally specified except for P0 and the time histories

of Q and R. These covariances not only account for actual

noises and disturbances in the physical system, but also are

a means of declaring how adequately the assumed model

represents the "real world" system. Once a particular

filter has been tuned, an error budget [19-24) can be

established. Essentially, this consists of repeated

covariance analyses in which the error sources in the truth











model are "tuned on" individually to determine the separate

effects of these sources.


Square Root Filtering

Measurement updating of the covariance matrix requires

a rather long wordlength to maintain acceptable numerical

accuracy. The difficulties encountered in converting a

Kalman filter tuned over a long wordlength, to an effective

algorithm on a small wordlength on line computer are well

documented [25]. To circumvent these problems in numerics

inherent to the Kalman filter algorithm, alternate recursion

relationships [25-27] have been developed to propagate and

update the covariance matrix in a square root sense. The

square root approach can yield twice the effective precision

of the conventional filter in ill-conditioned problems.

Another alternative is the U-D covariance factorization

filter which although is not actually a square root filter

is closely related to it.



Suboptimal Filtering

The dimensions of the matrices that are involved

sometimes become so large as to render their manipulation

virtually impossible in a particular digital computer.

Also, when the state vector is large, it is often

impracticable to compute estimates for every component and

it becomes necessary to consider only a subset of the











components. To reduce the dimensionality of the over-all

system and still retain a reasonable model, methods which are

not optimal in the strict sense have been devised [28].




Compensation of Linear Model Inadequacies

There is always some discrepancy between the

performance indication propagated by the filter and the

actual performance achieved in realistic applications,

because the model embodied in the filter cannot be exactly

correct. Such discrepancy is termed divergence

[18-24,29-32]. Compensation techniques have been used for

model inadequacy. The first of these methods is the

addition of pseudonoise to the assumed model and artificial

lower bounding of error covariance matrix elements

[29,30,32]. By adding such fictitious noise to the dynamic

model one "tells" the filter that it should decrease its

confidence in its own model. In the case in which a linear

model is indeed adequate, but only for a limited length of

time, then limiting of effective filter memory and

overweighting most recent data [30,33-35] and finite memory

filtering [30,36] methods are used. Finally extended Kalman

filtering [30,37] attempts to exploit linear models and

methods in the case in which a substantially better

depiction of the true system would be in the form of a

nonlinear model. The basic idea of the extended Kalman











filter is to relinearize about each state estimate once it

has been computed.



Linear Smoothing Problem

There are three classes of smoothing problems, namely

the fixed-interval, fixed-point and fixed-lag [8,38].

Fixed-interval technique is used for post-experiment data

reduction to obtain refined state estimates of better

quality than that provided by on-line filters, i.e.,

post-flight analysis of a missile. Fixed-point smoothing is

considered when there is a certain point (or points) in time

at which the value of the system state is considered

critical. For example conditions at engine burnout time are

critical to rocket booster problems. Finally the fixed-lag

technique is used when it is desirable to delay the

computation of a state estimate in order to take advantage

of additional information. It is particularly applicable to

communications and telemetry data reduction.



Observers

In some estimation problems, it may be desired to

reconstruct the state of a deterministic, linear dynamical

system -- based on exact observations of the system output.

For deterministic problems of this nature, stochastic

estimation concepts are not directly applicable. Luenberger

[39,40] formulated the motion of an observer for











reconstructing the state vector of an observable

deterministic linear system from exact measurements of the

output. The stochastic observer unifies the concepts of

deterministic Luenberger observer theory and stochastic

Kalman theory, for continuous linear systems. Analogous

results have been derived for discrete systems [41,42].


Adaptive Estimation

In order to improve the quality of the state estimates,

it would be desirable in many instances to estimate a number

of uncertain parameters in the dynamics of a measurement

model simultaneously in an online fashion. This is often

termed combined state estimation and system identification

[43-49]. In addition, one would like to readjust the

assumed noise strengths in the filter's internal model,

based upon information obtained in real time from the

measurements becoming available, so that the filter is

continually "tuned" as well as possible. Such an algorithm

is often termed an adaptive or self-tuning estimation

algorithm [45-48,50,51]. Several authors use the

information found in the innovation sequence to estimate the

unknown noise covariances [45-48,50].


Applications

The earliest applications of the Kalman filter dealt

with satellite orbit determination, tracking and navigation











problems. Tapley and Ingram [52] used an extended Kalman

filter to estimate the state and the unmodeled accelerations

acting on an orbiting space vehicles. Cambell et al. [53]

described the use of the filter for orbit determination for

the Voyager spacecraft during its Jupiter fly by. Dawn and

Fitzgerald [54] presented the application of a Kalman filter

in numerous phased array radars to track satellites, reentry

vehicles and missiles. Application of Kalman filtering to

spacecraft range residual prediction was studied by Madrid

and Bierman (55]. They developed a U-D covariance factored

Kalman filter to be able to evaluate the validity of range

measurements to a distant spacecraft in near-real time so as

to be able to detect receiving-station hardware problems as

they occur. The Kalman filter is admirably suited to the

problem of monitoring measurement consistency and the

related statistics. Additional study was done on several

versions of the extended Kalman filter which is used to

estimate the position, velocity and other key parameters

associated with maneuvering reentry vehicles [56]. The

Kalman filter is also well suited for application to the

problem of anti-aircraft gun fire control. Berg [57] made

use of the Kalman filter theory to develop an accurate,

numerically efficient scheme for estimating and predicting

the present and future position of maneuvering fixed-wing

aircraft. Application of modern estimation techniques to

multisensor navigation systems began in the mid-1960's.











Because the errors in a typical navigation system propagate

in essentially a linear manner and linear combinations of

these errors can be detected from external measurements, the

Kalman filter is ideally suited for their estimation

[58-61]. Application of Kalman filtering on ship motion

prediction and position control [62-64] have many features

in common with the tracking and navigation problems. The

modeling of ship/wave interactions requires the development

of simple, but effective, models using signal processing or

identification techniques. The authors are concerned with

the Kalman filter performance in the context of a control

system. For remote sensing applications [65,66] the Kalman

filter appears as a part of a signal processor that is used

for offline analysis. These applications are not concerned

with real-time processing. Nonetheless, simple models are

used to accomplish the processing, and the performance is

based on the achievement of the general goals of the system.

Estimation and control problems abound in industrial

processes and power systems [67-71]. These applications

impose the necessity for real-time operation, generally with

limited computational capability, using poorly defined

models of the process. Simplicity and adaptability are

basic concerns in these applications. Application of Kalman

filtering has been extended to the geophysics area. A

Kalman filter was developed as a white-noise estimator for

seismic data processing in oil exploration [72]. The Kalman











filtering approach was used to obtain optimal smoothed

estimates of the so-called reflection coefficient sequence.

The sequence contains important information about subsurface

geometry. Ruckebusch [73] described an application of

Kalman filtering to a geophysical subsurface estimation

problem. Other areas where the theory of optimal estimation

has applications are tracer studies in nuclear medicine,

statistical image enhancement, estimation of river flows,

network and load forecasting [74,75], classification of

vectorcardiograms and demographic models [76].



Logarithmic Number System

In 1971 Kingsbury and Rayner [77] first outlined logic

hardware and software function approximations to the

addition of two positive numbers in logarithmic arithmetic.

Swartzlander and Alexopoulos [78] concentrated on ROM based

hardware for the fastest addition but with a wordlength

limit of 12 bits. Later Lee and Edgar [79] developed

efficient algorithms to implement 8 and 16 bit logarithmic

arithmetic on microprocessors. They established that unlike

floating-point arithmetic, which provides accuracy at the

expense of speed, the LNS provides for both and is suitable

for advanced DSP applications. Kurokawa et al. [80]

applied LNS to the implementation of digital filters and

demonstrated that it gives filtering performance superior to











that of a floating-point system of equivalent wordlength and

range. Similar observation was made by Swartzlander et al.

[81] when applying LNS to Fast Fourier Transform processor.

The general conclusion of these studies is that addition,

subtraction, multiplication and division are very fast and

easy to implement in LNS. The advances in semiconductor

memory technology have renewed interest in the LNS over the

last few years. A brief survey of the literature in the LNS

will be presented.



Logarithmic Conversion

Multiplication and division operations in computers are

usually accomplished by a series of additions, subtractions

and shifts. Consequently, the time required to execute such

commands is much longer than the time required to execute an

add or subtract command. In an early work by Michell [82] a

method of computer multiplication and division is proposed

which uses binary logarithms. The logarithm of a binary

number is determined approximately from the number itself by

simple shifting and counting with no table lookups required.

Simple add (or subtract) and shift operations are all that

is required to multiply or divide. However, since the

logarithms used are approximations of the actual logarithms,

the operations yield errors. The paper discusses a method

of finding approximate base-two logarithms and an analysis

of maximum errors that may occur as a result of











approximations. A more vigorous analysis of errors in the

multiplication and division operations using binary

logarithms is presented by Hall et al. [83]. In this paper

the authors discuss algorithms for computing approximate

binary logarithms, antilogarithms and applications to

digital filtering computations and establish that

logarithmic techniques are well suited for parallel digital

filter banks and multiplicative digital filters. Marino

[84] presents a parabolic approximation method for

generation of binary logarithms with an increased precision

over Hall et al. [83], where a review base-two logarithm

and antilogarithm computation using a fixed number of

iterations is presented. Their method is claimed to yield

an improved accuracy and can be implemented as a subroutine

or using a hardware peripheral device.


Addition and Subtraction Algorithms

Kingsbury and Rayner [77] were the first to propose two

methods of adding or subtracting logarithmically encoded

numbers. In the direct method of addition and subtraction,

they use approximate evaluation of logarithms. However,

their second method is based on table lookups using ROM

which leads to potentially faster operations. The paper

also indicates a possible method of table reduction. In

1975, Swartzlander and Alexopoulos [78] unaware of the

previous work by Kingsbury and Rayner, proposed a











sign/logarithm system along with arithmetic algorithms

identical to those in [77]. Their paper suggests hardware

architectures for arithmetic units and includes comparison

of speeds with conventional arithmetic units. Ironically,

in a paper by Lee and Edger [79], the LNS was reinvented a

third time but enriched with ideas for implementation in 8

and 16 bit microcomputers. A more rigorous treatment of LNS

addition and subtraction algorithms with respect to storage

efficiency, accuracy, range, noise and array area is found

in Lee and Edger [85].


Multiplication and Division Algorithms

In LNS, multiplication and division are achieved by

simply adding or subtracting logarithms and are thus

straightforward operations. These algorithms along with

examples and hardware implementations may be found in

references 77-79 and 85.


Table Reduction Techniques

Current algorithms for addition and subtraction

operations in LNS require a table lookup operation. The

size of the table is governed by the wordlength size. A

16-bit arithmetic would require 64K of storage. The current

state of the art in high speed ROM and RAM is such that LNS

is useful for wordlengths up to 12-13 bits for high speed

operations and 16 bits for lower speeds. Due to the











discrete coding and nature of the functions filling lookup

tables, a large portion of the values is zero. The number

of zeros depends on the radix used in the number system and

the number of fractional bits allotted in the word format.

Kingsbury and Rayner [77] suggested two ways of reducing the

storage capacity by employing read and compare successive

approximation process and linear interpolation technique.

Lee and Edgar [85] determined a cutoff value for the address

below which the table output is zero.


Extended Precision LNS

The aforesaid address space limitations of contemporary

high speed memories limit the admissible logarithm wordlengh

to 12-13 bits. To overcome this problem, Taylor (86,87]

suggested a modified interpolation technique to be useful in

increasing the wordlength of practical logarithmic

arithmetic unit without adding any significant hardware

burden.


Digital Filtering Using LNS

An early work by Kingsbury and Rayner [77] implemented

a second order recursive lowpass filter with a 16 bit

logarithmic arithmetic and demonstrates the improvement in

dynamic range and performance compared to 16 bit fixed point

arithmetic. Kurokawa et al. [80] recently applied

logarithmic arithmetic to the implementation of digital











filters and present an analysis of roundoff error

accumulation in the direct realization of logarithmic

digital filter and establishes that LNS gives filtering

performance superior to that of a floating-point system of

equivalent wordlength and range. Statistical roundoff error

behavior of a single logarithmic arithmetic operation may be

found [80]. In another paper by Sicuranza [88], preliminary

results on two dimensional filters implemented with LNS are

presented.


FFT Implementation with LNS

Swartzlander et al. [81] presented a FFT

implementation scheme with LNS and established that this

approach resulted in an improved error performance for a

given wordlength over FFT's implemented with conventional

fixed or floating point arithmetic.


Other Applications with LNS

Swartzlander and Gilbert [89] concluded that

convolution units using LNS arithmetic have a figure of

merit twice that of two's complement merged arithmetic and

tens of orders higher than high speed modular array.

Besides, they claim that the advances in VLSI technology

offer a strong potential for LNS applications. Taylor

[86,87] presented an approach to a fast floating point

arithmetic unit employing LNS. His approach does not








20

require exponent alignment, supports high speed addition and

multiplication and admits a simple VLSI realization.


















CHAPTER THREE
KALMAN FILTERS






Statement of the Problem


Application of optimal estimation is predicated on the

description of a physical system under consideration by

means of mathematical models. Early work in control and

estimation theory involved system description and analysis

in the frequency domain. Recent works have involved system

descriptions in the time domain. The formulation used

employs state-space notation which offers the advantage of

mathematical and notational convenience. Moreover, this

approach to system description is closer to physical reality

than any of the frequency-oriented transform techniques. It

is particularly useful in providing statistical description

of system behavior. Initially, the work done in the time

domain was concerned with continuous systems, but recently

work was extended to the discrete case in which information

is available or designed only at specified time intervals.












With the use of digital computers or microprocessors,

discrete systems can be simulated easily and computation

burden is avoided considerably.

A discrete stochastic dynamic system can be represented

as


xi = Fxi- + wi (3.1)



Z. = Hxi. + v. (3.2)



Where

i -0,1,2 ......

x. = nxl state vector

F = nxn state transition matrix (constant)

w. nxl vector of Gaussian input white noise

z. mxl output vector

H = mxn output matrix (constant)

v. mxl vector of Gaussian measurement errors (white)

The above equations are illustrated in Figure (3-1).




Least Square Estimation


The linear least-squares problem involves using a set

of measurements z, which are linearly related to the unknown

quantities x, by the expression


z Hx + v



































DISCRETE SYSTEM


MEASUREMENT


Figure 3-1. A Discrete Stochastic Dynamic System.











where v is a vector of measurement "noise." The goal is to

find an estimate of the unknown, denoted by x. In

particular, given the vector difference


z Hx (3.3)

^
we wish to find the x that minimizes the sum of the squares

of the elements of z-Hx. The vector inner product generates

the sum of squares of a vector. Thus, we wish to minimize

the scalar cost function J, where


J = (z-Hx) (z-Hx) (3.4)



Minimization of a scalar, with respect to a vector, is

obtained when



S= 0 (3.5)
ax


and the Hessian of J is positive semidefinite.

Differentiating J and setting the result to zero yields



H Hx H z (3.6)



It is readily shown that the second derivative of J

with respect to x is positive semidefinite; and Equation
T
(3.5) does indeed define a minimum when H H possesses an

inverse (i.e. when it has a non-zero determinant); the











least-squares estimate is

T -1 T
x (HTH) H z (3.7)





Estimation of Parameters Using Weighted Least-Squares


A more optimal way to estimate the state vector x is

the weighted least-squares estimate. For this estimate we

choose x, as the value of x that minimizes the quadratic for


J=.5[(x-x') TM- (x-x')+(z-Hx)TR-l(z-Hx)] (3.8)


M- and R- are the inverse matrices of the expected values

of (x-x')(x-x')T and (z-Hx)(z-Hx)T respectively. Note that

x' is an estimate of the state before the measurements were

made. With this choice of weighting matrices, the weighted

least-squares estimate is identical to the

conditional-expected value estimate assuming Gaussian

distributions of x and v.

To determine x, consider the differential of Equation

(3.8)


dJ-dxT[M-1 (x-x')-HT R1 (z-Hx)] (3.9)


In order that dJ-0 for arbitrary dxT the coefficient of dxT

in Equation (3.9) must vanish











-1 T -1 -1 T -1
(M +H R H)x=M x'+H R z

-(M -I+HTR- H)x'+HTR- (z-Hx')

or
T -1
x x'+PH R (z-Hx')

where

P = M-I+HTR-1 H (3.10)



The quantity P in Equation (3.10) is the covariance matrix

of the error in the estimate x; that is, we have

T
P=E[(x-x)(x-x)1] (3.11)



Since M is the error covariance matrix before

measurement, it is apparent from Equation (3.10) that P, the

error covariance matrix after measurement, is never larger
T -1
than M, since H R -H is at least a positive semidefinite

matrix. Thus, the act of measurement, in the average

decreases (more precisely, it never increases) the

uncertainty in our knowledge of the state x.

Another noteworthy property of the estimate is the fact

that
T
E[(x-x)x] = 0

that is, the estimate and the error of the estimate are

uncorrelated. In this case where x and v are Gaussian, this

implies that x and (x-x) are independent; in other words no
A A
improvement in (x-x) can be obtained by knowledge of x or z.











Recursive Filters


A recursive filter is one in which there is no need to

store past measurements for the purpose of computing present

estimates.

Consider a discrete system, whose state at time tk is

denoted by x(tk) or simply Xk, where wk is a zero mean white

sequence of covariance Q.


xk FXk-1 + Wk-l (3.12)


Measurements are taken as linear combinations of the system

state variables, corrupted by uncorrelated noise. The

measurement equation is written in vector-matrix notation as


zk = Hxk + vk (3.13)



where, zk is the set of m measurements at time tk, namely,

Zlk,'z2k ...... Zmk arranged in a vector form. In addition, H

is the measurement matrix at time tk; it describes the

linear combinations of state variables which comprise z k in

the absence of noise. The dimension of the measurement

matrix is mxn, corresponding to m-dimensioned measurements

of a n-dimensioned state. The term vk is a vector of random

noise quantities (zero mean, covariance R) corrupting the

measurements.

Given a prior estimate of the system state at time tk,

denoted xk(-) we seek an update estimate, xk(+) based on use









of the measurement, zk. In order to avoid a growing memory

filter, the estimate is sought in the linear recursive form


xk(+)=KXk(-)+KkZk (3.14)


where K. and Kk are time-varying weighting matrices, as yet

unspecified. Although the following deviation is for an

assumed recursive, single stage filter, the result has been

shown to be the solution for a more general problem. If wk,

vk, are Gaussian, the filter we will find is the optimal

multi-stage filter; a nonlinear filter cannot do better [2].



Discrete Kalman Filters

It is possible to derive the Kalman filter by

optimizing the assumed form of the linear estimator. An

equation for the estimation error after incorporation of the

measurement can be obtained from Equation (3.14) through

substitution of the measurement Equation (3.13) and the

defining relations


Xk(+)-Xk+Xk(+) (3.15)


Xk(-)=Xk+Xk(-) (3.16)


The result is









Xk (+)-=(K+KkH-I)Xk+K(Xk(-)+KkVk (3.17)


By definition E(vk]=0. Also if E[(xk(-)]-0, this estimator

will be unbiased for any given state vector xk only if the

term K'+KkH-I is zero. Thus we require


K' = I-KkH

and the estimator takes the form


Xk(+)=[I-KkH]xk(-)+Kkzk


or alternatively
^ ^ ^
Xk(+)=Xk(-)+Kk[Zk-HXk(-)] (3.18)


The corresponding estimation error is from Equations (3.13),

(3.15), (3.16), and (3.18)


Xk(+)=[I-KkH]xk(-)+Kkvk (3.19)


Using Equation (3.19) the expression for the change in the

error covariance matrix when a measurement is employed can

be derived. From the definition

~ ~ T
Pk(+)=E[xk(+)xk(+)T]


Equation (3.19) gives

Pk+)EI-KH)Xk-)Xk-T T T T]
P k(+)=E[(I-K kH)x k(-)[x k(-) (I-K kH) +V k K k ]+









+KkV k[xk(-)T(I-KkH) +VkKk ]}

By definition

T
Pk(-)=E[xk(-)Xk(-)]

Rk=E[VkVk TI

and, as a result of measurement errors being uncorrelated

T T
E[xk(-)vk l=E[vkxk(-)T]=0

Thus

Pk(+)=[I-KkH]Pk(-)[I-KkH]T+KkRKkT (3.20)



The criterion for choosing Kk is to minimize a weight

scalar sum of the diagonal elements of the error covariance

matrix Pk(+). Thus, for the cost function we choose

Jk=E[xk(+)TSxk(+)] (3.21)


where S is any positive semidefinite matrix. The optimal
estimate is independent of S; hence, we may as well choose

S=I, yielding


Jk =Trace[Pk(+)]











This is equivalent to minimizing the length of the

estimation error vector. To find the value Kk which

provides a minimum, it is necessary to take the partial

derivative of Jk with respect to Kk and equate it to zero.

Use is made of the relation for the partial derivative of

the trace of the product of two matrices A and B (with B

symmetric),



[Trace(ABAT)] = 2AB
8A

From Equations (3.20) and (3.21) the result is

-2[I-KkH]Pk(-)H T+2KkR=0

Solving for Kk,


T T -lP-^ Rr (3.22)
k =Pk(-)H T[HPk(-)H T+R]-1 (3.22)


which is referred to as the Kalman gain matrix. Examination

of the Hessian of Jk reveals that this value of Kk does

indeed minimize Jk"

Substitution of Equation (3.22) into Equation (3.20)

gives, after some manipulation,

Pk(+)=Pk(-)-Pk(-)HT [HPk(-)H T+R]- HPk(-)


=[I-KkH]Pk(-)










which is the optimized value of the updated estimation error

covariance matrix.

Thus far we have described the discontinuous state

estimate and error covariance matrix behavior across a

measurement. The extrapolation of these quantities between

measurements is


xk(-)=Fxk-l(+)


Pk(-)=FPkl(I(+)FT +Q


The equations of the discrete Kalman filter are

summarized in Table (3-1). Nevertheless the application of

these equations requires knowledge of the input and output

noise covariances statistics, Q and R respectively. Figure

(3-2) illustrates these equations in block diagram form.

The Kalman filter to be implemented appears outside the

dashed-line box. In the linear discrete Kalman filter,

calculations of the covariance level ultimately serve to

provide Kk, which is then used in the calculation of mean
^
values (i.e. the estimate xk). There is no feedback from

the state equations to the covariance equations.











Table 3-1. Summary of Discrete Kalman Filter Equations.


System Model xk =Fxk-l+Wk-'I Wk~N(O,Q)

Measurement Model z k =Hx k+V, v~N(0,R)

Initial Conditions E[x0]=x0, E[(x0-x0)(x0-X0) T]-P0

Other Assumptions E[WkVj T=0 for all j,k

Extrapolation

State Estimate xk (-)=Fxkl(+)
T
Error Covariance Pk(-)=FPkl (+)F +Q


State Estimate Update xk(+)=xk(-)+Kk[zk-Hxk(-)]

Error Covariance Update Pk(+)=[I-KkH]Pk(-)

Kalman Gain Matrix Kk=Pk(-)HT[HPk(-)HT+R]-1















MATHEMATICAL MODEL
DISCRETE SYSTEM I _7 EASuFEMENT


DISCRETE KALMAN FILTER


xk


L--- __---------J


Figure 3-2. System Model and Discrete Kalman Filter.











Kalman Filters with Deterministic Inputs


When the system under observation is excited by a

deterministic time-varying input, u, whether due to a

control being intentionally applied or a deterministic

disturbance which occurs, these known inputs must be

accounted for by the optimal estimator. The dynamic system

containing deterministic inputs can be described as


xi Fx i-1 + wi_1 + u_1

Z. = Hx. + v.


In order for the estimator to be unbiased the state estimate

update is given as follows:


xk(+)=Fxkl(+)+ukl+Kk[Zk-HFx kl(+)-Hukl] (3.23)



By subtracting x from x, it is observed that the resultant

equations for x are precisely these obtained before. Hence,

the procedures for computing P,K, etc., remain unchanged.


















CHAPTER FOUR
LOGARITHMIC NUMBER SYSTEM





Various number systems have been in use for signal

processing. Traditionally, when a large dynamic range and

high precision are required, data are given a floating point

representation. The problem with floating point arithmetic

is that operations are slow compared to their fixed point

counterparts. On the other hand, fixed point arithmetic

offers high speed in the expense of dynamic range and

precision. Recent developments in the area of logarithmic

number systems (LNS) have proven it to offer some

significant advantages over other number weighted systems

like the fixed or floating point. The main advantages of

LNS are the extended dynamic range on which they operate and

the high speed of operations accompanied by a remarkable

regular data flow.









LNS Representation

In a floating point system, a real number X can be

approximated as

+e'
X-m r ; i/r

where r is the radix, mx is the unsigned M-bit mantissa and

e' is the unsigned E-bit exponent. In an LNS environment, a
x
real X is represented as
+e
x)
X-+r ; exzex+x; Ax=logr(mx)

where in practice X is coded as a (N+2)-bit word. The first

two bits represent the sign of X. The N-bit exponent ex is

represented as a N-bit word consisting of F fractional bits

and N-F=I integer bits.

By distributing the appropriate number of bits for the

integer and fractional parts, both large dynamic range and

high precision can be achieved. More specifically, this

system is characterized by


lxmax=r 2-(2 -l)=r2



-2- F(2N-I)~r-2I
lximin r N r



2-F+l(2N_1) 2 I1+1
RA=r -r









where Xl max and IXlmin are the largest and smallest

positive numbers in LNS, respectively. In addition, RA is

the range of this system (the ratio of the largest to the

smallest number).





LNS Arithmetic Operations


Following the previously adapted number representation,

arithmetic in LNS is described below. For the fastest

execution of operations like addition and subtraction, table

lookup techniques have been employed to support them.



Multiplication

In LNS, the logarithm of the product of two numbers is

given by the sum of their logarithms. Specifically,

multiplication in this system is performed in the following

manner.
e eb e
For A=r a, B=r C=r c, S=sign bit


C=AB; ec4ea+ea; Sc=SaOeb


where 9 denotes an "exclusive-or" operation.

Overflow occurs when the magnitude of the result

exceeds the largest number which can be represented for the

given word format. Overflow and underflow can be easily











detected, by a simple comparison with the largest and

smallest numbers representable in the system.



Addition and Subtraction

Addition and subtraction in this number system are

defined in terms of the identity A+B=A(1+B/A). In

particular, addition and subtraction in LNS are performed in

the following manner.
ea eb ec
For A-r B-r ,C-r S-sign bit

and assuming without loss of generality that A>B



C=A+B; ec<-ea+r (v); v=eb-ea; tv=logr(l+rV); Sc=Sa



C-A-B; ec-ea+er0 (v); v=eb-ea; ev=logr(l-rV); Sc=Sa


The realization of t and e is achieved through the use

of table lookup operations from high-speed, high-density

RAMs or ROMs. Here, A N-bit address would be presented to a

memory unit and the value of t and 9 would be obtained as

direct table lookup. Based on contemporary high-speed

memory chip limitations, the wordlength of v can be

estimated to be on the order of 12 bits. PLAs have been

also used in place of ROMs in order to reduce table

requirements.










A Hybrid Floating Point Logarithmic Arithmetic Unit


When a large dynamic range and high precision are

required, data must often be given a floating point

representation. Although standardization is on the way,

there remains today many floating point formats. The

problem with floating point arithmetic is that operations

are slow and complex. Furthermore, floating point addition

can cause a special problem. The time it takes to perform a

floating point addition can be markedly dependent upon the

relative values of the data to be added. As a result,

developing efficient real time code in floating point is

difficult. In addition, the multiplier and addition paths

are sufficiently different so as to demand two separate data

paths and hardware subsystems. As a result, the utilization

rate of the hardware floating point unit can be as low as

50%.

Below we present a variation of the floating point

arithmetic unit. It incorporates the precision and dynamic

range of the floating point system while eliminating the

disadvantages of high overhead and reduced throughput due to

the exponent alignment requirement. This new concept shall

be referred to as the Florida University Floating (Point)

Unit or (FU)2 [87].











Floating Point Format

In a floating point environment, a real number X can be

approximated by Equation (4.1). Arithmetic in the

floating-point system is performed in stages.

Multiplication is a multistep process given by

1. Multiply mantissas and add exponents.

2. Post-normalize mantissa and round result.

3. Adjust exponent if required.


In a commercial floating-point adder/subtractor unit,

up to one-third of an arithmetic cycle can be consumed in an

exponent alignment process. In addition, the length of time

required to complete an exponent alignment is data dependent

and therefore variable. For example, in a LSI-11, exponent

alignment can take up to 7 us with a basic mantissa add time

of 42 ps. Formally, in a floating-point add, the following

steps are taken:

1. Align exponents and shift mantissas accordingly.

2. Add mantissas.

3. Post-normalize mantissa and round result.

4. Adjust exponent if required.


2
(FU) Unit

The floating point logarithmic processor, shown in

Figure (4-1), consists of a floating point to LNS encoder, a

LNS arithmetic unit, and a LNS to floating point decoder.










2
Several key observation can be made about the (FU) 2 Unit.

They are

1. No exponent alignment is required for addition.

2. The multiplier is a subset of the adder.

3. Independent of the values of X and Y, all adds will
occupy the same amount of time.

2
As a result, both a fast and low overhead (FU) unit

can become a reality if the mapping denoted f and Y can be

mechanized efficiently as a table lookup operation. For

r-2, the computed address v-b-a can be presented to a

preprogrammed memory unit which responds with the value of

*(v)-log2(l+2V ). In order to be consistent with the
addressing space limitation found in modern memory units, v

is currently limited to about 12 bits of precision. Taylor

was able to extend this figure using a linear interpolation

scheme [86]. Lee and Edgar also introduced a data

compression policy which ignored all addresses which map,

under 4, into zero [79]. Along a similar line of reasoning,

it can be argued that the mapping Y and T~ can also be

performed as a table lookup operation.

The LNS arithmetic unit contains both an adder and

multiplier. This is unlike a conventional floating point

system where the adder and multiplier paths are separate.

Nested multiply and add operations can be interleaved among
2
several (FU) units. In a conventional floating-point unit

with distinct multipliers and adders, a utilization factor








FLOATING POINT TO LNS
ENCODER
7F-7 hJ


LNS ARITHMETIC UNIT
ADD


m
MULT
II i

x=M e' ADDITION
x j PATH I FLOATING POINT TO LNS
II DECODER
SI I ^ TADD


ADD -MULTI
. I ^* ~~MULTIPLICATION |'I --
UI MULTY'


y=m e
I |
I I
e'
Y L___ ------J



2
Figure 4-1. (FU) System Architecture.






44

2
as low as 50% can occur. The (FU) units can be made to

internally switch data paths so as to support both

operations (i.e. 100% utilization). This can be a very
2
important property if the (FU) is to be used to support

systolic or other data path determined architectures.


















CHAPTER FIVE
ERROR MODELS FOR LNS ARITHMETIC





The implementation of filters with digital devices

having finite wordlength introduces unavoidable quantization

errors. There are three main sources of quantization error

that can arise: input quantization, coefficient

quantization and quantization in arithmetic operations. The

development of quantization error models for LNS arithmetic

will be presented with special emphasis on input

quantization. Once a model to represent input quantization

error has been developed, models to represent the other two

types of error can be easily be obtained.




Input Quantization


A quantizer can be viewed as a nonlinear mapping from

the domain of continuous-amplitude inputs onto one of a

countable number of possible output levels. In this error

analysis the input is presented as an infinite precision

floating point number and it is mapped to LNS with the use









of a floating point to LNS encoder, illustrated in Figure

(4-1). The only source of error in forming ex is the one

resulting from the logarithmic mapping m -x logr mx which is

performed as a memory table lookup operation. The proposed

error analysis model is illustrated in Figure (5-1). Here,

the two paths associated with the formation of Ax and its

finite wordlength approximation are shown as parallel paths.

Upon receiving the real mantissa m x, the lower path provides

the ideal mapping Ax=log rmx while the upper path consists of

an input quantizer (QI, which provides a discrete value of
*
mx ), an ideal mapping of mx into L=logr mx and finally an

output quantizer (Qo, which provides the machine version of
*
L, namely Ax -RND(L).

The input and output quantization errors can be defined

as

m* m* *
EI=mx-mx ; E =log r(m x *)-[log r(m x )]


where E and E are uniform white noise sequences possessing

the following statistical properties:


E[EI]=E[Eo]=0; E[EoEoj]=(q 2/12)8



*2
E[mxE]=E[(log (m* ))Eo]=0; E[EiEij]=(q2/12)6k

-NI -NO
Here, qI=2 and qo=2N where N and N are the numbers

of bits available at the input and output of the logarithmic

























































Figure 5-1. Input Quantization Error Model.











mapping L. In addition we assume that mx is uniformly

distributed over [l/r,l). The final error metric E is

parameterized by D=EI/m x, as


E = Ax -Ax = logr (m x+E I)+Eo-logrm = logr(I+D)+EO


Through the application of the theory of functions of random

variables [90], the probability density function (p.d.f.)

fD(D) is shown to be given by

f(D)= mJf (Dmm) dm; -q < D < q
f D(D)= Im f ~Em (Dm x 'm x) dm X; -q1 D -wM


Direct evaluation of f D(D) results in


q, 1
4 4q,


3
4q1


for -qI < D < -qi/2


for -qi/2 < D < qi/2


q, 4 1
4D 4q,


for qi/2 < D < qI


Defining P-log r(l+D), then D=r P-1, and fp(P) satisfies

fD(D) p p
f p(P)=- =rp lnr fD(D)=rP lnr fD(r P-l); lnr-log r

IDID


fD(D) =











Assuming that P and E are statistically independent random

variables, the p.d.f. of the final error E-P+E0 can be

computed directly using


fE(E)-_ fp(P) fE (E-P) dP (5.1)


where the evaluation of (5.1) depends on the relative

magnitude of the input and output quantization step.

Consider the case where the input quantizer offers

lower precision than the output quantizer can provide (i.e.:

NIqo); then fE(E) is given by

logr(1-qi/2) log r(l+qi/2)

fE(E)= f A(P) fE (E-P) dP + f B(P) fE (E-P) dP +

logr(l-qI) logr(l-qi/2)


logr(l+qI)

+ f A(P) f E(E-P) dP (5.2)
"0 (5.2)
logr(1+q,/2)

where


A(P)-lnr rP( qI- 1 B(P)=(3/4q.) lnr
4(r P-1)2 4q


A direct evaluation of (5.2) results into a complicated

piecewise continuous expression consisting of seven

subintervals. The final error density function having a

trapezoidal shape is represented by the solid line in Figure











(5-2). Two distinct first order approximations, one for

each side, can be used to simplify the expression. The

resulting approximated p.d.f. of the final error is plotted

in Figure (5-2) outlined by the dotted line. In addition

experimental results, outlined by the dashed line, verify

the theoretical analysis.

Similar results can be achieved considering the other

case, where the input quantizer offers higher precision than

the output quantizer (i.e.: NI>N0 or q(
the fact that qo is the dominating factor, an investigation

was made concerning the necessary and sufficient conditions

under which the final error p.d.f. coincides with the

output quantization error p.d.f. with a minimum error. In

order words, how many more bits of precision in the input

quantizer will bring the two errors at comparable levels?

As a basic criterion of comparison the mean square error

criterion was used.

The mean square error of the output quantizer error is
2 2
o2qo2/12 and let e represent the maximum allowable
percentage error. Then a relationship between the number of

input and output bits can be found to satisfy the equation



E[[(log (m ))-log (m )]2] E[[log (m *)-log (m ) 2I -
r x r x r x r x 0


2 2
or E[P 2-Cao, where P has been defined previously as

P-logr (1+E I/m x). Using (5.2), then





































-- Theoretical

........ Approximation

--- Experimental


Figure 5-2. Final Error Probability Density Function.









55 4 4 2,3
,2,- A Y2Yl ) 2yI(Y2-yI) YIY2-Yl)
E[P ] 2 A 2 4 + 13- ]+

55 4 4 4 2 3 3
(x15x2) 2xl(xl-x2)4 +Xlx-X2)
+ 2 2 1- 2 + 1 1-2--+

+ r [z1(lnz1)2-z2(lnz2) -2z1lnz1+2z2lnz2+2zl-2z2]

where
3 3 (3 3
4( q-,"-7-- ) lnr(-q--S-) 3
A- 2 ; Q- -2-; r=- )24-I
(y2-Yl) (x2-x1) (lnr) 4q

Yl-logr(l-qi); xl1-logr(1+qi); z1-(l+qi/2)

Y2-log r(l-qi/2); x2-logr( l+qi/2); z2-(l-qi/2)


The resulting theoretical mean square error of the final
error is shown in Figure (5-3), designated by stars, for 4,
8, 12 output bits. In addition, experimental results
designated by circles verify the theoretical ones. Finally,
the straight lines represent the mean square error of the
uniform output quantization error. As it can been seen from
Figure (5-3), three extra bits of precision in the input
quantizer are required to make the final error density
function comparable to the output quantization error density
function. In addition, one bit of precision in the input
quantizer is gained since the leading bit in the mantissa is
always one, thus reducing the number of extra bits to two.
Using this assumption the effect of QI is ignored and each



































4 OUTPUT BITS


_-"-- 4.*-< ,l-4w)-- )-..- -

I I I I I
8 1O 12 14 16
NUMBER OF INPUT BITS


18


12 OUTPUT BITS


Mean Square Error of Final and Output Error.


10E-?-=

1OE-8-=-


- N i


J."J~~~ ~~~ 4- /- '',.. ,) I,- .- ..1 4 -.> -. -


10E-4 =


Figure 5-3.










infinite floating-point measurement can be interpreted as a

LNS word with a fractional exponent as follows:
*
ez =r ; ee+Ax ; Ax =Ax+w; Ax=log (ms)
I x r* x

where w is a uniform white sequence with the following

statistical properties:


E[W]=0; E [WkWj]=O; =(qo/12)kj; E[xkwk]=0


The simplified input quantization model is illustrated in

Figure (5-4).




Coefficient Quantization

ea
Any coefficient a=r is represented as a quantized

coefficient a when represented in a finite wordlength

register. The error Aea=e -e is inherently deterministic
a a a
in nature. Thus we model the quantized coefficient as

ea=ea+Aea. It should be noted that |&ea|
the quantization step.




Quantization in Arithmetic Operations


We assume throughout that no overflow occurs during any

LNS arithmetic operation. Multiplication does not introduce

quantization errors. However, during addition and

subtraction quantization errors are unavoidable due to











memory wordlength limitations. The addition/subtraction

quantization process can be adequately represented using the

model in Figure (5-5). Here e and e represent the finite
y
wordlength LNS exponents to be added. Let the symbol

represent logarithmic addition. For the operation of

addition the following sequence of events takes place. For

Z=X+Y and assuming w.l.o.g. that e >ey it follows that

e e e
Z=rz
Z .r ez = r +r Y

or


ez eey = exey x+(v); f(v)=logr (l+rV); v=ey-ex


The quantized addition e is obtained by rounding the

fractional part of e and can be modeled as
z


ez = (exOe )q = (exOe )+w = (ex+w)(e +w)


where w is a white sequence uncorrelated with ex and ey.





















/_A| I i/ i __ .JLe
m x O g l og ( m ) + / \ - ---\ e *
xv T x u\y^^z
--r x + T



ey
y




Figure 5-4. Simplified Input Quantization Model.











ex 0



v e


e 0 e*
y


Figure 5-5. Addition Quantization Model.


















CHAPTER SIX
LOGARITHMIC KALMAN FILTERS






LNS Kalman Filters


A discrete stochastic dynamic system

deterministic inputs can be represented as


containing


Xk+1 Fx k + uk + w k


(6.1)


Z Hxk + vk


Given the a priori estimate of the initial state x0 and

the state covariance P0(-) and given the a priori

statistical information of the input noise covariance Q and

output noise covariance R, an estimate of the state of the

system defined by (6.1) is obtained sequentially for

k-1,2,3... for the standard Kalman filter


Xk+1 Akk + uk + Ckzk, 0


(6.2)


where










Ak-(F-FKkH); Ck= FKk


Further Ak and Ck are assumed to be precomputed with a

high-precision computer. If the arithmetic operations in

(6.2) are implemented using theoretically ideal

infinite-precision arithmetic, then the mean and the

variance of the estimation error are well known to be

respectively,


E(Xk-xk] 0 (6.3)


T~^xXx )1 P()(6.4)
E((xkXk)(xk-Xk) P k(-) (6.4)



The LNS implementation of the Kalman filter (6.2) can

be presented as
^ ^ ^
exk+l=(eAkOexk)euk(eCkezk), ex0 (6.5)



where a and designate LNS matrix multiplication and LNS

vector addition, respectively.




Theoretical Error Analysis in LNS


We consider here the effect of LNS roundoff arithmetic

since in a digital environment real numbers must be replaced

by their finite wordlength equivalent. In this work, the

(N+2)-bit sign magnitude LNS data format described in

Chapter Four will be employed. We will assume that the











number of integer bits, I, is properly selected so no

overflow will occur during various arithmetic operations.

When the filter (6.5) is implemented on a finite-precision

machine, the standard performance measures (6.3) and (6.4)

will be altered because of input, coefficient, and

arithmetic quantizations. We predict these anticipated

degradations under the assumption that the models

representing various quantizations, as developed in Chapter

Five, are correct.
^ ,
Let ex k+l denote the LNS actual estimate of the state
*
Xk+1 given the quantized measurements ez0, ezl,...,ezk; this

is the estimate that is actually computed by the

finite-precision implementation of (6.5). By replacing each

finite precision operation in (6.5) by a quantizer model,

the following difference equation for the actual state

estimate is obtained.


exk+l[ (eAxkex) +eBkI{[eu k+CAk+eBk]



[(eC*Oez)*+eAk+cBk]}, ex0




ez k ez k + ez k


where cAk and cBk are zero-mean white noise sequences of

dimension n, which represent the LNS addition error

sequences. Furthermore, ez k is also zero-mean white noise

sequence of dimension m, which represents the input











quantization error. The effect of quantized coefficients

can be represented as


*
eAk = eAk + AeAk;



eCk = eCk + AeCk;


euk = euk + Aeuk



ex0 = ex0 + Aex0


Now define the error between the actual and the

theoretically ideal estimate as


ek = xk xk (6.6)



where


^*
exk+1 *
xk+1 = r = YlkAkY2kxk + Y3kUk + Y3kCkY4kZk


(6.7)


The LNS arithmetic errors are represented in the diagonal

matrices Ylk' Y2k' Y3k of dimension nxn and Y4k of dimension

mxm. Specifically


Bdiag(r ,lk B2k eBnk
Ylk = diag (r ,r ", .... r n )


niAx ni Ax nI cAx CA
Y2k=diag(r il lik ril 2ik r i=2 3ik,......r n(n-l)k)


Y3k = diag (r


Bk +Alk B 2k +A2k
r


.Bnk +Ank


-I ++mrEIn-i T1c
YZlk+i4l CZlik LZ2k i .1 CZ2ik
Y4k=diag(r ,r









m-i
rZ3k i=2 3ik


.r .z +cAm(m-l)k)


Subtracting (6.2) from (6.7) and use of (6.6) yields


ek+l 1 YlkAkY2kek + [YlkAkY2k-Ak]xk +


+ Y3kuk uk + [Y3kCkY4k-Ck]Zk


(6.8)


We define an auxiliary augmented state YkER n such that

y.(e k,x ,x ). Corresponding augmentation of (6.8), (6.2)

and (6.1) yields


(6.9)


Yk+ m G yk + dk + fk' y0


YlkAkY2k

0

0


Y3kUk-Uk
uk

uk


YlkAkY2k -Ak [Y3kCkY4k-CkJH

Ak CkH

0 F


Sfk


[Y3kCkY4k-Ck Vk

CkVk

wk


where


G k =


dk -


Axo
^0
xo

x0


I YO M









Mean Error Analysis

From the stochastic difference (6.9) for Yk' we can

write a recursive equation for the mean lk=E[yk] as follows:


(6.10)


k+1 -= Gklk + dk,


where


G = EY1 kAk EY2k EY1kAkEY2k-Ak [EY3kCkEY4k-Ck ]H

0 Ak CkH

0 0 F


10 = 0
^0
x0

x0


Furthermore


EYIk = E[Ylk]

EY2k E[Y2k]

EY3k = E[Y3k]

EY4k = E[Y4k]


= diag (p,p, .....,p)
,. ,n-i n-I n-2
- diag (n- ,n ,- ,n ....
2 2 2
= diag ( 2 ,2 .... ,1 )

= a m m m-l 2)
= diag (p ,/ ,p .... ,/ )


where
E rq/2 -q/2
S= E[r ] =r
q(lnr)


dk -


q=2








Variance Error Analysis
Computation of the variance is more complex since the
matrix Gk contains random variables. Decomposition of the
KT
covariance matrix, Pk=E[ (yk-lk)(Yk-lk) T and state matrix Gk
yields

P k EtYkYkT] 1 T (6.11)
k = E[ iklk

Gk Glk + G2k (6.12)

where

Gl k= 0 -Ak -CkH
0 Ak CkH
0 0 F


G2k = YlkAkY2k YlkAkY2k Y3kCkY4kH
0 0 0
0 0 0

From the stochastic difference (6.9) for Yk' and using
(6.12) we can write a recursive equation for the mean square
T
error Mk=E[kyk T] as follows:

Mkl=GkMkGlkT +'T G ~d'T G'MGT
Mk+1 = G1 lT + Gl kMG2T + Glklk k + G2kMkGlkT +
EG2kYkYkTG2 T [G2kYkdkT] d'llT
+ E[G2yyTG2TI + E[G2 ydT] + d+lkGlkT +

+ E[dkykTG2kT] + E(dkdkT] + E[fkfkT M









where


G2k = EY1kAkEY2k EYlkAkEY2k EY3kCkEY4kH

0 0 0

0 0 0


M = 0 0 0
0 0 0

0 0 P0(-)


We can now predict the mean and covariance of the

actual estimation error due to a finite precision LNS

implementation.



Simulation Studies

The purpose of this section is to investigate how well

analytical predictions match actual experimental results

obtained by simulation, and to compare the finite precision

conventional Kalman filter and the logarithmic Kalman filter

in terms of speed and accuracy.

For the purpose of comparisons during the simulation

studies, the double-precision operations of a VAX 750 with

64 bits are taken to be infinite-precision operations.

The performance of the Kalman filter (6.2) for a fourth

order elliptic lowpass filter was investigated for both LNS

and conventional fixed point implementations. In order to











avoid overflow operations four bits were used to form the

integer field of a N-bit sign-magnitude, fixed point

wordlength format, leaving F =N-4-1 fractional bits. In the

N-bit LNS data format two bits were used to cover the same

integer dynamic range resulting in F I-N-2-2 fractional bits.

A Monte Carlo simulation was carried out with 10,000 sample

paths for each register-length implementation for both cases

(i.e. LNS and conventional). The average degradation in

the LNS mean E[ek] and the LNS covariance cov[ek,ek] of the

actual estimate can be compared with the theoretical

estimate analytically predicted by (6.9) and (6.11). The

results of the simulation and the analytical prediction for

the degradation of LNS mean are shown in Figure (6-1), for

8, 12, and 16-bits while the corresponding LNS covariance

comparisons are shown in Figure (6-2). The simulation

results are indicated with marks while straight lines show

the analytical predictions. The corresponding results for

the conventional fixed point mean and covariance, using

Stripad's error model [91], are illustrated in Figures (6-3)

and (6-4).

There appears to be good agreement between the

analytical predictions and the averaged simulation results.

In addition, a comparison between LNS and conventional fixed

point results demonstrates the superiority of the LNS

architecture in terms of accuracy, as illustrated in Figure

(6-5). The straight lines represent the LNS error variance

















STATE: 1


Figure 6-1. Plot of LNS Error Mean.


8 BITS





12 BITS




16 BITS





















STATE: (1,1)


c7


-. -t.


U
I)


%V


,' \~I
/
/
/
/


r T T -i
1 2 3 4


S 1 T. .T -
5 6 7 U 9
DISCRETF TIME


Figure 6-2. Plot of LNS Error Covariance.


1-r -f


-4- 4.


I 0-


10-s


F. 10-"



U3
1c)-


w i0n-
Z

~ o-
i0-6




C2 1 fl
>
K10-
0
r-


10- -

!0-" -=


8 BITS


12 BITS


S16 BITS


1o-"


-.=
S


z























STATE: 1


H
--1

10-'j /-'


10- -=


7
10-' -
z





M -




1O:a -
I0-'-=



10-" -



I0-'-
1


8 BITS
*


12 BITS


16 BITS


(7 Y"\


Figure 6-3. Plot of Con. Fixed Point Error Mean.


I I I I I I I I I
2 3 4 5 6 7 8 9 10 11
DISCRETE T[ME




















STATE: (1,1)


1 ,,,


Tr -w


\ T


: 6 ( U
DISC)1RETE T1MTI;


9 10 11


Figure 6-4. Plot of Con. Fixed Point Error Covariance.


L0-

1.0-.


10 '


8 BITS



12 BITS


z




0


10-'




10-7 I_


10- __=
10' -



1 0-l" -
10-"


1


..---. -... 16 BITS


1~-.~ I
2 ~1


1, -







70












10-1^ 8 A J BITS

10 A A ,
0 ... 81 BITS



10
1Di 1- A A- A 2BT


S10- - A -


0- --= / A 16 BITS
-A _____----
S= / A ------





10-4 /"
-/
i0-1" --/ i
=o- --------------------------
101- I




1 2 3 4 5 6 7 8 9 1
DISCRETE TIME





Figure 6-5. Comparison between Con. Fixed Point and LNS
Error Covariances.










and marks denote the conventional fixed point error variance

for 8, 12 and 16 bits.

Throughput rates of the Kalman filter (6.2) can be

estimated in terms of the system's multiplication time

(TMULT) and add time (TADD). Specifically, the execution

time for a n-dimensional dynamic system with a measurement

vector of order m is


T n[(n+m)TMULT + (n+m)TADD] (6.13)


Using the Schottkey TTL (2500LS) design as a reference, and

2500LS parts with fixed point addition (TA) 37ns plus a 4K

35ns ROM parameters (TM), then TADD-2TA+TM-109nsand

TMULT-TA-37ns for LNS arithmetic. Using the same technology

for conventional fixed point arithmetic addition and

multiplication can be performed in 37 nsec and 150 nsec

respectively. It can be noted that the multiplication

operations in an LNS architecture are comparable to

conventional fixed point addition. However, there is an

advantage of LNS addition over fixed point multiplication

and considering the fact that there are as many

multiplications as additions in (6.13), the LNS superiority

becomes apparent. If pipeline configuration is used, where

throughput is a function of the slowest process element, the

throughput rate for LNS is 37ns and the throughput rate for

fixed point is 150ns.

















CHAPTER SEVEN
ADAPTIVE KALMAN FILTERS








The Effect of Erroneous Models
on the Kalman Filter Response


A well-known limitation in the application of the

Kalman-Bucy filter to real-world problems is the assumption

of known a priori statistics for the stochastic errors in

both the state and observation process. Wrong statistics

made a priori can lead to a "wrong" Kalman filter. A

derivation of the "wrong" Kalman filter is possible,

assuming wrong a priori statistics.

Given an a priori estimate of the "wrong" state

covariance P0 and given the "wrong" a priori statistical

information of the input (source) noise covariance QI and

output (measurement) noise covariance RI, a suboptimal state

of the system defined by (6.1) can be obtained sequentially

for k=1,2,3... similar to the standard Kalman filter having










a state estimation


Xk+1 F(I-K kH)Xk + uk + F k (7.1)


an error covariance


Pk+1 = F(I-KkH)P FT + QI (7.2)

and a Kalman Gain Matrix

T T+, -I
K PH (HPH +RI] (7.3)
k k k

Observe that for an optimal Kalman filter studied under the

case that QI=Q and RI=R, then P k=P k is the covariance of the

error in estimating the state. Now assume that the actual

gains employed in the filter design deviate from the optimal

gains, i.e., that


Kk Kk + Kk

In order to determine the effect of these errors we ask the

following question. How close do we come to the optimum by

using the suboptimal filter equation (7.1). The answer is

given by the "actual" covariance matrix of the estimation

error obtained when (7.1) is used instead of (6.2). The
"actual" covariance matrix denoted by PAk indicates the

quality of the suboptimal estimate and is defined as

Ak+l k+l-* Tlk+l- l
P~k1 =E[(k+l-Xk+l)(k+l-Xk+l)










Considering the error in the suboptimal case and using
both the system model (6.1) and suboptimal state estimation

(7.1), we obtain the recursive equation
(x ^, **
(Xk+l-Xk+) (F(I-KkH)(xk-xk) + wk +KkVk

and after some manipulation we have

TFT* *TFT
PAk+ F(I-KkH)PA(I-KkH)TFT + Q + FKkRK F
k+1 r~-kHPk .~kH k k

It should be noted that the "actual" steady state error
covariance PA depends upon both the a priori statistics of

the system model (i.e., "wrong" Kalman gain, K*) and the

true noise statistics R and Q.
In order to quantify the difference between the
"actual," optimal and "wrong" steady state error covariances

consider the single-input, single-output model (i.e., m-1

and wk a scalar time series) where Q is a diagonal matrix
with only one non-zero entry element and R is a scalar

matrix. Figures (6-1), (6-2), and (6-3) illustrate these

differences. Here, each iteration represents the Kalman

filtering of a 256 point time series with R-RI and QI

changing. The value of the true noise statistics R and Q
are shown in the upper right corner of the figures. The
trace of the steady state error covariances are plotted for
each iteration with the ratio (RI/QI) distributed
logarithmically from .01 to 100. The straight line

represents the optimal steady state error covariance. The









triangles represent the "wrong" steady state error

covariance while the stars represent the "actual" steady

state error covariance based on the "wrong" Kalman gains.

We can observe that the crossover point of these curves is

the optimal case in which we obtain the optimal Kalman gain.

In addition, the "actual" error covariance is increasing

after the crossover point. It can be explained by the
"wrong" Kalman filter gain which is decreasing in each

iteration, and as a result the innovation error is

magnified.





The Effect of Input and Output Noise Covariances
on the Kalman Gains


The performance of the Kalman filter depends on the

choice of the a priori noise statistics RI and QI which in

turn define the optimal steady state Kalman gain. A proof

will be presented showing that in order to obtain the

optimal Kalman gain the ratio of the a priori noise

statistics (RI/QI) must be equal to the ratio of the true

noise statistics (R/Q).

The steady state optimal Kalman gain and the steady

state optimal error covariance are given by


K = PHT [HPH T+R] -


(7.4)

























R=~AJ1


"C A


18


16


14
z

CilO
5 12


10

<6






4


2


0


I I 1111111

10-I


I I I 1111 + I i/QI
LOC 10'
RATIO (RI/Q[)


-r mTII1


Figure 7-1. Error Covariances for R-0.01 and Q-0.01


* A

(1 A
A
A
A^ -


10-2


* V A: fc -* 0 1. 1 . * *- .


I -


























R=3.1
iQ=0.01


180 -


160 -


140 -


A

A
4' A
* A
4' A


4' A


4' A
4' A
A


4' A
*4k4 ., -, ,,


I I I I I I I A- A, & it ....... ...


I III lii"'


10-'


I I I I I 1111


I FTTF1TiJ


10
RATIO (RI/Q[)


Figure 7-2. Error Covariances for R-0.1 and Q-0.01


100 -


80 -


60 -


40 -


20


0-


' ' '" 'I














































SA


*
A
*
*


- I I I I I. I I


I I I I I I I II


Figure 7-3. Error Covariances for R-0.01 and Q-0.1


18 -


16


14-


12 -


10 -


R=@J~1
Q=~. 1


V.


1 0-'


100
RATIO (RI/QD[


I T T I llI


6h >tl~


a i naL- aP ,"aa ,









P F[P-PHT(HPHT+R) HPIFT + Q


(7.5)


while the suboptimal steady state Kalman gain and suboptimal

steady state error covariances are shown to be


* HT [HP*H T ]-1
K P H (HP H +RI]


* *F[P T (HP*HT +RI)-1 H* T+ QI
P F[P -P H (HP H +RI) HP IF + Q1


(7.6)



(7.7)


Let the true noise statistics be equal to the ratio of the
"wrong" noise statistics, or


R a RI


Q a QI


(7.8)


(7.9)


Direct substitution of (7.8) and (7.9) in (7.5) yields


P = F(P-PHT (HPH T+aRI)- HP]FT + aQI


or


QI = -F(P/a)FT + FPHT (HPH T+aRI)- H(P/a)FT + (P/a)


(7.10)


Solving (7.7) for QI and equating the result to (7.10) we

have










0 -F[(P/a)-P*]FT + [(P/a)-P*l +


+ F[PH (HPHT +aRI) 1H(P/a)-P*HT(HP*HT +RI) 1HP* ]FT


The solution of the above equation is


P (P/a) or P a P (7.11)


Direct substitution of (7.11) and (7.8) in (7.4) results


K aP*H T[a(HP*H T+RI)]-I

T T -l
= P HT (HP H T+RI]- 1

*
= K

Therefore the optimal steady state Kalman gain can be

defined in terms of the ratio of the a priori noise

covariances (RI/QI). In other words knowledge of the

correct ratio of the noise covariances yields the optimal

steady state Kalman gain. This is illustrated in Figures

(7-4), (7-5) and (7-6). The straight line represents the

optimal steady state Kalman gain. The triangles represent

the "wrong" steady state Kalman gain while the stars

represent the "actual" steady state Kalman gain to be

defined in the next section. We can observe that the

crossover point of these curves is the optimal case and at

that point the ratio of the a priori noise statistics is

equal to the ratio of the true noise statistics.




























R=0.01
P =F-. 0 1


18 -


16 J


14 *


12-


10-


8-


A


A


* A


- ~ * *. 4 ~ 9~ 4 4' ~ ,~, .'


fA
d~a


I I I 1 i I j1 I I I t1 t i


I (J~'


I ()c
RAT[O RI/Q[)
RAT[O (RI/0[)


Figure 7-4. Kalman Gains for R=0.01 and Q-0.01


* A


10-2


I I i I I1i




























R=O.1
0Q=0.01


1 tg


16 A


14 1


A
A


* A
* A


I' *~ *1


.. I I I I I IIi I -I I TI *Ir Il


I I I 111111
10-I


I I I I 1 1 1


I Rc [
RATIO (RI/Q1i


Figure 7-5. Kalman Gains for R-0.1 and Q-0.01


12 -


10 -


8-


10-2


I i I I I I I


i I I I I I I I
























R=~J~11
D=~. 1


* A
l~ A
*


,g
A
A


.1 AAAA


I I 11111
iT'


I I II IIII1 I I I I1111
10C 10'
RATIO (RI/Q[)


I I I i I I I I


Figure 7-6. Kalman Gains for R-0.01 and Q-0.1


10-2


4* .











Adaptive Kalman Filtering


We have seen that to yield optimal performance for a

Kalman filter it is necessary to provide the correct a

priori descriptions of F, H, Q, R and P0. As a practical

fact, this is usually impossible; guesses of these

quantities must be advanced. Hopefully, the filter design

will be such that the penalty for misguesses is small. But

we may raise an interesting question; i.e., "Is it possible

to deduce non-optimal behavior during operation and thus

improve the quality af a priori information?" Within certain

limits, the answer is yes. Using the innovations property

of the Kalman filter an adaptive algorithm can be designed.

The innovations error process of a Kalman filter is

defined as


k M Zk Hxk


The innovations property states that, if K is the optimal

gain,


ElvkV 0; j 0


In other words the innovations process, v, is a white noise

process. Heuristically, there is no "information" left in
A
V k' if xk is an optimal estimate.

Let the autocorrelation function, for a lag variable k,

be denoted by A(k). Then the first two autocorrelation

measures of the innovation error process are









A(0) = E([k- ] T H PA H + R (7.12)


A(1) = E[kvkl = HF(PAHT K (HPAHT +R)]

HF[KA K* ][HPAHT+R] (7.13)



where KA is the steady state Kalman gain corresponding to

the "actual" error covariance. That is

Ai T T -1
KA = PAH [HPAH +R]- (7.14)


Note that when optimally configured, the optimal Kalman gain
K is given by KA in Equation (7.14) and A(1), in Equation

(7.13), diminishes to zero. This well known condition of

the innovations error being "white" will be used as an

adaptive filter design objective. Furthermore, this

observation suggests that if the ratio of the a priori

information, namely (RI/QI), is equal to the ratio of the

true noise statistics, namely (R/Q), then A(1) is always

zero. As the ratio varies up or down, A(1) moves positively

or negatively in harmony. Specifically if (R/Q)>(RI/QI),

the computed Kalman gain is greater than its "actual"

counterpart, as it can been seen in Figures (7-4), (7-5) and

(7-6), and forces A(1) to be negative. In the other case

where (RI/QI)>(R/Q), the results are opposite. This is

illustrated in Figures (7-7), (7-8) and (7-9). Here the











second term of the autocorrelated innovations error process,

namely A(l), is plotted as a function of the ratio RI/QI.

In this figure, the theoretical results are indicated with a

straight line while the stars denote the experimental

results.



Algorithm Implementation


To implement the proposed algorithm, several key

functional units need to implemented and they are

1. A shift-invariant Kalman filter.

2. An autocorrelator.

3. Adaptive search and decision rule.


They will be briefly developed as follows:



Kalman Filter

It is desired to realize a given four-tuple (F,H,QI,RI)

description of the nth order SISO Kalman filter as a high

throughput filter. Fast systolic and SIMO architectures

have been proposed to achieve this goal [92]. Pipelining

can also be used to increase the effective throughput of an

implemented Kalman filter. These methods will be presented

in Chapter 8. In addition using LNS Kalman filters higher

throughput and accuracy can be obtained.





































4.
4.


4..
4.

~


-0.5 -.



-1.0 -i


10'


1- RAT[O (RI0


I tl


I I I I


Figure 7-7.


Autocorrelation of Innovation Process
for R-0.01 and Q-0.01


r . . . L


R=0.01
Q=0.01





































*
I'
4' 4~,'
$'


I I II 11111
10-I


I 1 111 1 I I 1 1
10 l0'
RATIO (RI/Q1)


Figure 7-8.


Autocorrelation of Innovation Process
for R-0.1 and Q-=O.01


88















R=p. 1
Q= .D


-4 -



-1
t f i


-16 -


-20 -


-24 -


-28 -


V


-l .


10-3


I I








89












,* R= 01
|Q=0.1I


-4 1


4
V


4 -


- 'I


-1



2
.3
I /*'
-2 / ..



3 -i - i I I| - I I ll | - I i I ll l -- I 11 ill
to-" ia' too 10
I0-" 10-' 1T O l01
________________RATIO (RIiQ[),______


Figure 7-9.


Autocorrelation of Innovation Process
for R-0.01 and Q-0.1









Autocorrelation

A discrete autocorrelation function is given by


N
r(k) = Z x(n) x(n+k)
n-i


If r(k) is desired for a range of k on the order N/2, then

using the FFT to compute DFT- [x(f)x* (f)]-r(k) is

computationally efficient. For relatively short delay

analysis (i.e., k<
algorithm is optimum from a computationally complexity

viewpoint. The BP version of the autocorrelation algorithm

satisfies


p-i k
r(k) = Z x(2jk+i+k) [x(2jk+i) + x(2jk+i+2k)]
j=-0 i=1

k 1 2, ..., p; p-INT(N/2k)


For N>>p, the BP algorithm essentially halves the number of

multiplications normally associated with direct computation

and can also be considerably less than the associated with

FFT mechanizations.

The speed of the BP algorithm can be further improved

by eliminating the time consumed to compute data invariant

indices (i.e., (2jk+i+k), (2jk+i), (2jk+i+2k)). This can be

achieved through the use of so called inline code or knotted

code [93]. Therefore, in the context of the proposed

adaptive filter, it is apparent that the computation of the









autocorrelation function can be performed at real-time

speeds (i.e., t(autocorrelation) < t(Kalman filter)).
Another alternative to the BP algorithm is the

dedicated hardware unit, illustrated in Figure (7-10). It

consists of two multiply-add modules in parallel

configuration. Each module calculates each of the

autocorrelation coefficients, namely A(0) and A(l), at real

time speed.


Search Algorithm

In order to efficiently estimate the ratio R/Q, a

modified Fibonacci search has been developed, as illustrated

in Figure (7-11). Letting f(x) denote a real-valued

function of a single real variable x, assume that f(x) has a

invariant minimum at x in some real interval [a0,b0]. The

purpose of the search algorithm is to sequentially reduce

the search initial interval to [ann1<[an-1 bn-1] which also

contains x after the nth iteration. In the proposed

algorithm, the points x1 through x5, are chosen to divide

[a0,b0] into four subintervals of equal length and then

compute f(a0) and f(b0). If If(a0)|>If(b0)I, in the next

iteration a1=x2 and b1=b0. Otherwise, a1=a0 and b1-x4. As

a result, the new interval has been reduced by one

subinterval so that after two iterations the interval

[a0,b0] has been halved. This process is repeated until the

prespecified interval of uncertainty is reached.
























A(1)









A(O)


Figure 7-10. Hardware Autocorrelation Unit.