- Email: [email protected]
Inverse filtering based method for estimation of noisy autoregressive signals
Signal Processing 91 (2011) 1659–1664
Contents lists available at ScienceDirect
Signal Processing journal homepage: www.elsevier.com/locate/sigpro
Fast communication
Inverse filtering based method for estimation of noisy autoregressive signals Alimorad Mahmoudi a,n, Mahmood Karimi b a b
Electrical Engineering Department, Shahid Chamran University of Ahvaz, Khouzestan, Iran School of Electrical and Computer Engineering, Shiraz University, Iran
a r t i c l e i n f o
abstract
Article history: Received 7 June 2010 Received in revised form 21 December 2010 Accepted 12 January 2011 Available online 19 January 2011
In this paper we present a new method for estimating the parameters of an autoregressive (AR) signal from observations corrupted with white noise. The leastsquares (LS) estimate of the AR parameters is biased when the observation noise is added to the AR signal. This bias is related to observation noise variance. The proposed method uses inverse filtering technique and Yule–Walker equations for estimating observation noise variance to yield unbiased LS estimate of the AR parameters. The performance of the proposed unbiased algorithm is illustrated by simulation results and they show that the performance of the proposed method is better than the other estimation methods. & 2011 Elsevier B.V. All rights reserved.
Keywords: Autoregressive signals Least-squares method Yule–Walker equations
1. Introduction Modeling of random signals as autoregressive (AR) models are known to be a simple yet effective means for information invoked in many applications, e.g., geophysics, econometrics, speech processing, image processing and communications. The widespread use of AR modeling techniques is due to their excellent resolution performance and existence of a linear set of equations for computing unknown model parameters. It is known that the AR estimator usually shows a high sensitivity to the addition of observation noise. This is an important problem, which limits the utility of these parameter estimation techniques [1]. The noisy AR estimation techniques have been studied in [1–16]. The optimum solution for estimating the AR parameters based on mean squared error (MSE) criteria leads to low-order Yule–Walker equations. When the AR model is corrupted with white noise, due to observation
n
Corresponding author. E-mail addresses: [email protected] (A. Mahmoudi), [email protected] (M. Karimi). 0165-1684/$ - see front matter & 2011 Elsevier B.V. All rights reserved. doi:10.1016/j.sigpro.2011.01.008
noise variance, the zero lag autocorrelation is biased. Due to this, the solution of low-order Yule–Walker equation is biased. The noisy AR estimation techniques can be divided into two main categories. The techniques that belong to the first category avoid zero lag autocorrelations and use high-order Yule–Walker equations. In fact, in this category, the AR(p) signal is modeled by autoregressive moving average ARMA(p,p) model and AR parameters are estimated by ARMA estimation techniques. The modified Yule–Walker (MYW) method [1] is based on this idea. One problem concerning these methods is lack of enough data for computing higher autocorrelation lag estimates, which leads to error in the autocorrelation lags. The second category of common methods used to estimate noisy AR models is based on bias compensation principle. These methods try to remove bias from loworder YW equations. Based on this idea, in [2], Kay uses the Burg estimation method for noisy AR estimation and tries to remove bias from the Burg estimates. In this method, Kay assumed that the observation noise variance is known. This assumption is a limitation for the proposed method. Based on the idea proposed by Kay [2], Zheng [4–8] tried to remove bias from least-squares (LS) estimation of
1660
A. Mahmoudi, M. Karimi / Signal Processing 91 (2011) 1659–1664
noisy AR signals. These methods, which improve the least-squares estimate, are called ILS-based methods. Several ILS-based methods for estimating parameters of noisy AR processes have been proposed by Zheng [4–8]. The main difference between ILS-based methods is that observation noise variance estimation differs in each method. In [5], Zheng proposed an improved leastsquares (ILS) method with direct structure (called ILSD). In [7], he proposed a non-iterative ILS-based method (called ILSNI). In this paper, a new method for estimating AR parameters is proposed. This method is an iterative ILS-based method that uses inverse filtering technique and Yule–Walker equations to estimate the corrupting observation noise variance. In this method, the equations of the linear system are solved and an estimate of observation noise variance is obtained. Then, this estimate is used recursively to correct the bias in the AR parameters estimated by the LS method. Simulation results show that the performance of the proposed method is better than the other estimation methods.
where 2
rx ½0 6 6 rx ½1 Rx ¼ 6 6 ^ 4 rx ½p1
rx ½1
...
rx ½0
...
^ rx ½p2
& ...
rx ½1p
3
7 rx ½2p 7 7 7 ^ 5 rx ½0
ð6Þ
2
3 rx ½1 6 7 6 rx ½2 7 7 rx ¼ 6 6 ^ 7 4 5 rx ½p
ð7Þ
Using (3), (6) and (7), we get Ry ¼ Rx þ s2w Ip
ð8Þ
ry ¼ rx
ð9Þ
where the definitions of Ry and ry are similar to Rx and rx, respectively, and Ip is the p p identity matrix. Combining (5) and (8), and substituting (9), we obtain Ry as2w a ¼ ry
ð10Þ
Multiplying both sides of Eq. (10) by 2. Problem statement
ay ¼ as
ð1Þ
yðtÞ ¼ xðtÞ þ wðtÞ
ð2Þ
where e(t) is the zero hmean stationary driving noise T having variance s2e , a ¼ a1 a2 . . . ap is the vector of parameters of the AR model (T denotes the transpose operation), and w(t) is the zero mean stationary and white observation noise with variance s2w . It is assumed that e(t) and w(t) are mutually uncorrelated, i.e., E{w(t)e(n)}= 0 for all n and t, where E{ } is the expectation operator. The autocorrelation function of y(t) is given by 2 w,
rx ½k þ s
k¼0
rx ½k,
ka0
ð3Þ
rx[k] is the autocorrelation function of noise free AR process x(t). It can be seen from (3) that, due to the observation noise variance, the autocorrelation function of x(t) at zero lag is biased. The main concern of this paper is to estimate the AR coefficients ai (i= 1,2,y,p) using observation data {y(1),y(2),y,y(N)}. Consider the Yule–Walker equations as follows [1]: rx ½k ¼
p X
ay ¼ R1 y ry
kZ1
ð4Þ
i¼1
ð12Þ
ay is the LS estimate of a [4]. Eq. (11) shows that the bias in the LS estimate ay is caused by the observation noise variance s2w . By means of the bias compensation principle presented in [2], the ILS estimate of a may be achieved via a ¼ ay þ s2w R1 y a
ð13Þ
It is clear from Eq. (13) that an unbiased estimate of a is attainable if s2w can be estimated in some way. 3. The estimation algorithm 3.1. Estimation of the observation noise variance In this section, we derive the new method for estimating noisy AR parameters. We use an iterative algorithm for estimating s2w , s2e and a. This algorithm is named as the improved least-squares based on inverse filtering (IFILS). Consider the following vectors: T xt ¼ xðt1Þ xðt2Þ xðtpÞ ð14Þ wt ¼ wðt1Þ
ai rx ½ki,
ð11Þ
where
xðtÞ ¼ a1 xðt1Þ þa2 xðt2Þ þ þ ap xðtpÞ þ eðtÞ
ry ½k ¼ EfyðtÞyðtkÞg ¼
we get
2 1 w Ry a
The noisy AR model under study is represented by
(
R1 y ,
wðt2Þ
yt ¼ xt þ wt ¼ yðt1Þ
wðtpÞ
yðt2Þ
T yðtpÞ
ð15Þ T
ð16Þ
Evaluating (4) for k= 1,y, p gives the following linear system of equations:
If the noisy observation process is filtered by the AR coefficients, then the output z(t) is given by
Rx a ¼ rx
zðtÞ ¼ yðtÞyTt a
ð5Þ
ð17Þ
A. Mahmoudi, M. Karimi / Signal Processing 91 (2011) 1659–1664
Combining (1), (2), (14), and (16), and substituting into (17) gives zðtÞ ¼
wðtÞwTt a þ eðtÞ
ð18Þ
Now, we compute the autocorrelation function of z(t)from (17) and (18). For simplicity, we define a (p+ 1) 1 vector as follows: 1 b¼ ð19Þ a
1661
and 2
3 b b r ½ij i j y 6 7 6 i ¼ 0j ¼ 0 7 6 7 6X 7 p X p 6 7 6 7 b b r ½1 þij i j y 6 7 h ¼ 6 i ¼ 0j ¼ 0 7 6 7 6 7 ^ 6 7 6 p p 7 6 XX 7 4 bi bj ry ½q þij 5 p X p X
ð29Þ
i ¼ 0j ¼ 0
Using (17) and (19), z(t) is expressible as zðtÞ ¼
p X
bi yðtiÞ
ð20Þ
Minimizing the least-squares error of Eq. (26) with respect to h yields the LS estimate as follows:
i¼0
It follows from (20) and (3) that the autocorrelation function of z(t) is equal to rz ½k ¼ EfzðtÞzðtkÞg ¼
p X p X
bi bj ry ½k þ ij
ð21Þ
h ¼ ðHT HÞ1 HT h
ð30Þ
Note that we cannot obtain h from (30) because H and h are unknown. So, we will use a recursive algorithm for estimating h. This algorithm is discussed in Section 3.2.
i¼0j¼0
Using (18) and (19), z(t) is expressible as zðtÞ ¼
wðtÞwTt a þ eðtÞ ¼
p X
3.2. The IFILS algorithm bi wðtiÞ þ eðtÞ
ð22Þ
i¼0
From (22) and (3), the autocorrelation function of z(t), rz[k], can be written as follows: rz ½k ¼ EfzðtÞzðtkÞg ¼ s2w
p X
bi bi þ k þ s2e dðkÞ
ð23Þ
i¼0
where bi + k is taken equal to zero for the values of its index that are outside the interval 0 ri+ krp, and d(k) is the impulse function, which is defined as follows: ( 1, k¼0 dðkÞ ¼ ð24Þ 0, otherwise Using (21) and (23), we obtain
s2w
Xp i¼0
bi bi þ k þ s2e dðkÞ ¼
p X p X
bi bj ry ½k þ ij,
k ¼ 0,1,. . .,q
i¼0j¼0
ð25Þ where q is an arbitrary positive integer less than or equal to p. For illustration convenience, we rewrite Eq. (25) in matrix form as follows: Hh ¼ h
ð26Þ
where " 2# ð27Þ
p X
1 0 ^ 0
7 7 7 7 7 7 7 7 7 7 7 7 7 5
ð0Þ ^ a^ IFILS ¼ a^ y ¼ R^ 1 y ry
ð31Þ
ð32Þ
Step 3. Set i= i+ 1 and use r^ y ½k’s and a^ ði1Þ IFILS to compute ^ ðiÞ ^ ðiÞ , h the estimates H Step 4. Calculate the estimate of s2e and s2w in the following way: ^ ðiÞ ^ ðiÞ Þ1 H ^ ðiÞT h ^ ðiÞT H h^ ðiÞ ¼ ðH ^ ðiÞ s^ 2ðiÞ w ¼ h ð1Þ ^ ðiÞ s^ 2ðiÞ e ¼ h ð2Þ
ð33Þ
Step 5. Perform the bias correction ð34Þ
Step 6. If
:a^ ði1Þ IFILS :
ð28Þ
k ¼ 0,1,. . .,p þ q
^ y , r^ y . and use them to construct the estimates R Step 2. Set i= 0, where i denotes the iteration step. Then, compute the initial estimate:
^ ði1Þ :a^ ðiÞ IFILS aIFILS :
3
2
bi 6 6 i¼0 6 6X 6 p 6 bi bi þ 1 6 H¼6 i¼0 6 6 ^ 6 6 X p 6 4 bi bi þ q i¼0
N X 1 r^ y k ¼ yðtÞyðtkÞ, N t ¼ pþqþ1
^ ^ 2ðiÞ ^ 1 ^ ði1Þ a^ ðiÞ IFILS ¼ ay þ sw Ry aIFILS
sw h¼ s2e 2
The following steps illustrate our algorithm for estimating noisy AR parameters: Step 1. Compute autocovariance estimates using given samples {y(1),y(2),yy(N)}
rd
ð35Þ
where d is an appropriate small positive number and 99 99 is the Euclidean norm, the convergence is achieved and the iteration process must be terminated. Otherwise, go to step 3. The consistent convergence of the proposed algorithm can be established in a way similar to that of the previous ILS methods proposed by Zheng (see [5] for details).
1662
A. Mahmoudi, M. Karimi / Signal Processing 91 (2011) 1659–1664
4. Simulation results
^ ¼ mðaÞ
In this section, the performance of the proposed IFILS method is presented by simulation results and compared with MYW [1], LS, ILSD [5], and ILSNI [7] methods. We illustrate the performances of the proposed method using two numerical examples. In the first example, the poles of the AR model are close to the unit circle, and in the second example, the poles of the AR model are well inside the unit circle. In all simulations of this section, M is the number of trials set to 1000. The value of q in the IFILS methods is set to 2 and the value of d used in ILSD and IFILS methods for terminating the iterations is set to 0.001. Example 1. Consider a fourth order noisy AR process with 2 noise variance sw and noise free model parameters as follows: aT ¼ 0:5500 0:1550 0:5495 0:6241
s2e ¼ 1 Note that this AR model has a pair of complex conjugate poles located at z =0.3 70.8i and two real poles located at z = 0.95 and z =0.9, which are relatively close to the unit circle. In this example, the parameter q in the MYW method is set to 4 and N = 4000 data samples are used. In Table 1, the mean and standard deviation of M =1000 independent estimates of s2w , s2e , and elements of a are shown for MYW, LS, ILSD, ILSNI, and IFILS methods. The signal to noise ratio (SNR) is defined by SNR 10log10
ð36Þ
s2w
^ :mðaÞa: 99a99
ð37Þ
ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi r P ffi 2 M ^ =M : a a: m m¼1 RMSE ¼
ð39Þ
where a^ m stands for the estimate of a in the mth trial and M is the total number of trials. The results presented in Table 1 show that, for this example, the performance of IFILS method is much better than the other methods in terms of robustness in high noise. We can consider the power spectral density (PSD) obtained by using the AR estimation algorithms as another tool for comparing the performances of these algorithms. The PSD of an AR process is defined as follows:
s2
e Gx ðf Þ ¼ 1 Pp ai ej2pfi 2 i¼1
ð40Þ
In Fig. 1 the true PSD and the PSD of the average model obtained by each algorithm using M =1000 simulation runs are shown. It is clear from this figure that the PSD obtained by the IFILS method is closer to the true PSD than the other estimation algorithms. In Table 2 SNR is set to 10 dB and other simulation parameters are the same as in Table 1. In this case, it can be seen that the performance of the IFILS method is still better than the other methods. Example 2. In this example, we consider another fourth order noisy AR process with noise free model parameters as follows: aT ¼ 1 0:5 0:25 0:0625
s2e ¼ 1 2
mðxðtÞ2 Þ
where m( ) denotes the mean. In Table 1, the SNR is set to 5 dB. The parameter estimation methods are also compared in terms of relative error (RE) and the normalized root mean squared error (RMSE). The RE and RMSE are defined as follows: RE ¼
M 1 X a^ m Mm¼1
ð38Þ
99a99
The SNR is set to be 15 dB by choosing sw =0.1. In contrary to Example 1, this AR model has a pair of complex conjugate poles whose magnitudes are equal to 0.5 and thus, they are well inside the unit circle. Simulation results for N =4000 and 500 are presented in Tables 3 and 4, respectively. It can be seen that, in this case, the performance of the IFILS method is better than the other estimation methods. Note that in Table 4 the parameter q in the MYW method is set to 8. In general, convergence of the iterative ILS-based methods is not guaranteed [16]. The convergence of this type of methods may be affected by the number of samples, location of poles of the AR model, the initial estimate of AR parameters, and the Signal to Noise Ratio
Table 1 Simulation results (M= 1000, N = 4000, and SNR= 5 dB). True value
MYW
LS
ILSD
ILSNI
IFILS
a1 = 0.5500 a2 = 0.1550 a3 = 0.5495 a4 = 0.6241 s2w ¼ 7:3
0.5748 7 0.2648 0.1588 7 0.1663 0.5839 7 0.3200 0.6150 7 0.2234 –
0.0818 7 0.0156 0.0857 7 0.0164 0.0378 70.0181 0.1098 7 0.0167 –
0.5982 70.2800 0.1156 70.0912 0.5849 7 0.2893 0.65067 0.2326 7.08077 0.8765
0.5915 70.3709 0.1089 70.3161 0.5718 70.4524 0.6496 70.3279 7.0833 70.8165
0.5386 7 0.1107 0.1568 7 0.0445 0.5345 70.1253 0.6059 70.0999 7.35517 0.5099
s2e ¼ 1
–
–
1.18427 1.1520
1.1960 71.1277
1.0796 7 0.4668
RE (%) RMSE (%)
4.3176 49.7248
85.8518 85.9157
7.5695 47.5454
7.0065 73.7773
2.6004 19.9644
A. Mahmoudi, M. Karimi / Signal Processing 91 (2011) 1659–1664
1663
30
POWER SPECTRAL DENSITY (dB)
25
LS MYW ILSD IFILS True ILSNI
20 15 10 5 0 -5 -10 -15 0
0.05
0.1
0.15
0.2 0.25 0.3 FREQUENCY
0.35
0.4
0.45
0.5
Fig. 1. Comparison of various noisy AR estimation algorithms in terms of the estimated power spectral density.
Table 2 Simulation results (M = 1000, N = 4000, and SNR =10 dB). True value
MYW
LS
ILSD
ILSNI
IFILS
a1 = 0.5500 a2 = 0.1550 a3 = 0.5495 a4 = 0.6241 s2w ¼ 0:24
0.5521 70.0298 0.15407 0.0228 0.5505 7 0.0291 0.6229 70.0197 –
0.4421 70.0142 0.1772 70.0146 0.43127 0.0156 0.5143 70.0147 –
0.5478 7 0.0268 0.1548 7 0.0141 0.5465 7 0.0285 0.6207 7 0.0259 0.2342 7 0.0449
0.5493 7 0.0268 0.1546 7 0.0147 0.5481 7 0.0293 0.6220 7 0.0261 0.2363 7 0.0461
0.5495 7 0.0181 0.1543 7 0.0139 0.5484 70.0185 0.6227 7 0.0179 0.2390 70.0274
s2e ¼ 1
–
–
1.0013 7 0.0865
1.0036 7 0.0899
0.9968 7 0.0509
RE (%) RMSE (%)
0.2765 5.1039
19.3632 19.5834
0.5046 4.8809
0.2642 4.9273
0.1936 3.4175
Table 3 Simulation results (M = 1000, N = 4000, and SNR =15 dB). True value
MYW
LS
ILSD
ILSNI
IFILS
a1 = 1 a2 = 0.5 a3 = 0.25 a4 = 0.0625 s2w ¼ 0:0740
1.1001 7 13.3010 0.1589 7 8.2474 0.0921 7 13.5144 0.0043 73.7158 –
0.9002 70.0158 0.3642 7 0.0212 0.1564 7 0.0222 0.0236 7 0.0156 –
1.0266 7 1.2835 0.5262 7 1.2911 0.2627 7 0.5989 0.0675 7 0.1983 0.1216 72.1293
1.0378 7 0.3594 0.5638 7 0.6863 0.3053 7 0.6771 0.0888 7 0.3454 0.0841 7 0.0529
0.98087 0.0781 0.4772 70.1104 0.2367 70.0799 0.05747 0.0361 0.0551 70.0537
s2e ¼ 1
–
–
0.9255 73.3853
0.9591 7 0.2621
1.0317 7 0.1186
RE (%) RMSE (%)
43.2983 1830.7
17.1414 17.4560
3.4621 167.8856
8.3778 94.9148
2.8804 14.3346
(SNR). In the examples presented in this paper, the iterative algorithms ILSD and IFILS converged in all simulation runs in few iterations, but the convergence is not guaranteed in general. In all methods proposed in the
literature for estimation of the noisy AR parameters, a relatively high number of data samples is needed to result in a good convergence behavior. This is a shortage for these methods.
1664
A. Mahmoudi, M. Karimi / Signal Processing 91 (2011) 1659–1664
Table 4 Simulation results (M= 1000, N = 500, and SNR= 15 dB). True value
MYW q =8
LS
ILSD
ILSNI
IFILS
a1 = 1 a2 = 0.5 a3 = 0.25 a4 = 0.0625 s2w ¼ 0:0740
1.442470.4347 1.0554 70.6332 0.5425 70.5945 0.1365 70.2222 –
0.8966 70.0406 0.3642 70.0607 0.1584 70.0603 0.02647 0.0459 –
1.0629 7 0.1691 0.6133 7 0.2380 0.3499 7 0.1692 0.1125 7 0.0728 0.1276 7 3.7556
1.0333 71.1907 0.5525 72.3486 0.29107 2.3485 0.08107 1.1910 0.0875 70.0561
0.9835 70.3118 0.51047 0.5500 0.2889 70.5081 0.08997 0.2551 0.0193 70.0728
s2e ¼ 1
–
–
0.1276 7 3.7556
0.9573 7 0.7942
1.0545 7 0.2274
RE (%) RMSE (%)
67.2473 109.8006
17.1697 19.4637
14.9135 44.6114
6.6913 324.4726
4.4857 15.1227
5. Conclusion A new method for estimation of the parameters of noisy autoregressive processes was proposed. This method is based on least-squares estimation. Due to noise variance, the LS estimation of the AR parameters is biased. We used inverse filtering technique and Yule–Walker equations to derive a new method to estimate noise variance. The proposed unbiased parameter estimation algorithm is iterative. Using a simulation study, performance of the proposed method was evaluated and compared with other existing methods. Simulation results showed that the performance of the proposed method is better than the other methods in terms of smaller variance in estimation and better robustness against high observation noise. References [1] S.M. Kay, Modern Spectral Estimation, Prentice-Hall, Englewood Cliffs, NJ, 1988. [2] S.M. Kay, Noise compensation for autoregressive spectral estimates, IEEE Trans. Acoust., Speech, Signal Process. 28 (3) (1980) 292–303. [3] C.E. Davila, A subspace approach to estimation of autoregressive parameters from noisy measurements, IEEE Trans. Signal Process. 46 (2) (1998) 531–534. [4] W.X. Zheng, Unbiased identification of autoregressive signals observed in colored noise, Proc. ICASSP Conf. 4 (1998) 2329–2332. [5] W.X. Zheng, Autoregressive parameter estimation from noisy data, IEEE Trans. Circuits Syst. II: Analog Digital Signal Process. 47 (1) (2000) 71–75.
[6] W.X. Zheng, Fast identification of autoregressive signals from noisy observations, IEEE Trans. Circuits Syst. II 52 (1) (2005) 43–48. [7] W.X. Zheng, On estimation of autoregressive signals in the presence of noise, IEEE Trans. Circuits Syst. II 523 (12) (2006) 1471–1475. [8] A. Mahmoudi, M. Karimi, Estimation of the parameters of multichannel autoregressive signals from noisy observations, Signal Process. 88 (1) (2008) 2777–2783. [9] S.A. Fattah, W.P. Zhu, M.O. Ahmed, A novel technique for the identification of ARMA system under very low levels of SNR, IEEE Trans. Circuits Syst. I 55 (7) (2008) 1988–2001. [10] H. Tong, Autoregressive model fitting with noisy data by Akaikes information criterion, IEEE Trans. Inf. Theory IT-21 (1975) 476–480. [11] D. Aboutajdine, A. Adib, A. Meziane, Fast adaptive algorithms for AR parameter estimation using higher order statistics, IEEE Trans. Signal Process. 44 (1996) 1998–2009. [12] A. Nehorai, P. Stoica, Adaptive algorithms for constrained ARMA signals in the presence of noise, IEEE Trans. Acoust., Speech, Signal Process. 36 (1988) 1282–1291. [13] A. Mahmoudi, M. Karimi, Parameter estimation of autoregressive signals from observations corrupted with colored noise, Signal Process. 88 (1) (2010) 157–164. [14] D. Labarre, E. Grivel, Y. Berthoumieu, E. Todini, M. Najim, Consistent estimation of autoregressive parameters from noisy observations based on two interacting Kalman filters, Signal Process. 86 (10) (2006) 2863–2876. [15] W. Bobillet, R. Diversi, E. Grivel, R. Guidorzi, M. Najim, U. Soverini, Speech enhancement combining optimal smoothing and errors-invariables identification of noisy AR processes, IEEE Trans. Signal Process. 55 (12) (2007) 5564–5578. [16] C.Z. Jin, L.J. Jia, Z.J. Yang, K. Wada, On convergence of a BCLS algorithms for noisy autoregressive process estimation, in: Proceedings of the 41st IEEE Conference on Decision and Control, vol. 4, 2002, pp. 4252–4257.
Contents lists available at ScienceDirect
Signal Processing journal homepage: www.elsevier.com/locate/sigpro
Fast communication
Inverse filtering based method for estimation of noisy autoregressive signals Alimorad Mahmoudi a,n, Mahmood Karimi b a b
Electrical Engineering Department, Shahid Chamran University of Ahvaz, Khouzestan, Iran School of Electrical and Computer Engineering, Shiraz University, Iran
a r t i c l e i n f o
abstract
Article history: Received 7 June 2010 Received in revised form 21 December 2010 Accepted 12 January 2011 Available online 19 January 2011
In this paper we present a new method for estimating the parameters of an autoregressive (AR) signal from observations corrupted with white noise. The leastsquares (LS) estimate of the AR parameters is biased when the observation noise is added to the AR signal. This bias is related to observation noise variance. The proposed method uses inverse filtering technique and Yule–Walker equations for estimating observation noise variance to yield unbiased LS estimate of the AR parameters. The performance of the proposed unbiased algorithm is illustrated by simulation results and they show that the performance of the proposed method is better than the other estimation methods. & 2011 Elsevier B.V. All rights reserved.
Keywords: Autoregressive signals Least-squares method Yule–Walker equations
1. Introduction Modeling of random signals as autoregressive (AR) models are known to be a simple yet effective means for information invoked in many applications, e.g., geophysics, econometrics, speech processing, image processing and communications. The widespread use of AR modeling techniques is due to their excellent resolution performance and existence of a linear set of equations for computing unknown model parameters. It is known that the AR estimator usually shows a high sensitivity to the addition of observation noise. This is an important problem, which limits the utility of these parameter estimation techniques [1]. The noisy AR estimation techniques have been studied in [1–16]. The optimum solution for estimating the AR parameters based on mean squared error (MSE) criteria leads to low-order Yule–Walker equations. When the AR model is corrupted with white noise, due to observation
n
Corresponding author. E-mail addresses: [email protected] (A. Mahmoudi), [email protected] (M. Karimi). 0165-1684/$ - see front matter & 2011 Elsevier B.V. All rights reserved. doi:10.1016/j.sigpro.2011.01.008
noise variance, the zero lag autocorrelation is biased. Due to this, the solution of low-order Yule–Walker equation is biased. The noisy AR estimation techniques can be divided into two main categories. The techniques that belong to the first category avoid zero lag autocorrelations and use high-order Yule–Walker equations. In fact, in this category, the AR(p) signal is modeled by autoregressive moving average ARMA(p,p) model and AR parameters are estimated by ARMA estimation techniques. The modified Yule–Walker (MYW) method [1] is based on this idea. One problem concerning these methods is lack of enough data for computing higher autocorrelation lag estimates, which leads to error in the autocorrelation lags. The second category of common methods used to estimate noisy AR models is based on bias compensation principle. These methods try to remove bias from loworder YW equations. Based on this idea, in [2], Kay uses the Burg estimation method for noisy AR estimation and tries to remove bias from the Burg estimates. In this method, Kay assumed that the observation noise variance is known. This assumption is a limitation for the proposed method. Based on the idea proposed by Kay [2], Zheng [4–8] tried to remove bias from least-squares (LS) estimation of
1660
A. Mahmoudi, M. Karimi / Signal Processing 91 (2011) 1659–1664
noisy AR signals. These methods, which improve the least-squares estimate, are called ILS-based methods. Several ILS-based methods for estimating parameters of noisy AR processes have been proposed by Zheng [4–8]. The main difference between ILS-based methods is that observation noise variance estimation differs in each method. In [5], Zheng proposed an improved leastsquares (ILS) method with direct structure (called ILSD). In [7], he proposed a non-iterative ILS-based method (called ILSNI). In this paper, a new method for estimating AR parameters is proposed. This method is an iterative ILS-based method that uses inverse filtering technique and Yule–Walker equations to estimate the corrupting observation noise variance. In this method, the equations of the linear system are solved and an estimate of observation noise variance is obtained. Then, this estimate is used recursively to correct the bias in the AR parameters estimated by the LS method. Simulation results show that the performance of the proposed method is better than the other estimation methods.
where 2
rx ½0 6 6 rx ½1 Rx ¼ 6 6 ^ 4 rx ½p1
rx ½1
...
rx ½0
...
^ rx ½p2
& ...
rx ½1p
3
7 rx ½2p 7 7 7 ^ 5 rx ½0
ð6Þ
2
3 rx ½1 6 7 6 rx ½2 7 7 rx ¼ 6 6 ^ 7 4 5 rx ½p
ð7Þ
Using (3), (6) and (7), we get Ry ¼ Rx þ s2w Ip
ð8Þ
ry ¼ rx
ð9Þ
where the definitions of Ry and ry are similar to Rx and rx, respectively, and Ip is the p p identity matrix. Combining (5) and (8), and substituting (9), we obtain Ry as2w a ¼ ry
ð10Þ
Multiplying both sides of Eq. (10) by 2. Problem statement
ay ¼ as
ð1Þ
yðtÞ ¼ xðtÞ þ wðtÞ
ð2Þ
where e(t) is the zero hmean stationary driving noise T having variance s2e , a ¼ a1 a2 . . . ap is the vector of parameters of the AR model (T denotes the transpose operation), and w(t) is the zero mean stationary and white observation noise with variance s2w . It is assumed that e(t) and w(t) are mutually uncorrelated, i.e., E{w(t)e(n)}= 0 for all n and t, where E{ } is the expectation operator. The autocorrelation function of y(t) is given by 2 w,
rx ½k þ s
k¼0
rx ½k,
ka0
ð3Þ
rx[k] is the autocorrelation function of noise free AR process x(t). It can be seen from (3) that, due to the observation noise variance, the autocorrelation function of x(t) at zero lag is biased. The main concern of this paper is to estimate the AR coefficients ai (i= 1,2,y,p) using observation data {y(1),y(2),y,y(N)}. Consider the Yule–Walker equations as follows [1]: rx ½k ¼
p X
ay ¼ R1 y ry
kZ1
ð4Þ
i¼1
ð12Þ
ay is the LS estimate of a [4]. Eq. (11) shows that the bias in the LS estimate ay is caused by the observation noise variance s2w . By means of the bias compensation principle presented in [2], the ILS estimate of a may be achieved via a ¼ ay þ s2w R1 y a
ð13Þ
It is clear from Eq. (13) that an unbiased estimate of a is attainable if s2w can be estimated in some way. 3. The estimation algorithm 3.1. Estimation of the observation noise variance In this section, we derive the new method for estimating noisy AR parameters. We use an iterative algorithm for estimating s2w , s2e and a. This algorithm is named as the improved least-squares based on inverse filtering (IFILS). Consider the following vectors: T xt ¼ xðt1Þ xðt2Þ xðtpÞ ð14Þ wt ¼ wðt1Þ
ai rx ½ki,
ð11Þ
where
xðtÞ ¼ a1 xðt1Þ þa2 xðt2Þ þ þ ap xðtpÞ þ eðtÞ
ry ½k ¼ EfyðtÞyðtkÞg ¼
we get
2 1 w Ry a
The noisy AR model under study is represented by
(
R1 y ,
wðt2Þ
yt ¼ xt þ wt ¼ yðt1Þ
wðtpÞ
yðt2Þ
T yðtpÞ
ð15Þ T
ð16Þ
Evaluating (4) for k= 1,y, p gives the following linear system of equations:
If the noisy observation process is filtered by the AR coefficients, then the output z(t) is given by
Rx a ¼ rx
zðtÞ ¼ yðtÞyTt a
ð5Þ
ð17Þ
A. Mahmoudi, M. Karimi / Signal Processing 91 (2011) 1659–1664
Combining (1), (2), (14), and (16), and substituting into (17) gives zðtÞ ¼
wðtÞwTt a þ eðtÞ
ð18Þ
Now, we compute the autocorrelation function of z(t)from (17) and (18). For simplicity, we define a (p+ 1) 1 vector as follows: 1 b¼ ð19Þ a
1661
and 2
3 b b r ½ij i j y 6 7 6 i ¼ 0j ¼ 0 7 6 7 6X 7 p X p 6 7 6 7 b b r ½1 þij i j y 6 7 h ¼ 6 i ¼ 0j ¼ 0 7 6 7 6 7 ^ 6 7 6 p p 7 6 XX 7 4 bi bj ry ½q þij 5 p X p X
ð29Þ
i ¼ 0j ¼ 0
Using (17) and (19), z(t) is expressible as zðtÞ ¼
p X
bi yðtiÞ
ð20Þ
Minimizing the least-squares error of Eq. (26) with respect to h yields the LS estimate as follows:
i¼0
It follows from (20) and (3) that the autocorrelation function of z(t) is equal to rz ½k ¼ EfzðtÞzðtkÞg ¼
p X p X
bi bj ry ½k þ ij
ð21Þ
h ¼ ðHT HÞ1 HT h
ð30Þ
Note that we cannot obtain h from (30) because H and h are unknown. So, we will use a recursive algorithm for estimating h. This algorithm is discussed in Section 3.2.
i¼0j¼0
Using (18) and (19), z(t) is expressible as zðtÞ ¼
wðtÞwTt a þ eðtÞ ¼
p X
3.2. The IFILS algorithm bi wðtiÞ þ eðtÞ
ð22Þ
i¼0
From (22) and (3), the autocorrelation function of z(t), rz[k], can be written as follows: rz ½k ¼ EfzðtÞzðtkÞg ¼ s2w
p X
bi bi þ k þ s2e dðkÞ
ð23Þ
i¼0
where bi + k is taken equal to zero for the values of its index that are outside the interval 0 ri+ krp, and d(k) is the impulse function, which is defined as follows: ( 1, k¼0 dðkÞ ¼ ð24Þ 0, otherwise Using (21) and (23), we obtain
s2w
Xp i¼0
bi bi þ k þ s2e dðkÞ ¼
p X p X
bi bj ry ½k þ ij,
k ¼ 0,1,. . .,q
i¼0j¼0
ð25Þ where q is an arbitrary positive integer less than or equal to p. For illustration convenience, we rewrite Eq. (25) in matrix form as follows: Hh ¼ h
ð26Þ
where " 2# ð27Þ
p X
1 0 ^ 0
7 7 7 7 7 7 7 7 7 7 7 7 7 5
ð0Þ ^ a^ IFILS ¼ a^ y ¼ R^ 1 y ry
ð31Þ
ð32Þ
Step 3. Set i= i+ 1 and use r^ y ½k’s and a^ ði1Þ IFILS to compute ^ ðiÞ ^ ðiÞ , h the estimates H Step 4. Calculate the estimate of s2e and s2w in the following way: ^ ðiÞ ^ ðiÞ Þ1 H ^ ðiÞT h ^ ðiÞT H h^ ðiÞ ¼ ðH ^ ðiÞ s^ 2ðiÞ w ¼ h ð1Þ ^ ðiÞ s^ 2ðiÞ e ¼ h ð2Þ
ð33Þ
Step 5. Perform the bias correction ð34Þ
Step 6. If
:a^ ði1Þ IFILS :
ð28Þ
k ¼ 0,1,. . .,p þ q
^ y , r^ y . and use them to construct the estimates R Step 2. Set i= 0, where i denotes the iteration step. Then, compute the initial estimate:
^ ði1Þ :a^ ðiÞ IFILS aIFILS :
3
2
bi 6 6 i¼0 6 6X 6 p 6 bi bi þ 1 6 H¼6 i¼0 6 6 ^ 6 6 X p 6 4 bi bi þ q i¼0
N X 1 r^ y k ¼ yðtÞyðtkÞ, N t ¼ pþqþ1
^ ^ 2ðiÞ ^ 1 ^ ði1Þ a^ ðiÞ IFILS ¼ ay þ sw Ry aIFILS
sw h¼ s2e 2
The following steps illustrate our algorithm for estimating noisy AR parameters: Step 1. Compute autocovariance estimates using given samples {y(1),y(2),yy(N)}
rd
ð35Þ
where d is an appropriate small positive number and 99 99 is the Euclidean norm, the convergence is achieved and the iteration process must be terminated. Otherwise, go to step 3. The consistent convergence of the proposed algorithm can be established in a way similar to that of the previous ILS methods proposed by Zheng (see [5] for details).
1662
A. Mahmoudi, M. Karimi / Signal Processing 91 (2011) 1659–1664
4. Simulation results
^ ¼ mðaÞ
In this section, the performance of the proposed IFILS method is presented by simulation results and compared with MYW [1], LS, ILSD [5], and ILSNI [7] methods. We illustrate the performances of the proposed method using two numerical examples. In the first example, the poles of the AR model are close to the unit circle, and in the second example, the poles of the AR model are well inside the unit circle. In all simulations of this section, M is the number of trials set to 1000. The value of q in the IFILS methods is set to 2 and the value of d used in ILSD and IFILS methods for terminating the iterations is set to 0.001. Example 1. Consider a fourth order noisy AR process with 2 noise variance sw and noise free model parameters as follows: aT ¼ 0:5500 0:1550 0:5495 0:6241
s2e ¼ 1 Note that this AR model has a pair of complex conjugate poles located at z =0.3 70.8i and two real poles located at z = 0.95 and z =0.9, which are relatively close to the unit circle. In this example, the parameter q in the MYW method is set to 4 and N = 4000 data samples are used. In Table 1, the mean and standard deviation of M =1000 independent estimates of s2w , s2e , and elements of a are shown for MYW, LS, ILSD, ILSNI, and IFILS methods. The signal to noise ratio (SNR) is defined by SNR 10log10
ð36Þ
s2w
^ :mðaÞa: 99a99
ð37Þ
ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi r P ffi 2 M ^ =M : a a: m m¼1 RMSE ¼
ð39Þ
where a^ m stands for the estimate of a in the mth trial and M is the total number of trials. The results presented in Table 1 show that, for this example, the performance of IFILS method is much better than the other methods in terms of robustness in high noise. We can consider the power spectral density (PSD) obtained by using the AR estimation algorithms as another tool for comparing the performances of these algorithms. The PSD of an AR process is defined as follows:
s2
e Gx ðf Þ ¼ 1 Pp ai ej2pfi 2 i¼1
ð40Þ
In Fig. 1 the true PSD and the PSD of the average model obtained by each algorithm using M =1000 simulation runs are shown. It is clear from this figure that the PSD obtained by the IFILS method is closer to the true PSD than the other estimation algorithms. In Table 2 SNR is set to 10 dB and other simulation parameters are the same as in Table 1. In this case, it can be seen that the performance of the IFILS method is still better than the other methods. Example 2. In this example, we consider another fourth order noisy AR process with noise free model parameters as follows: aT ¼ 1 0:5 0:25 0:0625
s2e ¼ 1 2
mðxðtÞ2 Þ
where m( ) denotes the mean. In Table 1, the SNR is set to 5 dB. The parameter estimation methods are also compared in terms of relative error (RE) and the normalized root mean squared error (RMSE). The RE and RMSE are defined as follows: RE ¼
M 1 X a^ m Mm¼1
ð38Þ
99a99
The SNR is set to be 15 dB by choosing sw =0.1. In contrary to Example 1, this AR model has a pair of complex conjugate poles whose magnitudes are equal to 0.5 and thus, they are well inside the unit circle. Simulation results for N =4000 and 500 are presented in Tables 3 and 4, respectively. It can be seen that, in this case, the performance of the IFILS method is better than the other estimation methods. Note that in Table 4 the parameter q in the MYW method is set to 8. In general, convergence of the iterative ILS-based methods is not guaranteed [16]. The convergence of this type of methods may be affected by the number of samples, location of poles of the AR model, the initial estimate of AR parameters, and the Signal to Noise Ratio
Table 1 Simulation results (M= 1000, N = 4000, and SNR= 5 dB). True value
MYW
LS
ILSD
ILSNI
IFILS
a1 = 0.5500 a2 = 0.1550 a3 = 0.5495 a4 = 0.6241 s2w ¼ 7:3
0.5748 7 0.2648 0.1588 7 0.1663 0.5839 7 0.3200 0.6150 7 0.2234 –
0.0818 7 0.0156 0.0857 7 0.0164 0.0378 70.0181 0.1098 7 0.0167 –
0.5982 70.2800 0.1156 70.0912 0.5849 7 0.2893 0.65067 0.2326 7.08077 0.8765
0.5915 70.3709 0.1089 70.3161 0.5718 70.4524 0.6496 70.3279 7.0833 70.8165
0.5386 7 0.1107 0.1568 7 0.0445 0.5345 70.1253 0.6059 70.0999 7.35517 0.5099
s2e ¼ 1
–
–
1.18427 1.1520
1.1960 71.1277
1.0796 7 0.4668
RE (%) RMSE (%)
4.3176 49.7248
85.8518 85.9157
7.5695 47.5454
7.0065 73.7773
2.6004 19.9644
A. Mahmoudi, M. Karimi / Signal Processing 91 (2011) 1659–1664
1663
30
POWER SPECTRAL DENSITY (dB)
25
LS MYW ILSD IFILS True ILSNI
20 15 10 5 0 -5 -10 -15 0
0.05
0.1
0.15
0.2 0.25 0.3 FREQUENCY
0.35
0.4
0.45
0.5
Fig. 1. Comparison of various noisy AR estimation algorithms in terms of the estimated power spectral density.
Table 2 Simulation results (M = 1000, N = 4000, and SNR =10 dB). True value
MYW
LS
ILSD
ILSNI
IFILS
a1 = 0.5500 a2 = 0.1550 a3 = 0.5495 a4 = 0.6241 s2w ¼ 0:24
0.5521 70.0298 0.15407 0.0228 0.5505 7 0.0291 0.6229 70.0197 –
0.4421 70.0142 0.1772 70.0146 0.43127 0.0156 0.5143 70.0147 –
0.5478 7 0.0268 0.1548 7 0.0141 0.5465 7 0.0285 0.6207 7 0.0259 0.2342 7 0.0449
0.5493 7 0.0268 0.1546 7 0.0147 0.5481 7 0.0293 0.6220 7 0.0261 0.2363 7 0.0461
0.5495 7 0.0181 0.1543 7 0.0139 0.5484 70.0185 0.6227 7 0.0179 0.2390 70.0274
s2e ¼ 1
–
–
1.0013 7 0.0865
1.0036 7 0.0899
0.9968 7 0.0509
RE (%) RMSE (%)
0.2765 5.1039
19.3632 19.5834
0.5046 4.8809
0.2642 4.9273
0.1936 3.4175
Table 3 Simulation results (M = 1000, N = 4000, and SNR =15 dB). True value
MYW
LS
ILSD
ILSNI
IFILS
a1 = 1 a2 = 0.5 a3 = 0.25 a4 = 0.0625 s2w ¼ 0:0740
1.1001 7 13.3010 0.1589 7 8.2474 0.0921 7 13.5144 0.0043 73.7158 –
0.9002 70.0158 0.3642 7 0.0212 0.1564 7 0.0222 0.0236 7 0.0156 –
1.0266 7 1.2835 0.5262 7 1.2911 0.2627 7 0.5989 0.0675 7 0.1983 0.1216 72.1293
1.0378 7 0.3594 0.5638 7 0.6863 0.3053 7 0.6771 0.0888 7 0.3454 0.0841 7 0.0529
0.98087 0.0781 0.4772 70.1104 0.2367 70.0799 0.05747 0.0361 0.0551 70.0537
s2e ¼ 1
–
–
0.9255 73.3853
0.9591 7 0.2621
1.0317 7 0.1186
RE (%) RMSE (%)
43.2983 1830.7
17.1414 17.4560
3.4621 167.8856
8.3778 94.9148
2.8804 14.3346
(SNR). In the examples presented in this paper, the iterative algorithms ILSD and IFILS converged in all simulation runs in few iterations, but the convergence is not guaranteed in general. In all methods proposed in the
literature for estimation of the noisy AR parameters, a relatively high number of data samples is needed to result in a good convergence behavior. This is a shortage for these methods.
1664
A. Mahmoudi, M. Karimi / Signal Processing 91 (2011) 1659–1664
Table 4 Simulation results (M= 1000, N = 500, and SNR= 15 dB). True value
MYW q =8
LS
ILSD
ILSNI
IFILS
a1 = 1 a2 = 0.5 a3 = 0.25 a4 = 0.0625 s2w ¼ 0:0740
1.442470.4347 1.0554 70.6332 0.5425 70.5945 0.1365 70.2222 –
0.8966 70.0406 0.3642 70.0607 0.1584 70.0603 0.02647 0.0459 –
1.0629 7 0.1691 0.6133 7 0.2380 0.3499 7 0.1692 0.1125 7 0.0728 0.1276 7 3.7556
1.0333 71.1907 0.5525 72.3486 0.29107 2.3485 0.08107 1.1910 0.0875 70.0561
0.9835 70.3118 0.51047 0.5500 0.2889 70.5081 0.08997 0.2551 0.0193 70.0728
s2e ¼ 1
–
–
0.1276 7 3.7556
0.9573 7 0.7942
1.0545 7 0.2274
RE (%) RMSE (%)
67.2473 109.8006
17.1697 19.4637
14.9135 44.6114
6.6913 324.4726
4.4857 15.1227
5. Conclusion A new method for estimation of the parameters of noisy autoregressive processes was proposed. This method is based on least-squares estimation. Due to noise variance, the LS estimation of the AR parameters is biased. We used inverse filtering technique and Yule–Walker equations to derive a new method to estimate noise variance. The proposed unbiased parameter estimation algorithm is iterative. Using a simulation study, performance of the proposed method was evaluated and compared with other existing methods. Simulation results showed that the performance of the proposed method is better than the other methods in terms of smaller variance in estimation and better robustness against high observation noise. References [1] S.M. Kay, Modern Spectral Estimation, Prentice-Hall, Englewood Cliffs, NJ, 1988. [2] S.M. Kay, Noise compensation for autoregressive spectral estimates, IEEE Trans. Acoust., Speech, Signal Process. 28 (3) (1980) 292–303. [3] C.E. Davila, A subspace approach to estimation of autoregressive parameters from noisy measurements, IEEE Trans. Signal Process. 46 (2) (1998) 531–534. [4] W.X. Zheng, Unbiased identification of autoregressive signals observed in colored noise, Proc. ICASSP Conf. 4 (1998) 2329–2332. [5] W.X. Zheng, Autoregressive parameter estimation from noisy data, IEEE Trans. Circuits Syst. II: Analog Digital Signal Process. 47 (1) (2000) 71–75.
[6] W.X. Zheng, Fast identification of autoregressive signals from noisy observations, IEEE Trans. Circuits Syst. II 52 (1) (2005) 43–48. [7] W.X. Zheng, On estimation of autoregressive signals in the presence of noise, IEEE Trans. Circuits Syst. II 523 (12) (2006) 1471–1475. [8] A. Mahmoudi, M. Karimi, Estimation of the parameters of multichannel autoregressive signals from noisy observations, Signal Process. 88 (1) (2008) 2777–2783. [9] S.A. Fattah, W.P. Zhu, M.O. Ahmed, A novel technique for the identification of ARMA system under very low levels of SNR, IEEE Trans. Circuits Syst. I 55 (7) (2008) 1988–2001. [10] H. Tong, Autoregressive model fitting with noisy data by Akaikes information criterion, IEEE Trans. Inf. Theory IT-21 (1975) 476–480. [11] D. Aboutajdine, A. Adib, A. Meziane, Fast adaptive algorithms for AR parameter estimation using higher order statistics, IEEE Trans. Signal Process. 44 (1996) 1998–2009. [12] A. Nehorai, P. Stoica, Adaptive algorithms for constrained ARMA signals in the presence of noise, IEEE Trans. Acoust., Speech, Signal Process. 36 (1988) 1282–1291. [13] A. Mahmoudi, M. Karimi, Parameter estimation of autoregressive signals from observations corrupted with colored noise, Signal Process. 88 (1) (2010) 157–164. [14] D. Labarre, E. Grivel, Y. Berthoumieu, E. Todini, M. Najim, Consistent estimation of autoregressive parameters from noisy observations based on two interacting Kalman filters, Signal Process. 86 (10) (2006) 2863–2876. [15] W. Bobillet, R. Diversi, E. Grivel, R. Guidorzi, M. Najim, U. Soverini, Speech enhancement combining optimal smoothing and errors-invariables identification of noisy AR processes, IEEE Trans. Signal Process. 55 (12) (2007) 5564–5578. [16] C.Z. Jin, L.J. Jia, Z.J. Yang, K. Wada, On convergence of a BCLS algorithms for noisy autoregressive process estimation, in: Proceedings of the 41st IEEE Conference on Decision and Control, vol. 4, 2002, pp. 4252–4257.