Rapid High-Resolution Wavenumber Extraction from Ultrasonic Guided Waves Using Adaptive Array Signal Processing

Abstract

Quantitative ultrasound techniques for assessment of bone quality have been attracting significant research attention. The axial transmission technique, which involves analysis of ultrasonic guided waves propagating along cortical bone, has been proposed for assessment of cortical bone quality. Because the frequency-dependent wavenumbers reflect the elastic parameters of the medium, high-resolution estimation of the wavenumbers is required at each frequency with low computational cost. We use an adaptive array signal processing method and propose a technique that can be used to estimate the numbers of propagation modes that exist at each frequency without the need for time-consuming calculations. An experimental study of 4-mm-thick copper and bone-mimicking plates showed that the proposed method estimated the wavenumbers accurately with estimation errors of less than 4% and a calculation time of less than 0.5 s when using a laptop computer.

1. Introduction

The first quantitative ultrasound (QUS) technique for use in bone assessment was proposed by Langton et al. [1] in the 1980s. Since then, a number of different QUS techniques have been proposed and developed [2,3,4]. QUS techniques were developed for early osteoporosis detection and fracture risk evaluation screening applications. When compared with X-rays, QUS offers several advantages: it is a noninvasive process and can be performed using portable and low-cost equipment. Techniques that can evaluate the quality of cortical long bone structures have recently attracted significant research attention [2,3,5,6,7,8]. In this study, we have developed cortical long bone quality evaluation methods using the axial transmission (AT) technique [9,10,11,12,13].

A number of recent studies have used the AT technique, in which wide-band signals are emitted and the ultrasonic guided wave that propagates along the cortical bone is then analyzed. The cortical bone has previously been described as a transversely isotropic absorbing plate with finite thickness [14,15,16,17,18]. The ultrasonic guided wave that propagates in cortical bone consists of multiple propagation modes, and the frequency-dependent wavenumbers of these modes represent the elastic properties of the medium. Many techniques for frequency-wavenumber (f–k) analysis have therefore been proposed and reported.

Tran et al. proposed the Radon transform method, which estimates the phase velocity ($c_{\mathrm{p}}$) at each frequency [12]. The phase velocity is given by $c_{\mathrm{p}}=2\pi f/k$, and the basic premises of the f–k and f–$c_{\mathrm{p}}$ analysis methods are theoretically the same. The Radon transform method uses a single transmitter with multiple receivers, estimates the phase velocity using an iterative process, and produces high-resolution estimates. However, the computational cost of this method is not low because it requires inversion of a large-scale matrix. Sasso et al. [19] proposed a singular value decomposition (SVD)-based method and used a probe composed of multiple transmitters and emitters to analyze the guided waves. Minonzio et al. extended this method and produced accurate depictions of the phase velocities of the plate [10]. Use of SVD allows the method to detect the modes with low intensity. The method was effective for estimation of the elastic properties of the medium [20]. However, the measurement resolution is determined by the aperture size, and it is necessary to estimate the number of propagation modes in the received signal. Xu et al. proposed the sparse SVD (S-SVD) method, which combines the Radon transform method proposed by Tran et al. with the SVD method [21]. The resulting method acquires super-resolution estimates using a combination of SVD and an iterative process. However, this method also requires estimation of the number of propagation modes in the received signal and multiple calculations for inversion of a large matrix.

We recently developed an adaptive beamforming technique [22]. We used the estimation of signal parameters via rotational invariance technique (ESPRIT) algorithm. This algorithm uses eigenvalue decomposition to separate the desired signal from noise. The ESPRIT algorithm also requires estimation of the number of propagation modes and the conventional method includes an iterative process with eigenvalue decomposition [22]. Because the matrix inversion and eigenvalue decomposition processes incur large computational costs, reduction of these computational costs would be useful for practical applications. Therefore, in this study, we propose a method that estimates the number of propagation modes with a low computational cost.

