非線性回音消除之收斂性分析

全文

(1)非線性回音消除之收斂性分析. 學生：江冠諭. 指導教授：謝世福. 國立交通大學電信工程學系碩士班摘要. 為了補償非線性失真在免持聽筒或是視訊會議系統，回音的消除，通常使用一個無記憶性的多項式結合 NLMS 演算法的適應性濾波器。在傳統上，多項式採用次方級數展開的形式；而在本篇論文中，為了提昇非線性適應性濾波器的收斂速度我們使用正交多項式的非線性適應性濾波器。不論是級數型多項式或是正交多項式，我們都分析出它們的收斂理論值，並且從電腦模擬結果得知和我們的理論值是符合的，而模擬的結果也說明了我們提出的方式的確有較好的收斂性。除了使用適應性的方式之外，在實際的語音傳送之前，我們使用了訓練序列信號的方式來評估非線性濾波器的係數。我們分析它們的收斂理論值，並且由電腦模擬得到驗證，而且在訊雜比不佳的情況下，它的收斂值會比適應性濾波器來得佳。. i.

(2) Convergence Analyses of Nonlinear Acoustic Echo Cancellation Student : G. Y. Jiang. Advisor : S. F. Hsieh. Department of Communication Engineering National Chiao Tung University. Abstract In order to compensate the nonlinear distortion in the hands -free telephones or teleconferencing system, a memoryless polynomials NLMS adaptive filter can be used to cancel nonlinear acoustic echo. Conventional polynomials model employs a power-series expansion. In this thesis we propose an orthogonal polynomials adaptive filter and perform theoretical convergence analysis of residual echo power which proves its faster convergence rate owing to the reduced eigen spread of the input signal. Computer simulations justify our analysis and show the improved performance of the proposed nonlinear acoustic echo canceller. In addition to the adaptive method, the training sequence (TS) can be used to estimate the coefficients of the nonlinear acoustic echo cancellation. The convergence rate of the training method is derived analytically. Computer simulations show that the TS method performs better than the NLMS method at low SNR.. ii.

(3) Acknowledgments I am deeply grateful to my advisor Dr. S. F. Hsieh, for his support, patience, and encouragement throughout my master studies. My thanks also go to the members of my lab they give me many helpful comments. Last but by no means least; I would like to thank my parents, my family, and all of people who have helped me in the past two years.. iii.

(4) Contents 中文摘要……………………………………………………….…………………….i English Abstract……………………………………………………………………ii Acknowledgement………………………………………………..……….………iii Contents……………………………………………………….…………………….iv List of Figures……………………………………………….…………………….vii List of Tables……………………………………………….………………………ix 1. Introduction……………………………………………….……………………1. 2. Nonlinear Adaptive AEC……………………………….……………………7 2.1 Adaptive nonlinear NLMS AEC……………………………….…………7 2.2 Convergence analysis of nonlinear NLMS algorithm…………….……….11 2.2.1 Variance of nonlinear coefficients error………………………..……11 2.2.2 Residual echo power analysis……………………………….….……14 2.3 Convergence analysis of linear NLMS algorithm…………………………15 2.3.1 Variance of linear coefficients error…………………………...…….15 2.3.2 Residual echo power analysis……………………………….…….…16 2.4 Nonlinear processor with orthogonal polynomials………………………..17 2.4.1 Signal dependent orthogonal polynomials…………………………17 2.4.2 Eigenvalue spread analysis…………………………………………21 2.5 Application of dual loudspeakers system…………………………………28 2.6 Summary…………………………………………………………………..31. 3. Training Sequence and Coefficients Estimation………………………32 3.1 Correlation based nonlinear AEC………………………………………….33 3.2 Training sequence estimation algorithm…………………………………..35 iv.

(5) 3.2.1 Linear TS coefficients estimation algorithm………………………...35 3.2.2 Nonlinear TS coefficients estimation algorithm……………………38 3.3 Recursive training sequence algorithm……………………………………41 3.4 Linear convergence analysis………………………………………………44 3.4.1 Variance of linear coefficients error…………………………………44 3.4.2 Residual echo power analysis………………………………………..47 3.5 Nonlinear convergence analysis…………………………………………...48 3.5.1 Variance of nonlinear coefficients error……………………………48 3.5.2 Residual echo power analysis………………………………………49 3.6 TS for dual loudspeakers system…………………………………………..50 3.6 Summary…………………………………………………………………..51. 4. Computer Simulations………………………………………………………52 4.1 Parameters of simulations…………………………………………………52 4.2 Coefficients estimation based on NLMS adaptive algorithm……………57 4.2.1 Individual coefficients and residual echo power convergence………57 4.2.2 Convergence rate using orthogonal bases…………………………60 4.2.3 Joint convergence rate using orthogonal bases……………………62 4.2.4 Convergence rate of a speech input………………………………….65 4.2.5 Convergence rate of dual loudspeakers system……………………...66 4.3 TS-based coefficients estimation…………………..………………………70 4.3.1 Convergences of linear TS coefficients and residual TS echo power 70 4.3.2 Joint TS convergence rate…………………………………………73 4.3.3 Joint residual echo comparison of TS and NLMS algorithms………74 4.3.4 Training Sequences for dual loudspeakers system…………………..77 4.4 Experiments with a real echo path………………………………………..78. 5. Conclusions…………………………………………………………………84 v.

(6) Appendix……………………………………………………………………………85 Bibliography………………………………………………………………………90. vi.

(7) List of Figures 1.1 The simplified diagram of hands-free telephone system…………………………1 1.2 Nonlinear acoustic echo canceller………………………………………………2 2.1 Polynomial nonlinear acoustic echo canceller…………………………………7 2.2 Nonlinear processor……………………………………………………………17 2.3 Nonlinear AEC structure for dual loudspeakers system………………………28 3.1 Structure of nonlinear AEC based on training sequence………………………..34 4.1 Room impulse response of a typical office……………………………………..54 4.2 Room impulse response from the left loudspeaker……………………………54 4.3 Room impulse response from the right loudspeaker……………………………55 4.4 I/O mapping characteristic of the left nonlinear loudspeaker…………………55 4.5 I/O mapping characteristic of the right nonlinear loudspeaker…………………56 4.6 The speech signal………………………………………………………………56 4.7 Comparison of nonlinear coefficients misalignments (perfect h).……………58 4.8 Comparison of residual echo power (perfect h).………………………………..58 4.9 Comparison of linear coefficients misalignments (perfect a).………………….59 4.10 Comparison of residual echo power (perfect a).………………………………..59 4.11 Residual echo powers for uniform input (perfect h)…..………………………..60 4.12 Residual echo powers for WGN input (perfect h)...……………………………61 4.13 Residual echo powers for Laplacian input (perfect h).…………………………61 4.14 Joint-residual echo power for uniform input……………………………………63 4.15 Joint-residual echo power for Gaussian input…………………………………..64 4.16 Joint-residual echo power for Laplacian input………………………………….64 4.17 ERLE for a true speech input signal with perfect linear coeff. ………………65 4.18 Histogram of input speech………………………………………………………66 vii.

(8) 4.19 Nonlinear coefficients misalignments for dual loudspeakers…………………...68 4.20 Residual echo power for dual loudspeakers (perfect h’s) ……………………...68 4.21 Joint residual echo power for dual loudspeakers………………………………..69 4.22 Comparison of linear TS coeff. misal. (SNR=5 dB, perfect nonlinear coeff.)….70 4.23 Comparison of linear residual TS echo power (SNR=5 dB).…………………...71 4.24 Comparison of nonlinear TS coeff. misal. (SNR=5 dB, perfect linear coeff.)….72 4.25 Comparison of residual nonlinear TS echo power (SNR=5 dB)..………………72 4.26 Joint linear coefficients misalignment for training method (SNR=5 dB)………73 4.27 Joint residual echo power for training method (SNR=5 dB).…………………..74 4.28 Joint residual echo power comparison of NLMS and TS (SNR=10 dB)……….75 4.29 Joint residual echo power comparison of NLMS and TS (SNR=5 dB)………...76 4.30 Joint residual echo power comparison of NLMS and TS (SNR=0 dB)………...76 4.31 Comparison of residual echo powers for dual loudspeakers and training method (SNR=5 dB) …………………………………………………………………….77 4.32 Comparison of joint residual echo powers for dual loudspeakers for TS and NLMS (SNR=5 dB) ……………………………………………………………78 4.33 The speech input and its microphone output signal……………………………79 4.35 ERLE comparison between linear and nonlinear AEC for a true echo path……80 4.35 Average ERLE comparison of different nonlinear orders………………………82 4.36 Average ERLE comparison at different SNRs………………………………….83. viii.

(9) List of Tables 2.1 Comparison of computational cost, no. multiplication per iteration……………10 2.2 Signal-dependent orthogonal polynomials……………………………………...18 2.3 Orthogonal polynomials for a uniformly distributed signal…………………….19 2.4 Orthogonal polynomials for a Gaussian-distributed input……………………...20 2.5 Orthogonal polynomials for a Laplacian-distributed input……………………..21 2.6 Eigenvalue spread comparison…………………………….……………………27 3.1 Comparison of complexity computational……………………………………...43 4.1 Average joint-ERLE(dB) comparison for a speech input signal………………66 4.2 Average ERLE comparison between different linear filter lengths for true echo path……………………………………………………………………………81 4.3 Nonlinear coefficients for each AEC with different linear orders………………81 4.4 Average ERLE comparison of NLMS and TS methods for a true echo path…...82. ix.

