NTIA Report 02-395 

Estimation of System Gain and Bias 

Using Noisy Observations with 

Known Noise Power Ratio 

Stephen D. Voran 

U.S. DEPARTMENT OF COMMERCE 

Donald L. Evans, Secretary

Nancy J. Victory, Assistant Secretary 

for Communications and Information 

September 2002 

iii

ACKNOWLEDGEMENT 

The author gratefully acknowledges the contributions of Stephen Wolf to this work. He 

applied several of the algorithms presented in this report to objective and subjective video 

quality data and provided valuable feedback on the algorithms. 
iv
v

CONTENTS

FIGURES..................................................................................................................... vi

TABLES ...................................................................................................................... vi

ACRONYMS.............................................................................................................. vii

ABSTRACT...................................................................................................................1

1. INTRODUCTION .....................................................................................................1

2. LEAST SQUARES PLUS MANUAL SPLITTING (LSMS) ...................................6

3. TOTAL LEAST SQUARES (TLS)...........................................................................8

4. WEIGHTED JOINT MINIMIZATION (WJM)......................................................12 

5. CONSTRAINED JOINT MINIMIZATION (CJM)................................................14

6. DIRECT ESTIMATION (DE).................................................................................19

7. EXAMPLE SIMULATION RESULTS AND DISCUSSION ................................22

8. SUMMARY.............................................................................................................30

9. REFERENCES ........................................................................................................31

APPENDIX: SUMMARY OF ESTIMATION ALGORITHMS ...............................33 
vi

FIGURES 

Figure 1. Block diagram for noisy observations of a system with gain and bias ........3 

Figure 2. Alternative interpretation of Figure 1...........................................................4 

Figure 3. Example geometry of LSMS result ..............................................................8 

Figure 4. Example geometry of TLS result ...............................................................11 

Figure 5. Example geometry of WJM and CJM results.............................................14 

Figure 6. Block diagram for bias-free version of the estimation problem with mx

 and my removed..........................................................................................20 

Figure 7. Example results for five estimation algorithms; true gain is 0 dB .............23 

Figure 8. Example estimation errors for the case SNR=20 dB, r =-10 dB.................23 

Figure 9. Example estimation errors for the case SNR=20 dB, r =0 dB....................24 

Figure 10. Example estimation errors for the case SNR=20 dB, r =10 dB..................24 

Figure 11. Example estimation errors for the case SNR=10 dB, r =-10 dB.................25 

Figure 12. Example estimation errors for the case SNR=10 dB, r =0 dB....................25 

Figure 13. Example estimation errors for the case SNR=10 dB, r =10 dB..................26 

Figure 14. Example estimation errors for the case SNR=0 dB, r =-10 dB...................26 

Figure 15. Example estimation errors for the case SNR=0 dB, r =0 dB......................27 

Figure 16. Example estimation errors for the case SNR=0 dB, r =10 dB....................27 

Figure 17. Example estimation errors for correlated noise environment.....................29 

Figure 18. Example estimation errors for the case of inaccurate r values...................29 

TABLES 

Table A-1. Mathematical Constructions Leading to the Five Algorithms..................34 

Table A-2. Preparation Steps Common to the Algorithms.........................................35 

Table A-3. Basic Algorithm Results...........................................................................36 

Table A-4. Additional Algorithm Results...................................................................37 

Table A-5. Additional Algorithm Results...................................................................38 
vii

 ACRONYMS 

CJM constrained joint minimization 

DE direct estimation 

LSMS least squares manual splitting 

SNR signal-to-noise ratio 

SVD singular value decomposition 

TLS total least squares 

WJM weighted joint minimization 

XLS extended least squares 
viii

 
Estimation of System Gain and Bias Using Noisy Observations with 

Known Noise Power Ratio

Stephen Voran* 

The identification of linear systems from input and output observations is 

an important and well-studied topic. When both the input and output 

observations are noisy, the resulting problem is sometimes called the 

“errors in variables” problem. Existing work on this problem deals with 

the identification of multivariate systems and thus results in algorithms 

that are necessarily somewhat complex and often involve iteration. In this 

report we treat an important special case of the problem: estimation of a 

system bias and a system gain from noisy observations of system input 

and output. In addition, we invoke an input-output noise power ratio 

constraint. This constraint can also be interpreted as a parameter that 

moves the problem in a continuous fashion between two limiting cases, 

each of which is a conventional least-squares problem. We do not model 

the input signal, and we place minimal restrictions on the input and output 

observation noises. We develop five different low-complexity closedform solutions to the problem. The final two are the most satisfying and 

we explore these further through simulations. Our original motivation for 

working on this problem came from the need to calibrate objective and 

subjective estimates of perceived video or speech quality. We expect that 

our solutions may also find applications in remote sensing, active noise 

reduction, echo cancellation, channel estimation, and channel equalization. 

Key words: gain estimation, bias estimation, linear system identification, input noise, 

 errors in variables, total least squares, extended least squares, audio 

 quality estimation, speech quality estimation, video quality estimation 

1. INTRODUCTION 

ITS conducts research on the objective estimation of perceived video and speech quality. 

As part of this work it is often necessary to scale and shift subjective viewing or listening 

test results in order to compensate for the use of different laboratories, different 

populations of subjects, or different languages. The appropriate scaling (gain) and 

shifting (bias) factors are jointly estimated from subjective test results and objective 

estimates of perceived quality. Both of these can be considered to be noisy observations 

of a true underlying perceived quality value. This estimation problem is closely related 

to the problem of identifying linear systems from noisy input and output observations. 

________________________ 

* The author is with the Institute for Telecommunication Sciences, National 

Telecommunications and Information Administration, U.S. Department of Commerce, 

325 Broadway, Boulder, Colorado 80305.
2

The identification of linear systems from input and output observations is an important 

and well-studied topic. Much of this work has been based on noisy observations of 

system output and noise-free observations of system input. This is also described as the 

“errors in measurements” situation. While this approach is clearly useful in numerous 

situations, there are also plenty of situations where system identification must proceed 

from noisy output observations and noisy input observations. This could be described as 

the “errors in variables and measurements” situation, but it typically is identified by the 

shorter “errors in variables” moniker. A history of efforts at system identification using 

noisy output and input observations is provided in [1] which in turn points to [2] for 

coverage of earlier work. 

The classical solution to this problem is given by total least-squares (TLS) [3],[4] and 

numerous other innovative approaches have been offered as well. A minimum 

description length approach is based on the argument that the best approximation to the 

underlying linear system is the one that allows the observations to be encoded with the 

fewest bits [5]. This requires a joint probabilistic model for the observations and the 

noises. Examples of cumulant-based solutions can be found in [6] and [1]. In general, 

cumulant-based solutions require multiple steps, or may even be iterative. A frequencydomain approach utilizing bispectra or trispectra (subject to certain restrictions on the 

observation noises) is given in [7]. 

The solution in [8] stems from minimization of a residual norm (or output observation 

noise) subject to an absolute constraint on a coefficient perturbation matrix (or input 

observation noise). This solution requires a singular value decomposition (SVD) of the 

data matrix. In [9] this work is extended to include a family of absolute constraints 

where each applies a different portion of the coefficient matrix perturbation. An 

extended least squares (XLS) formulation of the problem is described in [10]. This 

formulation seeks to minimize a weighted combination of the model mismatch (or input 

observation noise) and measurement inaccuracies (or output observation noise) and two 

iterative algorithms result from this formulation. In certain specific cases, the XLS 

approach is equivalent to optimization of a joint maximum a posteriori – maximum 

