an effective approach to adaptive iir filtering - CiteSeerX
AN EFFECTIVE APPROACH TO ADAPTIVE IIR FILTERING M Harteneck, R W Stewart
J G McWhirter, I K Proudler
Signal Processing Division Department of EEE University of Strathclyde Glasgow, G1 1XW, U.K. email: [email protected]
DRA Malvern St Andrews Road Great Malvern Worchestershire, WR14 3PS, U.K. email: [email protected]
ABSTRACT
2. DERIVATION
In this paper an approach to adaptive IIR ltering based on a pseudo-linear regression and a QR matrix decomposition is developed. The algorithm has proved to be stable and has good convergence properties if the unknown system satis es the strictly positive real condition. The derivation of the algorithm is straightforward and the computational complexity is less than the computational complexity of the IIR-RPE algorithm. Simulation results of system identi cation with synthetic and real world data are shown comparing the algorithm with the IIR-RPE and the IIR-LMS algorithm.
The output of an adaptive lter at time k can be written as
1. INTRODUCTION
e(k) = y(k) ? y^(k);
Adaptive lters have a widespread use in many di erent signal processing applications, such as system modeling and inverse system modeling. Traditionally most of the current implementations are realized using FIR lters. However, in applications which require lters with very long impulse responses, the use of IIR lters would be advantageous. To date however, most of the adaptive IIR algorithms are based on gradient search methods such as Gauss-Newton [1, 2] to minimize an output-error performance criterion. These algorithms may converge to a local minimum of the multi-modal performance surface and usually require computationally expensive pole monitoring to ensure stability. Therefore, in this paper we investigate a pseudolinear regression (PLR) with a QR matrix decomposition algorithm in an attempt to solve a least squares criterion for an adaptive IIR lter.
where y(k) is the desired signal. In a least squares formulation, the aim of the weight update algorithm is to minimize the power of the error signal assuming the weights a^i (k) and ^bi (k) to be xed for all time prior to k. In the following equations this condition is denoted by subscript w^ . The performance criterion is therefore
y^(k) =
M ?1 X i=0
a^i (k)x(k ? i) +
LX ?1 i=1
^bi(k)^y (k ? i);
(1)
where x(k) is the driving input, y^(k) is the output sequence of the adaptive lter and fa^g and f^bg are the adaptive feedforward and feedback weights. The output error signal e(k) of the adaptive lter can then be written as
(k) =
k X i=0
i e2 (k ? i)jw^ ;
(2)
(3)
where is a forgetting factor smaller then 1. To present this performance criterion in a more structured way it is useful to de ne some vectors and write it in a matrix-
vector notation. The necessary de nitions are ^a(k) = [^a0(k) ^a1(k) : : : a^M ?1(k)]T b^(k) = [^b1(k) ^b2(k) : : : ^bL?1(k)]T x(k) = [x(k) x(k ? 1) : : : x(k ? M + 1)]T y^(k)jw^ = [^y(k ? 1)jw^ y^(k ? 2)jw^ : : : : : : y^(k ? L + 1)jw^ ]T y^(k) = [^y(k ? 1) y^(k ? 2) : : : y^(k ? L + 1)]T y(k) = [y(1) y(2) : : : y(k)]T e(k) = [e(1) e(2) : : : e(k)]T (k) =diagonal(k ; k?1; : : : ; ; 1)
(4a) (4b) (4c) (4d) (4e) (4f) (4g) (4h)
where \T" denotes matrix transpose. In particular note the di erence between (4d) and (4e). The elements of (4d) are calculated assuming the adaptive weights were xed for all time prior to k, whereas (4e) uses the time varying weights, i.e. (4e) is composed of the output signal y^(k) of the adaptive lter. Now (3) can be written as
2 3 2
3 62y(1) 3 2xT (1) y 7
T ^ (1) j 6 7
w ^ 66y(2) 7 6xT (2) y 7
T 7 ^ (2)jw^ 7 ^a(k) 77
1 66 7 6 ? (k) =
2 (k) 6 6 7 6 .. 75 b^ (k) 7 64 ... 5 4 .. . . 6
| {z }7 6 7
T T y ( k ) x ( k ) y ^ ( k ) j
4 w ^ w ^ (k ) 5
| {z } | {z }
y(k)
A(k)jw^
= 12 (k)y(k) ? 12 (k)A(k)jw^ w^ (k) 2
(5)
where k:k denotes Euclidean norm. The weight update algorithm has to minimize the norm of (5). This is a nonlinear minimization problem because of the dependency of the data matrix A(k)jw^ on the adaptive weights w^ . One possible approximation is to substitute y^(k)jw^ by y^(k) and so break the weight dependency. Eq. (5) can then be approximated as
22 3 2 T
y(1) xT (1) y^TT (1)3 3
2
1 66y(2) 7 6x (2) y ^ (2)77 ^a(k) 77
6 7 6 (k)
2 (k) 6 66 .. 7 ? 6 .. .. 75 b^(k) 75
44 . 5 4 . .
y(k) xT (k) y^T (k)
12 (k)y(k) ? 12 (k)A(k)w^ (k) 2
(6) This approximation is known in the literature as pseudo-linear regression (PLR) and its convergence criteria are reviewed in [3, 2] and studied in [4]. In [4] a convergence proof for the PLR is given if (for the
sucient-order case) the transfer function B(1z) of the unknown system is strictly positive real (SPR), i.e. 1 1 <e B (e|! ) ? 2 > 0; 8 ? < ! < : To achieve the minimization of (6), this equation is premultiplied by a k k unitary matrix Q(k) to produce a structured matrix equation which is easy to solve. The unitary matrix is chosen in such a way that 1 R ( k) 2 ; (7) Q(k) (k)A(k) =
0
where R(k) is an (M + L ? 1) (M + L ? 1) uppertriangular matrix and 0 is a (k ?M ?L+1) (M + L? 1) zero matrix. This matrix decomposition is referred to in the literature as QR decomposition [5, 6]. Substituting (7) in (6) gives (k) kQ(k) 2 (k)y(k) ? Q(k) 2 (k)A(k)w ^ (k)k2
2
p(k) ?vR(k()k)w^ (k)
1
1
(8)
where the vector Q(k) 21 (k)y(k) is partitioned into the vectors p(k) and v(k). R(k) is an upper-triangular matrix and so the weights which minimize this performance criterion can be extracted via backsubstitution. A time recursive and a direct residual extraction algorithm can be derived similar to the derivation shown in [5] for the FIR-QR lter. If standard Given's rotations are used for the calculation of the QR decomposition, then the overall complexity of the algorithm is 2N 2 + 5N multiply-accumulates (MACs), N square roots and N divides, where N is the overall number of adaptive weights (N = M + L ? 1). Note that the computational requirements for a simpli ed gradient RPE [2] is 5N 2 + 3N MACs and one division.
3. SIMULATIONS Simulations were done in a system identi cation setup where the unknown system was a 10th-order system where the feedforward polynomial was produced with arbitrary weights and the feedback polynomialis chosen as B (z ?1 ) = 1 ? 0:8z ?10: This model is often used to simulate acoustic transfer paths of teleconference systems [7] and produces poles close to the unit circle. More important is that it satis es the SPR condition. The pole-zero plot of the unknown system is shown in Fig. 1.
Smoothed Averaged Squared Error
40
2 1.5 1
Imag(z)
0.5 0
35 b
30
25 a 20
c
15 0
5000
10000
-0.5 -1 -1.5 -2 -2
-1.5
-1
-0.5
0 0.5 Real(z)
1
1.5
2
15000 Sample
20000
25000
30000
Figure 3: Smoothed Ensemble Squared Error Signals (over 50 Simulations) of System Identi cation with 15 dB observation noise, (a) FIR-LMS Algorithm, (b) IIRRPE Algorithm, (c) IIR-QR Algorithm
Figure 1: Pole-Zero Plot of the Unknown System
2 1.5
Forward Weights
1 0.5 0 -0.5 -1 -1.5 -2 -2.5 0
Averaged Squared Error (dB)
20
2000
4000
6000 8000 10000 12000 14000 Sample
Figure 4: Feedforward Weights with 15 dB observation noise.
0 a -20 b
-40 -60 -80
c -100 -120 0
5000
10000 Sample
15000
20000
Figure 2: Ensemble Squared Error Signals (over 50 Simulations) of System Identi cation with no observation noise, (a) FIR-LMS Algorithm, (b) IIR-RPE Algorithm, (c) IIR-QR Algorithm
The rst set of simulations was system identi cation in a noise free environment. In Fig. 2, the averaged squared error (from 50 simulations) against time are shown for (a) FIR-LMS algorithm [8] with 350 taps, (b) IIR-RPE algorithm [5] with 11 feedforward and 10 feedback weights and (c) IIR-QR algorithm with 11 feedforward and 10 feedback weights. The IIR-QR algorithm adapts almost instantaneously to a very low mean squared error (MSE). The level of this MSE is due to the forgetting factor of = 0:9999 and the xed point data used as input. A close investigation of the weights shows however, that the algorithm converged to the correct set of weights. The IIR-RPE algorithm converged to the same solution but it takes a much longer time. The FIR-LMS converges to a solution poorer than the minimum found by the IIR algorithms; however this FIR minimum could be lowered by increasing the number of weights used.
the IIR-QR needs considerably less taps than the FIRLMS algorithm.
1
Feedback Weights
0.8
4. CONCLUSION
0.6 0.4 0.2 0 -0.2 -0.4 0
2000
4000
6000 8000 10000 12000 14000 Sample
Figure 5: Feedback Weights with 15 dB observation noise. FIR-LMS 50 taps 100 taps 150 taps
ERLE IIR-QR 1.6 dB L = M = 9 2.9 dB L = M = 12 4.2 dB L = M = 21
ERLE 1.5 dB 2.7 dB 4.3 dB
Table 1: Real World Data: slightly reverberant environment Fig. 3 shows the same system identi cation but this time with an observation noise present so that 15 dB SNR in the desired signal y(k) is achieved. Fig. 4 and 5 show the convergence plots for the lter weights in the noisy case using the proposed IIR-QR algorithm. Inspection of the lter coecients in the above cases reveal that they are close to the true values. Table 1 and 2 show the results of the identi cation of a real world acoustic transfer path to compare the number of adaptive weights necessary to model the system. The 16bit/8kHz data was captured inside two acoustic enclosures which are slightly reverberant (length of impulse response about 0.06 seconds,Table 1) and less reverberant (length of impulse response about 0.02 seconds,Table 2). The criterion chosen for this comparison was the power of the desired signal y(k) divided by the power of the error signal e(k) after adaptation, which is called in echo cancelation systems echo-return-loss-enhancement (ERLE). The results show that to achieve the same amount of ERLE, FIR-LMS ERLE IIR-QR ERLE 50 taps 7.1 dB L = M = 19 7.1 dB 100 taps 10.2 dB L = M = 30 10.8 dB 150 taps 11.2 dB L = M = 35 11.5 dB Table 2: Real World Data: less reverberant environment
In this paper we proposed a straightforward and ef cient adaptive algorithm for IIR lters based on a pseudo-linear regression interpretation of the adaptive IIR ltering problem and on a QR matrix decomposition. To date oating point simulations with synthetic data and real acoustic data have demonstrated that the IIR-QR algorithm is a powerful and stable algorithm and superior to gradient search based algorithms in terms of adaptation speed, achievable minimum mean squared error and computational requirements.
5. REFERENCES [1] Feintuch P L. An Adaptive Recursive LMS Filter. Proceedings of the IEEE, pages 1622 { 1624, November 1976. [2] Shynk J J. Adaptive IIR Filtering. IEEE ASSP Magazine, pages 4 { 21, November 1989. [3] Ljung L. System Identi cation: Theory for the User. Prentice-Hall Information and Systems Sciences Series. Prentice-Hall, Inc., 1987. [4] Ljung L and Soderstrom T. Theory and Practice of Recursive Identi cation. MIT Press, Cambridge, Massachusetts, 1983. [5] Haykin S. Adaptive Filter Theory. Prentice Hall, 2nd edition, 1991. ISBN 0-13-012236-5. [6] Golub G H and Van Loan C F. Matrix Computations. John Hopkins University Press, 1989. [7] Kawabe S Chao J and Tsujii S. A New IIR Adaptive Echo Canceller : GIVE. IEEE Transactions on Selected Areas in Communications, 12(9):1530 { 1539, December 1994. [8] Widrow B and Stearns S D. Adaptive Signal Processing. Prentice Hall, Englewood Cli s, 1985.
J G McWhirter, I K Proudler
Signal Processing Division Department of EEE University of Strathclyde Glasgow, G1 1XW, U.K. email: [email protected]
DRA Malvern St Andrews Road Great Malvern Worchestershire, WR14 3PS, U.K. email: [email protected]
ABSTRACT
2. DERIVATION
In this paper an approach to adaptive IIR ltering based on a pseudo-linear regression and a QR matrix decomposition is developed. The algorithm has proved to be stable and has good convergence properties if the unknown system satis es the strictly positive real condition. The derivation of the algorithm is straightforward and the computational complexity is less than the computational complexity of the IIR-RPE algorithm. Simulation results of system identi cation with synthetic and real world data are shown comparing the algorithm with the IIR-RPE and the IIR-LMS algorithm.
The output of an adaptive lter at time k can be written as
1. INTRODUCTION
e(k) = y(k) ? y^(k);
Adaptive lters have a widespread use in many di erent signal processing applications, such as system modeling and inverse system modeling. Traditionally most of the current implementations are realized using FIR lters. However, in applications which require lters with very long impulse responses, the use of IIR lters would be advantageous. To date however, most of the adaptive IIR algorithms are based on gradient search methods such as Gauss-Newton [1, 2] to minimize an output-error performance criterion. These algorithms may converge to a local minimum of the multi-modal performance surface and usually require computationally expensive pole monitoring to ensure stability. Therefore, in this paper we investigate a pseudolinear regression (PLR) with a QR matrix decomposition algorithm in an attempt to solve a least squares criterion for an adaptive IIR lter.
where y(k) is the desired signal. In a least squares formulation, the aim of the weight update algorithm is to minimize the power of the error signal assuming the weights a^i (k) and ^bi (k) to be xed for all time prior to k. In the following equations this condition is denoted by subscript w^ . The performance criterion is therefore
y^(k) =
M ?1 X i=0
a^i (k)x(k ? i) +
LX ?1 i=1
^bi(k)^y (k ? i);
(1)
where x(k) is the driving input, y^(k) is the output sequence of the adaptive lter and fa^g and f^bg are the adaptive feedforward and feedback weights. The output error signal e(k) of the adaptive lter can then be written as
(k) =
k X i=0
i e2 (k ? i)jw^ ;
(2)
(3)
where is a forgetting factor smaller then 1. To present this performance criterion in a more structured way it is useful to de ne some vectors and write it in a matrix-
vector notation. The necessary de nitions are ^a(k) = [^a0(k) ^a1(k) : : : a^M ?1(k)]T b^(k) = [^b1(k) ^b2(k) : : : ^bL?1(k)]T x(k) = [x(k) x(k ? 1) : : : x(k ? M + 1)]T y^(k)jw^ = [^y(k ? 1)jw^ y^(k ? 2)jw^ : : : : : : y^(k ? L + 1)jw^ ]T y^(k) = [^y(k ? 1) y^(k ? 2) : : : y^(k ? L + 1)]T y(k) = [y(1) y(2) : : : y(k)]T e(k) = [e(1) e(2) : : : e(k)]T (k) =diagonal(k ; k?1; : : : ; ; 1)
(4a) (4b) (4c) (4d) (4e) (4f) (4g) (4h)
where \T" denotes matrix transpose. In particular note the di erence between (4d) and (4e). The elements of (4d) are calculated assuming the adaptive weights were xed for all time prior to k, whereas (4e) uses the time varying weights, i.e. (4e) is composed of the output signal y^(k) of the adaptive lter. Now (3) can be written as
2 3 2
3 62y(1) 3 2xT (1) y 7
T ^ (1) j 6 7
w ^ 66y(2) 7 6xT (2) y 7
T 7 ^ (2)jw^ 7 ^a(k) 77
1 66 7 6 ? (k) =
2 (k) 6 6 7 6 .. 75 b^ (k) 7 64 ... 5 4 .. . . 6
| {z }7 6 7
T T y ( k ) x ( k ) y ^ ( k ) j
4 w ^ w ^ (k ) 5
| {z } | {z }
y(k)
A(k)jw^
= 12 (k)y(k) ? 12 (k)A(k)jw^ w^ (k) 2
(5)
where k:k denotes Euclidean norm. The weight update algorithm has to minimize the norm of (5). This is a nonlinear minimization problem because of the dependency of the data matrix A(k)jw^ on the adaptive weights w^ . One possible approximation is to substitute y^(k)jw^ by y^(k) and so break the weight dependency. Eq. (5) can then be approximated as
22 3 2 T
y(1) xT (1) y^TT (1)3 3
2
1 66y(2) 7 6x (2) y ^ (2)77 ^a(k) 77
6 7 6 (k)
2 (k) 6 66 .. 7 ? 6 .. .. 75 b^(k) 75
44 . 5 4 . .
y(k) xT (k) y^T (k)
12 (k)y(k) ? 12 (k)A(k)w^ (k) 2
(6) This approximation is known in the literature as pseudo-linear regression (PLR) and its convergence criteria are reviewed in [3, 2] and studied in [4]. In [4] a convergence proof for the PLR is given if (for the
sucient-order case) the transfer function B(1z) of the unknown system is strictly positive real (SPR), i.e. 1 1 <e B (e|! ) ? 2 > 0; 8 ? < ! < : To achieve the minimization of (6), this equation is premultiplied by a k k unitary matrix Q(k) to produce a structured matrix equation which is easy to solve. The unitary matrix is chosen in such a way that 1 R ( k) 2 ; (7) Q(k) (k)A(k) =
0
where R(k) is an (M + L ? 1) (M + L ? 1) uppertriangular matrix and 0 is a (k ?M ?L+1) (M + L? 1) zero matrix. This matrix decomposition is referred to in the literature as QR decomposition [5, 6]. Substituting (7) in (6) gives (k) kQ(k) 2 (k)y(k) ? Q(k) 2 (k)A(k)w ^ (k)k2
2
p(k) ?vR(k()k)w^ (k)
1
1
(8)
where the vector Q(k) 21 (k)y(k) is partitioned into the vectors p(k) and v(k). R(k) is an upper-triangular matrix and so the weights which minimize this performance criterion can be extracted via backsubstitution. A time recursive and a direct residual extraction algorithm can be derived similar to the derivation shown in [5] for the FIR-QR lter. If standard Given's rotations are used for the calculation of the QR decomposition, then the overall complexity of the algorithm is 2N 2 + 5N multiply-accumulates (MACs), N square roots and N divides, where N is the overall number of adaptive weights (N = M + L ? 1). Note that the computational requirements for a simpli ed gradient RPE [2] is 5N 2 + 3N MACs and one division.
3. SIMULATIONS Simulations were done in a system identi cation setup where the unknown system was a 10th-order system where the feedforward polynomial was produced with arbitrary weights and the feedback polynomialis chosen as B (z ?1 ) = 1 ? 0:8z ?10: This model is often used to simulate acoustic transfer paths of teleconference systems [7] and produces poles close to the unit circle. More important is that it satis es the SPR condition. The pole-zero plot of the unknown system is shown in Fig. 1.
Smoothed Averaged Squared Error
40
2 1.5 1
Imag(z)
0.5 0
35 b
30
25 a 20
c
15 0
5000
10000
-0.5 -1 -1.5 -2 -2
-1.5
-1
-0.5
0 0.5 Real(z)
1
1.5
2
15000 Sample
20000
25000
30000
Figure 3: Smoothed Ensemble Squared Error Signals (over 50 Simulations) of System Identi cation with 15 dB observation noise, (a) FIR-LMS Algorithm, (b) IIRRPE Algorithm, (c) IIR-QR Algorithm
Figure 1: Pole-Zero Plot of the Unknown System
2 1.5
Forward Weights
1 0.5 0 -0.5 -1 -1.5 -2 -2.5 0
Averaged Squared Error (dB)
20
2000
4000
6000 8000 10000 12000 14000 Sample
Figure 4: Feedforward Weights with 15 dB observation noise.
0 a -20 b
-40 -60 -80
c -100 -120 0
5000
10000 Sample
15000
20000
Figure 2: Ensemble Squared Error Signals (over 50 Simulations) of System Identi cation with no observation noise, (a) FIR-LMS Algorithm, (b) IIR-RPE Algorithm, (c) IIR-QR Algorithm
The rst set of simulations was system identi cation in a noise free environment. In Fig. 2, the averaged squared error (from 50 simulations) against time are shown for (a) FIR-LMS algorithm [8] with 350 taps, (b) IIR-RPE algorithm [5] with 11 feedforward and 10 feedback weights and (c) IIR-QR algorithm with 11 feedforward and 10 feedback weights. The IIR-QR algorithm adapts almost instantaneously to a very low mean squared error (MSE). The level of this MSE is due to the forgetting factor of = 0:9999 and the xed point data used as input. A close investigation of the weights shows however, that the algorithm converged to the correct set of weights. The IIR-RPE algorithm converged to the same solution but it takes a much longer time. The FIR-LMS converges to a solution poorer than the minimum found by the IIR algorithms; however this FIR minimum could be lowered by increasing the number of weights used.
the IIR-QR needs considerably less taps than the FIRLMS algorithm.
1
Feedback Weights
0.8
4. CONCLUSION
0.6 0.4 0.2 0 -0.2 -0.4 0
2000
4000
6000 8000 10000 12000 14000 Sample
Figure 5: Feedback Weights with 15 dB observation noise. FIR-LMS 50 taps 100 taps 150 taps
ERLE IIR-QR 1.6 dB L = M = 9 2.9 dB L = M = 12 4.2 dB L = M = 21
ERLE 1.5 dB 2.7 dB 4.3 dB
Table 1: Real World Data: slightly reverberant environment Fig. 3 shows the same system identi cation but this time with an observation noise present so that 15 dB SNR in the desired signal y(k) is achieved. Fig. 4 and 5 show the convergence plots for the lter weights in the noisy case using the proposed IIR-QR algorithm. Inspection of the lter coecients in the above cases reveal that they are close to the true values. Table 1 and 2 show the results of the identi cation of a real world acoustic transfer path to compare the number of adaptive weights necessary to model the system. The 16bit/8kHz data was captured inside two acoustic enclosures which are slightly reverberant (length of impulse response about 0.06 seconds,Table 1) and less reverberant (length of impulse response about 0.02 seconds,Table 2). The criterion chosen for this comparison was the power of the desired signal y(k) divided by the power of the error signal e(k) after adaptation, which is called in echo cancelation systems echo-return-loss-enhancement (ERLE). The results show that to achieve the same amount of ERLE, FIR-LMS ERLE IIR-QR ERLE 50 taps 7.1 dB L = M = 19 7.1 dB 100 taps 10.2 dB L = M = 30 10.8 dB 150 taps 11.2 dB L = M = 35 11.5 dB Table 2: Real World Data: less reverberant environment
In this paper we proposed a straightforward and ef cient adaptive algorithm for IIR lters based on a pseudo-linear regression interpretation of the adaptive IIR ltering problem and on a QR matrix decomposition. To date oating point simulations with synthetic data and real acoustic data have demonstrated that the IIR-QR algorithm is a powerful and stable algorithm and superior to gradient search based algorithms in terms of adaptation speed, achievable minimum mean squared error and computational requirements.
5. REFERENCES [1] Feintuch P L. An Adaptive Recursive LMS Filter. Proceedings of the IEEE, pages 1622 { 1624, November 1976. [2] Shynk J J. Adaptive IIR Filtering. IEEE ASSP Magazine, pages 4 { 21, November 1989. [3] Ljung L. System Identi cation: Theory for the User. Prentice-Hall Information and Systems Sciences Series. Prentice-Hall, Inc., 1987. [4] Ljung L and Soderstrom T. Theory and Practice of Recursive Identi cation. MIT Press, Cambridge, Massachusetts, 1983. [5] Haykin S. Adaptive Filter Theory. Prentice Hall, 2nd edition, 1991. ISBN 0-13-012236-5. [6] Golub G H and Van Loan C F. Matrix Computations. John Hopkins University Press, 1989. [7] Kawabe S Chao J and Tsujii S. A New IIR Adaptive Echo Canceller : GIVE. IEEE Transactions on Selected Areas in Communications, 12(9):1530 { 1539, December 1994. [8] Widrow B and Stearns S D. Adaptive Signal Processing. Prentice Hall, Englewood Cli s, 1985.