The number of propagation modes in the received signal, which is denoted by M in this study, is usually estimated using a thresholding process. However, the eigenvalues that correspond to the signal are not equivalent to the signal intensities. Therefore, it is not easy to determine M when using a thresholding process that only uses the intensities of the eigenvalues or singular values. Therefore, in this study, we propose a new algorithm to estimate M using information theoretic criteria. This estimation procedure would be effective for both the SVD and ESPRIT methods. While many studies on estimation of M have been reported [23,24,25,26,27], the resulting methods were not applied to the analysis of ultrasonic guided waves when propagating along a transversely isotropic absorbing material such as cortical bone.

In the proposed method, we do not use an iterative process with high computational cost, such as the eigenvalue decomposition or matrix inversion methods, but use a diagonal loading (DL) technique that adds the diagonal loading matrix to the covariance matrix to determine M. We compare the computational cost of the proposed method with that of a conventional method and demonstrate the effectiveness of the proposed method via simple numerical simulations and experiments using a copper plate and a bone-mimicking plate.

2. Materials and Methods

2.1. System Model

In this study, we considered the system model that is shown in Figure 1. We placed a linear array probe on top of a free isotropic or transversely isotropic plate and analyzed the guided waves that propagate along the plate. When we consider multiple propagation signals with multiple phase velocities and ignore changes in amplitude, the received signals are given by

$$S_n\left(\omega \right)=S_1\left(\omega \right)\sum _{m=1}^{M\left(\omega \right)}exp\{-ȷ(n-1)k_m\left(\omega \right)x\}$$

(1)

where $S_n\left(\omega \right)$ is the received signal at the n-th receiver in the frequency domain with an angular frequency $\omega$ and x is the receiving array pitch. $k_m$ is the wavenumber of the m-th propagation mode, which is given by $k_m=\omega /c_m$, where $c_m$ is the phase velocity of the m-th mode. Note that M is dependent on the frequency.

2.2. Wavenumber Estimation Using the ESPRIT Algorithm

To realize super-resolution estimates, we used the ESPRIT algorithm. The basis of this algorithm is briefly explained here. When there is a single propagating wave, i.e., $M=1$ in Equation (1), we can then estimate the wavenumber by comparing the signals that are received at two receivers.

$$k_1\left(\omega \right)=\angle \left\{S_n\left(\omega \right)S_{n+1}^{*}\left(\omega \right)\right\}/x$$

(2)

where $\angle$ represents the angle of the complex signal.

As shown in Equation (2), the wavenumber is estimated directly without performing a peak search process. However, when $M>1$, we cannot estimate the wavenumber by simply comparing the phases of the two received signals alone. The ESPRIT algorithm estimates the wavenumbers of multiple modes via an eigenvalue decomposition of the covariance matrix to enable the separation of multiple signals.

The covariance matrix represents a correlation between the signals at each receiver. To estimate the covariance matrix, we use the sub-array averaging technique [28]. In the sub-array averaging technique, multiple sub-arrays are situated into the full-size array. Here, we define the signal matrix as follows:

$$\mathbf{A}\left(\omega \right)=\left(\begin{array}{ccc}{S}_{1,1}\left(\omega \right)& \cdots & S_{N,1}\left(\omega \right)\\ ⋮& \ddots & ⋮\\ S_{1,N_{\mathrm{sub}}}\left(\omega \right)& \cdots & S_{N,N_{\mathrm{sub}}}\left(\omega \right)\end{array}\right)$$

(3)

where $S_{i,j}$ is the received signal at the j-th receiver in the i-th sub-array, N is the number of sub-arrays contained within the whole array, $N_{\mathrm{sub}}$ is the number of receivers included within a single sub-array, i.e., the number of whole receivers is $N_{\mathrm{R}}=N+N_{\mathrm{sub}}-1$, and $\mathbf{A}$ is the signal matrix with size $N_{\mathrm{sub}}\times N$.

