## Abstract

In this paper, we present a heuristic method to simplify the liquid crystal adaptive optics system (LCAOS) into a single-input-single-output (SISO) system, then build the dynamic model of LCAOS based on subspace identification. Results show that the identified model could accurately describe the dynamical behavior of LCAOS (97% match), with extremely low complexity. The wonderful features of low complexity and high precision, make the identified model highly beneficial for model based controller design, system analysis and dynamical behavior simulation of liquid crystal adaptive optics systems.

© 2017 Optical Society of America

## 1. Introduction

Ever since Babcock first introduced the concept of adaptive optics in 1953 [1], adaptive optics (AO) systems have been widely used in astronomical telescopes for the past few decades, for real time compensation of the image degradation due to the Earth’s atmosphere turbulence. Up to now, most telescopes with an aperture larger than 1 meter have adopted AO systems to improve the image quality. To compensate wavefront aberration caused by atmosphere turbulence, traditional AO systems use deformable mirror (DM) as the correction device. However limited by manufacturing technology, the actuator number of DM could hardly satisfy the requirement of extreme adaptive optics systems (especially in visible wavelengths), which is the future of certain areas, like direct detection of extra solar planet.

With the advantages of low cost, high resolution, the ability of generating a desired wavefront with high-precision and no-hysteresis, liquid crystal wavefront corrector (LCWC) has become a very attractive substitute for deformable mirror (DM). Particularly, the high resolution feature makes LCWC very promising to be used in extreme AO systems. Indeed, LCWC also has some shortcomings, like low light energy efficiency and slow response speed, but they all have been subdued: the energy utilization efficiency could be improved to nearly 100% in the waveband of 0.4-0.9 um with an optimal energy-splitting open-loop system design [2]. And the response time of LCWC has been reduced to sub-millisecond with liquid crystal molecule design [3] and novel overdriving technique [4]. Up to now, LCWC has already been used in many adaptive optics systems, like retinal imaging systems [5,6], laser communication systems [7], microscopy systems [8] and large aperture telescope systems [9,10]. Z. L. Cao has reported an imaging result of 1.7 times the diffraction limit on a 1.2-meter telescope after the correction of a LCWC based AO system [2].

With certain hardware conditions, the performance of an AO system is mainly determined by the controller of the system. In order to achieve optimal control, a system model is indispensable [11–13]. And it is even more crucial for an open-loop AO system, since without a system model we would be blind to the corrector’s response to the driving signals we applied to it. And open-loop correction technology is not only adopted in LCAOS to improve energy efficiency, but also needed by DM based AO systems for extremely large telescopes (ELTs) [14] and multi-object adaptive optics (MOAO) [15]. Besides optimal controller design, during the process of appraising AO projects, a system model would be of great help for simulating the dynamical behavior of the system to make better decisions [16]. S. T. Wu’s has elegantly demonstrated the dynamic response characteristics of liquid crystal theoretically [17], however no study about modeling the whole liquid crystal adaptive optics system has been done up to date. As to DM based adaptive optics systems, there have been many studies on system modeling, both static (steady-state) [18,19] and dynamic [16,20].

The bottleneck of LCAOS modeling lies on the multi-input-multi-output (MIMO) nature of the system, which magnifies the complexity of the system modeling and reduces the precision of the acquired model. Hence we introduce a heuristic method to simplify LCAOS into a single-input-single-output (SISO) system, and build the dynamic model (SISO transfer function) of LCAOS based on subspace identification. In Section 2, we give a detailed illustration of the proposed method. Section 3 presents a real case of identifying the model of an actual LCAOS, and the conclusions are drawn in Section 4.

## 2. Principle of LCAOS system modeling

LCAOS mainly consists of a Shack-Hartmann wavefront sensor (WFS), a liquid crystal wavefront corrector (LCWC) and an imaging CCD, the principle of LCAOS is shown in Fig. 1 (where IM, CC and RQ are mere computational processes). Where ** Ф_{t}** represents the turbulence wavefront,

**is the slopes measured by WFS,**

*s***is the Zernike coefficients calculated with interaction matrix (IM) [21],**

*c*_{1}**is the output (Zernike coefficients) of the controller (CC),**

*c*_{2}**is the driving voltages acquired by wavefront reconstruction and quantification (RQ) with kinoform technique [22,23], and**