likelihood criterion [11]. The problem formulation of XLS is connected to our weighted 

joint minimization problem formulation as described later in this report. 

In general, the existing work allows for identification of multivariate systems and thus 

results in algorithms that are necessarily somewhat complex, and often involve iteration. 

In this report, we treat a relatively simple but important special case of the problem that 

yields low complexity closed-form solutions. This version of the problem is also unique 

in that a single piece of side information is used: the input-output noise power ratio r

2

. 

We simply wish to estimate a scalar gain (or equivalently a scalar attenuation) and bias 

using noisy observations of system input and output. We make no assumptions about the 

signal, we do not build a model for the signal, and we do not perform any spectral 

analysis. 
3

The situation is described in Figure 1 and by 

( ) , y T xT y x += + + ξ ξ a b 1 (1) 

where [ ] [ ] T T

123 12 3 , , , ..., and , , , ..., n n x y = = xxx x yy y y are the noisy 

observation vectors (corresponding to the input wT and the output Tz respectively), aT

and bT are the true scalar system gain and bias to be estimated, ξx and ξy are length n

column vectors of observation noise, and 1 is the length n column vector of ones. The 

additive noises in Figure 1 are defined to be -ξx and -ξy rather than the more conventional 

ξx and ξy. The only practical consequence of this definition choice is more readable 

mathematical derivations throughout this report. 

aT +

bT

+ +

wT

-ξx

wT = x + ξx

x = wT - ξx 

(noisy observation of wT)

y = zT – ξy = aT (x + ξx) + bT1 – ξy (noisy observation of zT) 

zT = aTwT + bT1 = aT (x + ξx) + bT1

Actual: y + ξy = aT (x + ξx) + bT1 

Estimated: y + εy = a (x + εx) + b1

Figure 1. Block diagram for noisy observations of a system with gain and bias. 

The absolute noise power levels are unknown, but they do conform to a cost-weighted 

noise power ratio constraint: 

2

2

2

E

0 ,

E

x

y

r

ξ

ξ

< =

C

C

 (2) 

where C is an n x n diagonal cost matrix populated with { } 1

n

i i

c = and 

2

1

0 < , i =1 to and 1.

n

i i

i

c nc =

∑ = (3) 

C can be used to appropriately weight the noisy observations if such weighting is 

indicated by the details of a particular problem under consideration. Otherwise C can be 

set to 1

n n×n I where n×n I is the n x n identity matrix. Given a pair of noisy observation 

 – ξy
4

vectors x and y, a cost matrix C, and the noise power ratio r

2

, we seek to find a≠0, b, εx, 

and εy such that 

y + ε y x= ++ a b ( x ε ) 1 (4) 

and 

2

2

2 = r . x

y

Cε

Cε

 (5) 

We will use a, b, w, εx, and εy as estimates of aT, bT, wT, ξx, and ξy respectively. While a

and b are the estimates sought, in some cases the estimation algorithms also generate w, 

εx, and εy as well. 

In Figure 1, the system input is wT and we have access only to noisy observations x and y. 

If we take zT to be the system output, then we have modeled a noise-free system that only 

scales and shifts the input wT. If we take y to be the system output, then we have modeled 

a noisy system that scales, shifts, and adds noise to the input wT. Figure 2 shows a 

mathematically equivalent interpretation of Figure 1. Here we consider x to be an 

observed noise-free input that is contaminated by the noise ξx before it enters the system 

that scales and shifts it. 

aT +

bT

+

x +

ξx

wT = x + ξx

-ξy

y = zT – ξy = aT (x + ξx) + bT1 – ξy

zT = aTwT + bT1 = aT (x + ξx) + bT1

Figure 2. Alternative interpretation of Figure 1. 

Note that the input-output noise power ratio r

2

 can be viewed as the parameter that lets us 

move in a continuous fashion between two conventional least squares problems. As r

goes to zero, the input noise vanishes (we have the errors in measurements case), and the 

problem that remains looks like a conventional least-squares problem. As r gets large, in 

relative terms the output noise vanishes (we have the errors in variables case) and the 

problem that remains looks like a closely related least-squares problem. In the case of the 

constrained joint minimization algorithm described later in this report (see Section 5), we 

find that the solutions are similarly parameterized by r, and they reduce to least-squares 

solutions in the limiting cases r = 0 and r → ∞. 
5

Note also that at this point we have not made any assumptions about the observation 

noises ξx and ξy. They could even be correlated to the signal (as in echoes or distortion 

products) but we will see later that the best estimates result from the case of uncorrelated 

observation noises. 

The problem we have described is a realistic one when environmental conditions or the 

physics of a situation yield prior knowledge of the ratio of the two observation noises, but 

not the absolute level of the observation noises. This might apply to multiple 

observations of signals by imperfect instruments or transducers (to perform 

interferometry or to form synthetic apertures), or multiple replications of an experiment 

at different times or places (as in meta-analysis [12]). Thus we expect that solutions to 

this problem may find applications in remote sensing, active noise reduction, echo 

cancellation, channel estimation, and channel equalization. 

In these applications, system physics often dictate that the bias is zero. Our original 

motivation for working on this problem comes from the objective estimation of perceived 

video or speech quality (see for example [13]-[16]) and in these applications the bias term 

is generally not zero. We maintain that subjective viewing or listening tests provide 

noisy observations (x and y) of an underlying true value of perceived quality (w). 

Similarly, objective quality estimates also give noisy observations of that true value. 

Subjective tests are not absolute; results may require scaling (gain) and shifting (bias) to 

compensate for the use of different laboratories, different populations of subjects, or 

different languages. But objective estimates are absolute. To advance the state of the art 

in objective estimation of perceived video and speech quality, it is necessary to relate 

objective estimator results to a number of disparate subjective test results, using a 

different gain and bias factor for each subjective test. These gain and bias factors must 

be estimated from noisy subjective (x) and objective (y) observations of the underlying 

true perceived video or speech quality. The absolute noise levels are not known. 

However, given our experience with objective estimators, we do have some confidence 

that the noise in the subjective test results is smaller than the noise in the objective 

estimates. In other words, r ≤1. This report provides tools for estimating the gain and 

bias factors that relate the results of a single subjective test with the corresponding 

objective estimates of perceived quality. 

In Sections 2-6, we develop five algorithms for estimating gain and bias in the situation 

described in Figure 1. We have adopted the following names for these algorithms: “least 

squares plus manual splitting” (LSMS), “total least squares” (TLS), “weighted joint 

minimization” (WJM), “constrained joint minimization” (CJM), and “direct estimation” 

(DE). The first four algorithms require no assumptions about the noises beyond their 

known noise power ratio given in (2). These algorithms are derived by minimizing four 

different noise quantities. The DE algorithm places additional restrictions on the noises, 

and does not involve the minimization of a noise. Rather, it follows directly from the 

algebraic manipulations of the noisy observations. 
6

We have generated simulation results for numerous situations including tone, chirp, 

noise, and speech signals and using both correlated and uncorrelated observation noises. 

These results and some discussion are presented in Section 7. The appendix summarizes 

the five algorithms through the use of several tables. 

2. LEAST SQUARES PLUS MANUAL SPLITTING (LSMS) 

The derivation of the five algorithms is facilitated by the following definitions: 

2 2

1 1

, ,

n n

x ii y ii

i i

m cx m cy = == = ∑ ∑ (6) 

x Cx y Cy ˆ ˆ =− =− ( ) m m x y1 1 , , ( ) (7) 

 T ˆ ˆ 1 when 0 , , sign( ) ˆ ˆ 1 otherwise,

ρ

ρ ρ

 < = = 

 −