The covariance matrix, $\mathbf{R}$, is then given by

$$\begin{array}{ccc}\hfill \mathbf{R}\left(\omega \right)& =& \mathbf{A}\left(\omega \right){\mathbf{A}}^{\mathrm{H}}\left(\omega \right)\hfill \end{array}$$

(4)

$$\begin{array}{ccc}\hfill & =& \sum _{i=1}^N{\mathbf{R}}_i\left(\omega \right)\hfill \end{array}$$

(5)

$$\begin{array}{ccc}\hfill & =& \sum _{i=1}^N{\mathbf{S}}_i\left(\omega \right){\mathbf{S}}_i^{\mathrm{H}}\left(\omega \right)\hfill \end{array}$$

(6)

where ${\mathbf{S}}_i\left(\omega \right)={[S_{i,1}\left(\omega \right),\cdots ,S_{i,N_{\mathrm{sub}}}\left(\omega \right)]}^{\mathrm{T}}$ is the signal vector at the i-th sub-array and ${}^{\mathrm{T}}$ and ${}^{\mathrm{H}}$ denote transpose and hermitian transpose of a matrix, respectively. As shown in Equations (5) and (6), the covariance matrix can be estimated by averaging the covariance matrices of the sub-arrays. Note that it is necessary to estimate M to enable estimation of the wavenumber.

The theoretical characteristics of the eigenvalues are expressed as follows:

$$l_1\left(\omega \right)\ge \cdots \ge l_{M\left(\omega \right)}\left(\omega \right)\ge l_{M\left(\omega \right)+1}\left(\omega \right)=\cdots =l_{N_{\mathrm{sub}}}\left(\omega \right)={\sigma }_{\mathrm{n}}^2$$

(7)

where $l_i\left(\omega \right)$ is i-th eigenvalue and ${\sigma }_{\mathrm{n}}^2$ is the noise intensity. Note that the absolute values of these eigenvalues do not match the signal intensity directly, i.e., when two waves have the same intensity, the eigenvalues that correspond to these signals do not have the same value. Therefore, the simple thresholding process is not suitable for accurate estimation of the number of signals.

In this section, we calculate the covariance matrix and apply the eigenvalue decomposition technique. When we apply SVD to the signal matrix, $\mathbf{A}$, we theoretically obtain the same result as that obtained when using N transmitters and $N_{\mathrm{sub}}$ receivers with the SVD technique.

2.3. Estimation of the Number of Signals

2.3.1. Overview of the Basic Theory

In this section, we propose an estimation technique for $M\left(\omega \right)$. A method to estimate $M\left(\omega \right)$ that uses information theoretic criteria called the minimum description length (MDL) principle has previously been reported [23,24]. The evaluation index $G\left(m\right)$ that is used to estimate M is given by

$$\begin{array}{c}\hfill G\left(m\right)=-log{\left({\displaystyle \frac{{\prod }_{i=m+1}^{N_{\mathrm{sub}}}{l}_i{\left(\omega \right)}^{{\displaystyle \frac{1}{N_{\mathrm{sub}}-m}}}}{{\sum }_{i=m+1}^{N_{\mathrm{sub}}}\frac{l_i\left(\omega \right)}{N_{\mathrm{sub}}-m}}}\right)}^{(N_{\mathrm{sub}}-m)N}+{\displaystyle \frac{1}{2}}m(2N_{\mathrm{sub}}-m)logN\end{array}$$

(8)

where $G\left(m\right)$ is the index, and the value of m that minimizes $G\left(m\right)$ represents the estimated $M\left(\omega \right)$. Here we define the estimated $M\left(\omega \right)$ as $M^{\prime }\left(\omega \right)$.

2.3.2. Diagonal Loading Technique

In Equation (8), the numerator and the denominator represent the geometric and arithmetic means of the eigenvalues, respectively. Smaller eigenvalues therefore make the estimates unstable. To stabilize the estimation procedure, we thus use a diagonal loading technique that adds a diagonal matrix to the covariance matrix [25,26,27].