*v***is the correcting wavefront generated by LCWC, the wavefront after correction (with residual error**

*Ф*_{c}**) is imaged with CCD. This is a typical open-loop control architecture, since the residual error is not measured by WFS and fed back to the controller.**

*e*Normally in most real time DM-based AO systems, the influence functions of the DM actuators are used to expand wavefront. However, LCWC is driven by applying voltage to pixels, which are too small (several microns) and too many (tens of thousands) to measure the influence function of every individual pixel. Hence LCAOS adopts Zernike coefficients for wavefront reconstruction, so generally the controller (CC) of LCAOS is applied to Zernike coefficients, seeing in Fig. 1. Then speaking of system modeling, what we are interested in is the transfer function from Zernike coefficients to Zernike coefficients. However usually the number of Zernike modes required for wavefront reconstruction of LCAOS could be larger than 100, and cross-talking normally exists between different modes, which makes LCAOS a complex multi-input-multi-output (MIMO) system. To identify such a system model is too difficult and the precision of the identified model is hard to guarantee. So system simplification, which would dramatically reduce the complexity and improve the accuracy of system modeling, is vitally needed.

Consider that besides applying controller (CC) to Zernike coefficients like in Fig. 1, the controller could also be applied to the driving voltages of LCWC, as shown in Fig. 2. Where ** c** is the Zernike coefficients calculated with interaction matrix (IM),

**is the driving voltages acquired with wavefront reconstruction and quantification (RQ),**

*v*_{1}**is the output (driving voltages) of the controller (CC).**

*v*_{2}Based on this kind of control architecture, what we want is the transfer function ** G_{o}** from driving voltages to driving voltages, as shown in Fig. 3.

**should include the dynamics of all the procedures of LCAOS: the dynamics of LCWC and WFS, and computational delays of IM and RQ.**

*G*_{o}Generally, the total pixel number of LCWC is very huge. For instance, the pixel number of a LCWC we currently use is 65536 (256 × 256), then the driving voltages of LCWC is a vector with 65536 elements. So, outwardly LCAOS becomes even more sophisticated if we try to modeling LCAOS based on driving voltages. However, considering that what makes LCWC modulating wavefront is the movement of liquid crystal molecules, which doesn’t involve any mechanical motion, so we could reasonably assume that every pixel of LCWC has the same dynamics. With this assumption, and the fact that every pixel of WFS has the same dynamics, consider all the pixels of LCWC as an entirety, we could use a SISO (single-input-single-output) model (voltage to voltage) to describe LCAOS, this would bring significant convenience to model-based controller design and dynamical behavior simulation of the whole system.

In order to identify this SISO system model, apply a sequence of random voltages to LCWC as the input of the system (*u**(k)* = [ ** v_{i}**(k)], k = 1, 2, … 100000). Noted here, at each time, same voltage is applied to all the pixels in LCWC, and try to measure the corresponding output of the system (

*y**(k)*= [

**(k)], k = 1, 2, … 100000) (with the assumption that every pixel of LCWC has the same dynamics,**

*v*_{o}**of all the pixels would also be the same). However the difficulty is: with same voltage applied to all pixels in LCWC at each time, LCWC would generate a piston wavefront, and by nature WFS is insensitive to piston. With a piston wavefront, the slopes**

*v*_{o}**measured by WFS would be**

*s***0**, then the Zernike coefficients

**calculated with interaction matrix (IM) would also be**

*c***0**, and finally the output voltage

**would be**

*v*_{o}**0**.

So the bottleneck lies on how to acquire the system output ** v_{o}** when the same voltage

**is applied to all the pixels of LCWC. Based on polarization optical theory, LCWC could not only modulate phase but also light intensity (just need to put a polarizer and an analyzer in front of LCWC). As shown in Fig. 4, where the polarizer plays both the role of polarizer and analyzer, then the intensity**

*v*_{i}**of the light reflected from LCWC is related to the phase**

*I***of the LCWC in the following relationship:**

*Ф*And ** Ф** is related to the applied voltage

**of LCWC as follows [24]:**

*v*Where ** Ф_{m}** is the maximal value of the phase retardance,

**and**

*n*_{‖}**are the reflective indices of LCWC for a light beam with polarization parallel and perpendicular to the initial orientation of the molecules, respectively.**

*n*_{⊥}**and**

