Issue 
Mechanics & Industry
Volume 18, Number 8, 2017
Experimental Vibration Analysis



Article Number  802  
Number of page(s)  10  
DOI  https://doi.org/10.1051/meca/2017044  
Published online  21 March 2018 
Regular Article
Modal parameter identification of a CMUT membrane using response data only
^{1}
FEMTO–ST Institute, DMA, BFC University,
25000
Besançon, France
^{2}
Université Clermont Auvergne, SigmaClermont, Institut Pascal,
BP 10448,
63000
ClermontFerrand, France
^{*} email: joseph.lardies@univfcomte.fr
Received:
12
December
2016
Accepted:
6
November
2017
Capacitive micromachined ultrasonic transducers (CMUTs) are microelectromechanical systems used for the generation of ultrasounds. The fundamental element of the transducer is a clamped thin metallized membrane that vibrates under voltage variations. To control such oscillations and to optimize its dynamic response it is necessary to know the modal parameters of the membrane such as resonance frequency, damping and stiffness coefficients. The purpose of this work is to identify these parameters using only the time data obtained from the membrane center displacement. Dynamic measurements are conducted in time domain and we use two methods to identify the modal parameters: a subspace method based on an innovation model of the statespace representation and the continuous wavelet transform method based on the use of the ridge of the wavelet transform of the displacement. Experimental results are presented showing the effectiveness of these two procedures in modal parameter identification.
Key words: Modal parameters / CMUT / subspace method / wavelet transform method / experimental identification
© AFM, EDP Sciences 2017
1 Introduction
Capacitive micromachined ultrasonic transducers (CMUTs) introduced in 1994 constitute a good alternative to conventional piezoelectric transducers in terms of broader bandwidth, better transduction efficiency and the possibility of onchip integration with electronics components. Basically, a CMUT cell is a capacitor structure composed of a thin flexible membrane, metallized or not metallized, and a bottom fixed plate.If an alternating voltage, applied between the membrane and the back plate, is superimposed on the bias voltage the modulation of the electrostatic force moves the membrane which generates ultrasonic waves. Figure 1 shows a schematic crosssection of a CMUT elementary cell.
A CMUT cell uses the deformation of a suspended thin membrane to transmit ultrasonic waves and the dynamic analysis of a CMUT has been studied in time domain using a FEM package [1,2]. They determined the normal displacement at the center of the membrane by finite element analysis and by experimental tests. The idea of generating acoustic waves by the electrostatic attraction force between the plates of a capacitor is very old and a tutorial can be obtained in [3]. Dynamic and acoustic modelization of capacitive micromachined ultrasonic transducers have been proposed in [4] where the acoustic pressure emitted by a CMUT is estimated by taking into account explicitly the dynamics of each membrane of the network.
CMUTs are fabricated using silicon processing routes, leading to a thin (0.5–2 μm thick) silicon nitride or polysilicon membrane above a small gap (typically about one micron in depth) to result in a capacitive structure. In general, such CMUTs tend to be less than 0.2 mm in diameter. We note R the membrane radius, t_{m} the membrane thickness and t_{g} the gap (or cavity) thickness.
A good design of a capacitive ultrasonic sensor requires a large displacement of the membrane and an optimum energy coupling between the membrane and the air that are achieved when the membrane vibrates at its mechanical resonance frequency. Stiffness is an important property and is vital to define both the sensitivity of CMUTs sensors and the dynamic behavior of the device that is strongly dependent on material characteristics and on clamping conditions of the membrane. The quality factor is also fundamental, as it provides immediate quantification of dissipations, which are generally addressable to air damping and internal friction in the material. In the paper the membrane displacement of a CMUT cell is reduced to a single degree of freedom and only the first vibration mode is considered. Higher order modes are neglected in our model and one objective of the paper is to identify the parameters of the first vibration mode from output data only.
Numerous papers have been presented on system parameter identification from measured data. Ho and Kalman [5] introduced the principle of minimal realization theory in which a statespace model has an alternative representation with a sequence of Markov parameters. Zeiger and Mc Ewem [6] proposed a concept of combining singular value decomposition (SVD) and minimal realization for parameter identification. Based on the developments of SVD and minimal realization algorithm, the eigen system realization algorithm (ERA) was proposed by Juang [7] to identify parameters of a system from the noisy inputoutput measurements. To reduce the noise effect Juang [7] proposed the eigen system realization algorithm with data correlation (ERA/DC). James et al. [8] developed the natural excitation technique (NExT) based on the principle that a correlation function measured under natural excitation can be expressed as a sum of decaying sinusoids. Each decaying sinusoid has an eigen frequency and a damping ratio that is identical to the one of the structural mode. A time domainbased stochastic subspace identification (SSI) method has been proposed by Van Overschee and De Moor [9] and is based on the statespace model of input/output signals and parameter identification is realized by extracting system matrices from input/output signals or output signals only. The key concept of SSI, which is very time consuming, is the projection of the row space of the future time data into the row space of the past time data. The system parameter identification is also possible in a timefrequency domain using the continuous wavelet transform (CWT) of signals [10]. The CWT is a time–frequency representation and applied to the free response of a system allows the identification of eigen frequencies, damping ratios and mode shapes. Indeed, the CWT of a signal results in a complex valued form, whose modulus and phase are related to the parameters of the analysed system [11–15].
In this paper, the dynamic characterization of the metallized membrane is aimed at defining its behavior in working conditions, as well as its performances in terms of resonance frequency, stiffness coefficient, damping factor, quality factor and halfpower bandwidth. Such coefficients are identified using two methods: the subspace method derived from a statespace innovation model and the continuous wavelet transform (CWT) method based on the use of the ridge of the CWT. These methods extract the modal parameters from structural response data only, using the time displacement of the membrane center. The procedures presented in the paper can identify closely eigen frequencies that cannot be identified by the traditional Fourier transform.
The paper is organized as follows: the second section is devoted to the modelization of a vibrating system and to the generation of a discrete time statespace model. An innovation representation with the estimation of the transition matrix is derived in the third section. The continuous wavelet transform is presented in Section 4. Validity of the modal parameter identification procedures is presented in the fifth section with simulated and an experimental test in laboratory. The paper is briefly concluded in Section 6.
Fig. 1 A conceptual crosssection of a capacitive micromachined ultrasonic transducer. 
2 Modeling a vibrating mechanical system
Consider a multivariate continuous time structural dynamic system described by a positive definite diagonal mass matrix M_{o}, a symmetric semidefinite viscous damping matrix C_{o} and a symmetric positive definite stiffness matrix K_{o}. These matrices are (n′ × n′). The system is represented by the linear differential equation and the observation equation [7] $${M}_{o}{\displaystyle \ddot{\xi}}(t)+{C}_{o}{\displaystyle \ddot{\xi}}(t)+{K}_{o}\xi (t)=\eta (t)$$(1) $$y(t)=L\xi (t)+v(t)$$(2) where ξ(t) is the (n′ × 1) vector of displacement of degrees of freedom; η(t) is the (n′ × 1) unmeasured excitation vector; y(t) is the (m × 1) output observation vector; the matrix L (m × n′) specifies which points of the system are observed, namely where the sensors are located and v(t) is an (m × 1) additive noise disturbance, a vector of measurement errors simulated by a Gaussian white noise. The modal characteristics (μ, Ψ_{μ}) of the system are solutions of: $$\left{M}_{o}{\mu}^{2}+{C}_{o}\mu +{K}_{o}\right=0;\text{\hspace{0.22em}}({M}_{o}{\mu}^{2}+{C}_{o}\mu +{K}_{o}){\mathrm{\Psi}}_{\mu}=0$$(3)where   represents a determinant.
The system described by equations (1) and (2) is equivalent to the following continuous time statespace model: $$\dot{x}}(t)={\displaystyle \tilde{A}}\text{\hspace{0.22em}}x(t)+{\displaystyle \tilde{B}}\text{\hspace{0.22em}}\eta (t)$$(4) $$y(t)=C\text{\hspace{0.22em}}x(t)+v(t)$$(5) where x(t) is the vector of dimension 2n′ = n $$x(t)=\left[\begin{array}{c}\hfill \xi (t)\hfill \\ \hfill {\displaystyle \dot{\xi}}(t)\hfill \end{array}\right]$$(6)and $\tilde{A}}\text{,\hspace{0.17em}}{\displaystyle \tilde{B}}\text{,\hspace{0.17em}}C$ are the (2n′ × 2n′), (2n′ × n′), (m × 2n′) matrices $$\tilde{A}}=\left[\begin{array}{cc}\hfill 0\hfill & \hfill I\hfill \\ \hfill {{M}_{o}}^{1}{K}_{o}\hfill & \hfill {{M}_{o}}^{1}{C}_{o}\hfill \end{array}\right];\text{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}}{\displaystyle \tilde{B}}=\left[\begin{array}{c}\hfill 0\hfill \\ \hfill {{M}_{o}}^{1}\hfill \end{array}\right];\text{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}}C=[\begin{array}{cc}\hfill L\hfill & \hfill 0\hfill \end{array}]$$(7)After sampling with constant period Δt and transformation of the 2n′ first order differential equations (4) into a discrete time equation, we obtain the following discrete time statespace model and the discrete time observation equation [7,15] $${x}_{k\text{+1}}=A\text{\hspace{0.22em}}{x}_{k}\text{+}{u}_{k}$$(8) $${y}_{k}=C\text{\hspace{0.22em}}{x}_{k}\text{+}{v}_{k}$$(9)where x_{k} represents the unobserved state vector of dimension n; $A={e}^{{\displaystyle \tilde{A}}\mathrm{\Delta}t}$ is the (n × n) discrete time statespace matrix or discrete time transition matrix and u_{k} is given by $${u}_{k}={\int}_{0}^{\mathrm{\Delta}t}{e}^{{\displaystyle \tilde{A}}s}\text{\hspace{0.22em}}{\displaystyle \tilde{B}}\eta (ts)\text{\hspace{0.22em}}ds$$(10)The excitation is not measured and can be considered (a priori) as a centred random process. Note that u_{k} and v_{k} are mutually uncorrelated. The matrix C is (m × n) and is called the observation matrix and y_{k} is an (m × 1) vector of observations.
The eigen values λ and eigen vectors Φ_{λ} of the discrete time transition matrix A are related to the modal characteristics (3) by [15] $$\lambda ={e}^{\mu \mathrm{\Delta}t}\text{\hspace{0.17em}and\hspace{0.17em}}C{\mathrm{\Phi}}_{\lambda}=L{\mathrm{\Psi}}_{\mu}$$(11) and the global modal parameters: the natural frequencies f_{i} and damping ratios ζ_{i} of the vibrating system can be determined by $${f}_{i}=\frac{1}{2\pi \mathrm{\Delta}t}\sqrt{\frac{{[\mathrm{ln}({\lambda}_{i}{{\lambda}_{i}}^{*})]}^{2}}{4}+{\left[{\mathrm{cos}}^{1}\left(\frac{{\lambda}_{i}+{{\lambda}_{i}}^{*}}{2\sqrt{{\lambda}_{i}{{\lambda}_{i}}^{*}}}\right)\right]}^{2}}$$(12) $${\zeta}_{i}=\sqrt{\frac{{[\mathrm{ln}({\lambda}_{i}{{\lambda}_{i}}^{*})]}^{2}}{{[\mathrm{ln}({\lambda}_{i}{{\lambda}_{i}}^{*})]}^{2}+4{\left[{\mathrm{cos}}^{1}\left(\frac{{\lambda}_{i}+{{\lambda}_{i}}^{*}}{2\sqrt{{\lambda}_{i}{{\lambda}_{i}}^{*}}}\right)\right]}^{2}}}$$(13)for i = 1, 2, … n′
Since the state vector x_{k} of dimension n = 2n′ is not observable we estimate it and transform equations (8) and (9) into a statespace innovation form, before we can obtain the transition matrix A.
3 Statespace innovation model
We estimate x_{k} by its orthogonal projection onto the subspace spanned by the stacked vector of past data vectors noted ${y}_{k1}^{}=[{{y}^{\prime}}_{k1},\text{\hspace{0.22em}}{{y}^{\prime}}_{k2},\text{\hspace{0.22em}}\dots {]}^{\prime}$. Here, the superscript ′ denotes the transpose operation and we use the notation ${z}_{k}=E[{x}_{k}\text{\hspace{0.22em}}\text{\hspace{0.22em}}{y}_{k1}^{}]$ to denote this orthogonal projection. We obtain z_{k}_{+1} using properties of orthogonal projections given by M. Aoki [16] $${z}_{k\text{+1}}=E[{x}_{k\text{+1}}\text{\hspace{0.22em}}\text{\hspace{0.22em}}{y}_{k}^{}]=E[{x}_{k\text{+1}}\text{\hspace{0.22em}}\text{\hspace{0.22em}}{y}_{k1}^{},{e}_{k}]=A\text{\hspace{0.22em}}E[{x}_{k}\text{\hspace{0.22em}}\text{\hspace{0.22em}}{y}_{k1}^{}]+B\text{\hspace{0.22em}}{e}_{k}=A\text{\hspace{0.22em}}{z}_{k}+B\text{\hspace{0.22em}}{e}_{k}$$(14) where ${e}_{k}={y}_{k}E[{x}_{k+1}\text{\hspace{0.22em}}\text{\hspace{0.22em}}{y}_{k1}^{}]$ is called the innovation component in the data vector y_{k}. It is an (m × 1) stochastic vector with zero mean and covariance matrix $Q=E[{e}_{k}\text{\hspace{0.22em}}{{e}^{\prime}}_{k}]$. Note that e_{k} is uncorrelated with ${y}_{k1}^{}$ and hence with z_{k} by construction, and that the elements of {e_{k}} are serially uncorrelated. z_{k} is an (n × 1) vector of unobservable states. The (n × m) matrix B is called the Kalman gain. Since $$E[{x}_{k+1}\text{\hspace{0.22em}}\text{\hspace{0.22em}}{e}_{k}]=E[{x}_{k\text{+1}}\text{\hspace{0.22em}}{{e}^{\prime}}_{k}](E[{e}_{k}\text{\hspace{0.22em}}{{e}^{\prime}}_{k}{])}^{1}{e}_{k}$$(15)we have: $$B=E[{x}_{k\text{+1}}\text{\hspace{0.22em}}{{e}^{\prime}}_{k}](E[{e}_{k}\text{\hspace{0.22em}}{{e}^{\prime}}_{k}{])}^{1}$$(16)Note that by assumption $E[{u}_{k}\text{\hspace{0.22em}}\text{\hspace{0.22em}}{y}_{k1}^{}]=0$ and $E[{v}_{k}\text{\hspace{0.22em}}\text{\hspace{0.22em}}{y}_{k1}^{}]=0$. Using the definition of e_{k} and equation (9), we obtain the discrete time observation equation which contains the innovation vector $${y}_{k}=E[C\text{\hspace{0.22em}}{x}_{k}\text{+}{v}_{k}\text{\hspace{0.22em}}\text{\hspace{0.22em}}{y}_{k1}^{}]+{e}_{k}=C\text{\hspace{0.22em}}{z}_{k}+{e}_{k}$$(17)The set of equations (8) and (9) is now replaced by $${z}_{k+1}=A\text{\hspace{0.22em}}{z}_{k}+B\text{\hspace{0.22em}}{e}_{k}$$(18) $${y}_{k}=C\text{\hspace{0.22em}}{z}_{k}+{e}_{k}$$(19)Such representation is called statespace innovation representation of the data generating process equations (8) and (9). The statespace innovation model in equations (18) and (19) is a suitable representation for the synthesis of temporal samples.
Denote the covariance matrix of the stacked vector ${y}_{k1}^{}$ by $R=E[{y}_{k1}^{}\text{\hspace{0.17em}}{y}_{k1}^{\prime}]$. This matrix is theoretically infinitedimensional, although in any actual implementation we cut off the stacked vector at y_{k}_{p} for some positive integer p. We use the hat to denote an estimate based on finite amount of data. Thus ${{\displaystyle \stackrel{\u02c6}{R}}}_{}$ is (mp × mp). We develop the state vector z_{k} $${z}_{k}=E[{x}_{k}\text{\hspace{0.17em}}\text{\hspace{0.17em}}{y}_{k1}^{}]=E[{x}_{k}\text{\hspace{0.17em}}{y}_{k1}^{\prime}]{(E[{y}_{k1}^{}\text{\hspace{0.17em}}{y}_{k1}^{\prime}])}^{1}{y}_{k1}^{}=K\text{\hspace{0.17em}}{R}_{}^{1}\text{\hspace{0.17em}}{y}_{k1}^{}$$(20) where $K=E[{x}_{k}\text{\hspace{0.17em}}{y}_{k1}^{\prime}]$ is called the controllability matrix. This is a reduced form expression for the state vector. At the moment equation (21) is not operational since we do not know the matrix K, which will be shown to be determined by a suitable factorisation of the block Hankel matrix defined later. Since $K=E[{x}_{k}\text{\hspace{0.17em}}{y}_{k1}^{\prime}]=E[{z}_{k}\text{\hspace{0.17em}}{y}_{k1}^{\prime}]$ it is also a matrix with a theoretically infinite number of columns $$K=[\begin{array}{cccc}\hfill G\hfill & \hfill AG\hfill & \hfill {A}^{2}G\hfill & \hfill \dots \hfill \end{array}]$$(21)where we define the covariance matrix between z_{k} and y_{k}_{1} as $G=E[{z}_{k}\text{\hspace{0.17em}}{y}_{k1}^{\prime}]$.
To estimate K it is necessary to introduce the stacked vector of future data vectors noted ${y}_{k}^{+}=[{{y}^{\prime}}_{k},\text{\hspace{0.17em}}{{y}^{\prime}}_{k\text{+1}},\text{\hspace{0.17em}}\dots {]}^{\prime}=[{{y}^{\prime}}_{k},\text{\hspace{0.17em}}{y}_{k+1}^{{\text{+}}^{\prime}}{]}^{\prime}$, then from equations (19) and (20) we obtain $${y}_{k}^{+}=O\text{\hspace{0.22em}}{z}_{k}+F\text{\hspace{0.22em}}{e}_{k}^{+}$$(22) where the matrix O is the observability matrix of the model (19),(20) and has the form $$O=\left[\begin{array}{c}\hfill C\hfill \\ \hfill CA\hfill \\ \hfill C{A}^{2}\hfill \\ \hfill \vdots \hfill \end{array}\right]$$(23)which theoretically has an infinite number of rows and F is a lower block triangular matrix formed by the matrices 0_{m}, I_{m}, A, B and C.
Then the covariance matrix between the future and the past is $$H=E[{y}_{k}^{\text{+}}\text{\hspace{0.17em}}{y}_{k1}^{\prime}]=O\text{\hspace{0.17em}}E[{z}_{k}\text{\hspace{0.17em}}{y}_{k1}^{{}^{\prime}}]=O\text{\hspace{0.17em}}K$$(24) The matrix H, called the block Hankel matrix, is formed of (m × m) covariance matrices ${R}_{i}=E[{y}_{k\text{+}i}\text{\hspace{0.17em}}{{y}^{\prime}}_{k}]$ as submatrices i = 1, 2, … and is factored as the product of two matrices O and K. Note that sample covariance matrices calculated from data y_{k}, k = 1, …, T $${{\displaystyle \stackrel{\u02c6}{R}}}_{i}={T}^{1}{\displaystyle \sum _{k=1}^{Ti}}{y}_{k+i}{{y}^{\prime}}_{k},\text{\hspace{1em}}i=1,\text{\hspace{0.17em}}2,\text{\hspace{0.17em}}\dots $$(25)are used to construct a sample block Hankel matrix $\stackrel{\u02c6}{H}$ which is factored into the estimated observability and controllability matrices $\stackrel{\u02c6}{O}$ and $\stackrel{\u02c6}{K}$. $$\stackrel{\u02c6}{H}}={\displaystyle \stackrel{\u02c6}{O}}\text{\hspace{0.22em}}{\displaystyle \stackrel{\u02c6}{K}$$(26)A second factorisation of the estimated block Hankel matrix is obtained using the singular value decomposition of $\stackrel{\u02c6}{H}$ $$\stackrel{\u02c6}{H}}={\displaystyle \stackrel{\u02c6}{U}}\text{\hspace{0.22em}}{\displaystyle \stackrel{\u02c6}{S}}\text{\hspace{0.22em}}{\displaystyle \stackrel{\u02c6}{V}$$(27)where $\stackrel{\u02c6}{U}$ and $\stackrel{\u02c6}{V}$ contain the singular vectors of $\stackrel{\u02c6}{H}$. The singular values of $\stackrel{\u02c6}{H}$ are contained in descending order in the diagonal matrix $\stackrel{\u02c6}{S}$. When estimating a model with n states, only the n largest singular values are used. Setting the smallest singular values to zero constitutes the set of exclusion restrictions used to identify the model. Denote the approximation to the block Hankel matrix created by enforcing these exclusion restrictions by ${{\displaystyle \stackrel{\u02c6}{H}}}_{n}$. The S.V.D. of ${{\displaystyle \stackrel{\u02c6}{H}}}_{n}$ is $${{\displaystyle \stackrel{\u02c6}{H}}}_{n}={\widehat{\mathit{U}}}_{\mathit{n}}\text{\hspace{0.17em}}{{\displaystyle \stackrel{\u02c6}{S}}}_{n}\text{\hspace{0.17em}}{{{\displaystyle \stackrel{\u02c6}{V}}}^{\prime}}_{n}$$(28)In the above the estimated singular values smaller than the cutoff value ${{\displaystyle \stackrel{\u02c6}{S}}}_{n}$ are all replaced by zero, and only the columns vectors of $\stackrel{\u02c6}{U}$ and $\stackrel{\u02c6}{V}$ corresponding to the included singular values are retained as ${{\displaystyle \widehat{\mathit{U}}}}_{n}$ and ${{\displaystyle \stackrel{\u02c6}{V}}}_{n}$. The model specification problem has been made an approximation problem. The two factorizations of the block matrix ${{\displaystyle \stackrel{\u02c6}{H}}}_{n}$ are $${{\displaystyle \stackrel{\u02c6}{H}}}_{n}={\widehat{\mathit{U}}}_{\mathit{n}}{{\displaystyle \stackrel{\u02c6}{S}}}_{n}\text{}{{{\displaystyle \stackrel{\u02c6}{V}}}^{\prime}}_{n}=({\widehat{\mathit{U}}}_{\mathit{n}}\text{}{{\displaystyle \stackrel{\u02c6}{S}}}_{n}^{\text{1/2}})({{\displaystyle \stackrel{\u02c6}{S}}}_{n}^{\text{1/2}}\text{}{{{\displaystyle \stackrel{\u02c6}{V}}}^{\prime}}_{n})$$(29) $$\text{and\hspace{0.17em}}{{\displaystyle \stackrel{\u02c6}{H}}}_{n}={{\displaystyle \stackrel{\u02c6}{O}}}_{n}\text{\hspace{0.22em}}{{\displaystyle \stackrel{\u02c6}{K}}}_{n}$$(30)We have formed the two factors $${{\displaystyle \stackrel{\u02c6}{O}}}_{n}={\widehat{\mathit{U}}}_{\mathit{n}}\text{\hspace{0.17em}}{{\displaystyle \stackrel{\u02c6}{S}}}_{n}^{\text{1/2}}$$(31) $${{\displaystyle \stackrel{\u02c6}{K}}}_{n}={{\displaystyle \stackrel{\u02c6}{S}}}_{n}^{\text{1/2}}\text{\hspace{0.22em}}{{{\displaystyle \stackrel{\u02c6}{V}}}^{\prime}}_{n}$$(32)with generalised inverses $${{\displaystyle \widehat{\mathit{O}}}}_{n}^{\#}={{\displaystyle \stackrel{\u02c6}{S}}}_{n}^{\text{1/2}}\text{\hspace{0.17em}}{{{\displaystyle \stackrel{\u02c6}{U}}}^{\prime}}_{n}$$(33) $${{\displaystyle \stackrel{\u02c6}{K}}}_{n}^{\#}={{\displaystyle \stackrel{\u02c6}{V}}}_{n}\text{\hspace{0.22em}}{{\displaystyle \stackrel{\u02c6}{S}}}_{n}^{\text{1/2}}$$(34)Other factorisations are possible, however the representation that corresponds to the choice of equations (31) and (32) is called balanced since the gramians are equal, we have ${\stackrel{\u02c6}{O}}_{\mathit{n}}^{\prime}\text{\hspace{0.17em}}{{\displaystyle \stackrel{\u02c6}{O}}}_{n}={{{\displaystyle \stackrel{\u02c6}{K}}}^{\prime}}_{n}\text{\hspace{0.17em}}{{\displaystyle \stackrel{\u02c6}{K}}}_{n}={{\displaystyle \stackrel{\u02c6}{S}}}_{n}$.
The unknown transition matrix A may now be found in terms of the known estimated covariance matrices ${{\displaystyle \stackrel{\u02c6}{R}}}_{i}$ with i = 0, 1, … and the singular value decomposition of ${{\displaystyle \stackrel{\u02c6}{H}}}_{n}$. The procedure is described below. Multiply equation (18) from the right by z′_{k} and take expectation to obtain $$E[{z}_{k+1}\text{\hspace{0.22em}}{{z}^{\prime}}_{k}]=E[(A\text{\hspace{0.22em}}{z}_{k}\text{+}B\text{\hspace{0.22em}}{e}_{k}){{z}^{\prime}}_{k}]=A\text{\hspace{0.22em}}E[{z}_{k}\text{\hspace{0.22em}}{{z}^{\prime}}_{k}]\text{+}E[B({e}_{k}\text{\hspace{0.22em}}{{z}^{\prime}}_{k})]$$(35) Dropping the last term since e_{k} and z_{k} are uncorrelated theoretically, we obtain $$E[{z}_{k+1}\text{\hspace{0.22em}}{{z}^{\prime}}_{k}]=A\text{\hspace{0.22em}}E[{z}_{k}\text{\hspace{0.22em}}{{z}^{\prime}}_{k}]$$(36)From equation (20) we have and we obtain $$E[{z}_{k+1}\text{\hspace{0.22em}}{{z}^{\prime}}_{k}]=K\text{\hspace{0.22em}}{R}_{}^{1}\text{\hspace{0.22em}}{{\displaystyle \overleftarrow{R}}}_{}\text{\hspace{0.22em}}{R}_{}^{1}\text{\hspace{0.22em}}{K}^{\prime}$$(37)where ${{\displaystyle \overleftarrow{R}}}_{}=E[{y}_{k}^{}\text{\hspace{0.22em}}{y}_{k1}^{{}^{\prime}}]$ denotes a shifted version of the block matrix R_{−} $$E[{z}_{k}\text{\hspace{0.22em}}{{z}^{\prime}}_{k}]=K\text{\hspace{0.22em}}{R}_{}^{1}E[{y}_{k1}^{}\text{\hspace{0.22em}}{y}_{k1}^{{}^{\prime}}]{R}_{}^{1}\text{\hspace{0.22em}}{K}^{\prime}=K\text{\hspace{0.22em}}{R}_{}^{1}\text{\hspace{0.22em}}{K}^{\prime}$$(38)From equation (36) we get $$K\text{\hspace{0.22em}}{R}_{}^{1}\text{\hspace{0.22em}}{{\displaystyle \overleftarrow{R}}}_{}\text{\hspace{0.22em}}{R}_{}^{1}\text{\hspace{0.22em}}{K}^{\prime}=A\text{\hspace{0.22em}}K\text{\hspace{0.22em}}{R}_{}^{1}{K}^{\prime}$$(39)and finally the transition matrix is given by $$A=(K\text{\hspace{0.22em}}{R}_{}^{1}\text{\hspace{0.22em}}{{\displaystyle \overleftarrow{R}}}_{}\text{\hspace{0.22em}}{R}_{}^{1}\text{\hspace{0.22em}}{K}^{\prime}){(K\text{\hspace{0.22em}}{R}_{}^{1}\text{\hspace{0.22em}}{K}^{\prime})}^{1}$$(40)Note that the controllability matrix K can be estimated from equation (26) using the generalised inverse of the observability matrix ${{\displaystyle \stackrel{\u02c6}{O}}}_{n}$ $$\stackrel{\u02c6}{K}}={{\displaystyle \widehat{\mathit{O}}}}_{n}^{\#}\text{}{\displaystyle \stackrel{\u02c6}{H}}={{\displaystyle \stackrel{\u02c6}{S}}}_{n}^{\text{1/2}}\text{}{{\displaystyle \stackrel{\u02c6}{U}}}_{n}^{\prime}\text{}{\displaystyle \stackrel{\u02c6}{H}$$(41)With estimates of the transition matrix A in hand we compute its eigen values and we can identify eigen frequencies and damping factors of the vibrating system following equations (12) and (13). However, all the subspace modal identification algorithms have a serious problem of model order determination. When extracting physical or structural modes, subspace algorithms always generate spurious or computational modes to account for unwanted effects such as noise, leakage, residuals, nonlinearity's For these reasons, the assumed number of modes, or model order, is incremented over a wide range of values and we plot the stability diagram [15]. The stability diagram involves tracking the estimates of eigen frequencies and damping factors as a function of model order. As the model order is increased, more and more modal frequencies and damping ratios are estimated, hopefully, the estimates of the physical modal parameters are stabilized using a criterion based on the modal coherence of measured modes and identified modes [15]. Using this criterion we detect and remove the spurious modes. A numerical example and an experimental tests in laboratory showing the effectiveness of the method in modal parameter identification are presented in Section 5.
4 The continuous wavelet transform
4.1 Definition and properties
The continuous wavelet transform (CWT) gives time and frequency information about the analyzed data. For all functions y(t) satisfying the condition ${\int}_{\mathrm{\infty}}^{+\mathrm{\infty}}\mathit{y}(\mathit{t}){}^{2}\mathit{d}\mathit{t}<\mathrm{\infty}$, which implies that y(t) decays to zero at ±∞, the wavelet transform of y(t) is defined as [10–15]: $$({W}_{\psi}y)(a,b)=\frac{1}{\sqrt{a}}{\int}_{\infty}^{+\infty}y(t){\psi}^{\ast}\left(\frac{tb}{a}\right)\hspace{0.17em}dt$$(42) where ψ(t) is an analyzing function called mother wavelet, a is the dilatation or scale parameter defining the analysing window stretching and b is the translation parameter localising the wavelet function in the time domain. The CWT represents the correlation between the signal y(t) and a scaled version of the function ψ(t) and the idea of the CWT is to decompose a signal y(t) into wavelet coefficients using the basis of wavelet functions. The decomposition is obtained locally at different time windows and frequency bands. Once the mother wavelet is chosen, the centre of the time window is controlled by the translation parameter b while the length of the frequency band is controlled by the dilatation parameter a. Hence, one can examine the signal at different time windows and frequency bands by controlling translation and dilatation.
4.2 Modal parameter identification by the CWT
The CWT gives time and frequency information about the analysed data. Any function ψ(t) can be used as an analysing wavelet when it satisfies the admissibility condition [10], this condition being necessary to obtain the inverse of the CWT. There are a number of different complex and real valued functions used as analysing wavelets, but the Morlet wavelet may be the most commonly considered and will be used in this communication. The Morlet wavelet is defined in the time domain as $$\psi (t)={e}^{j\text{\hspace{0.17em}}{\omega}_{o}t}{e}^{{t}^{2}\text{/2}}$$(43) where ω_{o} is the oscillating frequency of the analysing wavelet. The Morlet wavelet is basically a sinusoidal function, oscillating at the frequency ω_{o} and modulated by a Gaussian envelope of unit width. In the frequency domain the dilated Morlet wavelet becomes $$\psi (a\omega )=\sqrt{2\pi}{e}^{\frac{12}{}{(a\omega {\omega}_{o})}^{2}}$$(44)and the maximum is located at aω = ω_{o}.
Consider a modulated signal in amplitude and frequency: x(t) = A(t)cos(ϕ(t)). The CWT of x(t) is [10–13] $$({W}_{\psi}x)(a,b)=\frac{\sqrt{a}}{2}A(b){\mathrm{\Psi}}^{\ast}(a{\displaystyle \dot{\varphi}}(b)){e}^{j\varphi (b)}$$(45) where $\dot{\varphi}}(b)=d\varphi (b)/db$ is the instantaneous frequency of the signal. The dilatation parameter is obtained when (W_{ψ}x)(a, b) is maximum, that is $a=a(b)={\omega}_{o}/{\displaystyle \dot{\varphi}}(b)$ and we get then $$({W}_{\psi}x)(a(b),b)=\sqrt{\pi a(b)/2}A(b){e}^{j\varphi (b)}$$(46)Now, consider the impulse response of a viscously damped single degree of freedom system $$x(t)=B{e}^{\zeta {\omega}_{n}t}\mathrm{cos}({\omega}_{d}t+{\chi}_{o})$$(47)where ${\omega}_{d}={\omega}_{n}\sqrt{1{\zeta}^{2}}$ is the damped angular frequency and ξ the damping factor. It is easy to see that A(t) = Be^{−ζωnt} and ϕ(t) = ω_{d}t + χ_{o} so $\dot{\varphi}}(t)={\omega}_{d$. The amplitude of the wavelet transform is maximum for the constant scale $a={a}_{1}={\omega}_{o}/{\omega}_{d}$ and we obtain the ridge of the CWT. Equation (46) becomes $$({W}_{\psi}x)({a}_{1},b)=\sqrt{\pi {a}_{1}/2}B{e}^{\zeta {\omega}_{n}b}{e}^{j({\omega}_{d}b+{\chi}_{o})}$$(48)From equation (48) we can deduce $$\text{Arg}(({W}_{\psi}x)({a}_{1},b))={\omega}_{d}b+{\chi}_{o}$$(49) $$\mathrm{ln}(({W}_{\psi}x)({a}_{1},b))=\xi {\omega}_{n}b+\text{constant}$$(50)The phase function is a straight line whose slope represents the damped angular frequency. The log of the wavelet modulus crosssection is again a straight line whose slope is the decay rate ξω_{n}. Once the instantaneous frequency and the decay rate have been estimated we can identify the natural frequency and the damping ratio of the vibrating system.
5 Applications
5.1 Simulated data
To prove the effectiveness of the identification procedure based on the subspace analysis we consider a twoDOF system with very closely spaced modes. The parameters of the system are f_{1} = 2 MHz, f = 2.01 MHz, ξ_{1} = 0.01 and ξ_{2} = 0.02. Figure 2a shows the free response of the system where a Gaussian white noise has been added: the generated data were corrupted by a random noise. The sampling frequency is 5 MHz and only 151 time samples are used in the simulation. The signal duration is 30 μs.
We convert to the frequency domain this time response by taking the discrete Fourier transform of the noisy signal. It is impossible to identify the two frequencies components by using the FFT as shown in Figure 2b and c, where the power spectral density has been plotted.
In our identification procedure by subspaces, we plot the stabilization diagram on eigen frequencies and damping factors. Figure 3 shows stabilization diagrams with the modal coherence indicator: spurious modes have been eliminated and only physical modes are present.
Our procedure can separate closely spaced modes and the mean values on identified eigen frequencies and damping factors obtained by an average over the orders of the stabilization diagram are ${{\displaystyle \stackrel{\u02c6}{f}}}_{1}=2\text{\hspace{0.17em}MHz}$; ${{\displaystyle \stackrel{\u02c6}{f}}}_{2}=2.01\text{\hspace{0.17em}MHz}$; ${{\displaystyle \stackrel{\u02c6}{\xi}}}_{1}=0.01$; ${{\displaystyle \stackrel{\u02c6}{\xi}}}_{2}=0.02$. A very satisfactory estimation on eigen frequencies and damping factors has been obtained using simulated data.
We apply now the continuous wavelet transform to our numerical test. Figure 4 presents the continuous wavelet transform modulus and its ridge for the twoDOF system. The CWT is a timefrequency representation and the two modes appear after 18 μs of signal recording. The identified first mode is plot in red dash line and the identified second mode is plot in magenta dotdashline. The two vertical blue lines correspond to the limitations of edge effects. The identified modal parameters for the two dof system are ${{\displaystyle \stackrel{\u02c6}{f}}}_{1}=1.99\text{\hspace{0.17em}MHz}$; ${{\displaystyle \stackrel{\u02c6}{f}}}_{2}=2.01\text{\hspace{0.17em}MHz}$; ${{\displaystyle \stackrel{\u02c6}{\xi}}}_{1}=0.04$; ${{\displaystyle \stackrel{\u02c6}{\xi}}}_{2}=0.49$. A very satisfactory estimation on eigen frequencies are obtained, however damping factors are not well identified. This is due to the fact that the two modes are very close and only 151 time samples are used in the analysis. It is shown in [14] that it is necessary to have a long recording time and an important number of time samples to identify correctly closely spaced modes. Furthermore, the vibrating system is very quickly attenuated: the time response decreases rapidly and, it is shown in [14], the continuous wavelet transform gives good results if the time response is not attenuated abruptly. The limitations of the continuous wavelet transform concern the low number of time samples used in the identification procedure and the fast decay of the time response.
An experimental test in laboratory is presented in the following section.
Fig. 2 Time response for the analyzed twoDOF system (a) and its spectrum (b) and (c). 
Fig. 3 Stabilization diagram on eigen frequencies and damping factors for the twoDOF system. 
Fig. 4 Continuous wavelet transform modulus and its ridge for the twoDOF system. 
5.2 Modal parameter identification of a CMUT membrane
A CMUT cell is basically a metallized membrane suspended over an electrode and by applying an alternating voltage to the parallel plate capacitor, an electrostatic force will cause the deflection of the membrane and the generation of ultrasounds. The main purpose is to find the first resonance frequency of the membrane, the effectiveness of a CMUT cell being analyzed in terms of resonance frequency, stiffness and damping coefficients. It can be reminded that vibration amplitude of the membrane is maximized at the mechanical resonance frequency. A CMUT is commonly fabricated by means of the surface micromachining technology, using standard integrated circuits techniques. Several processes have been reported in the literature to fabricate CMUTs, using different materials and film deposition techniques. In our experimental test the CMUTs have been fabricated using the anodic bonding technology of a fixed thickness monocrystalline silicon layer of a SOI wafer (Silicon On Insulator technology) on a borosilicate glass substrate. The process evolution for each technological step and the fabrication runs are not discussed, but details can be found in Bellaredj [17]. Two photographs representing a circular CMUT cell and an array of CMUTs are presented in Figure 5. The bonding technology is a particularly promising technique for CMUTs fabrication because it is easy and reliable to control and reduces process cost and time.
The dynamic measurements of a metallized membrane are conducted in the time domain by means of a POLYTEC laser vibrometer. The vibrometer is composed of a fibercoupled vibrometer sensor head OFV534 and an OFV5000 controller. The controller has a DD300 decoder being able to measure displacements up to 150 nm peak to peak at 24 MHz. In our experimental test the external force acting in the membrane is a step force and Figure 6 shows the time response of the membrane center, where the sampling frequency of signal is 500 MHz and the duration is 18 μs. Only a part of this time response of the structure is used in the identification procedure (the damped response).
In this paper a single degree of freedom of a springmassdamper model is used to study the membrane behavior. The model is constituted by the following parameters: the metallized membrane mass m concentrated in its center, the damping coefficient c and the stiffness coefficient k. Table 1 shows the dimensions and material properties of the CMUT cell and the first theoretical resonance frequency of the metallized membrane can be computed by the relation given in Rossi [18] $$f=\frac{0.467\text{\hspace{0.17em}}{t}_{m}}{{R}^{2}}\sqrt{\frac{E}{\rho (1{v}^{2})}}$$(51) The metallized membrane mass is m = ρΠR^{2}t_{m} = 6.61 × 10^{‑11} kg. The membrane stiffness is given by k = m(2πf)^{2}, the damping coefficient is c = 4πmfξ, the quality factor is Q = 1/2ξ and the halfpower bandwidth of the system is ΔB = f/Q, where f is the resonance frequency of the membrane and ξ the damping factor. These two modal parameters are obtained by an average over the orders of the stabilization diagram presented in Figure 7. To our knowledge, it is the first time that the subspace algorithm is applied to a microelectromechanical system and in particular to a CMUT cell.
Figure 8 shows the modulus of the CWT and the identified ridge obtained from the membrane center displacement. From the ridge of the CWT we obtain Figure 9 and the instantaneous modal parameters of the membrane using the slope of straight lines. To our knowledge, it is the first time that the CWT algorithm is applied to the membrane of a CMUT cell.
Tables 2 and 3 show the identified metallized membrane parameters using respectively the subspace method and the CWT method. These two methods give very similar results.
Figure 10 shows the comparison between the measured time response of the metallized membrane and the reconstructed time response obtained using the identified modal parameters obtained by the subspace method. In this time response we have only taken in consideration the damped response portion between 6.5 and 20 μs.
An important part of the identification process is to assess the quality of the estimated model obtained from the identified parameters. A measure of model quality is obtained using the following expression [19] $$Mq=1\sqrt{\frac{{\displaystyle \sum _{k=1}^{N}}{(yk{\displaystyle \stackrel{\u02c6}{y}}k)}^{2}}{{\displaystyle \sum _{k=1}^{N}}{y}_{k}^{2}}}$$(52) which is a normalized estimation error of the used model. Increasing values of Mq indicate better fit and 1 indicates a perfect fit of the model to the measured data. In this expression y_{k} are the measured data and ${{\displaystyle \stackrel{\u02c6}{y}}}_{k}$ the estimated data, obtained from the model. For the metallized membrane we obtain Mq = 0.67 which is a satisfactory value. The decay function of the response is expressed by the relation g(t) = G_{m}exp(− βt), where G_{m} is obtained by initial conditions and β = 2πfξ is computed from the identified modal parameters. This decay function is plotted in Figure 7. Note that the coefficients G_{m} and β = c/(2m) can also be evaluated directly by interpolating the measured response of the membrane and in this case we can obtain the value of the damping coefficient c from β and m. In [20] nonlinearities are added in the model to analysis the performances of the CMUT for several DC voltages.
Fig. 5 A circular CMUT cell and an array of CMUTs. 
Fig. 6 Time response of the membrane center to a step force actuation. 
Dimensions and properties of the CMUT cell.
Fig. 7 Stabilization diagram on resonance frequency and damping factor for the experimental CMUT membrane. 
Fig. 8 Continuous wavelet transform modulus and its ridge for the experimental CMUT membrane center. 
Fig. 9 Amplitude and phase of the ridge of the CWT for the experimental CMUT membrane center. 
Identified parameters of the metallized membrane using the subspace method.
Identified parameters of the metallized membrane using the continuous wavelet transform method.
Fig. 10 Comparison between the measured time response of the membrane (in blue) and the reconstituted time response (in red). 
6 Conclusion
In modal parameter identification the key point is to determine the relationship between the system parameters and the measured dynamical data and a particular attention has been paid to modal parameter identification when only output measurements are used. Such parameters can then be used to improve the ultrasonic generating process by a CMUT, for model validation and for fault detection in CMUTs. In this paper, we have proposed two methods for modal parameter identification: the subspace method based in the statespace innovation model and the continuous wavelet transform method based in the use of the ridge of the CWT. Numerical and experimental results have shown the effectiveness of the procedures presented in the paper for modal parameter identification. However, the CWT has some limitations: if the number of data points is insufficient or if the vibrating system is quickly attenuated closely modes are not well identified in terms of damping factors.
References
 B. Bayram, O. Oralkan, A.S. Ergun, E. Hæggstrom, G.G. Yaralioglu, B.T. KhuriYakub, Capacitive micromachined ultrasonic transducer design for high power transmission, IEEE Trans.Ultrason. Ferroelectr. Freq. Control 52 (2005) 326–338 [CrossRef] [Google Scholar]
 G.G. Yaralioglu, B.T. KhuriYakub, Finite element analysis of capacitive micromachined ultrasonic transducers, IEEE Trans Ultrason. Ferroelectr. Freq. Control 52 (2005) 2185–2198 [CrossRef] [Google Scholar]
 F.V. Hunt, The analysis of transduction and its historical background, Harvard University Press, 1982 [Google Scholar]
 M. Berthillier, P. Le Moal, J. Lardiès, Dynamic and acoustic modelization of capacitive micromachined ultrasonic transducers, IEEE International Ultrasonics Symposium, Orlando, USA, 2011. [Google Scholar]
 B.L. Ho, R.E. Kalman, Effective construction of linear statevariable models from input/output data, Proceedings of the 3rd Annual Allerton Conference on Circuit and System Theory, University of Illinois, USA, 1965, pp. 449–459 [Google Scholar]
 H.P. Zeiger, A.J. McEwen, Approximate linear realizations of given dimension via Ho's algorithm, IEEE. Trans. Automatic Control 19, 153, 1974 [CrossRef] [Google Scholar]
 J.N. Juang, Applied system identification, Prentice Hall, 1994 [Google Scholar]
 G.H. James, T.G. Carne, J.P. Lauffer, The natural excitation technique for modal parameter extraction from operating structures, Int. J. Anal. Exp. Modal Anal. 10 (1995) 260–277 [Google Scholar]
 P. Van Overschee, B. De Moor, Subspace identification for linear systems: thearyimplementationapplications, Kluwer Academic Publishers, 1996 [CrossRef] [Google Scholar]
 B. Torrésani, «Analyse continue par ondelettes», Inter Editions/CNRS, 1995 [Google Scholar]
 M. Ruzzene, A. Fasana, L. Garibaldi, B. Piombo, «Natural frequencies and dampings identification using wavelet transform: application to real data», Mech. Syst. Signal Process. 11 (1997) 207–218 [CrossRef] [Google Scholar]
 J. Lardiès, S. Gouttebroze, Identification of modal parameters using the wavelet transform, Int. J. Mech. Sci. 44 (2002) 2263–2283 [CrossRef] [Google Scholar]
 P. Argoul, T.P. Le, Instantaneous indicators of structural behaviour based on the continuous Cauchy wavelet analysis, Mech. Syst. Signal Process. 17 (2003) 243–250 [Google Scholar]
 T.P. Le, P. Argoul, Continuous wavelet transform for modal identification using free decay response, J. Sound Vib. 277 (2004) 73–100 [CrossRef] [Google Scholar]
 J. Lardiès, M. Ta, Modal parameter identification of stay cables from outputonly measurements, Mech. Syst. Signal Process. 25 (2011) 133–150 [CrossRef] [Google Scholar]
 M. Aoki, State space modeling of time series, SpringerVerlag, Berlin, Heidelberg, 1987 [CrossRef] [Google Scholar]
 M.L.F. Bellaredj, Methods and tools for the fabrication of silicon micromachined ultrasonic transducers, PHD thesis, FrancheComté University, Besançon, France, 2013 [Google Scholar]
 M. Rossi, Audio Lausanne, Presses Polytechniques et Universitaires Romandes, 2007 [Google Scholar]
 J. Lardiès, Modal parameter identification of perforated microplates from output data only, The International Nanotech&Nano Science Conference, France, Paris, 2015 [Google Scholar]
 A. Jallouli, N. Kacem, G. Bourbon, P. Le Moal, V. Walter, J. Lardiès, Pullin instability tuning in imperfect nonlinear circular microplates under electrostatic actuation, Phys. Lett. A 380 (2016) 3886–3890 [CrossRef] [Google Scholar]
Cite this article as: J. Lardiès, G. Bourbon, P.L. Moal, N. Kacem, V. Walter, T.P. Le, Modal parameter identification of a CMUT membrane using response data only, Mechanics & Industry 18, 802 (2017)
All Tables
Identified parameters of the metallized membrane using the continuous wavelet transform method.
All Figures
Fig. 1 A conceptual crosssection of a capacitive micromachined ultrasonic transducer. 

In the text 
Fig. 2 Time response for the analyzed twoDOF system (a) and its spectrum (b) and (c). 

In the text 
Fig. 3 Stabilization diagram on eigen frequencies and damping factors for the twoDOF system. 

In the text 
Fig. 4 Continuous wavelet transform modulus and its ridge for the twoDOF system. 

In the text 
Fig. 5 A circular CMUT cell and an array of CMUTs. 

In the text 
Fig. 6 Time response of the membrane center to a step force actuation. 

In the text 
Fig. 7 Stabilization diagram on resonance frequency and damping factor for the experimental CMUT membrane. 

In the text 
Fig. 8 Continuous wavelet transform modulus and its ridge for the experimental CMUT membrane center. 

In the text 
Fig. 9 Amplitude and phase of the ridge of the CWT for the experimental CMUT membrane center. 

In the text 
Fig. 10 Comparison between the measured time response of the membrane (in blue) and the reconstituted time response (in red). 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.