The DL process is expressed as follows:

$$\begin{array}{c}\hfill {\mathbf{R}}^{\prime }\left(\omega \right)=\mathbf{R}\left(\omega \right)+\eta \left(\omega \right)\mathbf{I}\end{array}$$

(9)

where $\eta$ is a diagonal loading factor and $\mathbf{I}$ is the identity matrix. The process that was shown in the previous equation can then be rewritten as follows:

$$l_i^{\prime }\left(\omega \right)=l_i\left(\omega \right)+\eta \left(\omega \right)$$

(10)

where $l_i^{\prime }\left(\omega \right)$ is the i-th eigenvalue that is calculated from ${\mathbf{R}}^{\prime }\left(\omega \right)$. Therefore, even when a different DL factor is used with the same covariance matrix, we do not need to perform the eigenvalue decomposition process again. In other words, we use the DL technique to estimate M as an offset value and do not add the DL matrix to the covariance matrix directly. By replacing $l_i$ in Equation (8) with $l_i^{\prime }$, we can then obtain the modified estimates.

2.3.3. Determination of the DL Factor

Selection of the DL factor is an important aspect of the estimation of M. In this study, we used two DL factors that were dependent on the received signal intensity. A schematic illustration of these two factors is shown in Figure 2.

The DL factor is given by

$$\eta \left(\omega \right)={\eta }_1\left(\omega \right){\delta }_1\left(\omega \right)+{\eta }_2{\delta }_2$$

(11)

$${\delta }_1\left(\omega \right)=\mathrm{tr}\left\{\mathbf{R}\left(\omega \right)\right\}$$

(12)

$${\delta }_2=\frac{1}{{\omega }_2-{\omega }_1}{\int }_{{\omega }_1}^{{\omega }_2}\mathrm{tr}\left\{\mathbf{R}\left(\omega \right)\right\}\mathrm{d}\omega$$

(13)

where ${\eta }_1$ and ${\eta }_2$ are the two component DL factors. ${\eta }_2{\delta }_2$ is used for stabilization and is called the constant DL in Figure 2 because, at low signal intensities, ${\eta }_1\left(\omega \right){\delta }_1\left(\omega \right)$, which is called the frequency-dependent DL in Figure 2, approaches zero, and the estimates would thus be unstable.

To determine ${\eta }_1\left(\omega \right)$ for the wavenumber estimation procedure, we assumed that the optimal ${\eta }_1\left(\omega \right)$ is common within a specific frequency range. This assumption was made because, as mentioned below, $M^{\prime }$ shows a step-like change with respect to ${\eta }_1$. The optimal ${\eta }_1$ that gives the optimum M value has a range and is not a critical value.

We vary ${\eta }_1\left(\omega \right)$ and select the minimum value of ${\eta }_1$ that gives $M^{\prime }\left(\omega \right)$ such that it satisfies the following condition:

$$\begin{array}{c}\hfill max\left\{M^{\prime }\left(\omega \right)\right\}<M_{\mathrm{th}},\mathrm{with}(\omega -{\omega }_{\mathrm{w}}\le \omega \le \omega +{\omega }_{\mathrm{w}})\end{array}$$

(14)

where ${\omega }_{\mathrm{w}}$ is the width of the frequency window that is used for stable estimation and $M_{\mathrm{th}}$ is a threshold M value that is sufficiently large for wavenumber estimation.

2.4. Experimental Setup