*K*_{11}**represent the splay and bend elastic constants, respectively. And**

*K*_{33}**is the threshold voltage above which a reorientation takes place in the liquid crystal.**

*U*_{th}Combining Eqs. (2) and (3), we could derive the relationship between ** I** and

**:**

*v*Practically, we can measure this relationship ** f(v)** experimentally with optical layout shown in Fig. 5: apply different voltages to LCWC respectively (at each time, same voltage is applied to all the pixels of LCWC), then measure the reflected light intensity with WFS (integration of a certain area, like the red square frame in Fig. 5). Where MS is a narrow-band spectral filter, LCWC and WFS are well placed such that they are conjugated to each other. This process of measuring

**seems similar to the calibration procedure of LCWC [25]. However, the purposes are totally different. The calibration procedure in reference 25 is to calibrate the static gain between driving voltage (gray level) and induced phase change of LCWC. While our measurement of**

*f(v)***is to calibrate the relationship between driving voltage of LCWC and the light intensity change responded in WFS. This relationship is indispensable in the latter process of LCAOS system modeling, to be used to solve the problem of WFS’ incapacity of measuring LCWC’s response (piston) when the same voltage is applied to all the pixels of LCWC. Our purpose is to build the dynamic model of the whole system (which includes the response of LCWC and WFS, as well as the system delay), so WFS in Fig. 5 can’t be replaced with a camera like in reference 25.**

*f(v)*We have transferred the insensible piston into palpable light intensity for WFS when the same voltage ** v_{i}** is applied to all the pixels of LCWC, so the bottleneck of system modeling based on driving voltages has been sidestepped. As stated before, apply a sequence of random voltages to LCWC as the input of the system (

*u**(k)*= [

**(k)], k = 1, 2, … 100000). At each time, same voltage is applied to all the pixels in LCWC. And accordingly measure the light intensity sequence (**

*v*_{i}

*I**(k)*= [

**(k)], k = 1, 2, … 100000) responded in WFS (integration of the same area shown in Fig. 5), then we could calculate the output sequence (**

*i*

*y**(k)*= [

**(k)], k = 1, 2, … 100000) with**

*v*_{o}**, as shown in Fig. 6 (To guarantee the one-to-one mapping relationship of**

*f(v)***, we only choose those voltages where**

*f(v)***is a monotonic function). In order to make sure the identified system model**

*f(v)***containing the computation delays of IM and RQ, we still carry out those two proceedings during the procedure.**

*G*_{o}Since the high linearity of LCWC [26], a linear dynamic model of LCAOS could be identified from input sequence ** u** = [

**(k)] and output sequence**

*v*_{i}**= [**

*y***(k)] with predictor-based subspace identification (PBSID**

*v*_{o}_{opt}) algorithm [27,28]. According to PBSID

_{opt}, the linear dynamic model of LCAOS could be represented in discrete state space form as:

*x**(k)*∈ℝ

*is the state variable of the state space model at time instant*

^{n}*k*.

**∈ℝ**

*A**,*

^{n*n}**∈ℝ**

*B**and*

^{n*1}**∈ℝ**

*C**are the system matrices with*

^{1*n}*n*the order of the system.

**∈ℝ**

*K**is the Kalman gain,*

^{n*1}

*ɛ**(k)*is the one-step-ahead linear prediction error of

*y**(k)*.

The matrices ** A**,

**,**

*B***and**

*C***are identified in three steps [13]: (1) the Markov parameters are estimated from the block Hankel matrices by solving a least squares problem; (2) the states of the system are estimated from Markov parameters and the Hankel matrices; (3)**

*K***,**

*A***,**

*B***and**

*C***are identified based on the Markov parameters, the Hankel matrices and the states of the system by solving two least squares problems. With matrices**

*K***,**

*A***and**

*B***identified, the dynamic model of LCAOS could be constructed as (where**

*C**z*is the

*z*-transform shift operator,

**∈ℝ**

*I**is an identity matrix):*

^{n*n}And the quality of system model ** G_{o}** could be evaluated with the variance accounted for (VAF) [29] of actual system output

**and system model output**

*v*_{o}**(**

*v*^{’}_{o}**):**

*v*^{’}_{o}= G_{o}*v_{i}A VAF of 100% means there is no difference between model output ** v^{’}_{o}** and actual system output

**, which indicates that system model**