x y

x y

 (8) 

and 

[ ]

T

123 , , , ..., . w = www wn (9) 

Note that mx and my are cost-weighted means, xˆ and yˆ are cost-weighted, shifted versions 

of x and y such that Cxˆ andCyˆ are zero mean. Note also that ρ is the normalized crosscorrelation of xˆ and yˆ . We assume that 0 < xˆ and 0 < yˆ . If multiple pairs of 

observations are available, then quantities like ˆ ˆ ˆˆ , , , ,and T

m m x y x y xy are estimated by 

averaging over all available observations. 

The first algorithm, LSMS, comes from a conventional least-squares step followed by the 

manual splitting of the least-squares residual into two vectors. We can rewrite (4) as 

( ) . y x LSMS ab a x y + =+ − = + 1 ε ε y ε (10) 

Conventional weighted least squares allows us to minimize 2 2

LSMS LSMS ε = Cε with 

respect to a and b. The results are: 

ˆ , , ˆ LSMS LSMS y LSMS x a b ma m = =− ρ

y

x

 (11) 

( ) 2 2 2 ˆ ˆˆ ˆˆ ˆ , 1. ˆ LSMS LSMS LSMS = −= − = − a ρ ρ

y Cε xy xy Cε y

x

 (12) 
7

Next we manually decompose εLSMS into a pair of perfectly correlated noise estimates εx 

and εy that attain the desired noise power ratio r

2

 given in (5). Let 

1

LSMS

LSMS

a r

g

a r = +

 (13) 

and note that 0 < g < 1. We use g to decompose εLSMS into two portions: 

ε εε LSMS LSMS LSMS =− + (1 . g g ) (14) 

From (10) we see that 

. LSMS y LSMS x εε ε = − a (15) 

Equating the right hand sides of (14) and (15) leads us to define 

( ) 1 and . y LSMS x LSMS

LSMS

g g

a

εε ε ε = − =− (16) 

It follows that 

( ) ( ) ( ) ( )

2

2

2

sign 1 ˆˆ ˆˆ , , and . 1 1

x

x LSMS y LSMS

LSMS LSMS y

r

a ar

r a r a

− ρ = −= − = + +

Cε

Cε xy Cε x y

Cε

 (17) 

The total cost-weighted noise power is 

( ) ( )

2

2 2 2 2

2

1 ˆ 1 .

1 x y

LSMS

r

a r

ρ

+

+= −

+

Cε Cε y (18) 

When the noise power ratio is unity (r =1), and both cost-weighted shifted observations 

have unit power ( 2 2

x y ˆ ˆ = =1), (18) reduces to 

2 2 1

2 .

1 x y

ρ

ρ

− + =

+

Cε Cε (19) 
8

We can show that ˆ [ ] 0, T x CεLSMS = so the weighed error vector is orthogonal to the 

modified noisy observation vector xˆ . The two weighted noise estimates Cεx and Cεy are 

perfectly correlated since they are formed by scaling a single noise vector. Figure 3 

shows the geometry of the LSMS result in a simple two-dimensional case where mx = my 

= 0, r ≈ 2, and a ≈ 0.75. The figure shows the various additive relationships among the 

pertinent quantities in the LSMS algorithm, as well as scaling and orthogonality. The 

LSMS solution for this simple case can be described geometrically as follows. First 

C LSMS ε and a are found by projecting yˆ onto xˆ , then aCw is found by dividing C LSMS ε

into x −aCε and C y ε according to r. 

Figure 3. Example geometry of LSMS result. 

The LSMS algorithm satisfies the noise power ratio constraint given in (5). However, the 

estimates of a and b generated by the LSMS algorithm depend only on xˆ and yˆ . While 

xˆ and yˆ do implicitly provide the LSMS algorithm with information about r, the LSMS 

algorithm does not make explicit use of the r value that is available to the algorithm. 

This is inconsistent with our desire to use that side-information advantageously. An 

additional disadvantage is that the quantity minimized by the LSMS algorithm (15) is 

somewhat unnatural. Finally, the LSMS algorithm generates perfectly correlated noise 

estimates, but in real situations input and output noise sources are unlikely to be perfectly 

correlated. 

3. TOTAL LEAST SQUARES (TLS) 

The shortcomings of the LSMS algorithm motivate us to minimize a more meaningful 

quantity — the total cost-weighted noise power: 

2 2 2 . TLS x y ε = + Cε Cε (20) 

The problem is now to minimize 2

TLS ε with respect to w, a, and b and this results in a TLS 

algorithm. 

yˆ

xˆ

aCw

axˆ x aCε

C y ε

C LSMS ε
9

A necessary condition for minimizing 2

TLS ε is 

2 0, 1to TLS

i

i n

w

ε

∂ = = ∂ (21) 

and this results in 

( ) 2

1 . 1

a b

a

= +−     +

w xy 1 (22) 

An additional necessary condition for minimizing 2

TLS ε is 

2 0 TLS b

ε

∂ = ∂ (23) 

and this results in 

2

1

.

n

y ii

i

b m a cw == − ∑ (24) 

Using (22) in (24) and solving for b yields 

 . TLS y TLS x b m am = − (25) 

We can now use results (22) and (25) to rewrite the estimated noise terms and the total 

cost-weighted noise power: 

2 ( ) ( ) , 1 x xy

a

am m

a

− = − −−     +

ε x y 1 1 (26) 

2 ( ) ( ) 1 , 1 y xy am m

a = − −−     +

ε x y 1 1 (27) 

2 2 2 2

2

1 ˆ ˆ . 1 TLS x y a

a

ε =+= −

+

Cε Cε x y (28) 

Now (28) allows us to minimize 2

TLS ε with respect to the single variable a. The necessary 

condition 

2 0 TLS

a

ε

∂ = ∂ (29) 

yields 
10

2

2 ˆˆ ˆˆ 4 ˆˆ ˆˆ , 0, 2

0, 0.

TLS a

ρ

ρ

ρ

ρ

      −± − +  =     ≠ 



 =

yx yx

xy xy (30) 

To show that (30) is also a sufficient condition we must show that 

2

2

2 0 .

TLS

TLS

a a

ε

∂

<

∂ (31) 

Condition (31) is satisfied if we pick the correct sign in (30). Thus we arrive at the 

solution 

2

2 ˆˆ ˆˆ 4 ˆˆ ˆˆ , 0, 2

0, 0.

TLS a

ρ

ρ

ρ

ρ

      −+ − +  =     ≠ 



 =

yx yx

xy xy (32) 

Thus 2

TLS ε is minimized by (24) and (32), resulting in a total cost-weighted noise of 

2 2

2

1 ˆ ˆ . 1 TLS TLS

TLS

a

a

ε = − +x y (33) 

In the simple case where both cost-weighted shifted observations have unit power 

( 2 2

x y ˆ ˆ = =1), (32) reduces to aTLS=sign(ρ) when ρ≠0 and (33) reduces to 

2 2

1 . Cεx y + =− Cε ρ (34) 

From (26) and (27) it is clear that 

2

2

2 . x

y

= a

Cε

Cε

 (35) 

Note that (20) can also be minimized via the SVD-based TLS algorithm as described in 

[3] and [4]. By invoking the SVD, this approach can deal with the more general 
11

multivariate case. When applied to the present problem, the SVD approach would 

require that we extract the primary left singular vector u from the n x 2 observation matrix 

[ xˆ ˆy]. We would then form the projection operator ( ) 1 T TT

u

−

P = = u u u u uu . This 

operator projects the noisy observations onto the best fitting rank one subspace as defined 