The setup of the experimental study is the same as the system model shown at Figure 1, but we used a 128-element linear array probe with an element pitch of 0.375 mm that was manufactured by Japan Probe (Kanagawa, Japan). We selected a single transmitter and 28 receivers ($N_{\mathrm{R}}=28$) with an element pitch of 0.75 mm. The distance between the transmitter and the first receiver was 18.75 mm. The center frequency of the transmitted wave was 1.0 MHz. We used two test specimens: (1) a 4-mm-thick copper plate (with a shear wave velocity of 2260 m/s and a longitudinal wave velocity of 4650 m/s) and (2) a 4-mm-thick transversely isotropic bone-mimicking plate (Sawbones, Vashon, WA, USA). The elastic parameters of the bone-mimicking plate have been determined in previous studies [20,21]. The density, shear wave velocity, longitudinal wave velocity along fibers, and longitudinal wave velocity orthogonal to the fibers are 1.64 $\mathrm{g}·{\mathrm{cm}}^{-3}$, 1620 m/s, 3570 m/s, and 2910 m/s, respectively [20,21]. We put the probe parallel to the fibers.

In the proposed method, we set $M_{\mathrm{th}}=N_{\mathrm{sub}}-1$, ${\omega }_{\mathrm{w}}=0.5$ MHz, ${\eta }_2=-40$ dB, and $N_{\mathrm{sub}}=15$. We used ${\omega }_1=4.9$ MHz and ${\omega }_2=5.5$ MHz to select a frequency range that should include only noise.

We prepared the DL factor using values of ${\eta }_1$ within the $-100$ dB to 0 dB range at intervals of 10 dB. $M_{\mathrm{th}}=N_{\mathrm{sub}}-1$ represents the maximum value of M when using a sub-array size of $N_{\mathrm{sub}}$.

In the conventional SVD method, we used values of $N=5$ and $N_{\mathrm{sub}}=24$ to match those used in the conventional study [21]. In the previous study, multiple transmitters were used. The use of multiple transmitters and of sub-arrays is theoretically the same.

3. Results and Discussion

3.1. Evaluation of the Number of the Propagation Modes with DL

To evaluate the M estimation procedure using the DL technique, we performed a simple simulation because it is difficult to know ground truth of M in experimental study. We assumed that three waves with $k=1000,2000$, and 3000 rad/m were propagating. The signal-to-noise ratio (SNR) was 40 dB, $N_{\mathrm{R}}=28$, and $N_{\mathrm{sub}}=15$.

The estimation results showed step-like changes with increasing ${\eta }_1$. Within the $-100\mathrm{dB}\le {\eta }_1\le -90\mathrm{dB}$, $-80\mathrm{dB}\le {\eta }_1\le -20\mathrm{dB}$, and $-10\mathrm{dB}\le {\eta }_1\le 0\mathrm{dB}$ ranges, the corresponding estimated M values were 14, 3, and 0, respectively. The root mean square error (RMSE) of the wavenumber estimation process with ${\eta }_1$ that gives the correct M is 0.73 rad/m. An important point that should be noted here is that the true M value is estimated using a wide range for ${\eta }_1$, e.g., $-80\mathrm{dB}\le {\eta }_1\le -20\mathrm{dB}$.

3.2. Experimental Results

We first investigate the optimal size for the sub-array. Figure 3 shows the RMSEs when we use full-array sizes of $N_{\mathrm{R}}=28$ and 32 and change sub-array size. The measured specimen was bone-mimicking plate. The results show that the optimum sub-array size is $N=N_{\mathrm{R}}/2$. Therefore, with $N_{\mathrm{R}}=28$, the value of $N_{\mathrm{sub}}=15$ is the optimal size.

In the previous paragraph, we described the experimental investigation of the optimal size. Here, we describe the theoretical investigation of this size [29]. To recover the rank of the covariance matrix with M modes, a minimum of M times averaging, i.e., $N=M$, is required. In addition, to measure M modes, a minimum of $N_{\mathrm{sub}}=M+1$ is required. Under these two conditions, when we attempt to maximize M and $N_{\mathrm{sub}}$, the number of receivers $N_{\mathrm{R}}$ should be expressed using $N_{\mathrm{R}}=N_{\mathrm{sub}}+N-1=2N$. Theoretically, we can confirm that the selection $N=N_{\mathrm{R}}/2$ is reasonable.