(10) Chapter 1 Introduction Hands-free telephone or teleconferencing usually suffers from the annoying acoustic echo problem [1], which is the far end speech transmitted back to the far end user as a result of the coupling of the loudspeaker and microphone at the near end. A simplified diagram of hands-free telephone system is shown in Fig. 1.1. The main object of acoustic echo cancellation (AEC) is to estimate the unknown echo path and subtract the estimated echo components from the microphone output. Since the echo path may be time-variant due to objects moving around the room, an adaptive filter is commonly used for tracking the echo path to provide satisfactory speech communication quality [1]-[6]. Far-end room. Near-end room. Microphone. Loudspeaker. Far-end signal. echo AEC. Near-end Microphone. Loudspeaker. signal. Fig. 1.1 The simplified diagram of hands-free telephone system. The performance of AEC relies on its tracking capability. There are many recursive algorithms that have been proposed [6]. The least-mean-square (LMS) 1.

(11) algorithm is famous for its low computational cost and the recursive least-squares (RLS) algorithm has its advantage of fast convergence rate. Beside the tracking capability, the performance of AEC is restricted by noise, finite precision, truncation effects, under-model of the echo path and loudspeaker nonlinearities [7]. In this thesis, the issue of nonlinearities is our main focus. In [8], the nonlinearities of loudspeaker can be classified as nonlinearities with and without memory. Dynamic loudspeakers cause nonlinearity with memory when the power amplifier is not driven into saturation. The other is nonlinearity without memory when the power amplifier is overdriven. To compensate these kinds of nonlinearities distortions, several nonlinear AECs have been proposed recently [9]-[20]. The nonlinear AEC system is shown in Fig.1.2. The signal from the far end is passing through the nonlinear loudspeaker and the room impulse response and then is picked up by the microphone. The nonlinear AEC is supposed to cancel the echo signal. The echo can be cancelled perfectly if the nonlinear AEC filter is identical to the nonlinear loudspeaker and room impulse response. Different nonlinear structure has its own computation complexity, convergence speed and robustness. Nonlinear. Far end signal. Loudspeaker. Nonlinear AEC. Room. Echo Near end signal. Residual echo. Microphone. Fig.1.2 Nonlinear acoustic echo canceller 2.

(12) We summary some important nonlinear adaptive structures. z. Volterra model, which is attractive because it is a straightforward generalization of the linear system description and the behavior of many physical systems can be described with a Volterra filter. The Volterra series based filters have been proposed [10]-[12] etc. for line echo canceling. It can represent a large of class of nonlinear system. However, due to their high computational complexity they are limited use in practical systems.. z. Bilinear model, which is a parametric model that contains cross terms; it corresponds to a subclass of the NARMAX structure [13]. The NARMAX structure is a general parametric model but needs a pre-identification procedure.. z. Cascade model, which can be considered as a particular subclass of a Volterra series filter and its main advantage is to introduce fewer parameters for estimation. The neural network [14] with a cascade structure offers a new perspective but needs an extra reference microphone. Others cascade models include: . Hammerstein model [15], cascade of a memoryless polynomials filter and a FIR filter.. . Wiener model [4], cascade of a FIR filter and a memoryless polynomials filter.. . Wiener-Hammerstein model [26], cascade of a FIR filter, a memoryless polynomials filter and a FIR filter.. Among these nonlinear structures, the Hammerstein model will be used in this thesis for following reasons. First, the nonlinearity with memory only occurs in application of high quality loudspeakers [8]; so that for hands-free or power limited low cost application, the compensation for nonlinearity with memory is not necessary. 3.

(13) Second, considering that the far end signal passes through the loudspeaker first and then the room impulse response, the cascade model of a nonlinear processor and a FIR filter is a natural choice. Third, the Hammerstein model is widely used in nonlinear system identification such as nonlinear AEC, neural networks [15-17] etc. Its joint NLMS-type adaptation algorithm is well known [18]. Besides the polynomial function, a sigmoid function [6, 9] can also be used to model the nonlinear saturation. Similarly, a raised cosine function can also be used for nonlinear compensation [19]. The objective in [19] is to achieve a low computational complexity in implementing a nonlinear AEC. Both these two nonlinear models can use the NLMS adaptive algorithm to update its coefficients, except that its nonlinear component is generated by an exponential function. It only uses one parameter to control the nonlinearity therefore it has less freedom. By contrast, the polynomial type has more freedom in that each nonlinear order can be controlled with one coefficient. The power series polynomials are simple to implement but its high correlation among different polynomials orders leads to low convergence rate. To overcome this problem, recently some orthogonal structures have been developed. In [10], Mathews suggested to perform an orthogonalization procedure on the nonlinear bases outputs when the input signal is Gaussian distributed. Jenkins et al. in [25] proposed an orthogonal basis to represent the Volterra series thus orthogonalization procedure is not required but the input signal is also assumed to be have a Gaussian distribution and unity variance. Similarly, in [26] the Wiener-Hammerstein model is used and its nonlinearities is assumed to be expandable in a series of fixed orthogonal Hermite polynomials. The Hermite polynomials are a set of orthogonal polynomials on the 2. infinite interval with respect to the e x weight function therefore the works in [25, 26] are limited to the unit-variance Gaussian input signal only. In [20] Kuech et al. proposed an adaptive orthogonalized power filter to improve the convergence rate for 4.

(14) the input signal with any distribution, stationary or time-variant. The orthogonal basis is updated online at each iteration and the Gram Schmidt procedure is employed to find out the orthogonalization coefficients, as a result, computational complexity is increased. In this thesis we use fixed orthogonal polynomials to produce the nonlinear components. Its low computational complexity and fast convergence rate makes the orthogonal polynomials filter very promising for nonlinear AEC. For practical applications, the influence of the probability distribution of the input signal is insignificant because the performance of the fixed orthogonal polynomials remains relatively well. Unlike the previous researches, we perform the convergence analysis of the nonlinear AEC, from which we examine the effects of orthogonality due to various input probability distributions and conclude the superiority of the orthogonal polynomials. In addition to the orthogonal polynomial basis, sending a white sequence to train the coefficients of nonlinear AEC in advance of speech communication can also be used for speeding up the convergence rate. The training sequence, used in channel estimation, adaptive equalizer applications or echo path estimation, have been well studied in [21]-[24]. During the training mode, we can fast start up the adaptive filter, especially in noisy environments. The training sequence is generated by a training sequence generator. The estimation is done with the correlation method, where a portion of the training sequence is correlated with shifted versions of the received signal. Based on the difference between this known sequence and the received sequence, the coefficients of the unknown can be determined. Although the training method is to solve the Wiener-Hopf solution directly and often involves a matrix inverse, the solution is simple because the matrix is only a function of the known training sequence, and a pre-computed inverse of the matrix can be stored. For some special training sequences, an inverse matrix is not even required. 5.

(15) The other chapters of the thesis are organized as follows. z. Chapter 2, we perform convergence analyses of coefficients and residual echo power and we compare the convergence rates of orthogonal and non-orthogonal basis. Eigenvalue spread analysis is introduced for better illustration.. z. Chapter 3, we perform the linear coefficients convergence analysis based on the training sequence method. The recursive analytical form will also be introduced. We also show the analysis of dual loudspeakers system.. z. Chapter 4, we include many computer simulations that have been developed to illustrate the analyses in chapter 2 and 3. These simulations help us to compare the performances of different nonlinear AEC structure. Finally, for practical use, a true echo path experiment is performed.. z. Chapter 5, we give a conclusion of our work.. The main efforts in this thesis are: (1) For NLMS nonlinear algorithm, we derive individual convergence analyses of the linear and nonlinear coefficients and its residual echo power. (2) An orthogonal basis for Gaussian and nonGaussian input signal is used to accelerate the nonlinear coefficients convergence rate. (3) A training sequence algorithm is proposed for a nonlinear Hammerstein model. (4) Convergence analyses of linear coefficients based on training sequence algorithm is derived.. 6.