by the span of u. Thus the weighted noise estimates can be found by the orthogonal 

projection ( ) ˆ − =− C x u ε I Px and ( ) ˆ − =− C y u ε I P y . This approach ultimately yields the 

solution given above, but with significant extra complexity. 

We can show that [ ] [ ] [ ] 0, T T

x y = =     Cw Cε Cw Cε so both weighted noise estimates are 

orthogonal to the weighted input vector w. From (26) and (27) it is clear that the TLS 

algorithm generates weighted noise estimates Cεx and Cεy that are perfectly correlated. 

The geometry of the TLS result in a simple two-dimensional case with mx = m y= 0, a ≈

0.75, is shown in Figure 4. A geometric description of the TLS solution for this simple 

case follows. Given that x aCε and C y ε must be orthogonal to aCw as shown, we can 

swing aCw and slide a to find the unique solution that minimizes 2

TLS ε . 

Figure 4. Example geometry of TLS result. 

We have not included the noise power ratio constraint r

2

 in the derivation of the TLS 

algorithm, and consequently that noise power ratio will not be satisfied in general. 

Further, the TLS algorithm generates perfectly correlated noise estimates, but in real 

situations input and output noise sources are unlikely to be perfectly correlated. 

However, the result in (35) provides motivation for the WJM algorithm described in the 

next section. On the positive side, the TLS algorithm does minimize a natural and 

meaningful quantity. 

yˆ

xˆ

aCw

axˆ x aCε

C y ε
12

4. WEIGHTED JOINT MINIMIZATION (WJM) 

The TLS results suggest that minimizing a weighted version of (20) could lead to a 

weighted minimum noise solution that satisfies the noise power ratio constraint r

2

 given 

in (5). Equation (35) suggests the value of this weight. Thus for the WJM algorithm we 

define 

2 2 2

WJM x y

r

a

ε = + Cε Cε (36) 

and we seek to minimize 2

WJM ε with respect to w, a, and b. Note that (36) reproduces a 

special case of the XLS problem formulation in [10]. The notation of [10] defines a 

model matrix A, a parameter vector θ, a model weight matrix Wg, and a data weight 

matrix Wx. If we define these four XLS variables to be 

[ ] 2 2 , , , and ,

a r

b a

  = == =    

A -y w θ W C WC g x 1 (37) 

and note that the noise-free input x in [10] corresponds exactly to our input w, then the 

cost function C( ) x,θ given in equation (11) of [10] becomes our equation (36) exactly: 

 

() ( ) ( )

2 2 2

1

, 1

.

T

y x WJM

C

r

a

ε

 

= +− −        

= +=

T T

g x x θ θ AW A x x W x x

θ

Cε Cε

 

 (38) 

Because the WJM problem is constrained significantly relative to the general XLS 

problem, our procedure for minimizing our equation (36) is very different from the 

procedure given for minimizing (11) in [10]. 

Our solution steps parallel those in the TLS algorithm, but with some additional 

complexity introduced by the weighting factor. We first find the necessary conditions 

( ) 1 sign( ) 1

ra b

r a

=+ −     +

wx y 1 (39) 

and 

2

1

.

n

y ii

i

b m a cw == − ∑ (40) 
13

Using (39) in (40) and solving for b yields 

 . WJM y WJM x b mam = − (41) 

We can now use results (39) and (41) to rewrite the noise estimates and 2

WJM ε : 

( ) ( ) sign( ) , 1 x xy

r a

am m

r a

− = − −−     +

ε x y 1 1 (42) 

( ) ( ) 1 , 1 y xy am m

r a = − −−     +

ε x y 1 1 (43) 

( )

2 2 ˆ ˆ . 1 WJM

r

a

a ra

ε = − +x y (44) 

Equation (44) allows us to minimize 2

WJM ε with respect to the single variable a. The 

resulting necessary and sufficient condition is 

 

sign( ) ˆ . ˆ ˆˆ sign( ) WJM a

r r

ρ

ρ = + −

y

x y y

 (45) 

From (42) and (43) we see that the total cost-weighted noise is 

( )

2

2 2 2

2

1 ˆ ˆ

1 x y WJM

WJM

r

a

r a

+

+= −

+

Cε Cε x y (46) 

and the noise power ratio constraint is satisfied: 

2

2

2 . x

y

= r

Cε

Cε

 (47) 

In the simple case where the noise power ratio is unity (r = 1), and both cost-weighted 

shifted observations have unit power ( 2 2

x y ˆ ˆ = =1), (45) reduces to 

 

( )

sign( )

21 1

WJM a ρ

ρ

=

+ −

 (48) 

and (46) reduces to 
14

( ) 2 2

2 2 21 . x y ρ   + =− +   Cε Cε (49) 

Equations (42) and (43) make it clear that the WJM algorithm generates weighted noise 

estimates Cεx and Cεy that are perfectly correlated. On the other hand, the WJM 

algorithm does not result in any general orthogonality relations. The geometry of the 

WJM result in a simple two-dimensional case with mx = m y = 0, r ≈ 2, a ≈ 0.75, is 

shown in Figure 5. 

Figure 5. Example geometry of WJM and CJM results. 

The WJM algorithm minimizes a weighted version of the total noise power that is less 

contrived than the quantity minimized in the LSMS algorithm, but not as natural as the 

quantity minimized in the TLS algorithm. The WJM algorithm does satisfy the noise 

power ratio constraint. On the other hand, the WJM algorithm generates perfectly 

correlated noise estimates, but in real situations input and output noise sources are 

unlikely to be perfectly correlated. 

5. CONSTRAINED JOINT MINIMIZATION (CJM) 

In order to minimize a more natural noise quantity (as in TLS) and yet enforce the noise 

power ratio constraint in (5), we invoke a Lagrange multiplier λ and perform a 

constrained minimization on the total cost-weighted noise power 

2 2 2 2 2 2 ,0 . CJM x y x y ε λ λ r   = ++ − <     Cε Cε Cε Cε (50) 

In the following we seek to minimize 2

 CJM ε with respect to w, a, and b and this leads to 

the CJM algorithm. The solution steps partially parallel those in the TLS algorithm, but 

some significant differences are introduced by the Lagrange multiplier. 

yˆ

xˆ

aCw

axˆ x aCε

C y ε
15

The necessary condition 

2 0, 1to CJM

i

i n

w

ε

∂ = = ∂ (51) 

results in 

( )( ) ( ) ( ) ( ) 2 2 2 2 2 2 22 1 1 , 1to , iii i i i + − =− − − = λ λ c w x r a c aw y b i n (52) 

and (52) can be solved for w to get 

( ) ( ) ( ) ( ) ( ) 2

2 2

1 11 .

1 1

ra b

r a

λ λ

λ λ

= + +− −     + +−

w xy 1 (53) 

We can sum (52) over i and invoke the weighted noise power ratio constraint in (5) to get 

the relationship 

( ) ( ) 2 1

1 .

r a

r

λ

λ − + =± (54) 

Using (54) in (53) allows us to eliminate λ, giving 

( ) 1 1

1

1 br ar , when .

a r

− −

− = − ± ≠±     ± w yx 1 (55) 

The two sign choices in (55) come from the single sign choice in (54) and thus they must 

match. In other words, (55) describes two, not four possible solutions for w. 

Another necessary condition for minimizing 2

CJM ε is 

2 0 CJM b

ε

∂ = ∂ (56) 

and this results in 

2

1

.

n

y ii

i

b m a cw == − ∑ (57) 

Using (55) in (57) gives 

 . CJM y CJM x b mam = − (58) 
16