We determined that $N_{\mathrm{sub}}=15$ is the optimal size. Therefore, the following results employed $N_{\mathrm{sub}}$ of 15. Figure 4, Figure 5, Figure 6 and Figure 7 show the experimental results that were obtained for the copper plate and the bone-mimicking plate, respectively. The color maps show the results obtained using the conventional SVD method. The red dots represent estimated results obtained using the proposed method. The solid white lines show the theoretical curves. To remove any obvious false estimates, we removed all estimates with $c_{\mathrm{p}}<1200$ and $c_{\mathrm{p}}>50,000$ m/s. The proposed method thus successfully depicted the wavenumber without the need for peak search processes.

The RMSEs of the estimation results when using the proposed method with the copper plate and bone-mimicking plate were 108 rad/m and 121 rad/m, respectively. RMSEs are shown in

Table 1. Root mean square errors (RMSEs) of the experimental results with different settings. The units are rad/m.

Proposed Method with Copper Plate	Proposed Method with Bone-Mimicking Plate	ESPRIT with Fixed Threshold (−40 dB)	ESPRIT with Fixed Threshold (−30 dB)	Proposed Method with ${\mathit{\omega }}_{\mathbf{w}}=0.1$ MHz	Proposed Method with ${\mathit{\omega }}_{\mathbf{w}}=0.75$ MHz
108	121	184	143	144	118

. The wavenumbers of the shear waves that had larger effects on the theoretical curve than the longitudinal wave at the center frequency (1.0 MHz) of the copper plate and the bone-mimicking plate were 2780 rad/m and 3879 rad/m, respectively. When compared with the wavenumber of the shear wave, the RMSE was less than 4%.

The spectrum at the frequency that is indicated by the white dotted line in Figure 7 is shown in Figure 8. As Figure 8 illustrates, the proposed method has higher resolution than the conventional method because the method can depict two modes around 4000 rad/m with small error. The wavenumbers that were estimated using the proposed method were almost the same as those at the peak position in the SVD spectrum when the resolution of the SVD method was high enough or the resolution of the ESPRIT algorithm was insufficient to separate the signals. However, it is difficult to identify some of the peaks because of their small amplitudes and low prominence when using the SVD method, such as the peaks with wavenumbers of approximately 2800 and 4300 rad/m.

The wavenumber or phase velocity estimation steps that have high computational costs in the proposed method are eigenvalue decomposition of the covariance matrix with a size of $N_{\mathrm{sub}}\times N_{\mathrm{sub}}$, $\mathcal{O}\left(N_{\mathrm{sub}}^3\right)$ and the SVD of the matrix of size $N_{\mathrm{sub}}-1\times 2M$, $\mathcal{O}(min\left\{(N_{\mathrm{sub}}-1)\times 4M^2,{(N_{\mathrm{sub}}-1)}^2\times 2M\right\})$. The calculation steps are shown in detail in [30]. The corresponding steps in the conventional SVD method, the Radon transform method, and the S-SVD method are the SVD of the matrix of size $N\times N_{\mathrm{sub}}$, $\mathcal{O}(N^2\times N_{\mathrm{sub}})$ [30], the inversion of the matrix of size $N_{\mathrm{c}}\times N_{\mathrm{c}}$, $\mathcal{O}\left(N_{\mathrm{c}}^3\right)$, and the inversion of the matrix of size $N_{\mathrm{k}}\times N_{\mathrm{k}}$, $\mathcal{O}\left(N_{\mathrm{k}}^3\right)$, respectively, where $N_{\mathrm{c}}$ and $N_{\mathrm{k}}$ are the numbers of sampling points in the phase velocity and wavenumber axes, which is normally larger than 64. Therefore, the SVD method showed the smallest computational cost with the lowest resolution. The computational cost of the proposed method was significantly lower than the corresponding costs of the Radon and S-SVD methods.