*v*_{o}**could perfectly demonstrate the dynamics of LCAOS.**

*G*_{o}## 3. Experimental results and discussion

With the method described in section 2, we tried to identify a LCAOS with parameters as follows: the WFS has 225 (15 × 15) micro-lenses, 90 × 90 pixels and a 870 Hz frame rate; the LCWC has 256 × 256 pixels and a response time of 1.15ms; the wavefront reconstruction adopts 35 Zernike modes; and the computation delays of IM and RQ take 0.45 ms. The input voltage sequence ** v_{i}** used for system modeling and the measured corresponding output voltage sequence

**are shown in Fig. 7.**

*v*_{o}The complexity of a system model is mainly determined by system order, the larger the system order, the more complexity. And the order of system model ** G_{o}** could be estimated from the Hankel singular values [28] and VAF plots. As shown in Fig. 8(a), the third singular value is at least 50 times smaller than the first two singular values, so a model with order of 2 would be sufficient to describe LCAOS. Figure 8(b) (solid line) also confirms this conclusion, since the VAF value has already reached 97.3% with an order of 2, and hardly increases with larger orders. We have also tried to identify the system model of the same LCAOS based on coefficients of Zernike modes. As could be seen in Fig. 8(b) (dotted line), no matter how much we increase the system order, VAF of the model could only reach to 66.6%. Therefore system modeling of LCAOS based on driving voltages would dramatically reduce the complexity and increase the accuracy of the identified model, the system simplification is worth the trouble.

The identified system model ** G_{o}** is shown in Eq. (8), with step response shown in Fig. 9(a), and Bode plot in Fig. 9(b). Equation (8) is the system transfer function of LCAOS, which describes the relationship between system output and system input, and it also represents the

*z*-transform of system’s impulse response with all initial conditions assumed to be zero (where

*z*is the shift operator).

The actual system output ** v_{o}** and the system model output

**are compared in Fig. 10. As it can be seen, the identified model**

*v*^{’}_{o}**is able to predict the dynamic response of LCAOS with a very high precision (VAF = 97.3%).**

*G*_{o}In section 2, in order to simplify LCAOS into a SISO system, we made a strict assumption that the response characteristics of different areas of LCWC and WFS are identical. To verify this assumption, we applied a sequence of voltages to LCWC (at each time, same voltage was applied to every pixel of LCWC), then compared the output voltages of several different areas (as shown in Fig. 11, Area 1, 2, 3, 4 in WFS, corresponding to 4 different areas in LCWC).

The VAFs of any two different areas’ outputs are 98.8%, 98.9%, 99.0%, 98.7%, 98.8% and 99.6%, respectively. The largest variance takes place between area 2 and 3 with VAF = 98.7%, as shown in Fig. 12. This indicates that our assumption (the response characteristics between different areas of LCWC and WFS are identical) is solid, and our proposed system modeling method is well-grounded.

Even though we could obtain a high precision model of LCAOS with our proposed method, it would be time-consuming if system identification online is needed anytime we use LCAOS. So we hope LCAOS would be time-invariant, so that we could identify the system model beforehand. So as to evaluate the time-variance of the system, we compared the system outputs ** v_{o}** and

**under the same input voltage sequence**

*v*^{t}_{o}**between a time interval of 50 days. As shown in Fig. 13, the VAF of**

*v*_{i}**and**

*v*_{o}**is 99.6%. It suggests that the dynamics of LCAOS hardly change in nearly two months, LCAOS shows high time-invariance, which means we only need to do the system modeling once in at least several months.**

*v*^{t}_{o}As stated in section 1, such a high precision system model could be used for model based controller design. So we designed a simple first-order adaptive controller for the identified system, and simulated the performance of the controller. As could be seen in Fig. 14, compared to a traditional integral controller, a model based adaptive controller could reduce the RMS of residual error by 39%. Even though this is only simulation, the result already shows great potential in applying our proposed system modeling method to model based controller design.

## 4. Conclusion