We can now use results (55) and (58) to rewrite the noise estimates 

1 ( ) ( ) 1 , x xy am m

a r −

− = − −−     ±

ε x y 1 1 (59) 

( ) ( ) ( ) 1

1 . y xy am m

ra r −

± = − −−     ±

ε x y 1 1 (60) 

Note that all three of the sign choices in (59) and (60) must match. From (59) and (60) it 

is clear that the cost-weighted noise power ratio constraint in (5) is satisfied by the 

necessary conditions derived so far. Further, the total cost-weighted noise power is 

( ) ( )

( )

2

2 22 2 2

2 1

1

1 . ˆ ˆ xy x

r

r a

a r

−

−

−

+

+ =+ = −

±

Cε Cε Cε x y (61) 

We can minimize (61) by minimizing 2

C x ε with respect to the single variable a. The 

necessary condition 

2

0 x a

∂ = ∂

Cε (62) 

yields 

ˆ

ˆ

ˆ

ˆ

CJM

r

a

r

ρ

ρ

±

=

±

y

x

x

y

 (63) 

and the two sign choices must match. From (63) it follows that 

2

1

2

1 ˆ ˆ

ˆˆ ˆ CJM T

r

a r

r r

− ±

± =⋅

±

x y

x y x

 (64) 

and 

( )

( )

22 2 2

2

2 2

ˆˆ ˆ ˆ 1 ˆ ˆ

ˆˆ ˆ CJM

T

r

a

r

− ± ρ − =

±

x y x y

x y

xy x

 (65) 
17

with all sign choices matching. Thus we can rewrite (61) as 

( ) ( )( ) 2 2 2 2

2 2 2 2

2 2

ˆ ˆ 1 1

1 . ˆ ˆˆ ˆ 2 xy x T

r

r

r

ρ − − +

+ =+ =

± +

x y Cε Cε Cε

x xy y

 (66) 

From (66) it is clear that the total cost-weighted noise power will be minimized by 

picking the positive sign whenever 0 ˆ ˆ T < x y and picking the negative sign otherwise. 

This allows us to eliminate the sign choices from previous results (55), (59), (60), (63), 

(65), and (66) respectively: 

( ) 1 1

1

1 sign( ) , when sign( ) 0, sign( )

b r ar

a r

ρ ρ

ρ

− −

− = −+ + ≠     +

w yx 1 (67) 

1 ( ) ( ) 1 , sign( ) x xy am m

a r ρ −

− = − −−     +

ε x y 1 1 (68) 

( ) ( ) ( ) 1

sign( ) , sign( ) y xy am m

ra r

ρ

ρ − = − −−     +

ε x y 1 1 (69) 

ˆ

ˆ sign( ) , ˆ

ˆ

CJM

r

a

r

ρ

ρ

ρ

+

=

+

y

x

x

y

 (70) 

( ) ( )

( ) ( )

2 2 2 2

2

2 2

ˆˆ ˆ ˆ 1 sign ˆ ˆ ,

ˆˆ ˆ sign

CJM

T

r

a

r

ρ ρ

ρ

− + − =

+

x y x y

x y

xy x

 (71) 

( ) ( )( )

( )

2 2 2 2

2 2 2 2

2

ˆ ˆ 1 1

1 . ˆ ˆ sign xy x

r

r

r

ρ

ρ

− − +

+ =+ =

+

x y Cε Cε Cε

x y

 (72) 

Now we can use (70) to show that 

( ) 2

1

2

1 ˆ ˆ sign sign( ) sign( ) 0 ˆ ˆˆ CJM T

r

a r

r r

ρ

ρ ρ − +

+= ≠

+

x y

x xy

 (73) 

 

so the condition in (67) is satisfied. 
18

To show that (70) is a sufficient condition for the minimization of 2

C x ε we must show 

that 

2

2

2 0 .

CJM

x

a a

∂

<

∂

Cε (74) 

Using (73) and (71) we can show that 

( )

( )

( )

( ) ( )

2 2 2 1 2

2

2 4 1

22 2 2 2

2 4 2 1 2

2 sign( ) 2 ˆ ˆˆ

sign( )

2 sign( ) 2 ˆ ˆ ˆ ˆ ˆˆ ˆ

0

sign( ) ˆ ˆˆ

CJM

CJM CJM

x

CJM

T

T

CJM

a

a ra

a a r

r rr

a rr r

ρ

ρ

ρ ρ

ρ

−

−

−

∂ + −− = ∂ +

+ ++

= >

+ +

x x y Cε

x y x x xy y

x xy

 (75) 

and we conclude that (70) is both a necessary and sufficient condition for the 

minimization of 2

C x ε and hence for the minimization of 2 2Cεx y + Cε , resulting in the 

total cost-weighted noise level given in (72). 

Note that the CJM algorithm generalizes the two conventional least-squares approaches 

to fitting xˆ and yˆ and it reproduces the conventional least squares solutions in the 

limiting cases r = 0 and r → ∞. When r = 0, (70) gives 

ˆ

ˆ CJM a = ρ

y

x

 (76) 

which is the standard result for the least-squares problem 

aˆ ˆ = + . y x y ε (77) 

When r → ∞, (70) gives 

1 ˆ

ˆ CJM a

ρ = yx

 (78) 

which is the standard result for the least-squares problem 

1 ˆ ˆ . a = + x y x ε (79) 
19

Note however that (70) is not just a linear combination of (76) and (78). 

In the simple case where both cost-weighted shifted observations have unit power 

( 2 2

x y ˆ ˆ = =1) and r = 1, (70) reduces to 

 sign( ) CJM a = ρ (80) 

 and (72) reduces to 

2 2

1 Cεx y + Cε = − ρ (81) 

which are the same as the TLS results for this simple case. 

Equations (68) and (69) make it clear that the CJM algorithm generates weighted noise 

estimates Cεx and Cεy that are perfectly correlated. On the other hand, the CJM algorithm 

does not result in any general orthogonality relations. The geometry of the CJM result in 

a simple two-dimensional case with mx = m y = 0, r ≈ 2, a ≈ 0.75, is shown in Figure 

5. The CJM algorithm satisfies the noise power ratio constraint (5) and it minimizes a 

very natural quantity: the total cost-weighted noise power (72). On the negative side the 

CJM algorithm generates perfectly correlated noise estimates, but in real situations input 

and output noise sources are unlikely to be perfectly correlated. We consider it to be the 

most pleasing of the four algorithms presented so far. 

6. DIRECT ESTIMATION (DE)

A more direct estimation algorithm is possible when certain restrictions on the noises x ξ

and y ξ are satisfied. This DE algorithm results from algebraic manipulations of the noisy 

observations – it does not involve any minimization. It also differs from the other 

algorithms presented here in that the gain estimation is accomplished without explicit 

decomposition of the noisy observations xˆ and yˆ into signal and noise estimates. 

In analogy to all of the previous algorithms, we will eventually estimate the bias as in 

(58). This motivates us to start with the equivalent bias-free problem where mx and my

have been removed from the noisy observations and hence are no longer a part of the 

problem. This is described in Figure 6 and by 

x Cw y C w y C x C ˆ ˆ ˆˆ = − = − += + ( ) T T ξ x T y yT x ,, . (a a ξ ξξ ) ( ) (82) 
20

aT

+

wT

-ξx + -ξy

aTwT

x yˆ ˆ

C X X C

Figure 6. Block diagram for bias-free version of the estimation problem with mx and my 

 removed. 

We require that w and x ξ be independent, that w and y ξ be independent, and that x ξ

and y ξ also be independent. We also require that 2

x C ξ and 2y C ξ have zero mean. 

