© 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 on-chip 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 cross-section 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 state-space 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 input-output 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 domain-based stochastic subspace identification (SSI) method has been proposed by Van Overschee and De Moor [9] and is based on the state-space 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 time-frequency 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 half-power bandwidth. Such coefficients are identified using two methods: the subspace method derived from a state-space 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 state-space 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 cross-section 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 semi-definite 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;(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 state-space model: $${\displaystyle \dot{x}}(t)={\displaystyle \tilde{A}}x(t)+{\displaystyle \tilde{B}}\eta (t)$$(4) $$y(t)=Cx(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 ${\displaystyle \tilde{A}}\text{,\hspace{0.17em}}{\displaystyle \tilde{B}}\text{,\hspace{0.17em}}C$ are the (2n′ × 2n′), (2n′ × n′), (m × 2n′) matrices $${\displaystyle \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];{\displaystyle \tilde{B}}=\left[\begin{array}{c}\hfill 0\hfill \\ \hfill {M_o}^{-1}\hfill \end{array}\right];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 state-space model and the discrete time observation equation [7,15] $$x_{k\text{+1}}=A{x}_k\text{+}{u}_k$$(8) $$y_k=C{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 state-space matrix or discrete time transition matrix and u_k is given by $$u_k={\int }_0^{\mathrm{\Delta }t}{e}^{{\displaystyle \tilde{A}}s}{\displaystyle \tilde{B}}\eta (t-s)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{{[\ln ({\lambda }_i{{\lambda }_i}^{*})]}^2}{4}+{\left[{\cos }^{-1}\left(\frac{{\lambda }_i+{{\lambda }_i}^{*}}{2\sqrt{{\lambda }_i{{\lambda }_i}^{*}}}\right)\right]}^2}$$(12) $${\zeta }_i=\sqrt{\frac{{[\ln ({\lambda }_i{{\lambda }_i}^{*})]}^2}{{[\ln ({\lambda }_i{{\lambda }_i}^{*})]}^2+4{\left[{\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 state-space innovation form, before we can obtain the transition matrix A.

3 State-space innovation model

We estimate x_k by its orthogonal projection onto the subspace spanned by the stacked vector of past data vectors noted $y_{k-1}^{-}=[{y^{\prime }}_{k-1},{y^{\prime }}_{k-2},\dots {]}^{\prime }$. Here, the superscript ′ denotes the transpose operation and we use the notation $z_k=E[x_k|y_{k-1}^{-}]$ 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}}|y_k^{-}]=E[x_{k\text{+1}}|y_{k-1}^{-},e_k]=AE[x_k|y_{k-1}^{-}]+B{e}_k=A{z}_k+B{e}_k$$(14) where $e_k=y_{k}E[x_{k+1}|y_{k-1}^{-}]$ 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{e^{\prime }}_k]$. Note that e_k is uncorrelated with $y_{k-1}^{-}$ 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}|e_k]=E[x_{k\text{+1}}{e^{\prime }}_k](E[e_k{e^{\prime }}_k{])}^{-1}{e}_k$$(15)we have: $$B=E[x_{k\text{+1}}{e^{\prime }}_k](E[e_k{e^{\prime }}_k{])}^{-1}$$(16)Note that by assumption $E[u_k|y_{k-1}^{-}]=0$ and $E[v_k|y_{k-1}^{-}]=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{x}_k\text{+}{v}_k|y_{k-1}^{-}]+e_k=C{z}_k+e_k$$(17)The set of equations (8) and (9) is now replaced by $$z_{k+1}=A{z}_k+B{e}_k$$(18) $$y_k=C{z}_k+e_k$$(19)Such representation is called state-space innovation representation of the data generating process equations (8) and (9). The state-space 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_{k-1}^{-}$ by $R=E[y_{k-1}^{-}{y}_{k-1}^{-\prime }]$. This matrix is theoretically infinite-dimensional, 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 \hat{R}}}_{-}$ is (mp × mp). We develop the state vector z_k $$z_k=E[x_k|y_{k-1}^{-}]=E[x_k{y}_{k-1}^{-\prime }]{(E[y_{k-1}^{-}{y}_{k-1}^{-\prime }])}^{-1}{y}_{k-1}^{-}=K{R}_{-}^{-1}{y}_{k-1}^{-}$$(20) where $K=E[x_k{y}_{k-1}^{-\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{y}_{k-1}^{-\prime }]=E[z_k{y}_{k-1}^{-\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{y}_{k-1}^{-\prime }]$.

To estimate K it is necessary to introduce the stacked vector of future data vectors noted $y_k^{+}=[{y^{\prime }}_k,{y^{\prime }}_{k\text{+1}},\dots {]}^{\prime }=[{y^{\prime }}_k,y_{k+1}^{{\text{+}}^{\prime }}{]}^{\prime }$, then from equations (19) and (20) we obtain $$y_k^{+}=O{z}_k+F{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 ⋮\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{+}}{y}_{k-1}^{-\prime }]=OE[z_k{y}_{k-1}^{{-}^{\prime }}]=OK$$(24) The matrix H, called the block Hankel matrix, is formed of (m × m) covariance matrices $R_i=E[y_{k\text{+}i}{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 \hat{R}}}_i=T^{-1}{\displaystyle \sum _{k=1}^{T-i}}{y}_{k+i}{y^{\prime }}_k,i=1,2,\dots$$(25)are used to construct a sample block Hankel matrix ${\displaystyle \hat{H}}$ which is factored into the estimated observability and controllability matrices ${\displaystyle \hat{O}}$ and ${\displaystyle \hat{K}}$. $${\displaystyle \hat{H}}={\displaystyle \hat{O}}{\displaystyle \hat{K}}$$(26)A second factorisation of the estimated block Hankel matrix is obtained using the singular value decomposition of ${\displaystyle \hat{H}}$ $${\displaystyle \hat{H}}={\displaystyle \hat{U}}{\displaystyle \hat{S}}{\displaystyle \hat{V}}$$(27)where ${\displaystyle \hat{U}}$ and ${\displaystyle \hat{V}}$ contain the singular vectors of ${\displaystyle \hat{H}}$. The singular values of ${\displaystyle \hat{H}}$ are contained in descending order in the diagonal matrix ${\displaystyle \hat{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 \hat{H}}}_n$. The S.V.D. of ${{\displaystyle \hat{H}}}_n$ is $${{\displaystyle \hat{H}}}_n={\widehat{\mathit{U}}}_{\mathit{n}}{{\displaystyle \hat{S}}}_n{{{\displaystyle \hat{V}}}^{\prime }}_n$$(28)In the above the estimated singular values smaller than the cut-off value ${{\displaystyle \hat{S}}}_n$ are all replaced by zero, and only the columns vectors of ${\displaystyle \hat{U}}$ and ${\displaystyle \hat{V}}$ corresponding to the included singular values are retained as ${{\displaystyle \widehat{\mathit{U}}}}_n$ and ${{\displaystyle \hat{V}}}_n$. The model specification problem has been made an approximation problem. The two factorizations of the block matrix ${{\displaystyle \hat{H}}}_n$ are $${{\displaystyle \hat{H}}}_n={\widehat{\mathit{U}}}_{\mathit{n}}{{\displaystyle \hat{S}}}_n{{{\displaystyle \hat{V}}}^{\prime }}_n=({\widehat{\mathit{U}}}_{\mathit{n}}{{\displaystyle \hat{S}}}_n^{\text{1/2}})({{\displaystyle \hat{S}}}_n^{\text{1/2}}{{{\displaystyle \hat{V}}}^{\prime }}_n)$$(29) $$\text{and\hspace{0.17em}}{{\displaystyle \hat{H}}}_n={{\displaystyle \hat{O}}}_n{{\displaystyle \hat{K}}}_n$$(30)We have formed the two factors $${{\displaystyle \hat{O}}}_n={\widehat{\mathit{U}}}_{\mathit{n}}{{\displaystyle \hat{S}}}_n^{\text{1/2}}$$(31) $${{\displaystyle \hat{K}}}_n={{\displaystyle \hat{S}}}_n^{\text{1/2}}{{{\displaystyle \hat{V}}}^{\prime }}_n$$(32)with generalised inverses $${{\displaystyle \widehat{\mathit{O}}}}_n^{\#}={{\displaystyle \hat{S}}}_n^{-\text{1/2}}{{{\displaystyle \hat{U}}}^{\prime }}_n$$(33) $${{\displaystyle \hat{K}}}_n^{\#}={{\displaystyle \hat{V}}}_n{{\displaystyle \hat{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 ${\hat{O}}_{\mathit{n}}^{\prime }{{\displaystyle \hat{O}}}_n={{{\displaystyle \hat{K}}}^{\prime }}_n{{\displaystyle \hat{K}}}_n={{\displaystyle \hat{S}}}_n$.

The unknown transition matrix A may now be found in terms of the known estimated covariance matrices ${{\displaystyle \hat{R}}}_i$ with i = 0, 1, … and the singular value decomposition of ${{\displaystyle \hat{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}{z^{\prime }}_k]=E[(A{z}_k\text{+}B{e}_k){z^{\prime }}_k]=AE[z_k{z^{\prime }}_k]\text{+}E[B(e_k{z^{\prime }}_k)]$$(35) Dropping the last term since e_k and z_k are uncorrelated theoretically, we obtain $$E[z_{k+1}{z^{\prime }}_k]=AE[z_k{z^{\prime }}_k]$$(36)From equation (20) we have and we obtain $$E[z_{k+1}{z^{\prime }}_k]=K{R}_{-}^{-1}{{\displaystyle \overleftarrow{R}}}_{-}{R}_{-}^{-1}{K}^{\prime }$$(37)where ${{\displaystyle \overleftarrow{R}}}_{-}=E[y_k^{-}{y}_{k-1}^{{-}^{\prime }}]$ denotes a shifted version of the block matrix R₋ $$E[z_k{z^{\prime }}_k]=K{R}_{-}^{-1}E[y_{k-1}^{-}{y}_{k-1}^{{-}^{\prime }}]R_{-}^{-1}{K}^{\prime }=K{R}_{-}^{-1}{K}^{\prime }$$(38)From equation (36) we get $$K{R}_{-}^{-1}{{\displaystyle \overleftarrow{R}}}_{-}{R}_{-}^{-1}{K}^{\prime }=AK{R}_{-}^{-1}{K}^{\prime }$$(39)and finally the transition matrix is given by $$A=(K{R}_{-}^{-1}{{\displaystyle \overleftarrow{R}}}_{-}{R}_{-}^{-1}{K}^{\prime }){(K{R}_{-}^{-1}{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 \hat{O}}}_n$ $${\displaystyle \hat{K}}={{\displaystyle \widehat{\mathit{O}}}}_n^{\#}{\displaystyle \hat{H}}={{\displaystyle \hat{S}}}_n^{-\text{1/2}}{{\displaystyle \hat{U}}}_n^{\prime }{\displaystyle \hat{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{t-b}{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{\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 ${\displaystyle \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}\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) = ω_dt + χ_o so ${\displaystyle \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) $$\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 cross-section 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 two-DOF system with very closely spaced modes. The parameters of the system are f₁ = 2 MHz, f  = 2.01 MHz, ξ₁ = 0.01 and ξ₂ = 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 \hat{f}}}_1=2\text{\hspace{0.17em}MHz}$; ${{\displaystyle \hat{f}}}_2=2.01\text{\hspace{0.17em}MHz}$; ${{\displaystyle \hat{\xi }}}_1=0.01$; ${{\displaystyle \hat{\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 two-DOF system. The CWT is a time-frequency 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 dot-dash-line. The two vertical blue lines correspond to the limitations of edge effects. The identified modal parameters for the two dof system are ${{\displaystyle \hat{f}}}_1=1.99\text{\hspace{0.17em}MHz}$; ${{\displaystyle \hat{f}}}_2=2.01\text{\hspace{0.17em}MHz}$; ${{\displaystyle \hat{\xi }}}_1=0.04$; ${{\displaystyle \hat{\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 two-DOF system (a) and its spectrum (b) and (c).

Fig. 3 Stabilization diagram on eigen frequencies and damping factors for the two-DOF system.

Fig. 4 Continuous wavelet transform modulus and its ridge for the two-DOF 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 fiber-coupled vibrometer sensor head OFV-534 and an OFV-5000 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 spring-mass-damper 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.467t_m}{R^2}\sqrt{\frac{E}{\rho (1-v^2)}}$$(51) The metallized membrane mass is m = ρΠR²t_m = 6.61 × 10^‑11 kg. The membrane stiffness is given by k = m(2πf)², the damping coefficient is c = 4πmfξ, the quality factor is Q = 1/2ξ and the half-power 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 \hat{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 \hat{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_mexp(− β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.

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.

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 state-space 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

All Tables

All Figures

	Fig. 1 A conceptual cross-section of a capacitive micromachined ultrasonic transducer.
In the text

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

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

	Fig. 4 Continuous wavelet transform modulus and its ridge for the two-DOF 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