This paper presents a heuristic method to simplify liquid crystal adaptive optics system (LCAOS) into a single-input-single-output (SISO) system, and build the dynamic model (SISO transfer function) of LCAOS based on subspace identification. With system simplification, the system order could be reduce to 2, and the model precision could be improved from 66.6% to 97.3% (VAF). And we have also proved that LCAOS is highly time-invariant, which is a good feature that ensures us to identify the system offline. The low-complexity and high precision properties of the identified model make it highly beneficial for model based controller design and dynamical behavior simulation of liquid crystal adaptive optics systems. A simple stimulation suggests that a model based adaptive controller could reduce the RMS of LCAOS residual error by 39%. The proposed method could also be applied to build the dynamic models of any other liquid crystal devices (such as optoelectronic switch), so as to study the dynamic characteristics of those devices.

## Funding

National Natural Science Foundation of China (NSFC) (61475152).

## Acknowledgments

This work is supported by State Key Laboratory of Applied Optics, Changchun Institute of Optics, Fine Mechanics and Physics, Chinese Academy of Sciences.

## References and links

**1. **H. W. Babcock, “The possibility of compensating astronomical seeing,” Publ. Astron. Soc. Pac. **65**(386), 229–236 (1953). [CrossRef]

**2. **Z. Cao, Q. Mu, L. Hu, Y. Liu, Z. Peng, Q. Yang, H. Meng, L. Yao, and L. Xuan, “Optimal energy-splitting method for an open-loop liquid crystal adaptive optics system,” Opt. Express **20**(17), 19331–19342 (2012). [CrossRef] [PubMed]

**3. **Z. H. Peng, Q. D. Wang, Y. G. Liu, Q. Q. Mu, Z. L. Cao, H. Y. Xu, P. G. Zhang, C. L. Yang, L. S. Yao, L. Xuan, and Z. Y. Zhang, “Electrooptical properties of new type fluorinated phenyl-tolane isothiocyanate liquid crystal compounds,” Liq. Cryst. **43**(2), 276–284 (2016). [CrossRef]

**4. **H. Hu, L. Hu, Z. Peng, Q. Mu, X. Zhang, C. Liu, and L. Xuan, “Advanced single-frame overdriving for liquid-crystal spatial light modulators,” Opt. Lett. **37**(16), 3324–3326 (2012). [CrossRef] [PubMed]

**5. **K. Takeno and T. Shirai, “Chromatic aberration free liquid crystal adaptive optics for flood illuminated retinal camera,” Opt. Commun. **285**(12), 2967–2971 (2012). [CrossRef]

**6. **N. Kong, C. Li, M. Xia, D. Li, Y. Qi, and L. Xuan, “Optimization of the open-loop liquid crystal adaptive optics retinal imaging system,” J. Biomed. Opt. **17**(2), 026001 (2012). [CrossRef] [PubMed]

**7. **F. Feng, I. H. White, and T. D. Wilkinson, “Aberration Correction for Free Space Optical Communications Using Rectangular Zernike Modal Wavefront Sensing,” J. Lightwave Technol. **32**(6), 1239–1245 (2014). [CrossRef]

**8. **M. J. Booth, “Adaptive optical microscopy: the ongoing quest for a perfect image,” Light Sci. Appl. **3**(4), e165 (2014). [CrossRef]

**9. **D. Dayton, J. Gonglewski, S. Restaino, J. Martin, J. Phillips, M. Hartman, S. Browne, P. Kervin, J. Snodgrass, N. Heimann, M. Shilko, R. Pohle, B. Carrion, C. Smith, and D. Thiel, “Demonstration of new technology MEMS and liquid crystal adaptive optics on bright astronomical objects and satellites,” Opt. Express **10**(25), 1508–1519 (2002). [CrossRef] [PubMed]

**10. **Q. Q. Mu, Z. L. Cao, L. F. Hu, Y. G. Liu, Z. H. Peng, L. S. Yao, and L. Xuan, “Open loop adaptive optics testbed on 2.16 meter telescope with liquid crystal corrector,” Opt. Commun. **285**(6), 896–899 (2012). [CrossRef]

**11. **R. Muradore, E. Fedrigo, and C. Correia, “LQ control design for adaptive optics systems based on MIMO identified model,” Proc. SPIE **6272**, 62725H (2006). [CrossRef]

**12. **T. Ruppel, S. H. Dong, F. Rooms, W. Osten, and O. Sawodny, “Feedforward Control of Deformable Membrane Mirrors for Adaptive Optics,” IEEE T. Contr. Syst. T. **21**(3), 579–589 (2013). [CrossRef]