Under these conditions, we can calculate expectations over x ξ and y ξ : 

22 22 2 2 E E E, ˆ y =+ =+ ξ ξ r T xT x Cw C Cw C (83) 

2 2 2 2 E E , and ˆ T y = + a ξ T y Cw C (84) 

T 2

E( ) . ˆ ˆ T = a T x y Cw (85) 

Note also that (83) enforces the noise power ratio constraint in (2). 

We then multiply (84) by r

2

, subtract the result from (83), divide by (85), and solve for a

to find 

2 2 2 2 2

2 T

4 E E ˆ ˆ , where . 2 E( ) ˆ ˆ T

kk r r

a k

r

−± + − = = x y

x y

 (86) 

For a meaningful solution, we require sign(aT)=sign(ρ), and this requirement can be 

satisfied by 

( ) 2 2

2

sign 4 . 2 Tk kr

a

r

−+ + ρ = (87) 
21

Thus we have used the noisy observations xˆ and yˆ to calculate a value for aT. This 

calculated value of aT then naturally becomes our estimate of the system gain aDE. Thus 

the expansion of (87) yields 

2

2 2 22

2

ˆˆ ˆˆ 4 ˆˆ ˆˆ , 0, 2

0, 0.

DE

r rr

a

r

ρ

ρ

ρ

ρ

      −+ − +  =     ≠ 



 =

yx yx

xy xy (88) 

For notational clarity, (88) uses x y ˆ ˆ , , and ρ in place of 

T

2 2

2 2

E( ) ˆ ˆ E , E , and ˆ ˆ

E E ˆ ˆ

x y x y

x y

 respectively. In practice, each of these quantities 

will be estimated by conventional means using all available observations. 

We estimate the bias b as before: 

 . DE y DE x b m am = − (89) 

 

Note that when r =1, the result for aDE in (88) reduces to the TLS result for aTLS in (30). 

Thus for the special case of estimating a single system gain and a bias we might think of 

the DE algorithm as a generalization of the TLS algorithm. This is somewhat unexpected 

since the TLS algorithm is based on the minimization of an error expression but the DE 

algorithm is not. 

We have required wT and x ξ , wT and y ξ , and x ξ and y ξ , to be pairwise independent. 

From a geometric viewpoint these independence constraints correspond to orthogonality 

constraints, so no 2-dimensional example can be drawn. Figure 4 does show two of the 

three required orthogonality constraints, so we might consider this to be a 2-dimensional 

projection of a simple 3-dimensional DE example. In this figure we can view the DE 

algorithm operation geometrically as swinging aCw and sliding a to find the unique 

solution that solves the noise power ratio constraint given in (2). 

The DE algorithm generates estimates of a and b without any minimization. The noise 

power ratio constraint is directly incorporated into the DE algorithm derivation and thus 

it is always satisfied. We consider the DE and CJM algorithms to be the most pleasing 

algorithms for estimating gain and bias. We do note that the DE algorithm does not 

provide decompositions of noisy observations into signal and noise estimates. 
22

7. EXAMPLE SIMULATION RESULTS AND DISCUSSION 

This section contains the results of computer simulations of the five gain and bias 

estimation algorithms. The bias estimation result is identical in the five algorithms. 

Further, the bias issue is effectively removed from the problem before gain estimation. 

Thus we are able to simplify this section by treating only the case b=0, with no loss of 

generality in terms of the estimation of a. 

Figure 7 shows gain estimates for the five algorithms for the case a a dB = 20log 0 10 ( ) = , 

−≤ = ≤ 40 20log 40 r r dB 10 ( ) . Note that rdB controls the noise generators in our 

simulations and thus this is the true value of rdB. All algorithms use this true value of rdB

as an input as well. The input w is a chirp that sweeps from a normalized frequency of 

40

π

 to 

4

π radians/sample. The noise vectors x ξ and y ξ contain white Gaussian noise and 

are uncorrelated to w and to each other. For this problem we define an intuitive total 

signal-to-total noise ratio (SNR) as 

2 2

10 2 2 10log

x y

a

SNR

  + =  

  +  

w w

ε ε

 (90) 

and set it to 10 dB. All vectors have length n=4096 and we use 1

n n×n I for the cost 

function C. 

For the example shown in Figure 7, all five gain estimates show similar variances but 

they show very different means. As rdB gets large and negative, the LSMS mean 

estimation error vanishes. This is expected since this case corresponds to noise-free input 

observations and that is the situation that least-squares is meant to handle. The TLS 

algorithm has small mean estimation error only around rdB = 0. This is also expected 

since TLS tries to minimize both noise components equally. The WJM algorithm has 

small mean estimation error only when rdB is large and positive. The CJM algorithm has 

small mean estimation error near rdB = 0 and when |rdB| gets large, where it tends toward 

one of two different least-squares solutions. The DE algorithm has small mean 

estimation error across the full range shown. As expected, the CJM and DE algorithms 

give the best estimates across a wide range of r values. This result stands over a good 

range of SNR values, a values, N values, and input signal types including other chirps, 

tones, speech, and noise. Thus through the remainder of this section, we present 

simulation results for only the CJM and DE algorithms. 

Figures 8-16 show example errors in gain estimates for true gain values between -20 and 

+20 dB. The input and noises of Figure 7 are used again. The results are almost identical 

for the other signal types mentioned above. As before, n=4096. Figures 8-10 use 

SNR=20 dB, Figures 11-13 use SNR=10 dB, and Figures 14-16 use SNR=0 dB. Within 

each triple of figures we use three values of rdB :-10, 0, and +10 dB. 
23

Figure 7. Example results for five estimation algorithms; true gain is 0 dB. 

Figure 8. Example estimation errors for the case SNR=20 dB, r =-10 dB. 
24

Figure 9. Example estimation errors for the case SNR=20 dB, r =0 dB. 

Figure 10. Example estimation errors for the case SNR=20 dB, r =10 dB. 
25

Figure 11. Example estimation errors for the case SNR=10 dB, r =-10 dB. 

Figure 12. Example estimation errors for the case SNR=10 dB, r =0 dB. 
26

Figure 13. Example estimation errors for the case SNR=10 dB, r =10 dB. 

Figure 14. Example estimation errors for the case SNR=0 dB, r =-10 dB. 
27

Figure 15. Example estimation errors for the case SNR=0 dB, r =0 dB. 

Figure 16. Example estimation errors for the case SNR=0 dB, r =10 dB. 
28

As expected, Figures 8-16 show that the estimation errors are smaller at higher SNR

values. These figures also indicate that DE gain estimation errors have nearly zero mean. 

The CJM gain estimation errors have zero mean only when dB dB a r ≈ − . The CJM 

algorithm overestimates gain when dB dB a r < − and it underestimates gain when 

dB dB a r > − . Both algorithms show a smaller estimation error variance near 0 dB a = , and 

that variance increases as dB a increases. It also appears that the DE estimation error 

variance is smaller when dB dB a r ≈ − than when dB dB a r ≈ . Finally, we note that the CJM 

and DE estimation error variances are generally similar. The exception comes at lower 

SNR values (e.g. SNR=0 dB) where the CJM error variance is markedly smaller than the 

DE error variance when 0 dB a  . 

Figure 17 provides an example correlated noise case. Here the noise vectors x ξ and y ξ

contain a mixture of harmonic distortion (correlated to w) and white Gaussian noise 

(uncorrelated to w). For this example SNR = 20 dB and r = 0 dB. The chirp of Figure 7 is 

used for the signal w. Across most of the range of true gain values, the CJM algorithm 