The calculation time (not including loading time) for the experimental data was less than 0.5 s when using a commercial processor (Core(TM) i7-7500U; Intel, Santa Clara, CA, USA) on a laptop computer with MATLAB (R2017a, Mathworks, Natick, MA, USA) platform. The calculation time taken for the conventional Radon transform was less than 60 s, according to [12]. While noting here that we did not use the same processor or computer, and in fact used a faster processor, the reduction in the calculation time of 99% is interesting. In addition, the computational cost of the S-SVD method is greater than that of the Radon method. The calculation time of the conventional SVD method is around 0.1 s. Note that, as mentioned above, the conventional SVD method has considerably lower computational complexity than the proposed method. Because we ran the SVD method and the proposed method on MATLAB platform, the calculation time does not directly reflect the theoretical complexity. To examine the effectiveness of the proposed method in determination of M, we compared the performance of the proposed method with that of the ESPRIT method with a fixed threshold. Figure 9 and Figure 10 show the results obtained for a bone plate when we used eigenvalues that exceed $-40$ and $-30$ dB of the maximum eigenvalue, respectively. The red circles represent the results obtained using the proposed method. The blue cross marks show the estimates obtained using the ESPRIT algorithm with a fixed threshold. The solid gray lines show the theoretical curves. The RMSEs of the ESPRIT method with fixed thresholds of $-40$ and $-30$ dB were 184 and 143 rad/m, respectively. As shown in Figure 9, smaller thresholds tend to make the estimates unstable because the ESPRIT algorithm treated the noise as the desired signal. In contrast, larger thresholds can miss the weaker modes as shown in Figure 10. The threshold should thus be determined manually. The results are sensitive to the threshold value, so it is difficult to determine the threshold.

Figure 11 shows the wavenumbers for values of $N_{\mathrm{sub}}=5$ and 26. The red cross marks and blue circles show the results obtained for $N_{\mathrm{sub}}=5$ and 26, respectively. The black lines show the theoretical curves. The smaller sub-array had lower resolution and missed many of the modes, while the larger sub-array caused false estimates.

We used a frequency window, denoted by ${\omega }_{\mathrm{w}}$, of 0.5 MHz. When we used ${\omega }_{\mathrm{w}}$ values of 0.1 and 0.75 MHz in the experiments with the bone-mimicking plate, the RMSEs obtained were 144 and 118 rad/m, respectively. Therefore, the dependence on ${\omega }_{\mathrm{w}}$ was insignificant and a larger window width provided stable estimation.

We chose a DL factor of ${\eta }_2=-40$ dB. Because this value was added for stabilization when we analyzed the frequency range using a significantly smaller signal intensity, a smaller value was selected. When we used a larger value, more stable estimates were obtained as a result but the weaker modes were missed by the method.

We used the ESPRIT algorithm in this study. The novel advantage of the ESPRIT algorithm when used for wavenumber estimation of a guided wave propagating along the cortical bone is that it can estimate the wavenumber without the need for a peak search process. To extract the wavenumber from the f–k spectrum image when using conventional methods such as the Radon and SVD-based methods, multiple one-dimensional peak search processes are required [6,20]. While the computational cost is not high, parameters such as the intensity threshold and prominence must be set. Thus, the performance of the conventional method such as RMSE depends on the additional parameter setting.

4. Conclusions

In this study, we proposed a high-resolution and low-computational-cost technique for an AT device. We estimated the number of propagation modes using information theoretic criteria and the DL technique. We proposed a method to estimate the optimal DL value required for guided wave estimation. The proposed method did not involve processes of high computational cost. The proposed method was evaluated experimentally using 4-mm-thick copper and bone-mimicking plates. The estimation error was less than 4% and the calculation time when using the proposed method on a laptop computer was less than 0.5 s. While the processor is different and is in fact faster than that which was used in the previous conventional study, the calculation time is less than 1% of that of the conventional method. We believe that our proposed method thus has the potential to accurately characterize the elastic properties of cortical bone.