**13. **H. Song, R. Fraanje, G. Schitter, G. Vdovin, and M. Verhaegen, “Controller Design for a High-Sampling-Rate Closed-Loop Adaptive Optics System with Piezo-Driven Deformable Mirror,” Eur. J. Control **17**(3), 290–301 (2011). [CrossRef]

**14. **D. T. Gavel, “Adaptive optics control strategies for extremely large telescopes,” Proc. SPIE **4494**, 215–220 (2002). [CrossRef]

**15. **F. Vidal, E. Gendron, G. Rousset, T. Morris, A. Basden, R. Myers, M. Brangier, F. Chemla, N. Dipper, D. Gratadour, D. Henry, Z. Hubert, A. Longmore, O. Martin, G. Talbot, and E. Younger, “Analysis of on-sky MOAO performance of CANARY using natural guide stars,” Astron. Astrophys. **569**(3), 1123–1133 (2014).

**16. **A. Haber, A. Polo, S. Ravensbergen, H. P. Urbach, and M. Verhaegen, “Identification of a dynamical model of a thermally actuated deformable mirror,” Opt. Lett. **38**(16), 3061–3064 (2013). [CrossRef] [PubMed]

**17. **S.-T. Wu and C.-S. Wu, “Experimental confirmation of the Osipov-Terentjev theory on the viscosity of nematic liquid crystals,” Phys. Rev. A **42**(4), 2219–2227 (1990). [CrossRef] [PubMed]

**18. **C. R. Vogel, “Hysteresis Modeling and Go-To Control of Deformable Mirrors in Adaptive Optics,” in IEEE Aerospace Conference (2012). [CrossRef]

**19. **C. Vogel, G. Tyler, Y. Lu, T. Bifano, R. Conan, and C. Blain, “Modeling and parameter estimation for point-actuated continuous-facesheet deformable mirrors,” J. Opt. Soc. Am. A **27**(11), A56–A63 (2010). [CrossRef] [PubMed]

**20. **D. P. Looze, “Discrete-time model for an adaptive optics system with input delay,” Int. J. Control **83**(6), 1217–1231 (2010). [CrossRef]

**21. **X. Zhang, L. Hu, Z. Cao, Q. Mu, D. Li, and L. Xuan, “Improve the accuracy of interaction matrix measurement for liquid-crystal adaptive optics systems,” Opt. Express **22**(12), 14221–14228 (2014). [CrossRef] [PubMed]

**22. **J. A. Jordan Jr, P. M. Hirsch, L. B. Lesem, and D. L. Van Rooy, “Kinoform lenses,” Appl. Opt. **9**(8), 1883–1887 (1970). [CrossRef] [PubMed]

**23. **Y. J. Liu, Z. L. Cao, D. Y. Li, Q. Q. Mu, L. F. Hu, X. H. Lu, and L. Xuan, “Correction for large aberration with phase-only liquid-crystal wavefront corrector,” Opt. Eng. **45**(12), 128001 (2006). [CrossRef]

**24. **A. V. Kudryashov, J. Gonglewski, S. Browne, and R. Highland, “Liquid crystal phase modulator for adaptive optics. Temporal performance characterization,” Opt. Commun. **141**(5–6), 247–253 (1997). [CrossRef]

**25. **P. Prieto, E. Fernández, S. Manzanera, and P. Artal, “Adaptive optics with a programmable phase modulator: applications in the human eye,” Opt. Express **12**(17), 4059–4071 (2004). [CrossRef] [PubMed]

**26. **C. Li, M. Xia, Q. Mu, B. Jiang, L. Xuan, and Z. Cao, “High-precision open-loop adaptive optics system based on LC-SLM,” Opt. Express **17**(13), 10774–10781 (2009). [CrossRef] [PubMed]

**27. **A. Chiuso, R. Muradore, and E. Marchetti, “Dynamic Calibration of Adaptive Optics Systems: A System Identification Approach,” in Proceedings of IEEE Conference on Decision and Control (IEEE, 2008), pp. 750–755. [CrossRef]

**28. **A. Chiuso, R. Muradore, and E. Marchetti, “Dynamic Calibration of Adaptive Optics Systems: A System Identification Approach,” IEEE T. Contr. Syst. T. **18**(3), 705–713 (2010). [CrossRef]

**29. **M. Verhaegen and V. Verdult, *Filtering and System Identification: A Least Squares Approach* (Cambridge University, 2007), Chap. 10.