Experimental Evaluation of some Cross
Correlation Methods for Time-Delay Estimation
in Linear Systems
Svante Bj¨orklund Control & Communication Department of Electrical Engineering
Link¨opings universitet, SE-581 83 Link¨oping, Sweden WWW: http://www.control.isy.liu.se
E-mail: svabj@isy.liu.se 15th April 2003
AUTOMATIC CONTROL
COMMUNICATION SYSTEMS LINKÖPING
Report no.: LiTH-ISY-R-2513
Technical reports from the Control & Communication group in Link¨oping are available at http://www.control.isy.liu.se/publications.
Sammanfattning Abstract Nyckelord Keywords Rapporttyp Report: category Licentiatavhandling C-uppsats D-uppsats Övrig rapport Språk Language Svenska/Swedish Engelska/English ISBN
Serietitel och serienummer
Title of series, numbering
URL för elektronisk version
Titel Title Författare Author Datum Date Avdelning, Institution Division, department
Control & Communications
ISRN Examensarbete ISSN
X
LiTH-ISY-R-1400-3902
http://www.control.isy.liu.seDepartment of Electrical Engineering
2513
Experimental Evaluation of some Cross Correlation Methods for Time-Delay
Estima-tion in Linear Systems
In this report we study estimation of time-delays in linear dynamical systems with additive
noise. Estimating time-delays is a common engineering problem, e.g. in automatic control,
sys-tem identification and signal processing.
We have experimentally studied cross correlation methods for time-delay estimation in order to
find the "best" method or if none is best which method to use in different situations. We have
simulated systems with time-delay and studied internal quantities of the methods, the time-delay
estimates, RMS-errors of estimates, percentage of failed estimates, Analysis of Variance
(ANO-VA) and confidence intervals for different cases.
Our result is that no method is always best but one method is clearly inferior to the other. The
methods can be divided in two groups with similar properties: thresholding methods and area &
moment methods. The thresholding methods work better for fast systems than slow systems
whereas the area & moment methods works better for slow systems than fast. The thresholding
methods work better with random binary input signals than steps as input whereas the area &
moment methods works better with steps as input than random binary input.
time-delay, dead-time, estimation, system identification, linear dynamic systems, cross
correla-tion, CUSUM, impulse response, step response, area methods, moment methods, ANOVA,
con-X
2003-04-15
Contents
1 Introduction 1
2 Theory of cross correlation methods 2
2.1 Overview . . . 2
2.2 Impulse and step response estimation . . . 3
2.3 Thresholding methods . . . 4
2.3.1 Cross spectrum scale . . . 4
2.3.2 Test statistics . . . 4
2.3.3 Thresholding . . . 5
2.3.4 Change time estimation . . . 5
2.4 Phase methods . . . 5
2.5 Area and moment methods . . . 7
2.5.1 Area methods . . . 7
2.5.2 Moment methods . . . 8
2.5.3 An area and moment method with better noise properties . . . 10
2.6 Maximum likelihood methods . . . 10
3 Implementation of correlation methods 11 3.1 Impulse response estimation, estimpresp . . . 11
3.2 Step response estimation, eststepresp2 . . . 11
3.3 Thresholding of impulse response, idimp4 . . . 12
3.4 Thresholding of step response, idstep4 . . . 13
3.5 CUSUM detector on impulse response, idimpCusum3 . . . 13
3.6 CUSUM detector on step response, idstepCusum3 . . . 14
3.7 Area method, area1 . . . 14
4 Simulation setup 17
4.1 Fixed factors. . . 17
4.2 Varied factors. . . 17
5 Simple experiments 24 5.1 Impulse response estimation . . . 24
5.2 Step response estimation . . . 26
5.3 Thresholding impulse response . . . 26
5.4 Thresholding step response . . . 31
5.5 CUSUM detector on impulse response . . . 31
5.6 CUSUM detector on step response . . . 32
5.7 Area methods . . . 33
5.8 Moment methods . . . 35
6 Statistical analysis of simulations 36 6.1 Results with threshold on impulse response . . . 36
6.2 Results with threshold on step response . . . 36
6.3 Results with CUSUM detector on impulse response . . . 40
6.4 Results with CUSUM detector on step response . . . 41
6.5 Results with area methods . . . 46
6.6 Results with moment method . . . 46
7 Analysis with ANOVA and confidence intervals 52 7.1 Thresholding impulse response . . . 52
7.2 Thresholding step response . . . 54
7.3 CUSUM detector on impulse response . . . 54
7.4 CUSUM detector on step response . . . 56
7.5 Area methods . . . 56
7.6 Moment methods . . . 59
7.7 Summary of thresholding methods . . . 59
7.8 Summary of area and moment methods . . . 60
7.9 Summary of all methods . . . 65
7.10 Discussion . . . 65
8 Discussion and conclusions 71
8.1 Discussion . . . 71
8.1.1 Thresholding methods . . . 71
8.1.2 Area and moment methods . . . 71
8.1.3 Materials and methods . . . 72
8.2 Conclusions . . . 73
8.2.1 General conclusions . . . 73
8.2.2 Conclusions about specific methods . . . 73
8.3 Future work . . . 75
References 76 A ANOVA and confidence intervals 78 A.1 About ANOVA and confidence intervals . . . 78
A.2 ANOVA and confidence intervals to find the best method . . . 78
A.3 Validation of ANOVA of thresholding methods . . . 79
A.4 Validation of ANOVA of area and moment methods . . . 80
1
Introduction
The problem we address in this report is estimating time delays in linear dynamical systems with additive noise. A synonym for time delay is dead-time. Estimating time delays is a common engineering problem, e.g. in control performance monitoring of industrial processes [Hor00, Swa99], in design and tuning of controllers, in range estimation in radar [KQ92] and in direction estimation by time delay of arrival in signal intelligence [HR97, Wik02]. Dead-time estimation is also a necessary part in all system identification [Lju99]. The purpose with this work is to test and evaluate cross correlation methods for time-delay estimation, especially with automatic control applications in mind. Particularly interesting is it to determine the best method. Is one method best in all situations or should different methods be used for different situations?
In the next chapter, we group cross correlation methods in different groups and give a short theoretic description of them. A major part of this reports consists of simulations and analyzes of time-delay estimation methods. In Chapter 3, the used implementation of the methods is described in detail. Chapter 4 contains a description of the setup of the simulations. The methods have been analyzed with three different techniques. The first technique, used in Chapter 5, consists of simple experiments. The second is Monte Carlo simulations and studying of the RMS estimation error etc for different cases. The third technique, which is the most complicated is employed in Chapter 7. It consists of Analysis of Variance (ANOVA) and confidence intervals. Its predominant feature is that it can give conclusions with a degree of confidence. Following, Chapter 8 contains more discussion, conclusions and suggestions for further work. After a literature reference list, AppendixA contains a discussion about required prerequisites for the ANOVA and the confidence intervals to be valid.
2
Theory of cross correlation methods
All time-delay estimation methods we are studying in this report start by estimating the cross correlation between the output and input signals. If the input is white noise this cross correlation will be equal to the impulse response of the system. By integration we then can get the step response.
This chapter starts in Section 2.1 by grouping cross correlation methods into groups. Then in Section 2.2-2.6 we briefly give theory for different cross correlation methods.
2.1 Overview
Cross correlation methods can be divided in several types:
• Thresholding methods. These methods, here called thresholding methods, start with a cross correlation estimation, then uses three factors to modify the estimate and conclude with a thresholding. The factors are response type (impulse or step), cross spectrum scale (linear or logarithmic) and test statistics (direct or cumulative sum). When the factors are combined, eight combinations are obtained:
1. Cross correlation and thresholding.
2. Cross correlation, integration to step response and thresholding.
3. Cross correlation, FFT, logarithm, IFFT and thresholding. This resembles cepstrum methods for time delay estimation. Compare with [DC02, GLM01]. 4. Cross correlation, integration to step response, FFT, logarithm, IFFT and
thresholding.
5. Cross correlation and CUSUM [Gus00, GLM01] (inclusive thresholding). 6. Cross correlation, integration to step response and CUSUM (inclusive
thresh-olding).
7. Cross correlation, FFT, logarithm, IFFT and CUSUM (inclusive thresholding). 8. Cross correlation, integration to step response , FFT, logarithm, IFFT and
CUSUM (inclusive thresholding).
All these methods are different ways to discover the peak corresponding to the time-delay in the cross correlation function. Since the thresholding is conducted in the time domain, we could also call the methods time domain methods.
• Phase methods.
1. Cross correlation and the slope of FFT.
Since the estimation is conducted in the frequency domain, we could also call the methods frequency domain methods.
• Moment and area methods.
1. Cross correlation and moments.
• Maximum likelihood methods.
1. Cross correlation and maximum likelihood.
2. Cross correlation, integration to step response and maximum likelihood.
The moment and area methods and the maximum likelihood methods can also be used with impulse responses estimated by other methods than cross correlation.
2.2 Impulse and step response estimation
By correlation analysis the impulse response h(t) can be estimated by
h(τ ) = Ryu(τ )
λ (2.1)
when the input signal is white noise with auto correlation function Ru(τ ) = E{u(t + τ)u(t)} =
λδ(τ ) [Lju99, p. 170]. The quantity λ is the power of the input signal. An estimate of the cross correlation function Ryu(τ ) = E{y(t + τ)u(t)} can be obtained by
ˆ Ryu(τ ) = 1 N N X t=τ ytut−τ. (2.2)
An estimate of the input signal power is ˆ λ = 1 N N X t=1 u2t. (2.3)
The estimate of the impulse response is then
ˆ
h(τ ) = Rˆyu(τ ) ˆ
λ . (2.4)
Correlation analysis is used in matched filtering, in for example communication and radar. There the received signal is correlated with the sent signal. When the filter output is exceeding a threshold, the decision is made that this is indeed the sent signal.
A different way of estimating the impulse response is to estimate a parametric FIR model of the system. The model will be a linear regression and the estimate is given by a least-squares estimate:
ˆ θN = " 1 N X ϕtϕTt k | {z } # A −1 1 N X ϕtyt | {z } B (2.5)
with ϕt = [ut, ut−1, ...] [Lju99, p. 204]. This is essentially the same as equation (2.4). If
the input signal ut is white, then the matrix A in equation (2.5) will be A = ˆλ I and the
matrix B will be B = ˆ Ryu(0) ˆ Ryu(1) .. . (2.6)
and ˆθN = h0 h1 . . . is the same as ˆh(τ ) = h0 h1 . . . in equation 2.4.
The estimate ˆθN will be Gaussian distributed when the noise is Gaussian. Even if the
noise is not Gaussian, ˆθN will often be asymptotic Gaussian when N → ∞ [Lju99, p.
556].
When we have an estimate of the impulse response, we can get an estimate of the step response by simulating the estimated system (represented by its impulse response) with a step as input signal. This is the same as integrating the impulse response. The numeric integration of the impulse response estimate can be done in different ways. If the impulse response coefficients to integrate are Gaussian distributed then the step response will also be Gaussian since a sum of Gaussian random variables is Gaussian.
2.3 Thresholding methods
The principle of thresholding methods is to estimate the time point of the peak of the cross correlation or the start of the rise to the peak (which depends on if the output signal is only delayed or also transformed by a dynamic system). This is performed in the time domain.
2.3.1 Cross spectrum scale
Logarithmic cross spectrum Transform the estimated impulse or step response to the frequency domain and apply a logarithm and then transform back to the time domain. This is similar to Cepstrum methods for time-delay estimation [DC02, GLM01].
Linear cross spectrum Do not logarithm in the frequency domain.
2.3.2 Test statistics
Cumulative sum (CUSUM) Perform an averaging operation of the estimated esti-mated impulse or step response by CUSUM [Gus00] before the thresholding.
2.3.3 Thresholding
Fix threshold The threshold is fixed and data independent.
Threshold depending of the uncertainty of the response If we let the threshold be proportional to the uncertainty of the impulse or step response we will get a specified risk for false alarm (saying that the time-delay is over when it is not). With such a threshold, if the threshold never is crossed, then we can not assert that there is a connection between input and output.
Max threshold for detection Using a “max threshold for detection” is the same as finding the highest peak of the cross correlation function. This is a suitable approach when the output signal is a pure time delay of the input signal but not if there also is a dynamic system between input and output.
Constant False Alarm Ratio (CFAR) With CFAR the threshold is selected to give a constant false alarm rate by looking at the values of the estimated cross correlation function close to the time delay that is being tested. These close-by values are used in an estimate of the noise level. Compare with the Section “Threshold depending of the uncertainty of the response” above. This method is common in radar.
2.3.4 Change time estimation
When a change in the estimated estimated impulse or step response has been detected, the change time can be estimated by for example taking the peak of an interpolation of the estimated cross correlation function. See [Gus00] for estimating the change time when using CUSUM.
2.4 Phase methods
In the phase methods we estimate the slope of the phase of the cross spectrum in the frequency domain. This can be seen as dual methods to the thresholding methods where the estimation is conducted in the time domain.
If ¯G(s) is a linear continuous-time system then
d
dωarg ¯G(iω) = 0 if the system do not contain a time-delay and
d
dωarg ¯G(iω) =−Td if the system contains a time delay Td.
Since we have sampled data we would like to use the sampled frequency function G(eiωT)
instead of the continuous-time frequency function ¯G(iω). We are lucky as for low frequen-cies G(eiωT) does not differ much from the true frequency function ¯G(iω) [LG91, p. 74]. A rule of thumb is that the agreement is good for frequencies up to 1/10 of the sampling frequency. The sampled frequency function G(eiωT) can be estimated by nonparamet-ric methods like spectral analysis or by parametnonparamet-ric methods with arbitrary linear model structures [Lju99, LG91].
In time-delay estimation by spectral analysis, we utilize the formula [GLM01] Φyu(ω) =
G(eiωT)Φ
u(ω) or
G(eiωT) = Φyu(ω) Φu(ω)
, (2.7)
where the sampled frequency function G(eiωT) is connected with the cross spectrum Φyu(ω)
between output y(t) and input u(t) and the (auto) spectrum Φu(ω) of the input. A natural
estimate of G(eiωT) is then to use the estimated cross spectrum ˆΦ
yu(ω) the estimated (auto) spectrum ˆΦu(ω): ˆ G(eiωT) = Φˆyu(ω) ˆ Φu(ω) . (2.8)
Then the phase of ˆG(eiωT) is studied to give an estimate of the time-delay
ˆ Td=−
d
dωarg ˆG(e
iωT)
where the derivative is approximated in a suitable way.
Time-delay estimation by spectral analysis is equivalent to time-delay estimation by cross correlation analysis in the frequency domain. Since the cross spectrum Φyu(ω) is the
discrete-time Fourier transform of the cross correlation function Ryu(τ ) [GLM01, p. 54],
Φyu(ω) = Ts ∞
X
k=−∞
Ryu(k)eiωkTs,
and in the same way the auto spectrum the Fourier transform of the auto correlation function Ru(τ ), we see that equation (2.8) is the same as equation (2.1) in the frequency
domain. In the time domain (equation (2.1)) a time delay in the system will be a time delay in h(τ ). In the frequency domain (equation (2.8)) a time delay in the system will be a phase shift in ˆG(eiωT). For an example of time-delay estimation by nonparametric
methods see [Wik02].
If the model structure in a parametric method is a linear regression then the model can be estimated by cross and auto correlations. Compare section 2.2. For some examples of time-delay estimation by parametric methods see [Hor00, IHD01, Bj¨o02].
2.5 Area and moment methods
In [˚AH95] some methods for identifying simple process models of the kind
G(s) = K 1 + sTe sL (2.9) and G(s) = K (1− sT )2e sL. (2.10)
are described. Part of this identification is the estimation of the time-delay L, which is what we are interested in here.
2.5.1 Area methods
Area methods [˚AH95] use area calculations of a step response s(t). First the average residence time Tar is calculated as
Tar =
A0
K,
where K = s(∞) is the static gain and the area A0 is calculated from
A0=
Z ∞
0
(s(∞) − s(t))dt
The average residence time is a rough measure of the time for the step response to finish. Then we calculate
A1 =
Z Tar
0
For the first order model (2.9) the time constant is then given by
T = e
1A 1
K . (2.11)
The time delay is
L = Tar− T.
For the second order model (2.10) the time constant is given by
T = e
2A 1
4K and the time delay is
L = Tar− 2T.
See [˚AH95] for the derivation of these area methods. In [˚AH95] they use a real step response from a step response experiment when identifying the system with these methods. In this report we will use an arbitrary input signal and first estimate the step response before applying the methods.
2.5.2 Moment methods
In moment methods [˚AH95] we interpret the normalized impulse response
f (t) = R∞h(t)
0 h(t)dt
of the system as a probability density function. The quantity h(t) is the impulse response. Then we can call the quantities
mn=
Z ∞
0
tnf (t)dt (2.12)
K =
Z ∞
0
h(t)dt.
The average residence time will be
Tar = m1 (2.13)
and the time constant is solved from
T2= m2− Tar2.
The time delay is as for the corresponding area method:
L = Tar− T.
For the second order model (2.10) K and Tar is calculated as for the first order model.
The time constant is solved from
T2 = 1 2m2− 1 2T 2 ar
and the time-delay given by
L = Tar− 2T.
See [˚AH95] for the derivation of these moment methods. In [˚AH95] these methods require an impulse response experiment to be performed. In this report we will use arbitrary input signals and then estimate the impulse response before employing the moment methods above.
In [˚AH95] they also describe a method, by using moments of the input u(t) and output signals y(t), that can be used with any signal that decays enough quickly. They also describe a variant of this method that can give an approximative model when using input and output signals that do not decay but are bounded by eαt for large t, where α is a
2.5.3 An area and moment method with better noise properties
In [Ing03] they notice that in the calculation of moments, like equation (2.12), values of h(t), u(t) or y(t) at late times in the integration will have a higher weight tn. This will give a rapidly decreasing signal to noise ratio (SNR) for late times, since the noise level is about constant for all times but the signal is decaying. This results in inaccurate estimates of the moments. This is probably a major reason for the poor performance of the moment methods in section 5.8, 6.6 and 7.9 in this report.
In [Ing03] a modified area and moment method for the first order system (2.9) is described. By defining the signals yd(t) = dtdy(t) and ud(t) = dtdu(t) the factor t in the integral m1
(2.12) in equation (2.13) disappears giving integrals of only u(t) and y(t) for the calculation of Tar. By computing time constant T by equation (2.11) in the area method for the first
order system also the second order moment integral is avoided. The identification data are collected by a combination of a setpoint change experiment in closed loop and a step response experiment in open loop.
This method could be called a combined area and moment method. It could also be called a pure area method because that is what the result is.
2.6 Maximum likelihood methods
3
Implementation of correlation methods
In this chapter we describe the implementation of the time delay estimation methods used in the simulations in this report. The implementation of the methods has only been used for a sampling interval Ts = 1. No prewhitening of the input signal is performed in any
of the methods below. However, prewhitening can be executed before calling the methods if desired. In section 7.1-7.6 we give specific names to some special cases (special method parameters) of the methods in this chapter.
Since all time-delay methods in this report utilize estimated impulse and step responses, we start with the implementation of estimation these responses in Section 3.1 and 3.2. Then in Section 3.3-3.8 the description of the implementation of the time-delay estimation methods follows.
3.1 Impulse response estimation, estimpresp
The Matlab function estimpresp which estimates impulse response in this report is imple-mented in the following way:
• Estimate a FIR model (the impulse response) of length 71 by the arx command in [SIT] with the order specified as na= 0, nb = 71 and nk= 0 (number of time delays).
• The arx command also delivers an estimate of the standard deviation of the estimated impulse response coefficients.
• No prewhitening of the input signal is performed in this function. This implementation of the impulse response estimation was chosen because:
• Using the arx command instead of the impulse command in [SIT] gives a better understanding of what is done.
• Calculation of too few impulse response coefficients gives estimates of lower quality. The length of 71 coefficients is the same as in the impulse impulse command. • No prewhitening because want to test also without prewhitening.
Se Figures 5.1 and 5.2 for examples of impulse response estimation using this function.
3.2 Step response estimation, eststepresp2
1. Estimate the impulse response by the function estimpresp as in section 3.1. 2. Utilize the function sim in [SIT] to obtain the step response.
3. Let the estimated standard deviation of the estimated step response coefficients be either:
(a) The standard deviation delivered by the sim function (here called simvar ). The standard deviation increases with time. The reason is probably that the uncertainty accumulates as the impulse response is integrated to obtain the step response. This estimated standard deviation seems to be too large when manually studying some simulations (t132efunc.m), see Figures 5.5 and 5.6. (b) The standard deviation of the corresponding impulse response coefficients (here
called impulsevar ). The standard deviation appears to be independent of time. This estimated standard deviation seems to be too small when manually study-ing some simulations (t132efunc.m), see Figures 5.3 and 5.4.
The theoretically correct way to estimate the standard deviation should be 3a but it also gives an incorrect estimate. The choice between 3a or 3b depends on the intended use of the estimate.
This implementation of the step response estimation was chosen because:
• Using the commands arx and sim instead of the step command in [SIT] gives a better understanding of what is done.
• No prewhitening because we wanted to test also without prewhitening.
3.3 Thresholding of impulse response, idimp4
1. Estimate the impulse response with 71 coefficients by the function estimpresp as in section 3.1.
2. A threshold vector h(t) = hstd· ˆystd(t) is chosen where ˆystd(t) is the vector of
esti-mated standard deviations of the estiesti-mated impulse response coefficients and hstd
is a constant chosen by the user. There is thus a different threshold for each time point. The default value for hstd is 3. This default choice of threshold was reached
after manual inspection of about 100 simulation (t134dfunc.m) with systems G1 and
G6 (Section 4.2) at SNR=1 and SNR=100. Letting the threshold be as a constant
times the uncertainty was selected after manual inspection of estimated coefficients and estimated standard deviation in some simulations. This can also be theoreti-cally justified. If the estimated impulse response coefficients are Gaussian distributed N (y(t), ystd(t)), then ˆy(t)± 3ystd(t) covers the true impulse response coefficient y(t)
by a probability of 99.7%. Usually one says that the confidence interval has an confidence level of 99.7%. If this (confidence) interval contains the number zero it is possible that the time delay has not come to its end. The impulse response will be zero during the time delay. On the other hand, if the number zero is located outside the interval, the probability is large (99.7%) that the impulse response has started to rise from zero. In this way the risk for false alarm, i.e. the risk to say that the time delay has ended before is actually has, is low (0.3%). Because we choose a so low false alarm risk there will probably be some delay before we discover the start of the impulse response and hence the time-delay will be overestimated. We noted in Section 2.2 that if the noise is Gaussian then the estimated impulse response coefficients will also be Gaussian.
3. Return the first time when absolute value of the impulse response is higher than the threshold h(t). If the impulse response nowhere is higher than the threshold then return NaN (not a number), indicating a missed detection.
3.4 Thresholding of step response, idstep4
1. Estimate the step response with 71 coefficients by the function eststepresp2 as in section 3.2. Select the standard deviation of the corresponding impulse response coefficients as standard deviation of the step response estimate (choice 3b). This was selected since it appeared more true than choice 3a when examining some simulation graphs (t132efunc.m). See Figures 5.3-5.6 for some examples of such graphs.
2. A threshold vector h(t) = hstd·ˆystd(t) is chosen where ˆystd(t) is the vector of estimated
standard deviations of the estimated step response coefficients and hstd is a constant
chosen by the user. There is thus a different threshold for each time point. The default value for hstd is 5. This decision of the threshold was arrived at after about
125 simulations (t136dfunc.m) with G1, G2, G5 and G6 (Section 4.2) at SNR=1 and
SNR = 100. The choice of the threshold as a constant times the uncertainty can be motivated in the same way as for thresholding the impulse response in Section 3.3. If the estimated step response coefficients are Gaussian distributed N (y(t), ystd(t)),
then ˆy(t)±5ystd(t) covers the true impulse response coefficient y(t) by a probability of
99.99994%. We noted in Section 2.2 that if the noise is Gaussian then the estimated step response coefficients will also be Gaussian.
3. Return the first time when absolute value of the step response is higher than the threshold h(t). If the step response nowhere is higher than the threshold then return NaN (not a number).
This method is similar to idimp4. The large difference is that the step response is used instead of the impulse response. A complication with the step response is the choice of how the uncertainty of the response estimate is to be estimated.
It would probably be better to use the standard deviation ˆystd(t) delivered by the sim
function (choice 3a in section 3.2) and letting h(t) = hstd(t)· ˆystd(t), because the estimated
standard deviation delivered by this function increases with time. This seems correct when looking at estimated step responses (see Figures 5.3-5.6).Then a lower relative threshold would be required, which would correspond to a lower (and more correct) confidence level than above. The relative threshold would hopefully then be independent of the time.
3.5 CUSUM detector on impulse response, idimpCusum3
1. Estimate the impulse response with 71 coefficients by the function estimpresp as in section 3.1.
2. A single threshold h = hstd· ˆystd(0) is chosen where ˆystd(0) is the estimated standard
deviation of the estimated impulse response coefficients for the time point t = 0 and hstd is a constant chosen by the user. The same threshold is accordingly used for
all time points. Since the uncertainty of the impulse response coefficients seems to be about the same for all time points it should be enough with a single threshold. Default value for hstd is 2. This default choice of threshold was reached after some
hundred simulations (t130dfunc.m) of the systems G1, G2, G5 and G6(Section 4.2) at
SNR=1, looking at the estimated impulse response with estimated uncertainty and at the ”test statistics” g(t) in the CUSUM test (Section 2.3.2). A different method for impulse response estimation than in section 3.1 was employed when choosing the default threshold.
3. A drift ν = νstd· ˆystd(0) is chosen, where νstd is a constant chosen by the user. The
same drift is used for all time points. Default value for νstd is 1. This default drift
was selected in the same way as the threshold above.
4. Return the first time when the test statistic g(t) is higher than the threshold h. If the test statistic nowhere is higher than the threshold then return NaN (not a number), indicating a missed detection.
3.6 CUSUM detector on step response, idstepCusum3
1. Estimate the step response with 71 coefficients by the function eststepresp2 as in section 3.2. Select the standard deviation of the corresponding impulse response coefficients as standard deviation of the step response estimate (choice 3b). This was selected since it appeared more true than choice 3a when examining some simulation graphs (t132efunc.m). See Figures 5.3-5.6 for some examples. See also the discussion below.
2. A single threshold h = hstd· ˆystd(0) is chosen where ˆystd(0) is the estimated standard
deviation of the estimated step response coefficients for the time point t = 0 and hstd
is a constant chosen by the user. The same threshold is accordingly used for all time points. Since the used estimate of the uncertainty of the step response coefficients seems to be about the same for all time points it should be enough with a single threshold. Default value for hstd is 0.5. This default choice of threshold was reached
after some hundred simulations of the systems G1, G2, G5 and G6 (Section 4.2)
at SNR=1, looking at estimated step response with estimated uncertainty and at the ”test statistics” g(t) in the CUSUM test (Section 2.3.2). A different method for step response estimation than in section 3.2 was used when choosing the default threshold.
3. A drift ν = νstd· ˆystd(0) is chosen, where νstd is a constant chosen by the user. The
same drift is used for all time points. Default value for νstd is 4. This default drift
was selected in the same way as the threshold above.
4. Return the first time when the test statistic g(t) is higher than the threshold h. If the test statistic nowhere is higher than the threshold then return NaN (not a number). It would perhaps be better to use the standard deviation delivered by the sim function (choice 3a in section 3.2) and letting h = hstd· ˆystd, where ˆystd is the vector of estimated
standard deviations, because sometimes the estimated step response is very wrong for late times. The estimated standard deviation delivered by the sim function increases with time. It also would be more theoretically correct.
3.7 Area method, area1
This method is described in [˚AH95, p. 13-14, 24-27] and uses either the first order model (here called 1stord )
G(s) = K
1− sTe
or the second order model (here called 2ndord )
G(s) = K
(1− sT )2e
sL. (3.2)
The steps of the method are:
1. Estimate the step response ˆy(t) with 101 coefficients by the function eststepresp2 as in section 3.2. This length of the step response was selected after studying some plots of simulated step responses for systems G1-G6at SNR=100 and the input signal type
RBS 0-100% (Section 4.2) (t142dfunc.m). This length should be large enough for the step response to reach steady state.
2. Estimate the steady state gain of the system, ˆK, by taking the average of the 21 last step response coefficients. This estimate was selected in the same way as the length of the step response.
3. Get an estimate ˆA0 by numerically integrate ˆK− ˆy(t) from t = 0 to its end. This
integration can be conducted in two ways. Either by
(a) Using the trapezoidal method (here called trapez1 ). If the start or end points are not equal to one of the time points where the step response estimate is given, the integral from the start point to the next given time point and the integral from the previous given time point to the end point are also approximated by the trapezoidal method. The value of the step response at these start and end points are approximated by linear interpolation.
(b) A cumulative sum (here called sum1 ) of the step response coefficients. If the start or end point is not equal to one of the time points where the step response estimate is given, the integral is calculated from and to the nearest time point. 4. Calculate the average residence time: ˆTar= ˆA0/ ˆK.
5. Get an estimate ˆA1 by numerically integrate the estimated step response ˆy(t) from
t = 0 to t = ˆTar by the trapezoidal method.
6. The time constant ˆT is computed [˚AH95, p. 24]. 7. The time-delay estimate ˆL is computed [˚AH95, p. 25].
8. Return the time-delay estimate ˆTd = ˆL + 1. The one (+1) is added because this
gave the best estimates by manual inspection of some simulations. It can also be motivated by the fact that by sampling of a continuous-time system there will be an extra delay of one sampling period.
3.8 Moment method, moment1
This method is described in [˚AH95, p. 27-30] and also uses the first and second order models in equation (3.1) and (3.2).
1. Estimate the impulse response ˆh(t) with 101 coefficients by the function estimpresp as in section 3.1. This length of the step response was selected after looking at some estimated impulse responses.
2. Obtain approximations of the zeroth, first and second order moments of the impulse response h(t):
m0 =R0∞h(t)dt, m1 =R0∞t· h(t)dt, m2=R0∞t2· h(t)dt
The estimates ˆm0, ˆm1 and ˆm2 are acquired by numerical integration by the
trape-zoidal method (point 3a in section 3.7) of ˆh(t), t· ˆh(t) and t2· ˆh(t) from zero to the
end of the estimated impulse response.
3. The gain of the system is estimated by ˆK = ˆG(0) = ˆm0.
4. The average residence time is estimated as ˆTar = ˆm1/ ˆm0.
5. The (square of) the time constant is estimated by ˆT2 = ˆm
2/ ˆm0− ˆTar2.
6. The time-delay of the continuous-time system is computed as ˆL = ˆTar− ˆT .
7. As time-delay estimate for the sampled system we choose ˆTd = ˆL + 0.5. The 0.5 is
added because this gave the best estimates by manual inspection of some simulations and estimates. It could be argued that instead 1.0 should be added because of the fact that by sampling of a continuous-time system there will be an extra delay of one sampling period. However, there might be a bias in the estimate, perhaps due to the numerical integration, that makes 0.5 better.
4
Simulation setup
Simulations have been conducted in Matlab in order to experimentally study time delay estimation methods in this report. This section describes the setup of simulations. The output signal y(t) was simulated as
y(t) = G(s)u(t) + v(t) (4.1)
where y(t), u(t) and v(t) are the output, input and noise signals respectively. G(s) is a linear system with dead-time.
The factors of the simulations can be divided in fixed factors and varied factors. The different possible choices for the same factor are in this report called levels of the factor. The choice of level for all factors is called a factor level combination.
4.1 Fixed factors.
Fixed factors for the open loop simulations were
• The noise v(t) was white and Gaussian.
• The system G(s) was simulated by the function lsim in [CST] in continuous-time. • The sampling interval for the dead-time estimation was Ts = 1.
4.2 Varied factors.
The varied factors for the open loop simulations are given the names Method, Prewhite, threshold, drift, Sys, InType and SNR.
Method. Several dead-time estimation methods were tested:
• idimp4. See Section 3. • idstep4. See Section 3. • idimpCusum3. See Section 3. • idstepCusum3. See Section 3. • area1. See Section 3.
• moment1. See Section 3.
Prewhite. For each of the methods above the input and output signals u(t) and y(t) were either
1. Remove the mean value of u(t) and y(t). 2. Estimate a 10th order AR model for u(t). 3. Filter u(t) and y(t) through the AR model. • nopw. Not filtered through a prewhitening filter.
threshold. Relative threshold in thresholding methods. See Section 3. drift. Relative drift in CUSUM methods. See Section 3.
modelStruc. Model structure in area and moment methods. See Section 3.7. Possible choices are 1stord and 2ndord.
intMethod. Integration method in area and moment methods. See Section 3.7. Possible choices are trapez1 and sum1.
Sys. The following systems were simulated (see Figures 4.1 and 4.2 for impulse and step responses):
• A slow second order system (slow2 ): G1(s) = e−9s
1
(10s + 1)(s + 1) (4.2)
• A fast second order system (fast2 ): G2(s) = e−9s
1
(s + 1)(0.1s + 1) (4.3)
• A fourth order system with real poles (real4 ). The poles and zeros lie between the poles of system G1:
G5(s) = e−9s
0.05(s + 0.9)(s + 0.4)
(s + 1)(s + 0.6)(s + 0.3)(s + 0.1) (4.4)
• A fourth order system with complex poles (cplx4 ). The poles have the same distance to the origin as in system G1:
G6(s) = e−9s
1/36(s + 0.9)(s + 0.4)
(s− 0.1 · 2−1/2(−1 ± i))(s − 2−1/2(−1 ± i)) (4.5)
As can be seen in Figures 4.1 and 4.2 this system is very slow and the oscillations are weak.
• A slow first order system:
G7(s) = e−9s 1
(10s + 1) (4.6)
• A fast first order system:
G8(s) = e−9s
1
Time (sec.) Amplitude Impulse response G1 (t130g) 0 14 28 42 56 70 0 0.05 0.1 From: U(1) To: Y(1) Time (sec.) Amplitude Impulse response G2 (t130g) 0 14 28 42 56 70 0 0.5 1 From: U(1) To: Y(1) Time (sec.) Amplitude Impulse response G5 (t130g) 0 14 28 42 56 70 0 0.05 0.1 From: U(1) To: Y(1) Time (sec.) Amplitude Impulse response G6 (t130g) 0 14 28 42 56 70 −0.05 0 0.05 From: U(1) To: Y(1) Time (sec.) Amplitude Impulse response G7 (t130g) 0 14 28 42 56 70 0 0.05 0.1 From: U(1) To: Y(1) Time (sec.) Amplitude Impulse response G8 (t130g) 0 14 28 42 56 70 0 0.5 1 From: U(1) To: Y(1)
Time (sec.) Amplitude Step response G1 (t132g) 0 14 28 42 56 70 0 0.5 1 From: U(1) To: Y(1) Time (sec.) Amplitude Step response G2 (t132g) 0 14 28 42 56 70 0 0.5 1 From: U(1) To: Y(1) Time (sec.) Amplitude Step response G5 (t132g) 0 14 28 42 56 70 0 0.5 1 From: U(1) To: Y(1) Time (sec.) Amplitude Step response G6 (t132g) 0 14 28 42 56 70 0 0.5 1 1.5 From: U(1) To: Y(1) Time (sec.) Amplitude Step response G7 (t132g) 0 14 28 42 56 70 0 0.5 1 From: U(1) To: Y(1) Time (sec.) Amplitude Step response G8 (t132g) 0 14 28 42 56 70 0 0.5 1 From: U(1) To: Y(1)
Note that for all the systems above, the time delay will be 10 after the sampling. InType. The input signal was 500 samples long and could be of three different types:
• RBS 10-30%. (Random Binary Signal) with frequency contents between 10%-30% of the Nyquist frequency. It was generated by the function idinput in [SIT]. The same realization of the signal was used in all trials of the same Monte Carlo simulation. See Figure 4.3 for an example of this input signal type.
• RBS 0-100%. RBS with frequency contents between 0%-100% of the Nyquist fre-quency, i.e white noise. It was generated by the function idinput in [SIT]. The same realization of the signal was used in all trials of the same Monte Carlo simulation.See Figure 4.4 for an example of this input signal type.
• Steps. Three steps of the form (in Matlab code): [zeros(50,1);ones(150,1); -ones(150,1); zeros(150,1)]. See Figure 4.5 for an example of this input signal type. 0 50 100 150 200 250 300 350 400 450 500 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 samples Time signal RBS10−30% 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 10−2 10−1 100 101 Normalized frequency Power Power spectruml RBS10−30%
Figure 4.3: Time signal (left) and frequency spectrum (right) for a realization of the input signal type RBS 10-30%.
SNR. The signal-to-noise ratio (SNR or input SNR) could be chosen. The SNR was defined as the ratio of the input signal power and the output noise power. The used SNR has been either 1 or 100. See Figure 4.6 for SNR defined as the ratio of the output signal (without noise) power to the output noise power (output SNR) for different systems and different input signal types when the input SNR was fixed. As can be seen the two definitions of SNR differ at most with a factor of two.
Many of the factors and levels of the open loop simulations were chosen to be the same as in [Hor00, Bj¨o02], which makes comparisons easy. The names of the factors and factor levels in this section are also used in the following plots.
0 50 100 150 200 250 300 350 400 450 500 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 samples Time signal RBS0−100% 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 10−1 100 101 Normalized frequency Power Power spectruml RBS0−100%
Figure 4.4: Time signal (left) and frequency spectrum (right) for a realization of the input signal type RBS 0-100%. 0 50 100 150 200 250 300 350 400 450 500 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 samples Time signal steps
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 10−4 10−3 10−2 10−1 100 101 102 Normalized frequency Power
Power spectruml steps
Figure 4.5: Time signal (left) and frequency spectrum (right) for a realization of the input signal type Steps.
slow2 fast2 real4 cplx4 G_7 G_8 RBS10−30% RBS0−100 % steps 0 1 2 inType 1.91 1.04 1.99 1.04 1.41 1.89 1.64 1.04 1.92
t158a: mean, data(:,:)
1.02 system 1.04 1.89 1.02. MIN 1.05 1.99. MAX 1.04 1.46 1.66 slow2 fast2 real4 cplx4 G_7 G_8 RBS10−30% RBS0−100 % steps 0 0.1 0.2 inType 0.1 0.0209 0.119 0.015 0.0685 0.114 0.0825 0.0198 0.103 t158a: std, data(:,:) 0.0138. MIN system 0.0181 0.0949 0.0152 0.0211 0.122. MAX 0.0183 0.0687 0.0871
Figure 4.6: Estimated SNR defined as the ratio of the output signal (without noise) power to the output noise power for different systems and different input signal types. Mean value (left) and standard deviation (right). The SNR defined as the ratio of input signal power and output noise power was 1. (t158a.m)
5
Simple experiments
In this chapter we perform some simple experiments with different correlations methods for dead-time estimation. We have manually step-by-step inspected the calculations of time-delay estimates for different systems when using simulated input signals of different types and different SNRs.
Because all time-delay methods in this report utilize estimated impulse and step responses, we start by studying these in Section 5.1 and 5.2. Then in Section 5.3-5.8 the examinations of the time-delay estimation methods follows.
5.1 Impulse response estimation
In Figure 5.1 and 5.2 estimated impulse responses with the function estimpresp (Section 3.1) for systems G1and G2 (Section 4.2), different input signal types and different SNR are
depicted. 0 10 20 30 40 50 60 70 −0.02 0 0.02 0.04 0.06 0.08 Time [s] Impulse response ,10−30%, SNR=100, G1 0 10 20 30 40 50 60 70 −0.1 −0.05 0 0.05 0.1 0.15 Time [s] Impulse response ,10−30%, SNR=1, G1 0 10 20 30 40 50 60 70 −0.02 0 0.02 0.04 0.06 0.08 Time [s] Impulse response ,0−100%, SNR=100, G1 0 10 20 30 40 50 60 70 −0.05 0 0.05 0.1 0.15 Time [s] Impulse response ,0−100%, SNR=1, G1 0 10 20 30 40 50 60 70 −0.2 −0.1 0 0.1 0.2 0.3 Time [s]
Impulse response ,steps, SNR=100, G1
0 10 20 30 40 50 60 70 −2 −1 0 1 2 Time [s]
Impulse response ,steps, SNR=1, G1
Figure 5.1: Impulse response estimate of the system G1 (Section 4.2) by the function idimp4
(Section 3.3) for different input signal types (Section 4) and different SNRs. The solid line is the true impulse response. The circles are the estimated impulse response and the triangles mark ±two estimated standard deviations. Note the different ranges of the vertical axes. (t130f1.m)
0 10 20 30 40 50 60 70 −0.2 0 0.2 0.4 0.6 0.8 Time [s] Impulse response ,10−30%, SNR=100, G2 0 10 20 30 40 50 60 70 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 Time [s] Impulse response ,10−30%, SNR=1, G2 0 10 20 30 40 50 60 70 −0.2 0 0.2 0.4 0.6 0.8 Time [s] Impulse response ,0−100%, SNR=100, G2 0 10 20 30 40 50 60 70 −0.2 0 0.2 0.4 0.6 0.8 Time [s] Impulse response ,0−100%, SNR=1, G2 0 10 20 30 40 50 60 70 −0.2 0 0.2 0.4 0.6 0.8 Time [s]
Impulse response ,steps, SNR=100, G2
0 10 20 30 40 50 60 70 −2 −1 0 1 2 Time [s]
Impulse response ,steps, SNR=1, G2
Figure 5.2: Impulse response estimate of the system G2 (Section 4.2) by the function idimp4
(Section 3.3) for different input signal types (Section 4) and different SNRs. The solid line is the true impulse response. The circles are the estimated impulse response and the triangles mark ±two estimated standard deviations. Note the different ranges of the vertical axes. (t130f1.m)
We have manually inspected some estimated impulse responses (more than shown here) for different systems when using simulated input signals of different types and different SNRs. We make some observations from our simple experiments:
• Estimating more impulse response coefficients make the uncertainty in each coeffi-cient lower.
• The sampling misses the peak of the impulse response of the fast system.
• The estimated impulse response seems to be 1/2 sampling interval behind the true continuous-time impulse response (G1-G2, G5-G8). For fast rising impulse responses
the first high coefficient of the estimated impulse response is lower than the true one. This is because the true impulse response is sampled too seldom and therefore we miss the peak. This should give large errors in moment methods (Section 3.8, 5.8, 6.6).
• The input signal type Steps gives a larger uncertainty in the impulse response esti-mate than RBS 10-30% and RBS 0-100%.
• The uncertainty in the impulse response estimate is high for all input signal types at low SNR and for Steps also at high SNR.
5.2 Step response estimation
In Figures 5.3 and 5.4 estimated step responses with the function eststepresp2 using impul-sevar (Section 3.2) for systems G1and G2, different input signal types and different SNR
are displayed.
In Figures 5.5 and 5.6 estimated step responses with the function eststepresp2 using simvar (Section 3.2) for systems G1and G2, different input signal types and different SNR are
depicted.
We have manually inspected some estimated step responses (more than shown here) for different systems when using simulated input signals of different types and different SNRs. We make some observations from our simple experiments:
• For step responses, the estimated response is not behind the true one.
• The method impulsevar gives a too low and simvar a too high estimate of the uncertainty of the estimated step response.
5.3 Thresholding impulse response
We have manually inspected some graphs with estimated impulse response and threshold and the corresponding estimated time-delay for different systems when using simulated in-put signals of different types and different SNRs. See Figures 5.7 and 5.8 for two examples. We make some observations from our simple experiments:
0 10 20 30 40 50 60 70 −0.2 0 0.2 0.4 0.6 0.8 1 1.2 Time [s] Step response ,10−30%, SNR=100, G1 0 10 20 30 40 50 60 70 −0.2 0 0.2 0.4 0.6 0.8 1 1.2 Time [s] Step response ,10−30%, SNR=1, G1 0 10 20 30 40 50 60 70 0 0.2 0.4 0.6 0.8 1 Time [s] Step response ,0−100%, SNR=100, G1 0 10 20 30 40 50 60 70 −0.2 0 0.2 0.4 0.6 0.8 1 1.2 Time [s] Step response ,0−100%, SNR=1, G1 0 10 20 30 40 50 60 70 −0.2 0 0.2 0.4 0.6 0.8 1 1.2 Time [s]
Step response ,steps, SNR=100, G1
0 10 20 30 40 50 60 70 −2 −1 0 1 2 3 Time [s] Step response ,steps, SNR=1, G1
Figure 5.3: Step response estimate of the system G1 (Section 4.2) by the function idstep4
(Section 3.4) for different input signal types (Section 4) and different SNRs. The estimated standard deviations of the underlying impulse response estimate are taken as the estimated standard deviations of the step response estimate (choice 3b in Section 3.2). The solid line is the true step response. The circles are the estimated step response and the triangles mark ±two estimated standard deviations. Note the different ranges of the vertical axes. (t132f1.m)
0 10 20 30 40 50 60 70 −0.2 0 0.2 0.4 0.6 0.8 1 1.2 Time [s] Step response ,10−30%, SNR=100, G2 0 10 20 30 40 50 60 70 −0.5 0 0.5 1 1.5 2 Time [s] Step response ,10−30%, SNR=1, G2 0 10 20 30 40 50 60 70 0 0.2 0.4 0.6 0.8 1 1.2 1.4 Time [s] Step response ,0−100%, SNR=100, G2 0 10 20 30 40 50 60 70 −0.2 0 0.2 0.4 0.6 0.8 1 1.2 Time [s] Step response ,0−100%, SNR=1, G2 0 10 20 30 40 50 60 70 −0.2 0 0.2 0.4 0.6 0.8 1 1.2 Time [s]
Step response ,steps, SNR=100, G2
0 10 20 30 40 50 60 70 −2 −1 0 1 2 3 Time [s] Step response ,steps, SNR=1, G2
Figure 5.4: Step response estimate of the system G2 (Section 4.2) by the function idstep4
(Section 3.4) for different input signal types (Section 4) and different SNRs. The estimated standard deviations of the underlying impulse response estimate are taken as the estimated standard deviations of the step response estimate (choice 3b in Section 3.2). The solid line is the true step response. The circles are the estimated step response and the triangles mark ±two estimated standard deviations. Note the different ranges of the vertical axes. (t132f1.m)
0 10 20 30 40 50 60 70 −0.2 0 0.2 0.4 0.6 0.8 1 1.2 Time [s] Step response ,10−30%, SNR=100, G1 0 10 20 30 40 50 60 70 −0.5 0 0.5 1 1.5 Time [s] Step response ,10−30%, SNR=1, G1 0 10 20 30 40 50 60 70 −0.2 0 0.2 0.4 0.6 0.8 1 1.2 Time [s] Step response ,0−100%, SNR=100, G1 0 10 20 30 40 50 60 70 −0.2 0 0.2 0.4 0.6 0.8 1 1.2 Time [s] Step response ,0−100%, SNR=1, G1 0 10 20 30 40 50 60 70 −0.2 0 0.2 0.4 0.6 0.8 1 1.2 Time [s]
Step response ,steps, SNR=100, G1
0 10 20 30 40 50 60 70 −2 −1 0 1 2 3 Time [s] Step response ,steps, SNR=1, G1
Figure 5.5: Step response estimate of the system G1 (Section 4.2) by the function idstep4
(Section 3.4) for different input signal types (Section 4) and different SNRs. The standard deviations delivered from the sim are taken as the estimated standard deviations of the step response estimate (choice 3a in Section 3.2). The solid line is the true step response. The circles are the estimated step response and the triangles mark ±two estimated standard deviations. Note the different ranges of the vertical axes. (t132f1.m)
0 10 20 30 40 50 60 70 −0.5 0 0.5 1 1.5 Time [s] Step response ,10−30%, SNR=100, G2 0 10 20 30 40 50 60 70 −0.5 0 0.5 1 1.5 2 2.5 Time [s] Step response ,10−30%, SNR=1, G2 0 10 20 30 40 50 60 70 −0.2 0 0.2 0.4 0.6 0.8 1 1.2 Time [s] Step response ,0−100%, SNR=100, G2 0 10 20 30 40 50 60 70 −0.5 0 0.5 1 1.5 Time [s] Step response ,0−100%, SNR=1, G2 0 10 20 30 40 50 60 70 −0.2 0 0.2 0.4 0.6 0.8 1 1.2 Time [s]
Step response ,steps, SNR=100, G2
0 10 20 30 40 50 60 70 −1 −0.5 0 0.5 1 1.5 2 2.5 Time [s] Step response ,steps, SNR=1, G2
Figure 5.6: Step response estimate of the system G2 (Section 4.2) by the function idstep4
(Section 3.4) for different input signal types (Section 4) and different SNRs. The standard deviations delivered from the sim are taken as the estimated standard deviations of the step response estimate (choice 3a in Section 3.2). The solid line is the true step response. The circles are the estimated step response and the triangles mark ±two estimated standard deviations. Note the different ranges of the vertical axes. (t132f1.m)
0 10 20 30 40 50 60 70 −0.08 −0.06 −0.04 −0.02 0 0.02 0.04 0.06 0.08 0.1 0.12 Time [s] Impulse response 0 10 20 30 40 50 60 70 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 Threshold test Time [samples]
Figure 5.7: Left: Estimated impulse response with uncertainty. Right:Estimated impulse re-sponse and threshold. Simulated input signal of type RBS 10-30%. SNR was 1. System G1.
0 10 20 30 40 50 60 70 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 Time [s] Impulse response 0 10 20 30 40 50 60 70 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 Threshold test Time [samples]
Figure 5.8: Left: Estimated impulse response with uncertainty. Right:Estimated impulse re-sponse and threshold. Simulated input signal of type RBS 10-30%. SNR was 1. System G2.
The estimated time delay with idimp4 (hstd= 3) was ˆTd= 10. (t134d1.m)
5.4 Thresholding step response
We have manually inspected some graphs with estimated step response and threshold and the corresponding estimated time-delay for different systems when using simulated input signals of different types and different SNRs. See Figures 5.9 and 5.10 for two examples. We make some observations from our simple experiments:
• Sometimes the method misses to detect. The threshold is never crossed.
0 10 20 30 40 50 60 70 −0.2 0 0.2 0.4 0.6 0.8 1 1.2 Time [s] Step response 0 10 20 30 40 50 60 70 0 0.05 0.1 0.15 Threshold test Time [samples]
Figure 5.9: Left: Estimated step response with uncertainty. Right: Estimated step response and threshold. Simulated input signal of type RBS 10-30%. SNR was 1. System G1. The
estimated time delay with idstep4 (hstd = 5) was ˆTd= 12. (t136d1.m)
5.5 CUSUM detector on impulse response
We have manually inspected some graphs with test statistics g(t), threshold h and drift ν and the corresponding estimated time-delays for different systems when using simulated
0 10 20 30 40 50 60 70 −0.2 0 0.2 0.4 0.6 0.8 1 1.2 1.4 Time [s] Step response 0 10 20 30 40 50 60 70 0 0.1 0.2 0.3 0.4 0.5 0.6 Threshold test Time [samples]
Figure 5.10: Left: Estimated step response with uncertainty. Right: Estimated step response and threshold. Simulated input signal of type RBS 10-30%. SNR was 1. System G2. The
estimated time delay with idstep4 (hstd= 5) was ˆTd= 10. (t136d1.m)
input signals of different types and different SNRs. See Figures 5.11 and 5.12 for two examples. We make some observations from our simple experiments:
• Often the method misses to detect. The threshold is never crossed.
0 10 20 30 40 50 60 70 −0.06 −0.04 −0.02 0 0.02 0.04 0.06 0.08 0.1 0.12 Time [s] Impulse response G1 0 10 20 30 40 50 60 70 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07
CUSUM test statistics G1
Time [samples]
Test statistics threshold drift
Figure 5.11: Left: Impulse response with uncertainty. Right: Test statistics g(t), threshold h and drift ν for CUSUM on impulse response. Simulated input signal of type RBS 10-30%. SNR was 1. System G1. The estimated time delay with idimpCusum3 (hstd= 2 and νstd= 1)
was ˆTd= 11. (t146b1.m)
5.6 CUSUM detector on step response
We have manually inspected some graphs with test statistics g(t), threshold h and drift ν and the corresponding estimated time-delays for different systems when using simulated input signals of different types and different SNRs. See Figures 5.13 and 5.14 for two examples. We make some observations from our simple experiments:
0 10 20 30 40 50 60 70 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 Time [s] Impulse response G2 0 10 20 30 40 50 60 70 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16
CUSUM test statistics G2
Time [samples]
Test statistics threshold drift
Figure 5.12: Left: Impulse response with uncertainty. Right: Test statistics g(t), threshold h and drift ν for CUSUM on impulse response. Simulated input signal of type RBS 10-30%. SNR was 1. System G2. The estimated time delay with idimpCusum3 (hstd= 2 and νstd= 1)
was ˆTd= 7. (t146b1.m)
• Often the method misses to detect. The threshold is never crossed.
• The drift must be very high to counteract the cumulative summation which is per-formed to get the step response from the impulse response. Are CUSUM on impulse and step response equivalent?
0 10 20 30 40 50 60 70 −0.2 0 0.2 0.4 0.6 0.8 1 1.2 Time [s] Step response G1 0 10 20 30 40 50 60 70 0 0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04 0.045
CUSUM test statistics G1
Time [samples]
Test statistics threshold drift
Figure 5.13: Left: Step response with uncertainty. Right: Test statistics g(t), threshold h and drift ν for CUSUM on step response. Simulated input signal of type RBS 10-30%. SNR was 1. System G1. The estimated time delay with idimpCusum3 (hstd = 0.5 and νstd = 4) was
ˆ
Td= 11. (t146b2.m)
5.7 Area methods
We have manually step-by-step inspected some calculations of time-delay estimates with area methods for different systems when using simulated input signals of different types and different SNRs. We make some observations from our simple experiments:
0 10 20 30 40 50 60 70 −0.2 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 Time [s] Step response G2 0 10 20 30 40 50 60 70 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7
CUSUM test statistics G2
Time [samples]
Test statistics threshold drift
Figure 5.14: Left: Step response with uncertainty. Right: Test statistics g(t), threshold h and drift ν for CUSUM on step response. Simulated input signal of type RBS 10-30%. SNR was 1. System G2. The estimated time delay with idimpCusum3 (hstd = 0.5 and νstd = 4) was
ˆ
Td= 10. (t146b2.m)
• With the area method with the first order model (1stord) it is difficult to estimate the end value of the step because the uncertainty of the estimated step grows with time.
• Faster system are harder because the uncertainty in the step response is larger and the time span to integrate over is smaller. Despite simulated models of the correct model structure (G7 and G8) and a high SNR this method performs poorly. The
area method with the first order model (1stord ) and trapez1 quadrature: Slow step response is easier than fast. This is in contrast to many other methods, e.g. thresh-olding on impulse response. It is intuitive that fast system is harder because the A1
integral will be done over a very short time span where the step response is nonzero. Dead-time estimation of a slow system seems to work rather well even for low SNR and despite bad step response estimates with a high uncertainty.
• The area method with the second order model (2ndord) and trapez1 quadrature: Not as good as with the first order system for first order system (expected). Appears to be worse than the area method with the first order model for other systems also. Seems to give lower estimate than the area method with the first order model. • The uncertainty of the step response estimate in section 3.2 often increases with the
time (see Figures 5.5 and 5.6). Since the gain K, which is used in the area methods, of the system is estimated by the last step response coefficients this estimate often will be completely incorrect. The objective of using area and moment methods was to by integration avoid the sensitivity of a single or a few values of the step and impulse responses.
• For the area methods there are several things that can go wrong:
– The estimated gain ˆK can get negative if the estimated step response falls below zero at the end of the estimated step response.
– The estimated area ˆA0 can become negative if too much of the estimated step
response is above the estimated gain ˆK. This can happen if the estimated gain is too low.
– The time Tar can be estimated to a too high value. This can happen if the gain
K is estimated to a too low value.
– The time-delay estimate can be calculated to a negative value. This happens if either or both ˆTar is too low and ˆT is too large. The time constant ˆT may
become too large if the estimated gain ˆK is too low.
• The area methods seems to succeed (non-NaN as result) more often than the moment methods.
5.8 Moment methods
We have manually step-by-step inspected some calculations of time-delay estimates with moment methods for different systems when using simulated input signals of different types and different SNRs. We make some observations from our simple experiments:
• The area and moment methods with the second order model (2ndord) seem to be sensitive to the model structure being correct. These methods with the first order model (1stord ) works much better for G7, which is of the correct model structure
for 1stord.
• The area and moment methods with the first order model (1stord) appears to be more robust than these methods with the second order model (2ndord ).
• Moment methods seems to be more unreliable than area methods. Not even with a true simulated impulse response they give a good answer. The methods works worse for fast systems. Maybe an explanation of the problems with moment methods is that the integral around the knee will have a large error since the impulse response takes a large jump there. The area of the impulse response is smaller than the area of a step response. Therefore an error at the knee will have a large influence. In all our systems the time delay is a multiple of the sampling interval. Therefore the problem with the knee should be small.
• The moment methods (both estimated and true impulse response) give T2 < 0 for G 6
when SNR = 10000. Is it always so for oscillating impulse responses? The moment methods often give T2 < 0 for G2 when SNR = 10000 . Sometimes it is the same
for G8.
• The moment method with the second order model (2ndord) usually give a lower value than with the first order model (1stord ).
• For moment methods it is important how to integrate at the knee. The trapezoidal method tries to smooth the peak. Maybe a sum integration is better.
6
Statistical analysis of simulations
In this section we will study the time delay estimation methods with the help of statistical analysis of Monte Carlo simulations. We scrutinize RMS errors and percentages of failure for different cases.
6.1 Results with threshold on impulse response
In Figure 6.1 (top), the RMS error of the idimp4 (see Section 3.3) dead-time estimate is displayed for several values of the relative threshold hstd (1, 3 and 5). Estimation both
with and without prewhitening (Prewhite, Section 4.2) was used. The SNR was high (SNR ≈ 100) or low (SNR ≈ 1). The three different input signal types (InType) of Section 4.2 (RBS 10-30%, RBS 0-100% and Steps) were used. The four systems (Sys) G1, G2, G5
and G6 in Section 4.2 were employed. For each level combination of the factors threshold,
Prewhite, InType, SNR and Sys the time-delay was estimated from 1024 trials. However, some of these estimates failed. The reason was that, unfortunately, in many trials there was no detection, because the threshold was never crossed. The RMS error of the estimates was computed only for the successful trials. Then, for each level combination of the other factors, the system with the worst, i.e. highest, RMS error was plotted. Which systems that were chosen are not shown in the figure. Since it is important to know the number of failed estimates when interpreting the RMS error values, the percentage of failed estimates (missed detections) is depicted in Figure 6.1 (bottom).
We observe in Figure 6.1 that the prewhitening has no large influence. For the input signal types RBS 10-30% and RBS 0-100% at SNR=100, the thresholds hstd= 3 and 5 are
better than 1. For SNR=1 and for the input signal type Steps at SNR=100, the threshold hstd = 1 must be used. The thresholds hstd= 3 and 5 would give many missed detections
in this case. Unfortunately, there are different choices of method parameters (prewhitening and threshold) which are best for different environment factor level combinations (input signal type and SNR). The most robust threshold is hstd = 1 because it has almost no
missed detections. On the other side, its RMS error for RBS 10-30% and RBS 0-100% for SNR=100 is not so good.
Figure 6.2 displays the RMS estimation error and percentage of failures for the same simulation as in Figure 6.1 but now shows how the RMS error depends on the system. The maximum RMS error over the systems is thus not chosen. The graphs have the levels of the factors threshold, Prewhite and Sys on the axes. The mean value over the two SNRs and the three input signal types is shown.
In Figure 6.2 it can be seen that the time-delay estimation is easier with system G2
(fast second order system) than with the other used systems. It is also evident that the thresholds hstd= 3 and hstd= 5 give rise to many missed detections . The best choice of
the method parameter does not depend in a clear way on the system.
6.2 Results with threshold on step response
In Figure 6.3 (top), the RMS error of the idstep4 (see Section 3.4) dead-time estimate is displayed for several values of the relative threshold hstd (3, 5 and 7). The simulation
100*10−30% 100*0−100% 100*steps 1*10−30% 1*0−100% 1*steps nopw*1 nopw*3 nopw*5 pw*1 pw*3 pw*5 0 50 100 8.31 8.46 31.4 Prewhite*threshold 8.38 20.1 55. MAX 8.26 26.9 12.3 8.39 8.55 18.4 32.6 8.53 28.6 8.37 1.73 13.2
t138b2:030131 09:45 rms,max sys: −, data(:,:,:,:,m,m,m,m,&max,m,m,m,m,:)
8.63 21.2 1.77 1.23. MIN 8.37 25.5 11.6 2.6 8.46 18.1 27.8 8.44 1.71 32.1 SNR*InType 1.91 1.24 2.56 100*10−30% 100*0−100% 100*steps 1*10−30% 1*0−100% 1*steps nopw*1 nopw*3 nopw*5 pw*1 pw*3 pw*5 0 50 100 0 0 71.1 Prewhite*threshold 0 47 100 0 62.3 99.7 0 0 40 100. MAX 0.0977 64.6 0. MIN 0 99.8
t138b3:030131 09:46, %NaNs, max sys: −, data(:,:,:,:,m,m,m,m,&max,m,m,m,m,:)
0.0977 50 100 0 0 0 58.9 99.8 0 0 34.4 100 0 0 99.5 SNR*InType 0.0977 0 0.293
Figure 6.1: RMS error of time-delay estimates (top) and percentage failed estimations (bot-tom) for thresholding of impulse response ( idimp4) as a function of SNR, input signal type and relative threshold. For each combination of the other factors the system with the largest RMS error is used in the plot. The axis label “SNR*InType” means different (factor level) combinations of SNR and input signal type. The tick mark labels tell us what factor level com-binations there are in the different“rows”and“columns”. The level“10-30%”is an abbreviation of “RBS 10-30% and “0-100%” of “RBS 0-100%”. See Section 4.2 for further explanation.
slow2 fast2 real4 cplx4 nopw*1 nopw*3 nopw*5 pw*1 pw*3 pw*5 0 10 20 8.3 8.33 16.6. MAX 8.34 16.1 8.98 Prewhite*threshold 8.33 5.51 15.9 8.44 14.5 6.38
t138b2:030131 09:45 rms,max sys rms,max sys: −, data(:,m,:,m,m,m,m,m,:,m,m,m,m,:)
8.39 16.2 9.8 8.39 15.1 11.8 8.34 5.56 9.63 13.7 0.00784. MIN Sys 7.09 slow2 fast2 real4 cplx4 nopw*1 nopw*3 nopw*5 pw*1 pw*3 pw*5 0 50 100 0 0 36.1 0 35.9 66.6 Prewhite*threshold 0. MIN 10.5 66.5 0.0326 34.4 16.9
t138b3:030131 09:46, %NaNs, max sys, %NaNs: −, data(:,m,:,m,m,m,m,m,:,m,m,m,m,:)
0 34.5 66.5 0 32.5 66.6. MAX 0.0163 10.4 66.4 32.5 16.8 Sys 66.4
Figure 6.2: RMS error of time-delay estimates (left) and percentage failed estimations (right) for thresholding of impulse response ( idimp4) as a function of system and method parameters. Same simulation as in Figure 6.1.
setup and the processing of the estimates were the same as in Section 6.1. The percentage of failed estimates (missed detections) is depicted in Figure 6.3 (bottom).
We observe in Figure 6.3 that the prewhitening has no large influence except for the input signal type Steps at SNR=1, in which case prewhitening (pw ) clearly gives more missed detections. For the input signal types RBS 10-30% and RBS 0-100% at SNR=100, the thresholds hstd = 5 and 7 are better than 3. For SNR=1 and for the input signal type
Steps at SNR=100, the threshold hstd= 3 is the best. However, for the case with the input
signal type Steps at SNR=1 the threshold hstd = 3 should not be used. These results are
similar to those for thresholding on impulse response (Section 6.2). The number of failing estimates is lower when thresholding the step response than the impulse response. The RMS error for low SNR is lower with the step response except for Steps. For high SNR the RMS error for Steps is lower using the step response.
Figure 6.4 displays the RMS estimation error and percentage of failures for the same simulation as in Figure 6.3 but now shows how the RMS error depends on the system. The maximum RMS error over the systems is not chosen. The graphs have the levels of the factors threshold, Prewhite and Sys on the axes. The mean value over the two SNRs and the three input signal types is computed and shown. The best choice of the method parameter does not depend in a clear way on the system. The input signal Steps is the the hardest. It can be seen on the many missed detections.
By comparing Figure 6.1 and Figure 6.3 we notice that the choices of good method param-eters (especially threshold) are similar for thresholding on impulse and step response. For example, the lowest RMS errors are at the same places in the two bar plots. The factor level combinations “100*steps”, “1*10-30%” and “1*0-100%” have lower RMS error with thresholding on step response than on impulse response. The combination “1*steps” has a higher RMS error with step response. Otherwise, the best RMS values are comparable for thresholding on impulse and step response.