provides a smaller mean estimation error than the DE algorithm. We note here that the 

DE algorithm is explicitly predicated on an uncorrelated noise environment while the 

CJM algorithm is not explicitly predicated on an uncorrelated noise environment. In fact, 

while the CJM algorithm makes no assumptions about noise correlation, it does generate 

noise estimates that are perfectly correlated. Thus one might expect the CJM algorithm 

to be better suited to correlated noise environments. However, this expectation is not 

confirmed by our simulations. For some correlated noise environments, the DE algorithm 

performs better than the CJM algorithm. Results for correlated noise environments show 

modest dependence on signal type but strong dependence on noise type. 

Figure 18 shows example results for a case where the CJM and DE algorithms are given 

inaccurate values of r. This situation would occur when the actual r value in the noisy 

observations is different from an estimated, measured, or modeled r value. The signal 

and noises of Figure 7 are used, true rdB=0 dB, SNR=10 dB, and n=4096. These 

conditions would reproduce Figure 12, except that the estimation algorithms were given 

inaccurate values of r. The traces labeled “r 10 dB low” result from r =-10 dB and the 

traces labeled “r 10 dB high” result from r =10 dB. For these two situations, the CJM 

algorithm provides better estimates across most of the range of true gain values used in 

the simulation and this result is largely invariant to signal type. But this is not a general 

result. For some other r inaccuracies and true gain values, the DE algorithm will provide 

better estimates than the CJM algorithm. 
29

Figure 17. Example estimation errors for correlated noise environment. 

Figure 18. Example estimation errors for the case of inaccurate r values. 
30

From these results we conclude that the DE algorithm would generally be the best choice 

when noise is uncorrelated and r is well-known. When noise is correlated, or r is not 

well known, one should do simulations of the specific situation to determine if the CJM 

or DE algorithm is more likely to provide correct estimates. Note also that there are bias 

vs. variance trade-offs at work here. For example Figure 15 shows that when SNR = 0 

dB, and r = 0 dB, if the true gain is 15 dB, then the CJM algorithm gives a high bias 

(high mean error) low variance estimate of that true gain, while the DE algorithm gives a 

low bias, high variance estimate of that true gain. 

Finally, we mention here three additional algorithm properties. First, we note that we 

have been able to show that sign(a)=sign(ρ) for all five algorithms and this phase 

consistency is certainly a desirable property. Second we note a scaling property 

displayed by the LSMS, WJM, CJM and DE algorithms. For these four algorithms, we 

have been able to show that if the inputs x, y, and r yield the gain estimate a, then the 

inputs x, αy, and r/α yield the gain estimate αa. Third we note that the LSMS, TLS, 

WJM, and CJM algorithms decompose noisy observations into signal and noise estimates 

and thus have potential use in denoising applications. 

8. SUMMARY

We have provided background context and references for the general problem of 

identifying linear systems from noisy input and output observations and have noted that 

the existing work in this area allows for identification of multivariate systems and thus is 

fairly complex. We introduced the important special case of the estimation of just a 

system gain and bias from noisy input and output observations and added a single piece 

of side-information to the problem: the noise power ratio, r. We described the potential 

utility of solutions to this special case of the problem. We then derived five different 

solutions. 

The LSMS algorithm comes from a conventional least-squares step followed by the 

manual splitting of the least-squares residual into two vectors. This algorithm satisfies 

the noise power ratio constraint but the gain and bias estimates do not directly depend on 

r, and this is inconsistent with our desire to use that side-information advantageously. An 

additional disadvantage is that the quantity minimized by the LSMS algorithm is 

somewhat unnatural. 

The TLS algorithm minimizes a total noise power which is a more natural quantity, but it 

does not make use of r, and thus the noise power ratio constraint is not satisfied in 

general. We derived our TLS algorithm for this specific problem from first principles. 

Our algorithm has only a single step and does not require a SVD. The TLS algorithm 

does provide motivation for the weighting factor used in the WJM algorithm. 

The WJM algorithm minimizes a weighted version of the total noise power that is more 

natural than the quantity minimized in the LSMS algorithm, but not as natural as the 

quantity minimized in the TLS algorithm. The WJM algorithm does satisfy the noise 

power ratio constraint. 
31

The CJM algorithm derivation invokes a Lagrange multiplier which allows for the 

simultaneous minimization of a natural quantity while enforcing the noise power ratio 

constraint. This algorithm generalizes two conventional least-squares solutions and it 

reproduces conventional least squares solutions in the limiting cases r = 0 and r → ∞. 

The DE algorithm differs from the other four algorithms in that its derivation does not 

include the minimization of any noise quantity. Rather the DE algorithm results from 

algebraic manipulations of the noisy observations. The DE algorithm requires additional 

assumptions on the observation noises. Unlike the other algorithms, the DE algorithm 

does not provide decompositions of noisy observations into signal and noise components. 

The DE algorithm can be interpreted as a generalization of the TLS algorithm, even 

though they were derived in very different ways. 

From a mathematical perspective, we observed that the CJM and DE algorithms seemed 

most satisfying. A fundamental distinction between these two algorithms lies in the area 

of noise correlation. While the DE algorithm is explicitly predicated on an uncorrelated 

noise environment, the CJM algorithm actually produces noise estimates that are 

perfectly correlated. We conducted simulations that confirmed that these two algorithms 

would be most likely to give correct estimates of system gain. In our simulations we 

found that the CJM and DE estimation error variances are generally similar, but the 

estimation error means can be quite different. We concluded that when noise is 

uncorrelated and r is well-known, the DE algorithm would generally be the best choice. 

When noise is correlated, or r is not well known, further simulations of the specific 

situation are required to determine if the CJM or DE algorithm is more likely to provide 

correct estimates. 

We have also provided simplified 2-dimensional geometric interpretations for the LSMS, 

TLS, WJM, and CJM algorithms, noting the orthogonal and parallel vector relationships 

that would extend to higher dimensional examples. For convenience, the algorithm 

details are summarized in tabular form in the appendix. 

9. REFERENCES 

[1] J.K. Tugnait, “Stochastic system identification with noisy input using cumulant 

statistics,” IEEE Transactions on Automatic Control, vol. 37, no. 4, pp. 476-485, 

Apr. 1992. 

[2] T. Söderström, “Identification of stochastic linear systems in presence of input 

noise,” Automatica, vol. 17, no. 5, pp. 713-725, Sep. 1981. 

[3] G.H. Golub and C.F. Van Loan, Matrix Computations, Baltimore, MD: Johns 

Hopkins University Press, 1989. 

[4] L.L. Scharf, Statistical Signal Processing: Detection, Estimation, and Time Series 

Analysis, Reading, MA: Addison-Wesley, 1991. 

[5] J. Rissanen, “Estimation of errors-in-variables models,” Proceedings of the 27th

IEEE Conference on Decision and Control, Austin, Texas, Dec. 1988, vol. 3, pp. 

1828-1830. 
32

[6] J.M.M. Anderson and G.B. Giannakis, “Noisy input/output system identification 

using cumulants and the Steiglitz-McBride algorithm,” IEEE Transactions on 

Signal Processing, vol. 44, no. 2, pp. 1021-1024, Feb. 1996. 

[7] J. K. Tugnait, “Identification of multivariable stochastic linear systems using 

integrated polyspectrum given noisy input-output data,” Proceedings of the IEEE 

Signal Processing Workshop on Higher-Order Statistics, Banff, Alta, Canada, Jul. 

1997, pp. 239-244. 

[8] S. Chandrasekaran, G.H. Golub, M. Gu, and A.H. Sayed, “Parameter estimation 