(16) Chapter 2 Nonlinear Adaptive filter In Chapter 2 we will introduce the nonlinear AEC with cascade nonlinear processor and describe the joint NLMS-type adaptation algorithm to adaptive the both linear and nonlinear coefficients. In Section 2.2 we will analyses the nonlinear coefficients and residual echo power convergence rates under the assumption of linear coefficients have perfectly known. In Section 2.3 we will analyses the linear coefficients and its residual echo power when nonlinear coefficients have perfectly known. In Section 2.4, we discuss the eigenvalue spread of signal dependent orthogonal bases. Finally, in Section 2.5, we extend the analysis for dual loudspeakers system.. 2.1 Adaptive nonlinear NLMS AEC Nonlinear Loudspeaker. Far end signal x[n] (⋅) 2 ...... (⋅) N. (⋅)1 p1[n]. p2 [n]. p N [ n]. Room. a 1[ n] a 2 [ n] a N [ n]. h[n]. n] s[ NLMS algorithm. [n] FIR h n] y[. Residual echo e[n]. Echo. d [ n]. Near end signal v[n]. Microphone. Fig.2.1 Polynomial nonlinear acoustic echo canceller. 7.

(17) As shown in Fig 2.1, the signal x[n] from the far end is assumed to be nonlinearly distorted only in the power amplifier of loudspeaker. It is then passing through a room impulse response h[n]. Hence, the nonlinear processor is modeled as the loudspeaker and linear filter is modeled as a room impulse. The cascade filter structure is the same as the loudspeaker and the room impulse response. Let d [n] denote the desired signal. The nonlinear AEC output signal y[n] can be written as. y[n] = h T [n]s[n] h[n] = [h 0 [n], h 1[n],..., h M -1[n]]T. where h[n] represents the estimated coefficients vector of the linear FIR filter, M denotes the length of the FIR filter. s[n] is the output vector of the nonlinear filter s[n] = [ s[n], s[n -1],..., s[n - M + 1]]T .. For each s[n] is given by s[n] = [ x [n] x 2 [n]" x N [n]][al1[n] al2 [n]" amN [n]]T = x[n]T a [n] therefore, s[n] is given by s[n] = [ x[n]T a [n], x[n -1]T a [n -1],..., x[n - M + 1]T a [ n - M + 1] ]T .. pi is the polynomial basis of order i , for example p1[n] = x[n] and p2 [n] = x 2 [n] in case of a power series expansion basis. N is the order of the polynomials, and a [n] is the estimated coefficients vector of the nonlinear processor. The estimated error is e[n] = d [n] - y[n] T = d [n] - h [n]s[n]. The gradient of the error power e 2 [n] , as derived for linear transversal filter in [18]. 8.

(18) can be calculated according to:. ∂e2 [n] ∇h = = -2e[n]s[n] ∂h[n] ∇a =. ∂e 2 [n] = -2e[n]PT [n]h[n] ∂a[n]. where P[n] is nonlinear expanded matrix is defined by p2 [ n ] " pN [ n] ⎡ p1[n] ⎤ ⎢ ⎥ ⎢ p1[n -1] p2 [n -1] " pN [n -1] ⎥ ⎢# ⎥ # % # ⎢ ⎥ ⎣ p1[n - M + 1] p2 [n - M + 1] " pN [n - M + 1]⎦. The definition of error signal is different from [18], here, we only calculate the scalar s[n] of the vector s[n] at each iteration. If the coefficients vectors are updated with. step size μh and μa , a joint NLMS-type adaptive algorithm is given by h[n + 1] = h[n] +. a [n + 1] = a [n] +. μh s[n]e[n] 2. s[n]. (2.1.1). 2. μa 2. P [n]h[n] + δ T. PT [n]h[n]e[n]. (2.1.2). 2. At each iteration, the echo signal e[n] is the same for coefficients update in both (2.1.1) and (2.1.2). For computational complexity, we examine the number of multiplications required to make one complete iteration of the algorithm (2.1.1) and (2.1.2). s[n] in (2.1.1) and its 2-norm need N and M multiplications respectively thus the total requirement of (2.1.1) is about 2M + N. For (2.1.2), PT [n]h[n] and its 2-norm need MN and N multiplications respectively thus the total requirement of (2.1.2) is about MN+2N. It needs MN+2M+3N multiplications to find out the coefficients at each iteration. In addition to the adaptation algorithm, the nonlinear components need N-1 9.

(19) multiplications for power-series basis, x2[n] =x[n]x[n], x3[n] =x2[n]x[n],. …. , xN[n]. =xN-1[n]x[n]. When an orthogonal basis is used it needs more multiplications to produce the nonlinear components. According to Gram-Schmidt orthogonalization procedure it needs. 1 2. ( N 2 +N -2) . In Table 2.1, we show the total requirement of. multiplication at each iteration when Gram-Schmidt is used or not. The computational complexity only increases slightly when an orthogonal basis is used.. Table 2.1 Comparison of computational cost, no. multiplication per iteration Number of multiplication Without Gram-Schmidt procedure. MN+2M+4N-1. With Gram-Schmidt procedure. MN+2M+3.5N-1+0.5 N2. In the following sections, we assume that the nonlinear loudspeaker and room impulse response are time invariant. The near end signal v[n] only contains a white Gaussian noise (WGN) and double talk is not present.. 10.

(20) 2.2 Convergence analysis of nonlinear NLMS algorithm 2.2.1 Variance of nonlinear coefficients error. In this section we will derive the convergence rate of nonlinear coefficients under the assumption of perfect linear coefficients i.e., h[n] = h . The following analysis is similarly to [6]. First, the estimation error produced by the nonlinear AEC filter is expressed as T e[n] = d [n] - a [n]PT [n]h[n]. = aT PT [n]h + v[n] -[a + ε a [n]]T PT [n][h + ε h [n]] .. (2.2.1). Because the cascade structure we have more estimated error terms than [6], the joint error term produced by linear and nonlinear is difficult to perform its convergence analysis. For this reason, we assume that the linear coefficients are perfectly known, ε h [n] is equal to zeros. We express the estimated error as follows. e[n] = v[n] + εTa [n]PT [n]h .. (2.2.2). In (2.2.2) PT [n]h contains not only linear but also nonlinear order of input signal, this is different from [6] but the analyses procedures in [6] can still be used for here. We denote the nonlinear coefficients weight error by ε a [n + 1] = a-a [n + 1] 2. Using (2.1.2), (2.2.2) and let T denote PT [n]h + δ , we may rewrite ε a [n + 1] as 2. μ ε a [n + 1] = a-a [n]- a PT [n]h[n]e[n] T = ε a [ n] -. μa T. PT [n]h[n](v[n] + εTa [n]PT [n]h). μ ⎡ μ ⎤ = ⎢ I - a PT [n]hhT P[n]⎥ ε a [n] - a PT [n]hv[n] . T ⎣ T ⎦ According to the direct averaging method [6], when μa 1 , ε a [n + 1] can be approximated as follows: 11.

(21) ⎡ μ ⎤ ε a [n + 1] ≈ ⎢ I - a R PT [ n]h ⎥ ε a [n]+fa [n] ⎣ T ⎦. (2.2.3). where f a [ n] = -. μa. PT [n]hv[n]. T. and R PT [ n ]h is the correlation matrix of PT [n]h . By applying the unitary similarity transformation, R PT [ n ]h is transformed into a simpler form: QTa R PT [ n ]hQ a = Da where Q a is a unitary matrix and Da is a diagonal matrix consisting of the eigenvalues λai . Let K a [n] = QTa ε a [n] then we may transform (2.2.3) into the form ⎡ μ ⎤ QTa ε a [n + 1] = QTa ⎢ I - a R PT [ n ]h ⎥ ε a [n] + QTa fa [n] ⎣ T ⎦. ⎡ μ ⎤ K a [n + 1] = ⎢ I - a Da ⎥ K a [n] + Φ a [n] ⎣ T ⎦ where Φ a [n] = QTa f a [n] . The natural mode kai [n] , i-th entry of K a [n] , is stochastic with a mean and mean square value of its own. Let kai [0] denote the initial value of kai [n] and Φ ai [n] is denoted i-th entry of Φ a [n] . We may rewrite kai [n] as follows.. μa. kai [n] = (1-. T. μa. = (1-. T. λai )kai [n -1] + Φ ai [n -1] n -1. μa. j =0. T. λai ) n kai [0] + ∑ [1-. λai ]n-1- j Φ ai [ j ]. Hence, the first moment of kai [n] is given by E[kai [ n]] = (1-. = (1-. μa T. μa T. n -1. μa. j =0. T. λai ) n E[kai [0]] + ∑ [1-. λai ]n-1- j E[Φ ai [ j ]]. λai ) n kai [0]. where E[Φ a [n]] = QTa E[f a [n]] =-. μa. T = 0.. QTa hE[P[n]v[n]]. Since P[n] only contains the input signal and is independent of noise, we further 12.

(22) assume that the initial value of kai [0] is independent of Φ ai , therefore the second moment of kai [n] is given by. μa. 2. E[ kai [ n] ] = (1-. T. λai ) 2 n kai [0]. n -1 n -1. μa. g =0 j =0. T. + ∑∑ (1-. The kai [0]. 2. 2. μa. λai ) n-1- g (1-. T. λai ) n-1- j E[Φ ai [ g ]Φ ai [ j ]] .. (2.2.4). in the right-hand side of (2.2.4) is equal to the nonlinear coefficients. vector square 2-norm.The second term in right-hand side of (2.2.4) is zero when summation index g is not equal to j otherwise we can express E[Φ a [n]ΦTa [n]] as E[Φ a [n]Φ aT [n]] = ( =( =(. μa T. μa. ) 2 QTa hE[P[n]v[n]v[n]PT [n]]hT Q a ) 2 σ v2QTa R PT [ n ]hQ a. T. μa. ) 2 σ v2 Da .. T. (2.2.5). From (2.2.5), (2.2.4) can be written as 2. E[ kai [n] ] = (1-. μa T. μa =. T 2−. 2. λai ) 2 n a 2 + σ v2 (. σ v2. μa T. λai. μa T. μa 2. +( a 2 -. T 2−. ) 2 λai (1-. σ v2. μa T. λai. )(1-. μa T. μa T. n -1. μa. j =0. T. λai ) 2 n-2 ∑ (1-. λai ) 2 n .. λai )-2 j. (2.2.6). In (2.2.6) the error variance of nonlinear coefficients is given. Again, that (2.2.6) is under the assumption of perfect linear coefficients. The error variance of nonlinear coefficients can be determined with the knowledge of step size μa , noise power σ v2 , square 2-norm of nonlinear coefficients vector, eigenvalues λai of correlation matrix of s[n] and the sum of all eigenvalues. Because the step size and the eigenvalues are both positive, the second term of (2.2.6) will disappear when the iteration number 2. approaches to infinity. The steady state of E[ kai [n] ] is given in the first term of (2.2.6). 13.

(23) 2.2.2 Residual echo power analysis From (2.2.2), the mean square error (i.e., residual echo) due to an estimated error of nonlinear coefficients is given by 2 J a [n] = E ⎡ e[n] ⎤ ⎣ ⎦. = E[(v[n] + εTa [n]PT [n]h)(v[n] + hT P[n]ε a [n])] = σ v2 + E ⎡⎣εTa [n]PT [n]hhT P[n]ε a [n]⎤⎦ .. (2.2.7). Assume the variation of ε a [n] is slow compared with PT [n]h , hence. E ⎡⎣εTa [n]PT [n]hhT P[n]εa [n]⎤⎦ ≈ E ⎣⎡εTa [n]E ⎣⎡ PT [n]hhT P[n]⎦⎤ ε a [n]⎤⎦ = E ⎡εTa [n]R PT [ n ]h ε a [n]⎤ ⎣ ⎦. = E[K aT [n]QT R PT [ n ]h QK a [n]] = E[K Ta [n]Da K a [n]] N. 2 = ∑ λai E ⎡ kai [n] ⎤ . ⎣ ⎦ i =1. (2.2.8). From (2.2.6) and (2.2.8), the mean square error can be written as. 2 v. N. λai μa. μa. σ v2. μ + ∑ λai [ a 2 - T ](1- a λai )2 n . (2.2.9) J a [n] = σ + σ ∑ μ T T i =1 2 - λai i =1 2 - a λai T T 2 v. μa. N. 2. From (2.2.8) and (2.2.9), when the nonlinear coefficients error variance has been known the residual echo power can also be obtained. The nonlinear convergence rate depends on the values λai of PT [n]h , we may change the basis of P to have a faster convergence rate.. 14.

(24) 2.3 Convergence analysis of linear NLMS algorithm 2.3.1 Variance of linear coefficients error In this section we assume that the nonlinear coefficients are perfectly known, i.e. a [n] = a or ε a [n] = 0 . The estimation error produced by the nonlinear AEC filter is. similar to (2.2.1) but here ε a [n] = 0 and ε h [n] ≠ 0 T e[n] = d [n] - a [n]PT [n]h[n]. = aT PT [n]h + v[n] -[a + ε a [n]]T PT [n][h + ε h [n]] = v[n] + εTh [n]s[n]. (2.3.1). where s[n] = P[n]a . We denote the linear coefficients weight error by ε h [n + 1] = h-h [n + 1] . The following analysis procedures are similar to Section 2.2 thus we will omit the details. Using (2.1.1) and (2.3.1), we may rewrite ε h [n + 1] as ⎡ ⎤ μh ε h [n + 1] ≈ ⎢ I R s[ n ] ⎥ ε h [n]+fh [n] 2 ⎢⎣ s[n] 2 ⎥⎦. (2.3.2). where f h [n]= -. μh 2. s[n] 2. s[n]v[n] .. R s[ n ] is the correlation matrix of s[n] . The convergence of linear coefficients error variance can be obtained by using the same procedures as before in Section 2.2.1, the linear coefficients error variance is given by. μh 2. 2. E[ khi [n] ] =. s[n] 2 2−. μh. σ v2. 2. s[n] 2. 2. μh. 2. s[n] 2. λsi. +( h 2 -. 2−. 15. σ v2. μh 2. s[n] 2. λsi. )(1-. μh 2. s[n] 2. λsi ) 2 n. (2.3.3).

(25) where QTs R s Q s = Ds , Q s is a unitary matrix and Ds is a diagonal matrix consisting of the eigenvalues λsi , K h [n] = QTs ε h [n] and khi [n] is i-th entry of K h [n] . In (2.3.3) the error variance of linear coefficients is given under assumption of perfect nonlinear coefficients. It can be determined with the knowledge of step size μh , noise power σ v2 , square 2- norm of room impulse response, eigenvalues of correlation matrix of s[n] and the sum of all eigenvalues. Similarly, the steady state of 2. E[ khi [n] ] is given in the first term of (2.3.3). 2.3.2 Residual echo power analysis From (2.3.1), the mean square error (i.e., residual echo) due to estimated error of linear coefficients is given by 2 J h [n] = E ⎡ e[n] ⎤ ⎣ ⎦. = σ v2 + E ⎡⎣εTh [n]s[n]s[n]T ε h [n]⎤⎦ .. (2.3.4). Assume the variation of ε h [n] is slow compared with s[n] , hence the residual echo power can be obtained similarly by J h [n] = σ v2 +. μh. s[n] 2. i =1. 2-. μh M. + ∑ λsi [ h 2 i =1. λsi μh. M. σ v2 ∑ 2. 2. 2. s[n] 2 2−. 2. s[n] 2. λsi. σ v2. μh. 2. s[n] 2. λsi. ](1-. μh 2. s[n] 2. λsi ) 2 n. (2.3.5). Unlike Section 2.2, the nonlinear basis has no effect upon the eigenvalue in (2.3.5) when nonlinear coefficients are perfectly known. Therefore, in next section, we discuss the eigenvalues and nonlinear basis relationship only when linear coefficients are error free and nonlinear coefficients are not known.. 16.

(26) 2.4 Nonlinear processor with orthogonal polynomials 2.4.1 Signal dependent orthogonal polynomials In this section, we discuss the orthogonal and non-orthogonal basis of nonlinear processor. The nonlinear processor is shown in Fig. 2.2. To accelerate the convergence rate, we define a set of orthogonal sets { pi [n], i = 1, 2, ...} if their outputs are uncorrelated with each other E[ pm [ n] pn [n]] ∑ pm [n] pn [ x] f x [ x] x. = q δ[m - n]. (2.4.1). where f x [ x] is the probability density function (pdf) of the input signal x[n]. If q is equal to 1, then the polynomials are not only orthogonal but also orthonormal. Far end signal x[n] (⋅) 2 ...... (⋅) N. (⋅)1 p1[n]. p2 [n]. p N [ n]. a 1[ n] a 2 [ n] a N [ n] s[n]. Fig.2.2 Nonlinear processor The non-orthogonal polynomials set can be systematically modified to yield an orthogonal set by using Gram Schmidt orthonormalization in any interval with the weighting function f x [ x] . The weighting function depends on the probability distribution of the input signal, therefore we derive the general form first and then apply it to uniform, WGN, and Laplacian pdf’s. The orthogonal bases are described as follow p1[n] = x[n] j −1. p j [n] = x j [n] + ∑ C j ,i xi [n] 1< j ≤ N. i =1. 17.

(27) The orthogonalization coefficients C j ,i are chosen such that (2.4.1) is satisfactory, i.e., the coefficients C j ,i for the j-th order polynomial is obtained by solving. ⎡ m j +1 ⎤ ⎡ m2 " m j ⎤ ⎡ C j ,1 ⎤ ⎢ ⎥ ⎢ ⎥⎢ ⎥ ⎢ # ⎥ = -⎢ # % ⎥⎢ # ⎥ # ⎢ ⎥ ⎢ ⎥⎢ ⎥ ⎢m ⎥ ⎢ m " m ⎥ ⎢C ⎥ 2 j−2 ⎦⎥ ⎣⎢ j , j−1 ⎦⎥ ⎣⎢ 2 j−1 ⎦⎥ ⎣⎢ j where mi = E[ x i [ n]] is i-th moment of the input signal x[n] . We can find out the orthogonalization coefficients when the input signal pdf is known a priori. The orthogonalization coefficients are constant for stationary input otherwise the coefficients are time dependent. Here, we assume that the input signal is stationary and its pdf is pre-known therefore the ensemble average is performed to find out the orthogonalization coefficients. The moment estimation is performed when the input signal pdf is not known, in such case, the orthogonalization coefficients are time dependent as described in [12]. After the orthogonalization procedures, the first six orders of orthogonal polynomials are provided in Table 2.2. Here we have assumed pdf’s are symmetric with zero mean and all the odd-ordered moments are zero, mi = 0 , for i: odd.. Table 2.2 Signal-dependent orthogonal polynomials Polynomials. Coefficients. p1[n] = x[n] p2 [ x] = x 2 [n] - C2.1. C2.1 = m2. p3 [ x] = x 3 [n] - C3.1 x[n]. C3.1 =. m4 m2. p4 [ x] = x 4 [n] - C4,1 x 2 [n] - C4,2. C4,1 =. m6 - m2 m4 m4 − m22. C4,2 = -C4,1m2 + m4 18.

(28) p5 [ x] = x5 [n] - C5,1 x3 [n] - C5,2 x[n]. C5,1 =. m8 - mm42 m6 2. m6 - 2 mm42 m4 + mm42. C5,2 = -C5,1 p6 [ x] = x 6 [n] - C6,1 x 4 [n] - C6,2 x 2 [n] - C6 C6,1 =. m4 m6 + m2 m2 m10 - C4,1m8 - C4,2 m6. 2 m8 + C m4 + C4,2 - 2(C4,1m6 + C4,2 m4 + C4,1C4,2 m2 ) 2 4,1. C6,2 = -C6,1C4,1 +. m8 - m2 m6 m4 - m22. C6,3 = -C6,1C4,2 -. m8 - m2 m6 m2 + m6 m4 - m22. Next, we discuss three kinds of signal distribution models, uniform, Gaussian and Laplacian. A. Uniform input For uniform signal in the interval [-1,1], the moments are 1 1 1 1 1 m2 = , m4 = , m6 = , m8 = , m10 = . 3 5 7 9 11 According to Table 2.2, the orthogonal polynomials for the uniformly distributed input are listed in Table 2.3.. Table 2.3 Orthogonal polynomials for a uniformly distributed signal Order. polynomials. Order. polynomials. 1. p1[n] = x[n]. 4. 6 3 p4 [ x] = x 4 [n] - x 2 [n] 7 35. 2. p2 [ x ] = x 2 [ n ] -. 5. p5 [ x] = x 5 [n] -1.11x3 [n] + 0.24 x[n]. 6. p6 [ x] = x 6 [n] - 0.07 x 4 [n] - 0.65 x 2 [n] + 0.09. 3. 1 3. 3 p3[ x] = x 3[n] - x[n] 5. In Table 2.3, these polynomials are not like Legendre polynomial form. Because we do not set pn [1] = 1, n = 1 ~ 6. 19.

(29) B. Gaussian input For WGN input distribution signal we let m2 be equal to 1/9=0.1111. The white Gaussian signal x[n] with unit variance σ x2 = 1 is generated first and it is widely distributed over [−∞, ∞] . In practice we let the signal within the interval ±3σ x2 be normalized into[−1,1] , the signal outside this interval can be ignored. Therefore, we use the normalized. 1 x[n] as the Gaussian noise for comparison with other pdf’s. Its 3. 2nd moment is m2 = 0.1111 , the other moments are expressed as follows. m2 = 0.1111, m4 = 3m22 , m6 = 15m23 , m8 = 105m24 , m10 = 945m25 . Using these moments and Table 2.2, the orthogonal polynomials are listed in Table 2.4.. Table 2.4 Orthogonal polynomials for a Gaussian-distributed input Order. polynomials. Order. polynomials. 1. p1[n] = x[n]. 4. p4 [ x] = x 4 [n] - 0.666 x 2 [n] + 0.030. 2. p2 [ x] = x 2 [n] - 0.111. 5. p5 [ x] = x 5 [n] -1.11x3[ n] + 0.185 x[n]. 3. p3 [ x] = x 3 [n] - 0.333x[n]. 6. p6 [ x] = x 6 [n] - 0.416 x 4 [ n] - 0.277 x 2 [n] + 0.025. C. Laplacian input For Laplacian input distribution signal, in practice, the same reason as WGN signal we let m2=0.1111 and the other moments are given by m2 = 0.1111, m4 = 6m22 , m6 = 90m23 , m8 = 2520m24 , m10 = 113400m25 . The orthogonal polynomials for Laplacian distribution signal are listed in Table 2.5.. 20.

(30) Table 2.5 Orthogonal polynomials for a Laplacian-distributed input Order. polynomials. Order. polynomials. 1. p1[n] = x[n]. 4. p4 [ x ] = x 4 [n] -1.8667 x 2 [n] + 0.1333. 2. p2 [ x] = x 2 [n] - 0.1111. 5. p5 [ x] = x 5 [n] - 4.074 x 3[n] + 1.605 x[n]. 3. p3 [ x] = x 3[n] - 0.6666 x[n]. 6. p6 [ x] = x 6 [n] - 4.45 x 4 [n] + 2.307 x 2 [n] + 0.050. 2.4.2 Eigenvalue spread analysis In Section 2.2 we show that the convergence rate depends on the normalized eigenvalues. Equivalently, the eigenvalue spread is commonly used for convergence analysis, instead of normalized eigenvalues [6]. Here, we will examine the eigenvalue spreads for different input pdf’s versus different orthogonal polynomials. For simplicity we normalize h. 2. to be 1 and nonlinear distortion only contains. odd order harmonics and even orders are excluded, in case of a symmetric input/output nonlinear characteristic curve.. A. Power-series nonorthogonal polynomials basis I. Uniformly distributed input The correlation matrix, R PT [ n ]h , is expressed as follows and the orthogonal polynomials can be found in Table 2.3. R PT [ n ]h = E ⎡⎣ PT [n]hhT P[n]⎤⎦. 21.

(31) ⎡ r11 = ⎢⎢ r31 ⎢⎣ r51. with. r13 r33 r53. r15 ⎤ r35 ⎥⎥ r55 ⎥⎦. (2.4.2). r11 = E ⎡⎣( p1[n]h0 + p1[n -1]h1 + " + p1[n - M + 1]hM -1 ) 2 ⎤⎦ 2. = m2 h 2 . r33 = E ⎡⎣( p3[n]h0 + p3 [n -1]h1 + " + p3[n - M + 1]hM -1 ) 2 ⎤⎦ 2. = m6 h 2 . r55 = E ⎡⎣( p5 [n]h0 + p5 [n -1]h1 + " + p5 [n - M + 1]hM -1 ) 2 ⎤⎦ 2. = m10 h 2 . Unlike the orthogonal basis, the correlation term of R PT [ n ]h is not equal to zero as a result of the basis is not orthogonal r13 = r31 = E[( p1[n]h0 + p1[n -1]h1 + " + p1[n - M + 1]hM -1 ) ( p3[n]h0 + p3[n -1]h1 + " + p3 [n - M + 1]hM -1 )] ⎡ M -1 M -1 ⎤ ⎡ M -1 ⎤ = E ⎢ ∑ p1[n - i ] p3 [n - i ]hi2 ⎥ + E ⎢ ∑∑ p1[ n - i ] p3 [ n - j ]hi h j ⎥ ⎣ i =0 ⎦ ⎣ i =0 j ≠i ⎦ 2. = m4 h 2 . r15 = r51 = E[( p1[n]h0 + p1[n -1]h1 + " + p1[n - M + 1]hM -1 ) ( p5 [n]h0 + p5 [n -1]h1 + " + p5 [n - M + 1]hM -1 )] 2. = m6 h 2 . r35 = r53 = E[( p3 [n]h0 + p3[n -1]h1 + " + p3 [n - M + 1]hM -1 ) 22.

(32) ( p5 [n]h0 + p5 [n -1]h1 + " + p5 [n - M + 1]hM -1 )] 2. = m8 h 2 . Therefore, the correlation matrix becomes. R PT [ n ] h. 0.2 0.1428⎤ ⎡ 0.3333 ⎢ 0.1428 0.1111⎥⎥ , = ⎢ 0.2 ⎢⎣0.1428 0.1111 0.0909 ⎥⎦. which is no more diagonal. Its three eigenvalues are 0.0006, 0.0334, 0.5331 and its eigenvalue spread is as large as 888.5. II. Gaussian input The non-diagonal correlation matrix is ⎡ 0.1111 0.0370 0.0206 ⎤ R PT [ n ]h = ⎢⎢0.0370 0.0206 0.0160 ⎥⎥ . ⎢⎣0.0206 0.0160 0.0160 ⎥⎦. Three eigenvalues are 0.0008, 0.0170, 0.1297, thus the large eigenvalue spread is 162.1.. III. Laplacian input The correlation matrix for the Laplacian input is ⎡ 0.1111 0.0741 0.1235⎤ R PT [ n ]h = ⎢⎢ 0.0741 0.2344 0.3841⎥⎥ . ⎢⎣0.1235 0.3841 1.9204 ⎥⎦. Three eigenvalues are 0.0177, 0.1285, 2.0088, and its eigenvalue spread is 113.5.. B. Uniform orthogonal polynomials basis I. Uniform input In the first discussion, we use uniformly orthogonal polynomials which are found under. uniform. input. distribution.. When 23. the. input. signal. is. uniformly.

(33) distributed(matched to the polynomial bases), the nonlinear components between different orders have perfect orthogonality. The correlation matrix, R PT [ n ]h is expressed the same as (2.4.2) with r11 = E ⎡⎣( p1[n]h0 + p1[n -1]h1 + " + p1[n - M + 1]hM -1 ) 2 ⎤⎦ ⎡ M -1 M -1 ⎤ ⎡ M -1 ⎤ = E ⎢ ∑ p12 [n - i ]hi2 ⎥ + E ⎢ ∑∑ p1[n - i ] p1[ n - j ]hi h j ⎥ . ⎣ i =0 ⎦ ⎣ i =0 j ≠ i ⎦. Because the Input signal is zero mean and uncorrelated between different time index, (2.4.2) can be expressed as 2. r11 = m2 h 2 . Similarly, r33 and r55 can be written as r33 = E ⎡⎣( p3[n]h0 + p3 [n -1]h1 + " + p3[n - M + 1]hM -1 ) 2 ⎤⎦ 6 9 2 = (m6 - m4 + m2 ) h 2 . 5 25 r55 = E ⎡⎣( p5 [n]h0 + p5 [n -1]h1 + " + p5 [n - M + 1]hM -1 ) 2 ⎤⎦ 2. = (m10 - 2.22m8 + 1.71m6 - 0.53m4 + 0.057m2 ) h 2 . The other terms of the correlation matrix are zeros as a result of orthogonal property. Therefore, we can write the correlation matrix as. R PT [ n ] h. 0 0 ⎤ ⎡0.3333 0.0029 0 ⎥⎥ , = ⎢⎢ 0 ⎢⎣ 0 0 0.0015⎥⎦. whose eigenvalue spread is 222.22.. II. Gaussian input The correlation matrix becomes. 24.

(34) ⎡ 0.1111 −0.0296 0.0061 ⎤ R PT [ n ]h = ⎢⎢ −0.0296 0.0161 −0.0016 ⎥⎥ ⎢⎣ 0.0061 −0.0016 0.0024 ⎥⎦. with an eigenvalue spread of 59.1. III. Laplacian input The correlation matrix becomes ⎡ 0.1111 0.0074 0.0679 ⎤ R PT [ n ]h = ⎢⎢0.0074 0.0746 0.2241⎥⎥ ⎢⎣0.0679 0.2241 1.2461 ⎥⎦. with an eigenvalue spread of 39.3.. C. Gaussian orthogonal polynomials basis I. Gaussian input For its matched input signal with Gaussian distribution, the correlation matrix can be found by the same procedure with the orthogonal polynomials given in Table 2.4. 0 0 ⎤ ⎡ 0.0020 0.0082 0 ⎥⎥ , R PT [ n ]h = ⎢⎢ 0 ⎢⎣ 0 0 0.1111⎥⎦. Three eigenvalues are 0.0020, 0.0082, 0.1111, thus the eigenvalue spread is 55.6. II. Uniformly distributed input The correlation matrix becomes ⎡ 0.3333 0.0889 −0.0175⎤ R PT [ n ]h = ⎢⎢ 0.0889 0.0466 −0.0046 ⎥⎥ ⎣⎢ −0.0175 −0.0046 0.0024 ⎦⎥. with an eigenvalue spread of 246. III. Laplacian input The correlation matrix becomes. 25.

(35) ⎡ 0.1111 0.0370 0.0618⎤ R PT [ n ]h = ⎢⎢0.0370 0.0864 0.2402 ⎥⎥ ⎢⎣ 0.0618 0.2402 1.2389 ⎥⎦. with an eigenvalue spread of 42.. D. Laplacian orthogonal polynomials basis I. Uniformly distributed input The correlation matrix becomes ⎡ 0.3333 −0.0222 −0.1369 ⎤ R PT [ n ]h = ⎢⎢ −0.0222 0.0243 −0.0586 ⎥⎥ ⎢⎣ −0.1369 −0.0586 0.2584 ⎥⎦. with three eigenvalues of 0.0001, 0.1773, and 0.4386, and an eigenvalue spread of 3000.. II. Gaussian input The correlation matrix becomes ⎡ 0.1111 −0.0370 0.0480 ⎤ R PT [ n ]h = ⎢⎢ −0.0370 0.0206 −0.0404 ⎥⎥ ⎢⎣ 0.0480 −0.0404 0.0950 ⎥⎦. with three eigenvalues of 0.0002, 0.05512, and 0.1714, and an eigenvalue spread of 871.4.. III. Laplacian input. The correlation matrix becomes 0 0 ⎤ ⎡0.1111 ⎢ 0.0741 0 ⎥⎥ R PT [ n ] h = ⎢ 0 ⎢⎣ 0 0 0.5538⎥⎦. which has a small eigenvalue spread of 7.5.. 26.

(36) E. General comparison Finally, we discuss the eigenvalue spread under non-perfect orthogonality, a case of mismatch between the pdf of the input signal and the nonlinear polynomial basis. Table 2.6 compares the eigenvalue spreads for different input pdf’s versus different orthogonal polynomial bases. If the nonlinear components have perfect orthogonality, the minimum eigenvalue spread can be achieved, as can be seen from the diagonal entries in Table 2.6. When the input signal’s pdf is unknown, according to Table 2.6 the uniformly orthogonal polynomial basis is recommended.. Table 2.6 Eigenvalue spread comparison PDF. Basis. Uniform. Gaussian. Laplacian. Power-series. 222.2. 246. 3000★. 888.5. 55.5. ◆. Uniform Gaussian. 59.1. Laplacian 39.3 42 eigenvalues are 0.0001, 0.17728, 0.43863 ◆ eigenvalues are 0.00019671, 0.055117, 0.17141. 871.4. 162. 7.473. 113.5. ★. We recall the residual echo power in (2.2.9). J a [ n] = σ + 2 v. μa T. σ. 2 v. N. ∑ i =0. λi μa. μa. σ v2. μ + ∑ λi [ a 2 - T ](1- a λi ) 2 n . μ T i =1 2 - λi 2 - a λi T T N. 2. The convergence rate is dependent on the eigenvalue spread; a smaller eigenvalue has faster convergence rate. The Laplacian polynomial basis has larger eigenvalue spread than the others because it has an eigenvalue closely to zero. But, it is worthy to note that a smaller eigenvalue has less contribution to the residual echo power. Therefore, we may calculate the eigenvalue spread excluding zero eigenvalue in such case the Laplacian type basis has a smaller eigenvalue than the power series basis.. 27.

(37) 2.5 Application of dual loudspeakers system In this section we consider the hands-free system with dual loudspeakers. 3D sound is an essential element of the new services (mobile, multimedia, etc.). It enriches sound playback more vividly. In fact the 3D sound effect can be generated by dual loudspeakers. The AEC is also essential to achieve satisfactory speech quality in such system. We use two nonlinear AEC with a cascade structure to cancel the left and right channel echo signal. The nonlinear AEC structure for dual loudspeakers system is shown in Fig. 2.3. Let xL [n] , xR [n] denote the left-channel and right-channel signal and uncorrelated with each other. Nonlinear AEC1 and AEC2 are used to cancel the echo components from the left and right loudspeaker, respectively. SP1. xL [n] SP2. xR [n]. Nonlinear. Nonlinear. AEC1. AEC2. y [ n] 1. y [n] 2. h R [ n]. h L [ n] N ea r-en d. d [ n]. e[n]. s ig n a l v [ n ] Microphone. Fig. 2.3 Nonlinear AEC structure for dual loudspeakers system. The signals, passing through the echo path h L [n] and h R [n] , respectively, are picked up by the microphone. The microphone output signal is expressed as. 28.

(38) ⎡ hL 0 ⎤ ⎡ xL ,1[n] xL ,1[n -1] " xL ,1[n - M + 1] ⎤ ⎢⎢ hL1 ⎥⎥ d [n] = [ a1a3 ] ⎢ ⎥ ⎣ xL ,3 [n] xL ,3[ n -1] " xL ,3[n - M + 1] ⎦ ⎢ # ⎥ ⎢ ⎥ ⎣ hLM -1 ⎦ ⎡ hR 0 ⎤ ⎡ xR ,1[n] xR ,1[n -1] " xR ,1[n - M + 1] ⎤ ⎢ hR1 ⎥ ⎥ + v[n] . +[b1b3 ] ⎢ ⎥⎢ ⎢ ⎥ + " [ ] [ -1] [ 1] x n x n x n M # R ,3 R ,3 ⎣ R ,3 ⎦ ⎢ ⎥ ⎣ hRM -1 ⎦. We rewrite d [n] as vector inner products: d [n] = v[n] + aT PxTL [n]h L + bT PxTR [n]h R .. (2.5.1). In this section we derive the residual echo power under the assumption of perfect linear coefficients. We derive the nonlinear coefficients convergence rate first. It is easy to see that (2.5.1) is not the same with the case of a single loudspeaker, here, we have two kinds of nonlinear coefficients in the right-hand side of (2.5.1). We can cascade up the nonlinear coefficients a and b together into a new vector c . By the same procedure, we merge PxTL [n] and PxTR [n] into a block-diagonal matrix PsT [n] . The desired signal can be rewritten as ⎡ PxTL [n]h L ⎤ d [n] = v[n] + [a b ] ⎢ T ⎥ ⎣⎢ PxR [n]h R ⎦⎥ T. T. T. T. ⎡ PxTL [n] 0 ⎤ ⎡h L ⎤ = v[n] + [a1[n]a3 [n]b1[n]b3 [n]] ⎢ ⎥ ⎢ ⎥ PxTR [n]⎥⎦ ⎣h R ⎦ ⎢⎣ 0. = v[n] + cT [n]PsT [n]h s where h s denotes [h L. (2.5.2). hR ]. T. (2.5.2) is almost the same as the single loudspeaker case. The nonlinear AEC output and residual echo become. 29.

(39) y[n] = c T [n]P T [n]h s s e[n] = d [n] - y[ n] T = d [n] - c [n]PsT [n]h s .. The NLMS adaptation algorithm becomes ∂e 2 [n] ∇c = = -2e[n]PsT [n]h s ∂c[n]. PsT [n]h s. c[n + 1] = c[n] + μc. 2. Ps [n]h s + δ T. e[n] .. (2.5.3). 2. The residual echo convergence rate can be found following the same procedures in Section 2.2.2. The nonlinear coefficients error variance and residual echo power are described in (2.5.4) and (2.5.5).. μa. 2. E[ kai [n] ] =. μa. σ v2. σ v2. μ +( c 2 - T )(1- a λi ) 2 n μ μ T 2 − a λi 2 − a λi T T T. 2 v. N. λi μa. (2.5.4). μa. σ v2. μ J a [ n] = σ + σ ∑ + ∑ λi [ c 2 - T ](1- a λi ) 2 n (2.5.5) μ T T i =0 i =1 2 - λi 2 - a λi T T 2 v. μa. 2. N. where 2. T = PsT [n]h s + δ 2. T. ⎡ PxTL [n] 0 ⎤ ⎡h L ⎤ T Ps [n]h s = ⎢ ⎥ ⎢ ⎥ PxTR [n]⎥⎦ ⎣h R ⎦ ⎢⎣ 0 2. 2. c 2 = [a b] 2 .. 30. 2.

(40) 2.6 Summary In this chapter we show a nonlinear AEC with a signal dependent orthogonal basis to decrease the correlation among different polynomial orders and increase the convergence rate. We have presented convergence analyses of linear and nonlinear NLMS algorithm, under the assumption of the other part of coefficients are perfectly known. Because the linear and nonlinear coefficients errors affect each other in the cascade structure, it is difficult to perform the joint error analysis theoretically. For an input with unknown pdf, the orthogonality of the polynomial bases may not be perfect, the eigenvalue spread analyses in Section 2.4 shows that we also have faster convergence rate than conventional power series basis. For dual loudspeakers case, we cascade the coefficients together therefore the analyses in single loudspeaker case can be easily extended to the dual loudspeakers case.. 31.

(41) Chapter 3 Training Sequence and Coefficients Estimation The training sequence (TS) can also be used in adaptive equalizer or echo path estimation. Here we will apply it to nonlinear AEC with a cascade structure. There are applications where it is necessary to compare one reference signal with one or more signals to determine the similarity between these two and to determine additional information based on the similarity. For example, in digital communications, a set of data symbols are represented by a set of unique discrete time sequences. If one of these sequences is transmitted, the receiver has to determine which particular sequence has been received by comparing the received signal with every member of possible sequences from the set. Similarly, in radar and sonar applications, the received signal reflected from the target is the delayed version of the transmitted signal and by measuring the delay; one can determine the location of the target. In Section 3.2, we will derive the TS algorithm under expectation operator. In practice, we replace expectation with sample mean. We will show the recursive form of TS algorithm in Section 3.3. In Section 3.4, we will perform the convergence analysis of linear TS coefficients and nonlinear convergence analysis in Section 3.5. In Section 3.6, we extend applications to dual loudspeakers system where the linear convergence will be given.. 32.

(42) 3.1 Correlation based nonlinear AEC In Chapter 2 we have known that the NLMS adaptive algorithm is based on a simple stochastic gradient. The kind of adaptive algorithm has good performance when the background noise or double talk is not present. In Fig. 2.1 the error signal e[n] is used to update both the linear filter and the nonlinear processor. However, when the background noise is present and/or larger than the echo signal y[n] , the desired signal contains not only the echo signal but also the background noise. The coefficients may diverge when the error signal includes a significant near-end speech signal. In this chapter we use a white sequence (non polar signal) to train both linear and nonlinear coefficients to overcome the problem of adaptive filter under low SNR condition. Because of the cascade structure we have to modify the Wiener-Hopf equation to fulfill our requirement. In the case of our structure, the TS runs through a multipath nonlinear processor and weighted by nonlinear coefficients, the sum of multipath signals passes through a linear filter. Therefore, the vector direction of the estimated coefficients is parallel to the vector of the room impulse response and its vector length is composed of nonlinear coefficients. The other problem is how to separate each nonlinear coefficient from the length of estimated coefficients. Here we only have one cross correlation length between the input and the microphone signal, in order to get more information we can use the nonlinear order of input signal and the microphone output signal to generate the other information about the nonlinear loudspeaker. The nonlinear coefficients can be found by solving these equation sets. According to this concept, the structure of a nonlinear AEC based on TS is shown in Fig. 3.1.. 33.

(43) Nonlinear Loudspeaker. non polar white signal x[n] (⋅) 2 ...... (⋅) N. (⋅)1 p1[n]. p2 [n]. p N [ n]. Room. a 1[ n] a 2 [ n] a N [n]. h[n]. Correlator s[n] Echo. [n] FIR h. y[n]. d [ n]. Near end signal v[n]. Residual echo e[n]. Microphone. Fig.3.1 Structure of a nonlinear AEC based on training sequence In Fig. 3.1, the correlator produces the correlation between p1[n] and d [n] . To find the nonlinear coefficients, the correlator will also create the correlation between p j [n] and d [n] , for 2 ≤ j ≤ N. 34.

(44) 3.2 Training sequence estimation algorithm In this section we will derive the linear TS coefficients algorithm first in Section 3.2.1. In Section 3.2.2 we find the nonlinear TS coefficients algorithm. We also give the TS algorithm for nonlinear coefficients when an orthogonal basis is used.. 3.2.1 Linear TS coefficients estimation algorithm. In the section we will derive how to use the TS to find out the linear coefficients. First, we generate a white and zero mean sequence x[n] by TS generator at the near end. This sequence cannot be a polar signal in order to have nontrivial higher order moments, which will be explained later. The sequence is injected to the loudspeaker and through the near end room it is picked up by the microphone. The desired signal is expressed as follows: M -1. N. d [n] = ∑ at ∑ hi xt [n - i ] + v[n] t =1. i =0. M -1. M -1. M -1. i =0. i =0. i =0. = a1 ∑ hi x[n - i ] + a2 ∑ hi x 2 [n - i ] + " + aN ∑ hi x N [n - i ] + v[n] .(3.2.1) M is the length of the linear filter, N is the nonlinear order. We multiply (3.2.1) by x[n - k ] and take the expectation value to get M -1. N. E[ x[n - k ]d [n]] = ∑ at ∑ hi E[ x[n - k ]xt [n - i]] t =1. i =0. M -1. M -1. i =0. i =0. = a1 ∑ hi E[ x[n - k ]x[n - i ]] + a2 ∑ hi E[ x[n - k ]x 2 [n - i ]] + " M -1. + aN ∑ hi E[ x[n - k ]x N [n - i ]] + E[ x[ n - k ]v[n]] , k = 0,1.". (3.2.2). i =0. The two expectations in (3.2.2) may be written as M -1. N. rxd [-k ] = ∑ at ∑ hi rxxt [i - k ] t =1. i =0. M -1. M -1. M -1. i =0. i =0. i =0. = a1 ∑ hi rxx [i - k ] + a2 ∑ hi rxx2 [i - k ] + " + aN ∑ hi rxx N [i - k ]. where. 35. (3.2.3).

(45) rxd [-k ] = E[ x[n - k ]d [n]] is defined as the cross correlation function between the input signal x[n] and the microphone output signal d [n] for a lag of -k ; and rxxt [i - k ] = E[ x[n - k ]xt [n - i ]] ,. rxxt [i - k ] is defined as the autocorrelation function of the input signal for a lag of i - k . The last term of (3.2.2), E[ x[n - k ]v[n]] , is equal to zero due to the input signal is. independent of the background noise. We define rxx = [rxx [0 - k ], rxx [1- k ]," , rxx [ M -1- k ]]. and h = [h0 , h1 ," , hM -1 ]T ,. then we rewrite (3.2.3) into the form of inner products: N. rxd [-k ]= ∑ at rxxt h t=1. = a1rxxh + a 2rxx2 h + " + a N rxx N h .. (3.2.4). According to (3.2.4), it is similar to the Wiener-Hopf equation. The left-hand side of (3.2.4) contains the cross correlation, but its right-hand side contains not only linear but also nonlinear coefficients. Therefore, we need to take some procedures to estimate the linear and nonlinear coefficients. Before the procedure, we extend (3.2.4) into matrix form. The impulse response is defined by the finite set of tap weights, i.e., h = [h0 , h1 ," , hM -1 ]T , hence we let the lag index k go from 0 to M-1. We extend (3.2.4) in matrix form as follows:. R xd. ⎡ rxd [0] ⎤ ⎢ ⎥ rxd [-1] ⎥ ⎢ = E [ x[n]d [n]] = ⎢ # ⎥ ⎢ ⎥ ⎣ rxd [1- M ]⎦. 36.

(46) N. = ∑ at R xxt h t =1. =a1R xxh + a2 R xx2 h + " + aN R xx N h. (3.2.5). x[n] denotes the M-by-1 vector of the tap inputs x[n], x[n -1],..., x[n - M + 1] . R xd denotes the M-by-1 cross correlation vector between the input signal and the microphone output signal. R xx denotes the M-by-M correlation matrix of x[n] ⎡ x[ n ] x[ n ] ⎢ x[ n - 1] x[ n ] R xx = E ⎢ # ⎢ ⎢⎣ x[ n - M + 1] x[ n ]. x[ n ] x[ n - 1]. ". %. #. % ". ⎤ ⎥ ⎥ # ⎥ ⎥ x[ n - M + 1] x[ n - M + 1] ⎦ x[ n ] x[ n - M + 1]. ". = m2 I M ×M . R xx3 = m4 I M ×M #. R xxt = mt +1I M ×M for 1 ≤ t ≤ N and t is odd where the moment is defined as mi = E[ x i [n]] . We note that R xx2 , R xx4 ," , R xxt are zero matrices for 2 ≤ t ≤ N , t is even , since the input signal is white zero mean, a symmetric pdf, and uncorrelated between different orders, E[ x[n]x 2 [n]] = 0 . Here we let N be an odd integer and rewrite (3.2.5) as R xd =. N. ∑. t =1 t is odd. at mt +1h. = a1m2h + a3m4h + " + aN mN +1h .. (3.2.6). We may then pre-multiply both sides of (a1m2 + a3 m4 + " + aN mN +1 )-1 I M ×M and solve (3.2.6) for h , but the nonlinear coefficients are unknown yet. We assume that h. 2 2. equals to 1 thus the direction of R xd is equal to h , then the room impulse can be found by. 37.

(47) h=. R xd R xd. .. (3.2.7). 2. 3.2.2 Nonlinear TS coefficients estimation algorithm We have to solve the nonlinear coefficients but we only have one equation in (3.2.6). Because the input signal pdf is symmetric, the odd-order moments are zeros. We divide the nonlinear TS algorithm into two groups, odd and even-ordered. First, we solve the odd-ordered nonlinear coefficients. Again, we multiply (3.2.1) by x3 [n - k ] ,…, xt [n - k ] for 1 ≤ t ≤ N , t is odd , respectively and take the expectation value. Similar to (3.2.6), we can get the other equations: R x3d = (a1m4 + a3m6 + " + aN mN +3 )h. (3.2.8). #. R x N d = (a1mN +1 + a3 mN +3 + " + aN m2 N )h .. (3.2.9). We multiply (3.2.6), (3.2.8), and (3.2.9) by h T respectively and then we can solve the following equations to find the odd-ordered nonlinear coefficients a1 , a3 ," , aN ⎡hT R xd ⎤ ⎡ m m4 " mN +1 ⎤ ⎡ a1 ⎤ ⎢ T ⎥ ⎢ 2 m6 " mN +3 ⎥⎥ ⎢⎢ a3 ⎥⎥ ⎢h R x3d ⎥ ⎢ m4 . ⎢ ⎥=⎢ # # % # ⎥ ⎢# ⎥ ⎢ # ⎥ ⎢ ⎥⎢ ⎥ ⎢hT R N ⎥ ⎣ mN +1 mN +3 " m2 N ⎦ ⎣ aN ⎦ x d⎦ ⎣ . G odd. Hence, the odd-ordered nonlinear coefficients can be found as m4 " mN +1 ⎤ ⎡ a1 ⎤ ⎡ m2 ⎢a ⎥ ⎢ m m6 " mN +3 ⎥⎥ ⎢ 3 ⎥=⎢ 4 ⎢# ⎥ ⎢ # # % # ⎥ ⎢ ⎥ ⎢ ⎥ mN +1 mN +3 " m2 N ⎦ ⎣ aN ⎦ ⎣ . G odd. We note that the matrix. G odd. -1. ⎡hT R xd ⎤ ⎢ T ⎥ ⎢h R x3 d ⎥ ⎢ ⎥. # ⎢ ⎥ ⎢ hT R N ⎥ x d⎦ ⎣. (3.2.10). in (3.2.10) has to be nonsingular in order to have a. unique solution for the nonlinear coefficients a1 , a2 ,..., aN . As mentioned earlier, the input training sequence x[n] cannot be a polar ±1 signal. Up to now, we have found 38.

(48) the. linear. coefficients. h. and. odd-ordered. nonlinear. coefficients at for. 1 ≤ t ≤ N , t is odd . Next we derive the even-ordered nonlinear coefficients. Similarly, we start with (3.2.2) but there are some differences from previous procedure. Although the input signal is white and zero mean but the even-ordered moments are not equal to zero. In order to get a diagonal correlation matrix we multiply (3.2.1) by x 2 [n - k ] - m2 and take expectation to get R ( x2 -m ) d = a2 (m4 - m22 )h + a4 (m6 - m2 m4 )h + " + a4 (mN + 2 - m2 mN )h 2. (3.2.11). Again, we multiply (3.2.1) by x 4 [n - k ] - m4 ," , x N -1[ n - k ] - mN −1 and take expectation value #. R ( x N -1 -m. N -1 ) d. = a2 (mN +1 - m2 mN −1 )h + a4 (mN +3 - m4 mN −1 )h + " + aN −1 (m2 N − 2 - mN2 −1 )h. (3.2.12). By multiplying (3.2.11) and (3.2.12) by hT respectively, we can find out the nonlinear coefficients at for 1 ≤ t ≤ N , t is even , by solving the following equations ⎡ hT R ( x 2 - m ) d ⎤ 2 2 m6 - m2 m4 " mN + 2 - m2 mN ⎤ ⎡ a2 ⎤ ⎢ ⎥ ⎡ m4 - m2 ⎢ ⎥⎢ T ⎥ ⎢h R ( x 4 -m ) d ⎥ m6 - m2 m4 m8 - m42 " mN +3 - m4 mN -1 ⎥ ⎢ a4 ⎥ ⎢ 4 . ⎢ ⎥= ⎥⎢ # ⎥ # # % # ⎢ # ⎥ ⎢ ⎥⎢ ⎥ ⎢ T ⎥ ⎢⎢ m - m m mN +3 - m4 mN -1 " m2 N -2 - mN2 -1 ⎦⎥ ⎣ aN -1 ⎦ N +1 2 N -1 ⎣ h R N − 1 ⎢⎣ ( x - mN −1 ) d ⎥ ⎦ . G even. Hence, the even-ordered nonlinear coefficients can be found as m6 - m2 m4 " mN + 2 - m2 mN ⎤ ⎡ a2 ⎤ ⎡ m4 - m22 ⎥ ⎢a ⎥ ⎢ 2 m8 - m4 " mN +3 - m4 mN -1 ⎥ ⎢ 4 ⎥ = ⎢ m6 - m2 m4 ⎥ ⎢ # ⎥ ⎢ # # % # ⎥ ⎢ ⎥ ⎢ 2 mN +1 - m2 mN -1 mN +3 - m4 mN -1 " m2 N -2 - mN -1 ⎥⎦ ⎣ aN -1 ⎦ ⎢⎣ . G even. -1. ⎡ hT R ( x 2 - m ) d ⎤ 2 ⎥ ⎢ T ⎢h R ( x 4 - m ) d ⎥ 4 ⎥. ⎢ ⎥ ⎢ # ⎥ ⎢ T h R ( x N −1 -m ) d ⎥ N −1 ⎦ ⎣⎢ (3.2.13). 39.

(49) Again, we have assumed that the training white sequence x[n] is non polar signal to avoid the singularity of the matrix. G even .. Finally, in Chapter 2 we mentioned a. nonlinear processor with orthogonal basis. When the orthogonal basis is used the same procedure can also be used to find out the coefficients. The desired signal can be written as M -1. N. d [n] = ∑ at ∑ hi pt [n - i ] t =1. i =0. M -1. M -1. M -1. i =0. i =0. i =0. = a1 ∑ hi p1[n - i ] + a2 ∑ hi p2 [n - i ] + " + aN ∑ hi pN [n - i ] + v[n]. where pi. becomes the nonlinear orthogonal polynomial. Following the same. procedure from (3.2.8) to (3.2.13) , we can find the linear and nonlinear coefficients 0 " 0 ⎤ ⎡ a1 ⎤ ⎡ mP1 ⎢ ⎥ ⎢a ⎥ 0 mP2 " 0 ⎥ ⎢ 2 ⎢ ⎥= ⎢# ⎥ ⎢ # # % # ⎥ ⎢ ⎥ ⎢ ⎥ 0 " mPN ⎥⎦ ⎣ aN ⎦ ⎢⎣ 0 . -1. ⎡ hT R p1d ⎤ ⎢ T ⎥ ⎢ h R p2 d ⎥ ⎢ ⎥. # ⎢ ⎥ ⎢hT R p d ⎥ N ⎦ ⎣. (3.2.14). G orthogoanl. (3.2.14) shows the nonlinear coefficients, where m p is equal to the second moment i of the basis, pi. The matrix inversion in (3.2.14) is simpler than that of non-orthogonal basis because it is a diagonal matrix. Because of the orthogonality (3.2.16) can represent the odd and even-ordered coefficients simultaneously.. 40.

(50) 3.3 Recursive Training Sequence algorithm. In order to compute the linear and nonlinear coefficients in (3.2.7), (3.2.10), and (3.2.13), we need to compute the autocorrelation function and cross correlation function. In practical system, we replace the expectation operation with sample mean. The sample mean of a set x[1], x[2],..., x[ n] of n observations is defined by sample mean =. 1 n ∑ x[i] n i =1. In Section 3.2, the expectation operator is performed to compute coefficients thus the true coefficients notation h is used. In this section the sample mean method is used, therefore we use the estimation’s notation as that in Chapter 2. At first, we use the sample mean as follows. l xd [n] = R. 1 n ∑ x[i]d[i] n i=1. =. 1 ⎡ n -1 n-1 ⎤ x[i ] y[i ] + x[n] y[n]⎥ ∑ ⎢ n ⎣ n -1 i =1 ⎦. =. n -1 l 1 R xd [n -1] + x[n]d [ n] . n n. (3.3.1). According to (3.3.1), we can rewrite (3.2.7) as follows. l xd [n] R h[n] = l xd [n] R =. =. 2. 1. l xd [n] R. (. n -1 l 1 R xd [n -1] + x[n]d [n]) n n. 2. n -1 1 h[n -1] + x[n]d [n] . n l xd [n] n R. (3.3.2). 2. In (3.3.2),. l xd [n] is equal to (a 1[n -1]m + a 3 [n -1]m + " + a N [n -1]m ) , the R 2 4 N +1 2. nonlinear coefficients can be obtained from the last iteration of (3.3.4). (3.3.2) is similar to (2.1.1), the second term on the right hand side of (3.3.2) represents the 41.

(51) adjustment that is applied to the current estimate of linear coefficients. But the error signal does not appear in (3.3.2), it is different from the NLMS-type adaptation algorithm with an error feedback structure. The TS estimation of nonlinear coefficients in (3.2.10) can be written as a recursive form too. First, we have to replace the correlation matrix of (3.2.10) with the sample mean. We can rewrite the cross correlation as l x d [n] ⎤ l x d [n -1] ⎤ ⎡R ⎡R ⎡ x [n]d [n] ⎤ ⎢ ⎥ ⎢ ⎥ l x3d [n] ⎥ n -1 ⎢ R l x3d [n -1] ⎥ 1 ⎢ x3[n]d [n] ⎥ ⎢R ⎢ ⎥. ⎢ ⎥= n ⎢ ⎥+ n ⎢ ⎥ # ⎢ # ⎥ ⎢ # ⎥ ⎢ ⎥ N ⎢R ⎢R l x N d [ n ]⎥ l x N d [n -1]⎥ ⎢ ⎥⎦ [ ] [ ] x n d n ⎣ ⎣ ⎦ ⎣ ⎦. (3.3.3). According to (3.3.3), (3.2.10) can be written into a recursive form ⎡ a 1[n] ⎤ ⎢ ⎥ ⎢ a 3 [n] ⎥ -1 ⎢ ⎥ = G odd ⎢ # ⎥ ⎢ a N [n]⎥ ⎣ ⎦. ⎡ ⎡h T [n]R ⎡h T [n]x [n] ⎤ ⎤ l x d [n -1]⎤ ⎢ ⎢ ⎥ ⎢ ⎥⎥ T ⎢ n -1 ⎢ T l 3 ⎥ 3 ⎥ ⎢ ⎥ h [n]R x d [n -1] 1 h [n]x [n] ⎥ ⎢ + d [ n ] ⎢ ⎥ ⎢ ⎥ ⎢ n ⎢ # ⎥ n ⎢ # ⎥⎥ ⎢ ⎢ T ⎥ ⎢ T ⎥⎥ l x N d [n-1] ⎥ ⎢⎣ ⎢⎣h [n]R ⎢⎣h [n]x N [n]⎥⎦ ⎥⎦ ⎦. ⎡ a 1[n -1] ⎤ ⎢ ⎥ n -1 ⎢ a 3 [n -1] ⎥ 1 = + d [n]G -1odd ⎢ ⎥ n n ⎢ # ⎥ ⎢ a N [n -1]⎥ ⎣ ⎦. ⎡h T [n]x [n] ⎤ ⎢ ⎥ ⎢h T [n]x3 [n] ⎥ ⎢ ⎥. # ⎢ ⎥ ⎢ T ⎥ ⎢⎣h [n]x N [n]⎥⎦. (3.3.4). −1 in the right-hand side of (3.3.4) can be pre-computed so The matrix inverse G odd. long as each moment of input signal in (3.2.10) are known. The same procedure can be applied to (3.2.13); the recursive form of (3.2.13) becomes ⎡h T [n](x 2 [n] - m ) ⎤ ⎡ a 2 [n] ⎤ ⎡ a 2 [n -1] ⎤ 2 ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ T 4 ⎢ ⎢ a 4 [n] ⎥ n -1 ⎢ a 4 [n -1] ⎥ 1 h [n](x [n] - m4 ) ⎥ -1 G d [ n ] = + ⎢ ⎥ even ⎢ ⎥ ⎥ n n ⎢ # # ⎢ ⎥ # ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ T ⎢ a N -1[n]⎥ ⎢ a N -1[n -1]⎥ ⎢⎣h [n](x N -1[n] - mN -1 ) ⎥⎦ ⎣ ⎦ ⎣ ⎦ where the moment matrix G even are given in (3.2.13). 42. (3.3.5).