in the presence of bounded modeling errors,” IEEE Signal Processing Letters, 

vol. 4, no. 7, pp. 195-197, Jul. 1997. 

[9] A.H. Sayed and S.Chandrasekaran, “Parameter estimation with multiple sources 

and levels of uncertainties,” IEEE Transactions on Signal Processing, vol. 48, no. 

3, pp. 680-692, Mar. 2000. 

[10] A. Yeredor, “The extended least squares criterion: Minimization algorithms and 

applications,” IEEE Transactions on Signal Processing, vol. 49, no. 1, pp. 74-86, 

Jan. 2001. 

[11] A. Yeredor, “The joint MAP-ML criterion and its relation to ML and to extended 

least squares,” IEEE Transactions on Signal Processing, vol. 48, no. 12, pp. 

3484-3492, Dec. 2000. 

[12] L.V. Hedges and I. Olkin, Statistical Methods for Meta-Analysis, Orlando, 

Florida: Academic Press, 1985. 

[13] S. Wolf, “Measuring the end-to-end performance of digital video systems,” IEEE 

Transactions on Broadcasting, vol. 43, no. 3, pp. 320-328, Sep., 1997. 

[14] C.J. van den Branden Lambrecht, D.M. Costantini, G.L. Sicuranza, and M. Kunt, 

"Quality assessment of motion rendition in video coding," IEEE Transactions on 

Circuits and Systems for Video Technology, vol. 9, no. 5, pp. 766-782, Aug. 1999. 

[15] S.D. Voran, “Objective estimation of perceived speech quality, part II: Evaluation 

of the measuring normalizing block technique,” IEEE Transactions on Speech 

and Audio Processing, vol. 7, no. 4, pp. 383-390, Jul. 1999. 

[16] A.W. Rix, J.G. Beerends, M.P. Hollier, and A.P. Hekstra, “Perceptual evaluation 

of speech quality (PESQ) – A new method for speech quality assessment of 

telephone networks and codecs,” Proceedings of the IEEE International 

Conference on Acoustics, Speech and Signal Processing 2001, Salt Lake City, 

May 2001, vol. 2, pp. 749-752.
33

APPENDIX: SUMMARY OF ESTIMATION ALGORITHMS 

 In this appendix we summarize the five gain and bias estimation algorithms in 

tabular form. Table A-1 summarizes the mathematical constructions that lead to each of 

the five algorithms. Table A-2 gives the preparation steps that are common to all five 

algorithms. Table A-3 shows how each of the five algorithms calculates estimates of a

and b. This table also gives the resulting noise power ratio. Tables A-4 and A-5 provide 

additional algorithm results that may be of interest. 
Table A-1. Mathematical Constructions Leading to the Five Algorithms 

Algorithm Assumptions beyond Figure 1 Quantity Minimized Additiona

Least squares plus manual 

splitting (LSMS) 

None ( ) 2 2LSMS y x ε = − C ε εa 

x

y

Cε

Cε

Total least squares (TLS) None 2 2 2

TLS x y ε = + Cε Cε N

Weighted joint 

minimization (WJM) 

None 2 2 2

WJM x y

r

a

ε = + Cε Cε N

Constrained joint 

minimization (CJM) 

None 2 2 2

CJM x y ε = + Cε Cε

x

y

Cε

Cε

Direct estimation (DE) w and x ξ are independent 

w and y ξ are independent 

x ξ and y ξ are independent 

2

x C ξ and 2 y C ξ are zero mean 

No Minimization 

E

E

ξ

ξ

C

C

34 
Table A-2. Preparation Steps Common to the Algorithms 

Step Equations

1. Calculate cost-weighted means. 2 2

1 1

,

n n

x ii y ii

i i

m cx m cy

= =

= = ∑ ∑

2. Calculate cost-weighted shifted versions of x and y. xx yy

ˆ ˆ

=− =− Cm Cm ( x y 1 1 ) , ( )

3. Calculate normalized cross correlation. 

T ˆ ˆ

ˆ ˆ ρ =

x y

x y

4. Calculate correlation sign. 1 when 0 sign( ) 1 otherwise ρ ρ  < =  −

35 
Table A-3. Basic Algorithm Results 

Algorithm Estimated Bias Estimated Gain Noise 

LSMS LSMS y LSMS x b ma m = − ˆ

ˆ LSMS

a

=

ρ

y

x

C

C

TLS TLS y TLS x b m am = − 2

2

ˆˆ ˆˆ

4

ˆˆ ˆˆ

, 0

2

0, 0

TLS

a

ρ

ρ

ρ

ρ

        −+ − +    

=     ≠ 

 =

yx yx

xy xy

C

C

WJM WJM y WJM x b mam = − sign( ) ˆ

ˆ ˆˆ sign( ) WJM

a

r r

ρ

ρ

=

+ −

y

x y y

C

C

CJM CJM y CJM x b mam = − ˆ

ˆ

sign( ) ˆ

ˆ

CJM

r

a

r

ρ

ρ

ρ

+

=

+

y

x

x

y

C

C

DE 

DE y DE x b m am = − 2

2 2 22

2

ˆˆ ˆˆ

4

ˆˆ ˆˆ

, 0,

2

0, 0

DE

r rr

a

r

ρ

ρ

ρ

ρ

      −+ − +   

=     ≠ 

 =

yx yx

xy xy

E

E

36 
Table A-4. Additional Algorithm Results 

Algorithm Weighted Input Noise Estimate

C x ε

Weighted Output Noise Estimate 

C y ε

Total Weighted No

2 

Cεx 

+ Cε

LSMS ( )[ ] sign ˆ ˆ 1 LSMS LSMS r a r a − ρ − + x y ( ) [ ] 1 ˆ ˆ 1 LSMS LSMS a r a − + x y ( ) 

2 

2 1 ˆ 

1 LSMS ra r 

+ 

+ y

TLS [ ] 2 ˆ ˆ 1 TLS TLS TLS a a a − − + x y [ ] 2 1 ˆ ˆ 1 TLS TLS a a − + x y 

2 1 ˆ 

1 TLS TLS a a 

+ x 

WJM [ ] sign( ) ˆ ˆ 1 WJM WJM r a r a − ρ − + x y [ ] 1 ˆ ˆ 1 WJM WJM a r a − + x y ( ) 2 

2 1 

1 W

WJM r a r a+ 

+ 

CJM [ ] sign( ) ˆ ˆ 1 CJM CJM r a r a − ρ − + x y [ ] 1 ˆ ˆ 1 CJM CJM a r a − + x y ( )

( ) 2 2 2 

ˆ ˆ 1 

ˆ 

sign 

ρ ρ − 

+ x yx 

DE Not Calculated Not Calculated Not Calcula

37 
Table A-5. Additional Algorithm Results 

Special Case 2 2 Algorithm Estimated Input w x y ˆ ˆ = = = 1, 1 r

Estimated Gain Total Weighted Nois

2 

Cεx y + Cε

LSMS ( ) ( )( ) 1 sign 1 LSMS r b r a + − ρ + x y 1 ρ 121 ρρ −+

TLS ( ) 2 1 1 TLS TLS a b a  + −    + x y 1 sign ( ρ ) 1− ρ

WJM ( ) 1 sign( ) 1 WJM r b r a  + − ρ    + x y 1 ( ) sign( ) 21 1 ρ+ − ρ 2 2 21( 

ρ  

− +  

CJM ( ) 1 sign( ) 1 CJM r b r a  + − ρ    + x y 1 sign ( ρ ) 1− ρ

DE No Estimate sign ( ρ ) 1− ρ

38