Four 8103 hydrophones were mounted on a bracket, and so the sound pressure at four different radial positions could be measured at the same time, as shown in Figure 2. The hydrophone bracket was made from polyvinyl chloride (PVC), which has characteristic impedance similar to that of water, greatly reducing the scattering of sound waves as they passed through the bracket. Three other hydrophones measured the near-field signal at 0, 0.05, 0.1, and 0.15 m from the axial center of the transducer. After the pipeline was filled with water, it was allowed to stand for more than 30 h to eliminate the effect of bubbles on the experiment.

After the fluid column (water) was excited by the transducer, the hydrophones were moved along the axial slots away from the source, and recordings were taken at 10 cm intervals; thus, each slot had 9 or 10 recording points. The use of the hydrophone bracket ensured that the radial positions of the hydrophones remained unchanged when the axial position was changed. In this way, the sound pressure distributions of the fluid column along both the axial and radial directions were measured.

The analysis bandwidth of the collector was set to 0–25 600 Hz, in accordance with the theoretical normal frequencies, and the corresponding sampling frequency was set to 65,536 Hz by the Nyquist law, with a sampling time of 10 s. The power spectrum of the corresponding working condition was obtained through fast Fourier transform (FFT) processing of the time-domain signal. To determine the propagation path of the sound in the experimental system and compare it with the acoustic signal of the fluid in the pipeline at the same time, the vibration of the wall was also measured in this experiment. There were three rows of vibration sensors encircling the outside of the wall, each with four sensors, with a separation of 0.33 m between each row.

The deployment of the hydrophones and vibration sensors and the corresponding labels are shown in Figure 5.

**Figure 5.** Deployment and labeling of hydrophones and vibration sensors.

#### **4. Results and Discussion**

#### *4.1. Behavior of Normal Waves in the Pipeline*

For working condition 1, the transmitting voltage level response of the transducer is shown in Figure 6. The measurement results from the hydrophones at different positions along the axial direction in the pipeline are shown in Figure 7. It can be seen that the acoustic signal in the pipeline exhibits a significant cutoff phenomenon, with a cutoff frequency of 4.7 kHz, which is close to the theoretically determined cutoff frequency of 4.6 kHz.

**Figure 6.** Transmitting voltage level response of transducer.

**Figure 7.** Measurement results from hydrophones at different positions along the axial direction under white noise excitation: (**a**) 0.05, 0.1, and 0.15 m from the source; (**b**) 1.1 m from the source; (**c**) 5.95 m from the source; (**d**) 0.1, 5, 6, and 7 m from the source.

The sound energy of the far field can be divided into four intervals from the spectrum: (1) Below 4.7 kHz; (2) 4.7–10.6 kHz; (3) 10.6–16 kHz; (4) 16 kHz and above. The boundary points of these intervals are the frequencies of the normal wave in the pipeline, which are basically consistent with the calculated frequencies of the corresponding orders, as shown in Table 4.

**Table 4.** Comparison of calculated and experimental normal frequencies.

As mentioned before, the experimental pipeline is slotted. PE is not a very rigid material and, as a result of the slotting process, the tube is deformed radially. This is the main reason for the error between the actual measured frequency and the theoretically calculated frequency. In addition, the parameters in Table 1 are the material elastic parameters of standard high-density PE, which are not necessarily exactly the same as those of the wall of the experimental pipeline, which also leads to an error between the theoretical and experimental results.

In interval (1), the curve of the power spectrum is very close to the background noise. Close to 4.7 kHz, however, the curve suddenly rises, exhibiting a cutoff phenomenon. Then, in interval (2), the curve changes relatively gently, which is basically consistent with the corresponding frequency band of the transmitting transducer frequency response curve. As the distance between the measurement

points and the source becomes greater, it can be seen that the curve decreases monotonically; in interval (3), the curve changes more sharply, and many resonance peaks appear. As the frequency increases, the distance between adjacent peaks also increases, and the appearance of the curve in interval (4) is similar to that in interval (3).

In terms of the behavior of the normal wave, interval (1) is below the cutoff frequency, and the normal wave is attenuated exponentially in the axial direction. The curve of the power spectrum in this frequency band is basically the same as the background, and it can be seen that the power spectrum decays exponentially with frequency in the near field.

Interval (2) lies between the first-order and second-order normal wave frequencies; only the first-order normal wave can propagate, and the curve of the power spectrum hardly changes with frequency. The trend in this section of the curve is related to the transmitting response of the source, with the curve monotonically decreasing along the axial direction as a result of absorption of sound waves by the tube wall.

The first and second orders of the normal wave propagate simultaneously in interval (3), and the first, second, and third orders propagate in interval (4), where there is strong interference leading to large fluctuations in the curve and to the appearance of many resonant peaks.

The results of the vibration signal measurements are shown in Figure 8. The frequency response of the vibration signal is similar to that of the acoustic signal in its overall trend, with a cutoff phenomenon, and the division of the modal frequencies of each order is obvious. The vibration signal in the frequency band above the cutoff frequency is transmitted from the acoustic signal in the water.

**Figure 8.** Measurement results from vibration sensors on the outside of the pipeline wall: (**a**) 2 m from the source; (**b**) 3 m from the source.

#### *4.2. Variation of the Sound Field along the Axial Direction*

The variation of the sound field along the axial direction can be seen more clearly from analysis of the response to a single-frequency source compared with the response to the white noise in working condition 1. In working conditions 2 and 3, two representative single-frequency signals were used, 4.2 and 5.2 kHz, which are respectively below and above the cutoff frequency. The main frequency power spectrum was obtained as an average of the measurements by the four hydrophones on the bracket, and its variation along the axial distance is shown in Figure 9. It can be seen from Figure 9a that for a source frequency below the cutoff frequency, the power spectrum of the main frequency is attenuated very rapidly, indeed exponentially, as the distance increases. Acoustic signals below the cutoff frequency cannot propagate axially over long distances in the pipeline. For a source frequency above the cutoff frequency, as shown in Figure 9b, the power spectrum of the main frequency hardly changes with distance. There is only 4 dB attenuation from 3 to 10 m, and this attenuation is a result of acoustic absorption by the pipeline wall.

**Figure 9.** Variation of the main frequency average power spectrum along the axial direction: (**a**) 4.2 kHz source; (**b**) 5.2 kHz source.

#### *4.3. Variation of the Sound Field along the Radial Direction*

To explore the variation of the sound field in the radial direction, measurements were performed before and after the hydrophone bracket was raised by 1.5 cm, as shown schematically in Figure 10. The results of these measurements are shown in Figure 11.

**Figure 10.** Raising of hydrophone bracket.

**Figure 11.** Results of hydrophone measurements at different depths 5.05 m from the source: (**a**) Before and (**b**) after lifting of hydrophone bracket.

The frequency band between the first- and second-order normal frequencies was analyzed. Before lifting, the 2# and 3# hydrophones were at the same distance from the axis, and similarly for the 1# and 4# hydrophones. Therefore, in Figure 11a, the curves of the power spectra from the 2# and 3# hydrophones are the same, as are those from the 1# and 4# hydrophones. After lifting, the 3# hydrophone is nearest to the axis, followed in order by the 2#, 4#, and 1# hydrophones, which is consistent with the increasing strengths of the respective power spectra. This is in accordance with the theoretical radial dependence on the Bessel function *J*0(*krr*) in Equation (8).

The following is a quantitative analysis of the radial distribution. The distance between each pair of hydrophones is 50 mm, and the hydrophone bracket is initially at the radial center position. After lifting, the distances of the 3#, 2#, 4#, and 1# hydrophones from the axis are 10, 40, 60, and 90 mm, respectively. At the cutoff frequency, *kr* = *k*0, depending on the distance *r* from the axis, the theoretical differences between the sound pressure self-spectrum measured by the 3# hydrophone and those measured by the other three hydrophones can be calculated, and the results are compared with the experimental measurements in Table 5. For convenience of exposition, the distance between the 3# hydrophone and the axis is set as *r*0, and the differences between the sound pressure power spectrum of the 1#, 2#, and 4# hydrophones and the 3# hydrophone as *X*1, *X*2, and *X*3, respectively.

**Table 5.** Comparison of hydrophone measurements and theoretical values at different radial locations.

There is good agreement between theory and experiment, and the radial distribution of the normal wave in the tube is quantitatively confirmed to follow the Bessel function behavior.

The reasons for the error are as follows. In the experiment, the magnitude of the lifting was controlled manually, and not very accurately, which is the main source of error: If the lifting range were slightly larger, and the 3# and 4# hydrophones closer to the axis, the amplitude would be higher. If the 1# and 2# hydrophones were further away from the axis, the amplitude would be smaller, which would lead to an increase in *X*2; according to the properties of the Bessel function, the closer the value of the function is to the axis, the slower is its rate of change. Therefore, the power spectrum at the 4# hydrophone would increase more if it were lifted, so the value of *X*3 would be reduced. If the 1# hydrophone were closer to the upper slot and the outside medium were air (which can be regarded as an absolutely soft boundary), the sound would be totally reflected, so the amplitude would become higher, leading to a decrease in *X*1. In addition, the radial deformation of the pipeline due to slotting, scattering by the hydrophone bracket, and the fact that the orientation of the hydrophone was not strictly in the axial direction are all possible sources of error.

#### *4.4. Measurements under Mechanical Excitation*

Working conditions 4 and 5 use mechanical force excitation from outside the pipeline wall in a position directly under that of the sound source in the previous working conditions, as shown in Figure 12. Corresponding to working conditions 2 and 3, the single frequencies of excitation applied in working conditions 4 and 5 are 4.2 and 5.2 kHz, respectively. The measurement results from the hydrophones and vibration sensors are shown in Figures 13–16. In contrast to the response of the acoustic signal, the excitation of the exciter to the pipe wall was a single-point excitation. Therefore, when the sound propagated mainly along the wall, the sound source excitation and the exciter excitation had completely different radial distribution laws. When the sound propagated mainly along the liquid in the tube, both had the same radial distribution. This is an important basis for judging the propagation path of sound in a pipeline.

**Figure 12.** Exciter and excitation point in working conditions 4 and 5.

**Figure 13.** Mechanical excitation at a single frequency of 4.2 kHz. Measurement results from hydrophones at different distances from the excitation point: (**a**) 1.1 m; (**b**) 5.1 m.

Even if the pipeline wall is excited by a mechanical force, an acoustic signal below the cutoff frequency cannot propagate a long distance in the case of weak excitation. It can be seen from Figure 13b that the hydrophones 5.1 m from the excitation point have difficulty in picking up a signal. In Figure 13a, the acoustic signal at the main frequency exhibits a radial dependence that is different from the Bessel function: The hydrophone near the lower outer wall close to the excitation point receives a stronger signal.

**Figure 14.** Mechanical excitation at a single frequency of 4.2 kHz. Measurement results from vibration sensors at different distances from the excitation point: (**a**) 1.3 m; (**b**) 5.3 m.

From the vibration measurement results in Figure 14, it can be seen that a vibration signal at the main frequency can still be measured on the wall at 5.3 m, but it is attenuated compared with the measurement from the sensor at 1.3 m, and the hydrophones in the pipeline are unable to detect any signal at 5.3 m. It can be deduced that the signal at 4.2 kHz is propagated mainly through the pipe wall in the form of vibrations. The second and third peaks in Figure 14 result from frequency doubling. The exciter generates frequency doubling when transmitting a single-frequency signal, which is a consequence of its own physical structure and has no effect on the results of this experiment.

**Figure 15.** Mechanical excitation at a single frequency of 5.2 kHz. Measurement results from hydrophones at different distances from the excitation point: (**a**) 1.1 m; (**b**) 1.1 m (local); (**c**) 5.1 m; (**d**) 5.1 m (local).

**Figure 16.** Mechanical excitation at a single frequency of 5.2 kHz. Measurement results from vibration sensors at different distances from the excitation point: (**a**) 1.3 m; (**b**) 5.1 m.

For excitation at 5.2 kHz, which is above the cutoff frequency, Figure 15a,b shows the measurement results from hydrophones 1.1 m from the excitation point, and Figure 15c,d those from hydrophones 5.1 m from the excitation point, with Figure 15b,d being partial displays of the frequency band near the main frequency shown in Figure 15a,c, respectively. Figure 16 shows the measurement results from the wall vibration sensors 1.3 and 5.3 m from the excitation point.

The sound power at the main frequency suffers almost no attenuation in the axial direction from 1.1 to 5.1 m, and conforms to the Bessel function dependence in the radial direction. It can be deduced that the signal at 5.2 kHz is propagated mainly in the form of sound through the fluid in the pipeline.

#### **5. Conclusions**

The first four orders of normal frequencies in a fluid-filled PE pipeline were calculated, and the distributions of sound in the axial and radial directions were analyzed. The acoustic propagation characteristics of such a pipeline were also studied in an experimental system.

Both the theoretical and experimental investigations have revealed the following:

**Author Contributions:** Conceptualization, Q.L. and J.S.; methodology, Q.L.; software, J.S.; data validation, J.S.; formal analysis, D.S.; writing—original draft preparation, J.S.; writing—review and editing J.S. and D.S.; supervision and project administration, Q.L.

**Funding:** This research was funded by Acoustic Science and Technology Laboratory, Harbin Engineering University (SSJSWDZC2018010) and by the National Science Foundation of China (11874131).

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

### *Article* **Regularization Factor Selection Method for** *l***1-Regularized RLS and Its Modification against Uncertainty in the Regularization Factor**

#### **Junseok Lim 1,\* and Seokjin Lee 2**

Received: 5 December 2018; Accepted: 4 January 2019; Published: 8 January 2019

#### **Featured Application: This algorithm can be applied to various kinds of sparse channel estimations, e.g., room impulse response, early reflection, and underwater channel response.**

**Abstract:** This paper presents a new *l*1-RLS method to estimate a sparse impulse response estimation. A new regularization factor calculation method is proposed for *l*1-RLS that requires no information of the true channel response in advance. In addition, we also derive a new model to compensate for uncertainty in the regularization factor. The results of the estimation for many different kinds of sparse impulse responses show that the proposed method without a priori channel information is comparable to the conventional method with a priori channel information.

**Keywords:** *l*1-regularized RLS; sparsity; room impulse response; total least squares; regularization factor

#### **1. Introduction**

Room impulse response (RIR) estimation is a problem in many applications that use acoustic signal processing. The RIR identification [1] is fundamental for various applications such as room geometry related spatial audio applications [2–5], acoustic echo cancellation (AEC) [6], speech enhancement [7], and dereverberation [8]. In [9], the RIR has relatively large magnitude values during the early part of the reverberation and fades to smaller values during the later part. This indicates that most RIR entries have values close to zero. Therefore, the RIR has a sparse structure. The sparse RIR model is useful for estimating RIRs in real acoustic environments when the source is given a priori [10]. There has been recent interest in adaptive algorithms for sparsity in various signals and systems [11–22]. Many adaptive algorithms based on least mean square (LMS) [11,12] and recursive least squares (RLS) [14–17] have been reported with different penalty functions. Sparse estimation research, such as that done by Eksioglu and Tanc [17], has proposed a sparse RLS algorithm, *l*1-RLS, which is fully recursive like the plain RLS algorithm. The algorithm of *l*1-RLS in [17] proposed a proper calculation method for the regularization factor. These recursive algorithms have the potential for sparse RIR estimation; however, the regularization factor should be established prior to applying these algorithms. The regularization factor calculation method requires information about a true sparse channel response for a good performance. The authors in [18,19] have also proposed recursive regularization factor selection methods; however, these methods still need the true impulse response in advance.

In this paper, we propose a new regularization factor calculation method for *l*1-RLS algorithm in [17]. The new regularization factor calculation needs no information for the true channel response in advance. This makes it possible to apply *l*1-RLS algorithm in various room environments. In addition, we derive a new model equation for *l*1-RLS in [17] with uncertainty in the regularization factor and show that the new model is similar to the total least squares (TLS) model that compensates for uncertainty in the calculated regularization factor without the true channel response. For the performance evaluation, we simulate four different sparse channels and compare channel estimation performances. We show that, without any information of the true channel impulse response, the performance of the proposed algorithm is comparable to that of *l*1-RLS with the information of the true channel impulse response.

This paper is organized as follows. In Section 2, we summarize *l*1-RLS in [17]. In Section 3, we summarize the measure of sparsity. In Section 4, we propose a new method for the regularization calculation. In Section 5, we show that *l*1-RLS with uncertainty in the regularization factor can be modeled as the TLS model. In Section 6, we summarize *l*1-RTLS (recursive total least squares) algorithm as a solution for *l*1-RLS with uncertainty in the regularization factor. In Section 7, we present simulation results to show the performance of the proposed algorithm. Finally, we give the conclusion in Section 8.

#### **2. Summarize** *l***1-RLS**

In the sparse channel estimation problem of interest, the system observes a signal represented by an *M* × 1 vector **x**(*k*)=[*xk*, ··· , *xk*−*M*+1] *T* at time instant n, performs filtering, and obtains the output *y*(*i*) = **x***T*(*k*)**w***o*(*k*), where **w***o*(*k*)=[*wk*, ··· , *wk*−*M*+1] *T* is the M dimensional actual system with finite impulse response (FIR) type. For system estimation, an adaptive filter system applies with M dimensional vector **w**(*k*) to the same signal vector **x**(*k*) and produces an estimated output *y*ˆ(*k*) = **x***T*(*k*)**w**(*k*), and calculates the error signal *e*(*k*) = *y*(*k*) + *n*(*k*) − *y*ˆ(*k*) = *y*6(*k*) − *y*ˆ(*k*), where *n*(*k*) is the measurement noise, *y*(*k*) is the output of the actual system, and *y*ˆ(*k*) is the estimated output. In order to estimate the channel impulse response, an adaptive algorithm minimizes the cost function defined by

$$\mathbf{w} = \operatorname\*{argmin}\_{\mathbf{w}} \frac{1}{2} \sum\_{m=0}^{k} \lambda^{k-m} (e(m))^2. \tag{1}$$

From the gradient based minimization, Equation (1) becomes

$$\mathbf{R}(k)\mathbf{w}(k) = \mathbf{r}(k),\tag{2}$$

where **R**(*k*) = *k* ∑ *m*=0 *λk*−*m***x**(*m*)**x***T*(*m*) and **r**(*k*) = *k* ∑ *m*=0 *λk*−*my*6(*m*)**x**(*m*). This equation is the normal equation for the least squares solution. Especially, **w***o*(*k*) is considered as a sparse system when the number of nonzero coefficients *K* is less than the system order of *M*. In order to estimate the sparse system, most estimation algorithms exploit non-zero coefficients in the system [11–17]. In [17], Eksioglu proposed a full recursive *l*1-regularized algorithm by the minimization of the object function as shown in Equation (3).

$$J\_k = \frac{1}{2}\varepsilon\_k + \gamma\_k \|\mathbf{w}\|\_{1^\prime} \tag{3}$$

where *εk* = *k* ∑ *m*=0 *λk*−*m*(*e*(*m*))2 . From the minimization of Equation (3), a modified normal equation was derived as shown in Equation (4).

$$\mathbf{R}(k)\mathbf{w}(k) = \mathbf{r}(k) - \gamma\_k \nabla^s \|\mathbf{w}(k)\|\_1 = \mathbf{\hat{p}}(k). \tag{4}$$

When we solve Equation (4), we should select the regularization factor as shown in Equation (5).

$$\gamma\_k = \frac{2\frac{tr\{\mathbf{R}^{-1}(k)\}}{M}}{\left\|\mathbf{R}^{-1}(k)\nabla^s f(\mathbf{w}(k))\right\|\_2^2} \times \left[ (f(\mathbf{w}(k)) - \rho) + \nabla^s f(\mathbf{w}(k))\mathbf{R}^{-1}(k)\varepsilon(k) \right],\tag{5}$$

where *f*(**w**(*k*)) = **w**(*k*)1 and the subgradient of *f*(**w**) is ∇*s* **w**1 = sgn(**w**). In Equation (5), the regularization factor has the parameter, *ρ*, which should be set beforehand. In [17], the parameter was set as *ρ* = *f*(**w***true*) = **w***true*1 with **w***true* indicating the impulse response of the true channel. There was no further discussion about how to set *ρ*. However, it is not practical to know the true channel in advance.

#### **3. Measure of Sparseness**

In [20], the sparseness of a channel impulse response is measured by Equation (6).

$$\chi = \frac{L}{L - \sqrt{L}} \left( 1 - \frac{||\hat{\mathbf{w}}||\_1}{\sqrt{L} ||\hat{\mathbf{w}}||\_2} \right) \tag{6}$$

where **w**ˆ *p* is the p-norm of **w**ˆ and *L* is the dimension of **w**ˆ . The range of *χ* is 0 ≤ *χ* ≤ 1. That is dependent on the sparseness of **w**ˆ . As **w**ˆ becomes sparser, the sparsity, *χ*, comes close to 1, and as **w**ˆ becomes denser, *χ* comes close to 0. We often have small and none-zero value of *χ*, even in a dense channel. For example, Figure 1 shows the relation of the value of *χ* and the percentage of none-zero components in **w**ˆ with *L* = 215. In Figure 1, we consider all possible cases of none-zero components in **w**ˆ .

**Figure 1.** Sparsity (*χ*) vs. the percentage of none zero coefficients in the channel impulse response.

#### **4. New** *ρ* **Selection Method in the Sparsity Regularization Constant** *γk*

Section 2 shows that the regularization constant *γk* in Equation (5) needs *ρ* to be set as *ρ* = true system impulse response1 = **w***true*1. However, we need a new method in the constant selection because Equation (5) is not practical. Therefore, Section 4 proposes a new method to set this constant.

For a practical method for the constant selection, we can consider using the estimated vector **w**ˆ instead of using the true vector **w***true* because **w**ˆ , the solution with *l*1-norm, will be closer to the sparse true vector than the solution of the conventional RLS. The more iteration is repeated, the more **w**ˆ converges to the true value. Conventional RLS also converges to the true value; however, the solution with *l*1-norm, is closer to the sparse true value. Therefore, we can use sparse estimate **w**ˆ instead of **w***true* when we set *ρ*, and the uncertainty arising from this is compensated through a TLS solution in the next section. When we determine *ρ* using the estimated **w**ˆ , we choose between the average *ρ* and the current estimate **w**ˆ 1. Table 1 summarizes the *ρ* selection steps.

The determination method for *ρ* value shown in Table 1 is as follows. In Step 1, the sparsity of the estimated **w**ˆ is calculated. The sparsity represents the sparseness of **w**ˆ as a number [23]. In Step 2, *l*1-norm of the estimated **w**ˆ is scaled and the value is averaged with the previous *ρ* value. The scaling value approaches 1 as the sparsity, *χ*, gets close to 1. However, the scaling value gets close to *e*−1 0.37 as the sparsity, *χ*, gets close to 0. Therefore, the scaling does not change *l*1-norm of **w**ˆ for the sparse **w**ˆ . Instead the scaling changes the *l*1-norm smaller for the dense **w**ˆ . In Step 3, the smaller one between the averaged *ρ* and the *l*1-norm of the estimated **w**ˆ is selected as the new *ρ* value. In this case, the *ρ* value becomes completely new if the *l*1-norm of the estimated **w**ˆ is selected, otherwise the previous trend is maintained. In Figure 1, the reference value 0.75 used in Step 3 means that less than 16% of all the impulse response taps are not zero.

**Table 1.** *ρ* selection method in the sparsity regularization constant *γk*.

#### **5. New Modeling for** *l***1-RLS with Uncertainty in the Regularization Factor**

If we set *ρ* = constant, the regularization factor becomes

$$\begin{array}{lcl}\widetilde{\gamma}\_{k} &= \frac{2\frac{\mathsf{tr}\left(\mathbf{k}^{-1}(k)\right)}{M}}{\left\|\mathbf{R}^{-1}(k)\nabla^{\mathsf{c}}f(\mathbf{w}(k))\right\|\_{2}^{2}} \times \left[\left(f(\mathbf{w}(k))-\text{constant}\right)+\nabla^{\mathsf{s}}f(\mathbf{w}(k))\times\mathbf{R}^{-1}(k)\varepsilon(k)\right]} \\ &= \frac{2\frac{\mathsf{tr}\left(\mathbf{k}^{-1}(k)\right)}{M}}{\left\|\mathbf{R}^{-1}(k)\nabla^{\mathsf{c}}f(\mathbf{w}(k))\right\|\_{2}^{2}} \times \left(f(\mathbf{w}(k))-\left\|\mathbf{h}\right\|\_{1}+\left\|\mathbf{h}\right\|\_{1}-\text{constant}\right) +\frac{2\frac{\mathsf{tr}\left(\mathbf{k}^{-1}(k)\right)}{M}\nabla^{\mathsf{c}}f(\mathbf{w}(k))\mathbf{R}^{-1}(k)\varepsilon(k)}{\left\|\mathbf{R}^{-1}(k)\nabla^{\mathsf{c}}f(\mathbf{w}(k))\right\|\_{2}^{2}}.\end{array} \tag{7}$$

Then,

$$\tilde{\gamma}\_k = \gamma\_k + \frac{2\frac{tr\{\mathbf{R}^{-1}(k)\}}{M}(\|\mathbf{h}\|\_1 - \text{constant})}{\left\|\mathbf{R}^{-1}(k)\nabla^s f(\mathbf{w}(k))\right\|\_2^2} = \gamma\_k + \Delta\gamma. \tag{8}$$

Using Equation (8), Equation (4) becomes

$$\mathbf{R}(k)\mathbf{w}(k) = \mathbf{r}(k) - (\gamma\_k + \Delta\gamma)\nabla^s ||\mathbf{w}(k)||\_1 \quad . \tag{9}$$

∇*s* **w**1 = sgn(**w**) is represented as

$$\nabla^s \|\mathbf{w}(k)\|\_1 = \begin{bmatrix} \ddots & & & \\ & \ddots & & \\ & & \frac{1}{|\mathbf{w}\_i|} & & \\ & & & \ddots & \\ & & & & \ddots \end{bmatrix} \mathbf{w}(k). \tag{10}$$

By applying Equation (10) to Equation (9),

$$\left(\mathbf{R}(k) + \Delta\gamma \begin{bmatrix} \ddots & & & \\ & \frac{1}{|\mathbf{w}\_{i}|} & & \\ & & \ddots & \\ & & & \ddots & \\ \end{bmatrix}\right) \mathbf{w}(k) = \mathbf{r}(k) - \gamma\_{k} \|\mathbf{w}(k)\|\_{1'} \tag{11}$$

where **w***i* is *i*-th element of **w**(*k*). Then it is simplified as

$$\left(\mathbf{R}(k) + \Delta\gamma \begin{bmatrix} \ddots & & & \\ & \frac{1}{|\mathbf{w}\_i|} & & \\ & & \ddots & \\ & & & \ddots \end{bmatrix}\right) \mathbf{w}(k) = \mathbf{\hat{p}}(k) \tag{12}$$

Equation (12) is very similar to the system model in Figure 2 that is contaminated by noise both in input and in output. Suppose that an example of the system in Figure 2 is represented as

$$\begin{bmatrix} \mathbf{x}\_{k} + \boldsymbol{n}\_{i,k} & \cdots & \mathbf{x}\_{k-N+1} + \boldsymbol{n}\_{i,k-N+1} \\ \mathbf{x}\_{k-1} + \boldsymbol{n}\_{i,k-1} & \cdots & \mathbf{x}\_{k-N} + \boldsymbol{n}\_{i,k-N} \\ \vdots & \ddots & \vdots \\ \mathbf{x}\_{k-N+1} + \boldsymbol{n}\_{i,k-N+1} & \cdots & \mathbf{x}\_{k-2N+2} + \boldsymbol{n}\_{i,k-2N+2} \end{bmatrix} \times \mathbf{w}(k) = \begin{bmatrix} y\_{k} + \boldsymbol{n}\_{o,k} \\ y\_{k-1} + \boldsymbol{n}\_{o,k-1} \\ \vdots \\ y\_{k-N+1} + \boldsymbol{n}\_{o,k-N+1} \end{bmatrix},\tag{13}$$

where *xk* is *x*(*k*), *ni*,*k* is *ni*(*k*), and *no*,*k* is *no*(*k*). Equation (13) is simplified as

$$\mathbf{A}\mathbf{w}(k) = \mathbf{b}.\tag{14}$$

**Figure 2.** The model of a noisy input and noisy output system.

If we multiply Equation (14) by **A***H* and average it, we get

$$E\left(\mathbf{A}^H \mathbf{A}\right)\mathbf{w}(k) = E\left(\mathbf{A}^H \mathbf{b}\right). \tag{15}$$

We can rewrite Equation (15) as follows

$$\begin{bmatrix} r\_{xx}(0) + \sigma\_{\eta}^{2} & r\_{xx}(1) & \cdots & r\_{xx}(N-1) \\ r\_{xx}(1) & r\_{xx}(0) + \sigma\_{\eta}^{2} & \cdots & r\_{xx}(N-2) \\ \vdots & \vdots & \ddots & \vdots \\ r\_{xx}(N-1) & r\_{xx}(N-2) & \cdots & r\_{xx}(0) + \sigma\_{\eta}^{2} \end{bmatrix} \mathbf{w}(k) = \begin{bmatrix} r\_{xy}(0) \\ r\_{xy}(1) \\ \vdots \\ r\_{xy}(N-1) \end{bmatrix} \tag{16}$$

Then, it can be represented as

$$\mathbf{w}\left(\mathbf{R} + \sigma\_n^2 \mathbf{I}\right)\mathbf{w}(k) = \breve{\mathbf{p}}(k). \tag{17}$$

When we compare Equation (12) with Equation (17), the two system models have almost the same form. Therefore, it is feasible that the TLS method can be applied to Equation (12) [24–30]. Therefore, we expect to obtain almost the same performance as *l*1-RLS with the true channel response if we apply the TLS method by the regularization factor with the new *ρ* in Table 1. In the next section, we summarize *l*1-RTLS (recursive total least squares) algorithm in [29].

#### **6. Summarize** *l***1-RTLS for the Solution of** *l***1-RLS with Uncertainty in the Regularization Factor**

Lim, one of the authors of this paper, has proposed the TLS solution for *l*1-RLS known as *l*1-RTLS [30]. In this section, we summarize *l*1-RTLS in [30] for the solution of Equation (11).

The TLS system model assumes that both input and output are contaminated by additive noise as Figure 2. The output is given by

$$
\tilde{y}(k) = \tilde{\mathbf{x}}^T(k)\mathbf{w}\_o + n\_o(k), \tag{18}
$$

where the output noise *no*(*k*) is the Gaussian white noise with variance *σ*2 *o* . The noisy input vector in the system is modeled by

$$
\widetilde{\mathbf{x}}(k) = \mathbf{x}(k) + \mathbf{n}\_i(k) \in \mathbb{C}^{M \times 1}, \tag{19}
$$

where **n***i*(*k*) = [*ni*(*k*), *ni*(*k* − 1), ··· *ni*(*k* − *M* + 1)] *T* and the input noise *ni*(*k*) is the Gaussian white noise with variance *σ*2 *i* . For the TLS solution, we set the augmented data vector as

$$\overline{\mathbf{x}}(k) = \left[\widetilde{\mathbf{x}}^T(k), \widetilde{\mathbf{y}}(k)\right]^T \in \mathbb{R}^{(M+1)\times 1}.\tag{20}$$

The correlation matrix is represented as

$$
\overline{\mathbf{R}} = \begin{bmatrix}
\ \overset{\scriptstyle \mathbf{\tilde{R}}}{\mathbf{p}} & \mathbf{p} \\
\mathbf{p}^T & c
\end{bmatrix},
\tag{21}
$$

where **p** = *E*{6**x**(*k*)*y*(*k*)}, *c* = *E*{*y*(*k*)*y*(*k*)}, **R** = *E* 7 **x**(*k*)**x***T*(*k*) 8 and **R**6 = *E* 3 6**x**(*k*)6**x***T*(*k*) 4 = **R** + *σ*2 *i* **I**. In [27,28], the TLS problem becomes to find the eigenvector associated with the smallest eigenvalue of **R**. Equation (22) is the typical cost function to find the eigenvector associated with the smallest eigenvalue of **R**.

$$J(k) = \frac{1}{2}\tilde{\mathbf{w}}^T(k)\overline{\mathbf{R}}(k)\tilde{\mathbf{w}}(k),\tag{22}$$

where **R**(*k*) is a sample correlation matrix at *k*-th instant, and **w**6 (*k*) = - **w**ˆ *T*(*k*), −1 *T* in which **w**ˆ (k) is the estimation result for the unknown system at *k*-th instant. We modify the cost function by adding a penalty function in order to reflect prior knowledge about the true sparsity system.

$$J(k) = \frac{1}{2}\breve{\mathbf{w}}^T(k)\overline{\mathbf{R}}(k)\widetilde{\mathbf{w}}(k) + \lambda\left(\widetilde{\mathbf{w}}^T(k)\widetilde{\mathbf{w}}(k-1) - 1\right) + \gamma\_k f(\breve{\mathbf{w}}(k)),\tag{23}$$

where *λ* is the Lagrange multiplier and *γk* is the regularized parameter in [13]. We solve the equations by ∇**w**ˆ *J*(*k*) = 0 and ∇*λ J*(*k*) = 0 simultaneously. ∇**w**ˆ *J*(*k*) = 0 :

$$2\overline{\mathbf{R}}(k)\widetilde{\mathbf{w}}(k) + \lambda\widetilde{\mathbf{w}}(k-1) + \gamma\_k \nabla^s f(\widetilde{\mathbf{w}}(k)) = 0,\tag{24}$$

$$
\nabla\_{\lambda} I(k) = 0: \quad \tilde{\mathbf{w}}^{\top}(k)\tilde{\mathbf{w}}(k-1) = 1,\tag{25}
$$

where the subgradient of *f*(**w**6 ) = **w**6 1 is ∇*s* **w**6 **w**6 1 = sgn(**w**6 ). From (24), we obtain

$$
\widetilde{\mathbf{w}}(k) = -\frac{\lambda}{2} \overline{\mathbf{R}}^{-1}(k) \widetilde{\mathbf{w}}(k-1) - \gamma\_k \overline{\mathbf{R}}^{-1}(k) \nabla^s f(\widetilde{\mathbf{w}}(k)).\tag{26}
$$

Substituting Equation (26) in Equation (25), we get

$$\left(-\frac{\lambda}{2}\overline{\mathbf{R}}^{-1}(k)\widetilde{\mathbf{w}}(k-1) - \gamma\_k \overline{\mathbf{R}}^{-1}(k)\nabla^s f(\widetilde{\mathbf{w}}(k))\right)^T \times \widetilde{\mathbf{w}}(k-1) = 1,\tag{27}$$

or

$$\lambda = -2 \frac{1 + \gamma\_k \nabla^s f(\tilde{\mathbf{w}}(k))^T \mathbf{\tilde{R}}^{-1}(k) \tilde{\mathbf{w}}(k-1)}{\tilde{\mathbf{w}}^T(k-1) \mathbf{\tilde{R}}^{-1}(k) \tilde{\mathbf{w}}(k-1)}. \tag{28}$$

Substituting *λ* in Equation (26) by Equation (28) leads to

$$\check{\mathbf{w}}(k) = \frac{1 + \gamma\_k \nabla^s f(\check{\mathbf{w}}(k))^T \overline{\mathbf{R}}^{-1}(k) \check{\mathbf{w}}(k-1)}{\check{\mathbf{w}}^T(k-1) \overline{\mathbf{R}}^{-1}(k) \check{\mathbf{w}}(k-1)} \times \overline{\mathbf{R}}^{-1}(k) \check{\mathbf{w}}(k-1) - \gamma\_k \overline{\mathbf{R}}^{-1}(k) \nabla\_{\check{\mathbf{w}}} f(\check{\mathbf{w}}(k)). \tag{29}$$

Equation (29) can be expressed in a simple form as

$$
\widetilde{\mathbf{w}}(k) = a \overline{\mathbf{R}}^{-1}(k) \widetilde{\mathbf{w}}(k-1) - \gamma\_k \overline{\mathbf{R}}^{-1}(k) \nabla^s f(\widetilde{\mathbf{w}}(k)), \tag{30}
$$

where *α* = 1+*γk*∇*s f*(**w**6 (*k*)) *T***R**−1 (*k*)**w**6 (*k*−1) **w**6 *T*(*k*−1)**R**−1 (*k*)**w**6 (*k*−1) . Because asymptotically **w**6 (*k*) = 1 as *k* → ∞, Equation (29) can be approximated as the following two equations.

$$\tilde{\mathbf{w}}(k) \simeq \overline{\mathbf{R}}^{-1}(k)\tilde{\mathbf{w}}(k-1) - \gamma\_{\tilde{k}}\Big(\tilde{\mathbf{w}}^{T}(k-1)\overline{\mathbf{R}}^{-1}(k-1)\tilde{\mathbf{w}}(k-1)\Big)\overline{\mathbf{R}}^{-1}(k)\nabla^{s}f(\tilde{\mathbf{w}}(k-1)).\tag{31}$$

$$
\tilde{\mathbf{w}}(k) = \tilde{\mathbf{w}}(k) / \|\tilde{\mathbf{w}}(k)\|. \tag{32}
$$

Finally, we obtain the estimated parameter of the unknown system as

$$
\hat{\mathbf{w}}(k) = -\hat{\mathbf{w}}\_{1:M}(k) / \hat{\mathbf{w}}\_{M+1}(k). \tag{33}
$$

For Equation (23), we can use the modified regularization factor *γk* in [30]

$$\gamma\_k = \frac{2\frac{tr\left(\overline{\mathbf{R}}^{-1}(k)\right)}{M}}{\left\|\overline{\mathbf{R}}^{-1}(k)\nabla^s f\left(\mathfrak{w}\_{\text{aug}}(k)\right)\right\|\_2^2} \times \left[\left(f\left(\mathfrak{w}\_{\text{aug}}(k)\right) - \rho\right) + \nabla^s f\left(\mathfrak{w}\_{\text{aug}}(k)\right)\overline{\mathbf{R}}^{-1}(k)\varepsilon(k)\right],\tag{34}$$

where **w**ˆ *aug*(*k*) = - **w**ˆ *T*(*k*), −1 *T* , **w**ˆ *aug*,*RLS*(*k*) = - **w**ˆ *T RLS*(*k*), −1 *T* , *ε*(*k*) = **w**ˆ *aug*(*k*) − **w**ˆ *aug*,*RLS*(*k*), and **w**ˆ *RLS*(*k*) is the estimated parameter by recursive least squares (RLS). As *f*(**w**ˆ ) = **w**ˆ 1, the subgradient of *f*(**w**ˆ *aug*(*k*)) is

$$\left. \nabla^{s} \right| \left| \mathfrak{w}\_{aug}(k) \right| \big|\_{1} = \text{sgn}(\mathfrak{w}\_{aug}(k)). \tag{35}$$

As mentioned in Section 4, we apply new constant *ρ* in Table 1, to the regularization factor *γk* in Equation (34) instead of **w***true*1, where **w***true* is the true system impulse response.

#### **7. Simulation Results**

This section confirms the performance of the proposed algorithm in sparse channel estimation. In the first experiment, the channel estimation performance is compared with other algorithms using randomly generated sparse channels. In this simulation, we follow the same scenario in the experiments as [17]. The true system vector **w***true* is 64 dimensions. In order to generate the sparse channel, we set the number of the nonzero coefficients, S, in the 64 coefficients and randomly position the nonzero coefficients. The values of the coefficients are taken from an *N*(0, 1/S) distribution, where *N*( ) is the normal distribution. In the simulation, we estimate the channel impulse response by the proposed algorithms that are *l*1-RLS using the *ρ* in Table 1 and *l*1-RTLS using the *ρ* in Table 1. For the comparison, we estimate the channel impulse response by *l*1-RLS using the true channel response; in addition, we also execute the regular RLS algorithm in an oracle setting (oracle-RLS) where the positions of the true nonzero system parameters are assumed to be known. For the estimated channel results, we calculate the mean standard deviation (MSD), where MSD = *E* |**w**ˆ − **w***true*| 2 , **w**ˆ is the

estimated channel response and **w***true* is the true channel response. For the performance evaluation, we simulate the algorithms in the sparse channels for S = 4, 8, 16, and 32.

Figure 3 illustrates the MSD curves. For S = 4, Figure 3a shows that the estimation performance of *l*1-RTLS using the regularization factor with the *ρ* in Table 1 is almost the same as the *l*1-RLS using regularization with a true channel impulse response. However, the performance of *l*1-RLS using the regularization factor with the *ρ* in Table 1 is gradually degraded and shows a kind of uncertainty accumulation effect. In the other cases of S, we can observe the same trend in the MSD curves. Therefore, we can confirm that the new regularization factor selection method and the new modeling for *l*1-RLS can estimate the sparse channel as good as *l*1-RLS using the regularization with the true channel impulse response. In all the simulation scenarios, oracle RLS algorithm produces the lowest MSD as expected.

**Figure 3.** Steady-state MSD for S = 4, 8, 16, and 32 when applying the new *ρ* method to the regularization factor (-o-: *l*1-RLS with the true channel response, -×-: *l*1-RLS with the new *ρ* method, -\*-: proposed *l*1-RTLS with the new *ρ* method, --: oracle-RLS).

Table 2 summarizes the steady-state MSD values as varying S from 4 to 32. The results show that the proposed *l*1-RTLS with the new *ρ* is comparable to *l*1-RLS with the true channel.

In the second experiment, we compare channel estimation performance using room impulse response. The size of the room is (7.49, 6.24, 3.88 m). The position of the sound source is (1.53, 0.96, 1.12 m) and the position of the receiver is (1.81, 5.17, 0.71 m), respectively. T60 is set to 100 ms and 400 ms. The impulse response of the room is generated using the program in [31]. We focus on the direct reflection part and the early reflection part in the RIR because the direct reflection and early reflection part of the RIR has a sparse property. This is the part that is estimated in the AEC applications [32]. This part is also related to localization and clarity in room acoustics [33–35]. Comparing the impulse response (IR) generated by setting T60 = 100 ms to the channel with 65 coefficients used in the first experiment, it is equivalent to S = 4 in the channel with 65 coefficients. In the same manner, the IR generated by setting T60 = 400 ms is equivalent to S = 10.

**Table 2.** MSD (mean square deviation) comparison.

Table 3 summarizes the steady-state MSD values. The results also show the same trend as Table 2. In RIR estimation, the proposed *l*1-RTLS with the new *ρ* is also comparable to *l*1-RLS with the true channel.

**Table 3.** MSD (mean square deviation) comparison in sparse RIR estimations.

#### **8. Conclusions**

In this paper, we have proposed the regularization factor for recursive adaptive estimation. The regularization factor needs no prior knowledge of the true channel impulse response. We have also reformulated the recursive estimation algorithm as *l*1-RTLS type. This formulation is robust to the uncertainty in the regularization factor without a priori knowledge of the true channel impulse response. Simulations show that the proposed regularization factor and *l*1-RTLS algorithm provide good performance comparable to *l*1- RLS with the knowledge of the true channel impulse response.

**Author Contributions:** Conceptualization, J.L.; Methodology, J.L.; Validation, J.L. and S.L.; Formal analysis, J.L.; Investigation, J.L. and S.L.; Writing—original draft preparation, J.L.; Writing—review and editing, S.L.; Visualization, J.L.; Project administration, J.L.; Funding acquisition, J.L. and S.L.

**Funding:** This research received no external funding.

**Acknowledgments:** This research was supported by Agency for Defense Development (ADD) in Korea (UD160015DD).

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

### *Article* **Evaluation of Cracks in Metallic Material Using a Self-Organized Data-Driven Model of Acoustic Echo-Signal**

### **Xudong Teng 1,2,†, Xin Zhang 3,†, Yuantao Fan 4 and Dong Zhang 1,\***

Received: 24 November 2018; Accepted: 22 December 2018; Published: 28 December 2018

**Abstract:** Non-linear acoustic technique is an attractive approach in evaluating early fatigue as well as cracks in material. However, its accuracy is greatly restricted by external non-linearities of ultra-sonic measurement systems. In this work, an acoustical data-driven deviation detection method, called the consensus self-organizing models (COSMO) based on statistical probability models, was introduced to study the evolution of localized crack growth. By using pitch-catch technique, frequency spectra of acoustic echoes collected from different locations of a specimen were compared, resulting in a Hellinger distance matrix to construct statistical parameters such as *z*-score, *p*-value and T-value. It is shown that statistical significance *p*-value of COSMO method has a strong relationship with the crack growth. Particularly, T-values, logarithm transformed *p*-value, increases proportionally with the growth of cracks, which thus can be applied to locate the position of cracks and monitor the deterioration of materials.

**Keywords:** crack growth; acoustic echo; COSMO; *p*-value

#### **1. Introduction**

Nonlinear ultrasonic behaviors, such as harmonics, mix frequencies, and the resonance frequency shift, have been proven to be sensitive to structure imperfections and early degradation of materials [1–6]. In the early stage of damage, material fatigue can induce a number of micro-cracks with a typical length of 1–100 μm by continuous loading cycles, then the micro-cracks further grow, coalesce with other micro-cracks and eventually form macro-cracks [7]. Since the fatigue cracks are localized [8], not uniformly distributed in the structure, the generated nonlinear response is basically dependent on the configuration of the crack area related to the localized hysteretic deformation. Clapping between the contacting surface and dissipative mechanism due to frictional sliding and so on [7,8], leads to a much stronger nonlinearity than the surrounding material behaving linearly [8]. However, non-linear effects induced by the localized cracks in the materials are not obvious enough to conveniently be measured and analyzed [9,10]. In addition, the use of power amplifier, transducers, and coupling media in ultrasonic testing system also bring about external non-linear change. Since it is difficult to separate structure-induced non-linearity in materials from external non-linearity, non-linear ultrasonic technology is not applied widely to accurately evaluate and locate the structure imperfections in practical applications [11].

In terms of the uncertainties in real-life testing conditions, "big data" sets, collected from acoustic echo-signals of vast amounts of damaged material, show different but distinguishable statistical characteristics compared with intact material [12–17]. Some statistical models, also called data-driven models, relate the degree of damage to the probability of detection (PoD). They were first introduced by the National Aeronautics and Space Administration (NASA) in 1973 and were soon accepted as a standard method [18,19]. Later Lu and Meeker further developed statistical methods to estimate a time–to–failure distribution for a broad class of degradation structures [16]. Gebraeel employed Bayesian to update a method to develop a closed-form residual-life distribution for the monitored device [20]. Gang Qi et al. proposed a framework to meet the challenge by systematically evaluating material damage based on large data sets collected by using acoustic emission (AE) [12]. Zhou et al. investigated AE relative energy, amplitude distribution as well as amplitude spectrum to discern the delamination damage mechanism of the composites [21]. K*u*˚s et al. employed the model-based Clustering (Hellinger divergence) method to classify certain attributes of the original pure data obtained directly from the acoustic emission signals and form normed frequency spectra to perform physical separation tasks of AE random signals [22].

However, most data-driven methods are probabilistic and obtained historical degradation data or empirical knowledge [23,24], which require a large amount of experimental data to construct the reference curves or the preset feature threshold. Hence, in our previous work [25], we first introduced consensus self-organizing model (COSMO), neither any domain knowledge nor supervision to extract useful features, detecting a single flaw located at a fixed positon of steel specimens based on acoustical echo-signals. Nevertheless, as the contacting surfaces of the flaw produced by an electrical discharge machine were totally separated in [25], evaluating the actual fatigue cracks by applying COMSO method would potentially be problematic. In this paper, four cracks distributed along the length direction of a steel specimen were investigated, and the crack's growth (produced using fatigue testing [26]) was further discussed in details, which is probably accord with the mechanisms involved in contact acoustic nonlinearity and hysteretic nonlinearity. Both Numerical simulation and experimental measurements showed that COSMO models are effective in NDT inspection, as well as health monitoring for regular metallic structures.

#### **2. The Consensus Self-Organizing Models (COSMO)**

The COSMO method identifies the typical variability within a group of systems and to evaluate the likelihood of any individual being significantly different from the majority. In an ultrasonic testing system, a group of acoustic echo-signals were collected using pitch-catch method from N different locations (i.e., the total amount of samples) on a measured object. The spectral density of the acoustic signals were then obtained so that the difference between spectral density of two acoustic signals were compared [27,28], and Hellinger distance *di,j* was employed [22–24]:

$$d\_{i,j} = \frac{1}{\sqrt{2}} \sqrt{\sum\_{i,j=1}^{m} \left(\sqrt{p\_i} - \sqrt{q\_j}\right)^2},\tag{1}$$

where *pi* and *qj* are normalized histograms, representing spectral density of two acoustic signals. *m* is the number of sampled points of each acoustic signal. Two acoustic signals with different spectral density will yield a greater value of Hellinger distance than two similar ones. From the perspective of clustering analysis, samples with larger sum of distances to all peers are prone to be outliers.

For all pairs of these histograms, Hellinger distance was computed, resulting in a symmetric distance matrix D:

$$D = \begin{pmatrix} 0 & d\_{1,2} & \cdots & d\_{1,\mathbb{N}} \\ d\_{2,1} & 0 & \cdots & d\_{2,\mathbb{N}} \\ \vdots & \vdots & \vdots & \vdots \\ d\_{\mathbb{N},1} & \cdots & d\_{\mathbb{N},\mathbb{N}-1} & 0 \end{pmatrix}.$$

Within D, the row with minimum sum of distances is chosen as the most central pattern *c*, representing the most typical testing sample within this group. Based on the most central pattern *c*, the *z*-score is computed for each sample *k*, representing the percentage of all samples that are further away from *c* than sample *k*, that is,

$$z(k) = \left| \frac{\{i = 1, 2, \ldots, \text{N:} d\_{i, \varepsilon} < d\_{k, \varepsilon} \}}{\text{N}} \right| \quad (k = 1, 2, \ldots, \text{N:}) \; , \tag{2}$$

where *dk,c* is the distance between sample *k* and *c*; *di,c* denotes the distance between sample *i* and *c*; N is the total amount of samples. A sample with *z*-score close to 0 indicates a large distance to its peers. If there is no micro-crack at a particular position, the *z*-scores of a set of testing samples from this position should be uniformly distributed between (0, 1). We approximate the distribution of the average based on *n* samples using a normal distribution with mean 0.5 and variance 1/12n, i.e., [23,24]

$$z \sim N[\frac{1}{2}, \left(1/12n\right)^{1/2}].\tag{3}$$

Therefore, the one-sided *p*-value, which is the probability that a single observation *z* picked from a normal distribution with parameters (1/2, (1/12n)1/2) will fall in the interval (−∞, *z*], can be computed as:

$$p = \frac{2}{\sqrt{2\pi}} \int\_{-\infty}^{\infty} e^{-\frac{\left(z - 1/2\right)^2}{2\left(1/\sqrt{12n}\right)^2}} dz. \tag{4}$$

Based on a uniformity test of *z*-scores over an area, the resulting *p*-value of this test is obtained to estimate whether the inspected contains any micro-flaws.

#### **3. Crack Identification Based on COSMO**

In this section, we employed the COSMO model to identify and locate the cracks to steel specimens by numerical simulation and experimental tests.

#### *3.1. Numerical Simulation*

A two-dimensional model of a steel board embedded with a single crack was simulated by using a commercial software (Comsol Multiphysics V4.3a. COMSOL, Inc., Palo Alto, CA, USA). Figure 1 shows the schematic illustration of an ultrasonic measurement system in the simulation. Two longitudinal transducers with 60◦ wedge as the transmitter were typically used to carry out the ultrasonic inspection on the top surface of the steel board with the length of 220 mm and the height of 30 mm. A fixed distance between the two transducers was kept at a constant interval to make sure the first back wall echo was fully collected by the receiver. The single crack with a depth of *d* was located at the middle of the steel board along the x direction, i.e., x = 110 mm.

**Figure 1.** Schematic illustration of scanning on a simulated specimen with a single crack.

A 0.5 MHz continuous sinusoidal signal with signal-noise ratio (SNR) of 15 dB was applied as the exciting signal. Both the transmitter and the receiver were moved simultaneously to scan the steel board along the x direction. The receiving waveforms, at eight different positions, were spaced by 20 mm on the top of the simulated steel board, recorded, and the corresponding spectrums were then analyzed. Figure 2a depicts spectrums at eight positions of the simulated steel board with crack length *d* = 2 mm.

**Figure 2.** (**a**) Spectrums and (**b**) *z*-score distribution of echoes at eight monitoring points (depth of crack d = 2 mm) along x direction.

According to the COSMO algorithm model, the corresponding spectrums, at eight different observation positions, were saved then a group of *z*-scores were calculated by Equation (2) after every scanning process, finally 30 groups of *z*-score were obtained by scanning repeatedly 30 times. It is clearly shown in Figure 2b that *z*-score of observation points are distributed almost evenly between 0.3 and 1 except *x* = 110 mm, while most of *z*-score at *x* = 110 mm are mainly distributed below 0.4, just right at the crack's position. It is shown that the distribution of *z*-scores could be used to locate and identify cracks or defects in materials, i.e., *z*-score of damaged regions might be below 0.4. However, the conclusions need to be subjected to hypothesis tests to reach statistical significance, which determines whether a null hypothesis can be rejected or retained.

Figure 3a shows the calculated level of significance testing for crack depth of 2, 4, and 5 mm by Equation (4). It could be seen from Figure 3a that the *p*-value has a much smaller value than 0.1 around the crack region from 90 to 150 mm, which suggests that the imperfect structure of this region is significant. To make the comparisons and analysis clearly, an indicator called deviation level is defined as,

$$\mathbf{T} = -\lg(\mathbf{p}),\tag{5}$$

i.e., logarithm transformed *p*-value, obviously the small T-value indicates little significant probability of crack. The T-value curve of significance testing is shown in Figure 3b. Obviously, the maximum T-value occurs around the position of crack (*x* = 110 mm) for crack depth of *d* = 2, 4 and 5 mm, respectively. Furthermore, the maximum T-value increases simultaneously as crack growth, e.g., the maximum T-value is close to 4 when crack depth is 2 mm, and the peak of T-value up to 10 when crack depth equal to 5. This result indicates that the higher level T-value is strongly correlated with the crack depth, which could become an index to exhibit the evolution of crack growth inside materials.

**Figure 3.** The calculated level of significance testing (**a**) *p*-value curve and (**b**) deviation level T-value curve can be plotted for crack depth of 2, 4 and 5 mm by Equations (4) and (5) along *x* direction on the simulated board.

When the crack depth changed from 0 mm to 5 mm, the peaks of T-value around the cracks were obtained and thus the relationship curve between the maximum T-value and crack depth d was plotted, as shown Figure 4. It can be seen that as the crack depth increased from 0 to 1, the slope of the curve sharply increased. When crack depth was less than 1 mm, the T-value was not larger than 5, which is basically considered the formation stage of crack, due to the relative small change of crack depth, thus this phase is called stage I. As the crack depth gradually expanded form 1.5 mm to 3.7 mm, the maximum T-value increased slowly from 5 to 7, at stage II. When crack depth was larger than 4 mm, at stage III, the value of the curve increased rapidly up to 10, and as high as 2 times than that of stage I, which means that the small cracks had already expanded to macro-cracks. Therefore, the peaks of T-value might track the progression of damage and evaluate the evolution of crack growth.

**Figure 4.** Deviation level T-value grows monotonically as depth of crack increases.

#### *3.2. Experimental Measurement*

A specimen made of Q235 (See Table 1) with dimensions 800 mm × 250 mm × 20 mm was used in the experimental measurement, as shown in Figure 5. Four sections embedded with cracks with average depth of 6 mm, 2 mm, 1 mm, and 0.5 mm were manufactured in the specimen, mainly located at 150 mm, 300 mm, 450 mm, and 600 mm, respectively, denoted by B, C, D and E. Additionally, Section A and F represented as undamaged regions located in the two ends of the specimens. A portable TOFD ultrasonic detector (PXUT-920, Nantong Union Digital Tech., China) was used to excite a narrow-pulse

acoustic signal with 200 ns in width and stored the echo signals from inspected cracked region. Scanning was manually carried out by a scanner unit with one pair of 5 MHz normal transducer, i.e., the transmitter and the receiver, with 60◦ wedges for longitudinal waves. Two transducers, spaced 62 mm apart, were located at equidistant over the crack region center, and scanning was done by moving the scanner in the length direction of steel plate parallel to the crack region. The echo signal sampled by the detector, containing 1496 points, was acquired every 0.5 mm along the length direction. After a scanning, a total of 1600 echo signals (A-scan) were obtained and stored in the ultrasonic detector.

**Table 1.** Mechanical characteristics of Q235 carbon steel. (Provided by HBO Windpower Equipment Co., Ltd, Nantong, China).

**Figure 5.** (**a**) A specimen with four cracks with average depth of (**b**) 6 mm, (**c**) 2 mm, (**d**) 1 mm, and 0.5 mm located at Section B, C, D, E, respectively. Section A and F without cracks.

COSMO method was applied to analyze this dataset of echo-signal recorded by TOFD ultrasonic detector. Firstly, the Hellinger distance matrix D was constructed using Equation (1), and the row with minimum sum was chosen in the metric D so that the *z*-score could be determined using Equation (2). Figure 6a shows the *z*-score distribution for 30 scanning. It is shown that the *z*-score in undamaged sections is much larger than those in the region with cracks. For example, the *z*-score for positions A and F are about 0.4~1, while those for positions B, C, D and E are below 0.4. The results suggest that the *z*-score is closely related to cracks of specimen, just as the simulated results.

**Figure 6.** (**a**) *z*-score distribution (**b**) the curve of T-value of specimen with cracks and along *x* direction.

By using Equations (1)~(4) and (5), the deviation level T-Value is calculated to make a significant analysis. It can be observed from Figure 6b that the T-value at B, C, D and E are of high level compared to those in uncracks region. For instance, the T-value at x = 150 up to 13, and T- value at x = 600 sharply increasing from 0 to 7, but T-value of sections without cracks almost equal to 0, far less than T-value at crack region. In addition, T-value increase almost linearly with the depth of cracks. The relationship curve between the peaks of T-value and crack depth is plotted in Figure 7. It is shown that the peak of T value increases with crack growth from 0 to1.5 mm quickly up to 5, which is exactly in stage I. When the crack depth is larger than 2, the slope of curve go slow but still faster than the simulated results. It is worth noting that the change of slope is not distinct enough to easily recognize stage II or III when depth of crack exceeding 1.5 mm, different from the simulated curve, which might attributed to the result of multi-physical mechanisms.

**Figure 7.** The relationship curve of specimen between T-value peaks and crack depth.

#### **4. Conclusions**

Defects or cracks can significantly increase acoustic non-linearity, and the nonlinear acoustical parameters, and thereby can be exploited to evaluate the state of material damage. However, the harmonics are usually too weak to be detected in early fatigue. Therefore, the non-linear ultrasonic technique is rarely used to qualify crack growth. In this work, the COSMO method was applied to compare the spectrum of different positions by ultra-sonic scanning in order to obtain the distribution

of *z*-scores as well as the corresponding significance level in every scanning position. The results show that: (1) the *z*-scores in the location with cracks are distributed below 0.4, while the *z*-scores in the location without cracks are above 0.4; (2) the deviation level T - value at locations with crack are much larger than those at locations without cracks, and the T-value would get larger with the increase of crack depth; (3) based on the quantitative relations between the T values and the crack depth, we can evaluate and monitor the online state of the structural health by COSMO model. However, it is noted that the COSMO model is still a simple model that does not consider some other factors, such as the shape and mechanical properties of structure, as well as the requirements on data size. Therefore, the reliability of COSMO needs to be further optimized to reach a solution for non-destructive evaluation in future.

**Author Contributions:** D.Z., conceived and designed the experiments. X.T., and X.Z., performed the experiments. X.T., X.Z., and Y.F., analyzed the data. X.T., X.Z., and D.Z., wrote the manuscript.

**Funding:** This research was funded by the National Natural Science Foundation of China (Grant no's., 81627802, 11674173 and 11874216), QingLan Project, and the Fundamental Research Funds for the Central Universities.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

© 2018 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

### *Article* **Influence of Piano Key Vibration Level on Players' Perception and Performance in Piano Playing**

#### **Matthias Flückiger \* , Tobias Grosshauser and Gerhard Tröster**

Electronics Laboratory, ETH Zurich, Gloriastrasse 35, 8092 Zurich, Switzerland<EMAIL_ADDRESS>(T.G<EMAIL_ADDRESS>(G.T.)

**\*** Correspondence<EMAIL_ADDRESS>
Received: 29 October 2018; Accepted: 16 December 2018; Published: 19 December 2018

**Abstract:** In this study, the influence of piano key vibration levels on players' personal judgment of the instrument quality and on the dynamics and timing of the players' performance of a music piece excerpt is examined. In an experiment four vibration levels were presented to eleven pianists playing on a digital grand piano with grand piano-like key action. By evaluating the players' judgment of the instrument quality, strong integration effects of auditory and tactile information were observed. Differences in the sound of the instrument were perceived by the players, when the vibration level in the keys was changed and the results indicate a sound-dependent optimum of the vibration levels. By analyzing the influence of the vibration levels on the timing and dynamics accuracy of the pianists' musical performances, we could not observe systematic differences that depend on the vibration level.

**Keywords:** piano playing; vibrotactile feedback; interaction; musical performance; auditory perception; sensors; actuators

#### **1. Introduction**

Playing the piano is a complex multi-modal task, where the pianist controls the instrument through his or her intention and perceived instrument feedback. There are four main musician-musical instrument interaction modalities: visual feedback, auditory feedback, force feedback, and vibrotactile feedback. The interaction with a musical instrument can be modeled as a feedback controller [1], where the musician's brain controls his or her body, arms, and fingers to modify the instrument's behavior based on changes in sensory inputs. This closed-loop model implies that if the instrument's feedback is altered, the pianist will adapt his or her playing to compensate for and retain the desired instrument behavior. Vibrotactile feedback can support the precise control of finger force, as shown by Ahmaniemi [2] with a basic force repetition experiment on a rigid sensor box. Furthermore, Goebl and Palmer [3] demonstrated that tactile sensations from the finger-key surface interaction support some pianists to improve timing accuracy and precision of finger movements.

In piano playing, vibrotactile feedback is perceived through the fingers in contact with the keys and the feet in contact with the pedals. The keybed and soundboard vibrations excite the piano keys and pedals [4]. Askenfelt and Jansson [5] measured the vibrations of a depressed piano key and a depressed piano pedal. Piano key vibrations comprise broadband and tonal parts [4,6]. The tonal parts come from the string vibrations, and the broadband parts come from mechanical impacts (e.g., hammer-string impact and key–keybed impact) of the piano action when the piano key is played [4,6].

The levels of the vibrations' tonal part can rise to the micrometer range [7]; the vibrations are close to the limits of human vibration perception and are often sensed subconsciously [8]. However, piano key vibrations can be detected up to the middle octave of the keyboard [9]. Further, the vibration levels vary considerably among different pianos, but it remains unclear if pianists can perceive these differences [7].

Keane and Dodd [8] found that the ratio of broadband and tonal parts of piano key vibrations influenced the instrument's perceived sound. An upright piano was mechanically modified to reduce the broadband parts' amplitude, with the expectation to improve the instrument's quality. The pianists preferred the modification with regard to tone and loudness in an evaluation study. Interestingly, the participants did not report differences with regard to touch or vibrations.

Fontana et al. [10] showed by ratings of evaluation criteria that realistic piano key vibrations rendered on a digital keyboard are preferred to a no-vibration condition. In the same study, key vibrations did not show a significant effect on pianists' timing and dynamics accuracy during a scale playing task.

In addition to the state-of-the-art, the influence of four piano key vibration levels on pianists' personal judgment of an instrument's sound, control, and feel is investigated in this study. We designed the experiment, such that the control of vibration levels was independent of the sound of the instrument and aimed to explore connections between vibrotactile feedback and the perceived quality of the instrument. To test if pianists adapt their playing to vibrotactile feedback and to analyze if the vibrations support the control of finger forces, the effect of the vibration levels on timing and dynamics accuracy in pianists' performances is also studied in this paper.

#### **2. Methods**

#### *2.1. Equipment*

The pianists played on an AvantGrand N3X, a digital hybrid grand piano from Yamaha. This instrument was chosen because it simulates piano key vibrations, features state-of-the-art grand piano sound-rendering algorithms, and has a piano action resembling that of acoustic grand pianos. Musical instrument digital interface (MIDI) messages and the headphone audio output of the instrument were recorded. The pianists played with closed-back headphones to block the small amount of sound that vibrating keys radiate.

#### *2.2. Experiment Design*

The target was to control the key vibration levels independent of the sound and to cover the level range of piano key vibrations of acoustic concert grand pianos. Independent control could not be achieved with the built-in vibrotactile feedback rendering system of the AvantGrand N3X; therefore, it was extended as illustrated in Figure 1.

The mono audio output signal of the AvantGrand N3X was processed with a digital signal processor (DSP). Through a combined approach of vibrometer measurements and subjective evaluation by playing on the instrument, the DSP's filter stage was tuned to create vibration level *V*3 (see Figure 1), which approached the maximum vibration levels previously measured on acoustic grand pianos [7]. After implementation of *V*3, vibration levels *V*2 and *V*1 were created by attenuating the signal in steps of 6 dB. The no-vibration condition *V*0 completed the levels of the experiment. As shown in Figure 1, the vibration levels cover the range of acoustic grand pianos for notes A2, A3, and A4. For notes A0 and A1, the levels are more than 10 dB lower. The chosen music piece avoided the lowest notes. The deviations of the vibration level curves in Figure 1 are due to non-idealities of the excitation system.

The experiment was created to study the influence of four vibration levels on players' personal judgment of the instrument quality and on the dynamics and timing of the players' performance of a music piece excerpt. The experiment was designed so that the participants were unaware of the independent variable, and the session was split into three parts to steer the pianists' attention to different instrument properties. Free verbalizations were used to assess the players' judgment, allowing for unrestrained and possibly unexpected answers. Since a small influence of the piano key vibrations on musical performance and the players' judgments was assumed—as natural levels are close to the threshold of vibration sensation [7–9]—numerous repetitions were included in the protocol, and participants with high levels of playing experience were selected.

(**a**) Block diagram of the extended vibrotactile feedback rendering system

(**b**) Comparison of the measured tonal part of the vibration levels to levels measured on acoustic grand pianos

**Figure 1.** (**a**) Block diagram of the vibrotactile feedback rendering system to generate the key vibrations; the mono audio output of the N3X was filtered and attenuated with a DSP. Thereafter the signal was power amplified to drive the transducer of the built-in vibrotactile rendering system of the N3X. (**b**) Comparison of the tonal part of the vibration levels *V*1, *V*2, *V*3 to vibration levels of four acoustic concert grand pianos. The comparison is based on vibrometer measurements of *forte* keystrokes [7]. (*V*0 is not shown because it corresponds to no vibrations.)

#### *2.3. Participants*

Eleven pianists participated in the study: seven piano students, 22–26 years of age, with an average playing experience of 17 years; and four professional pianists, 31–40 years of age, with an average playing experience of 26 years. None of them reported having auditory or tactile impairments.

#### *2.4. Procedure*

The session for each participant lasted around 1.5 h. The participants were asked to prepare an interpretation of a music piece excerpt. The excerpt was 15 bars long and the participants were instructed to adhere to the tempo, dynamics, accents, and pedaling information. The excerpt was taken from Klage by Gretchaninov [11] and was edited to cover the dynamic range from *pianissimo* to *fortissimo* (see Figure A1 in the Appendix A). The participants were not informed about the purpose of the experiment beforehand. The participants were only told that their judgments of various settings of the instrument will be evaluated.

The experiment comprised a warm-up (with a duration of around 5–10 min), questionnaires, and three parts (A, B, and C) with a duration of roughly 20 min each.

During the warm-up, the pianists were free to play whatever they wanted and were instructed to evaluate the instrument in a way comparable to choosing an instrument for a concert or for purchase. After this familiarization, the pianists were asked to express their first impression by answering a set of questions about the sound and touch of the AvantGrand N3X and by comparing the instrument to their main instrument.

During the main parts of the experiment (A, B, and C) the pianists were asked to repeat the excerpted music piece accurately in 12 direct comparisons of different instrument settings. After each trial the pianists were asked to indicate a personal preference ("better", "worse", or "similar") of the current setting relative to the previous setting and to describe their impression in a few words. The participants were told that the differences between the comparisons can be small and that some might be perceptually irrelevant. In part A the pianists were told that a slight adjustment of the instrument (not further specified in order to not suggest an answer or category) was made between each repetition, in part B it was claimed that a small adjustment to the sound was made between each

trial, and in part C the pianists' attention was directed to the keyboard by asking about the instrument's control and feel.

In fact, the only independent variable throughout the experiment was the key vibration level (*V*0, *V*1, *V*2, *V*3). The sequence of parts (A, B, and C) was the same for all participants. Each part had a different randomized sequence of vibration levels. All pair of levels were compared twice in each part—once for each order. In total, nine trials per vibration level and per participant were recorded.

At the very end of each experiment session, personal information was collected and the participant was asked about his or her experience and preference of piano key vibrations in piano playing, before we disclosed and explained the purpose of the experiment.

#### **3. Results**

#### *3.1. Influence of Vibration Levels on Perceived Instrument Sound, Control, and Feel*

In the analysis of the preference ratings ("better", "worse", or "similar"), a high variance and for some participants also controversial ratings were observed. We decided to present the ratings across all parts and participants here, because it was not possible to draw conclusions from the ratings per participant or per part. The result is presented in Figure 2.

**Figure 2.** Analysis of the preference ratings ("better" (1), "worse" (−1), or "similar" (0)) of the vibration levels across all parts and all participants. The ratings are based on direct comparisons between all pairs of levels. The ratings are relative. For example, a positive value for *V*0:*V*1 indicates a preference of *V*0 over *V*1 and a negative value a preference of *V*1 over *V*0. The shaded area marks the standard deviation of the ratings.

The high variance of the ratings in Figure 2 reflects the closeness of the key vibration levels to the limits of human perception. Additional factors that can have disturbed the ratings are the mood and fatigue of the player. Also a self-evaluation of the playing, the difficulty of the task, or the imposed expectation of a difference between the settings might have disturbed the ratings. However, a visual comparison of the mean values in Figure 2 indicates a tendency in the preference of the players toward vibration level *V*2. The preference of vibration level *V*2 is confirmed by the evaluation of the player's verbal self-reports presented hereafter.

The free verbalizations were analyzed with an approach presented by Pate et al. [12], where concepts by Dubois [13] were applied to musical instrument evaluation.

Based on the context and for each participant, the meaning, category, and preference of each statement was identified via linguistic tools such as reformulations, oppositions, and comparatives. Thereafter, the statements were classified into positive and negative statements for three categories: sound, control, and feel. Statements covering multiple categories were split before classification. Statements not indicating a preference or a perceived difference were counted as "no difference". The results are summarized in Table 1. Significance was evaluated with Pearson's *χ*2-tests at a confidence level of 95%. The test was performed on all statements (positive and negative) per category.

**Table 1.** Evaluation per vibration level derived from the free verbalizations of the pianists. The number of positive and negative statements per category was counted for each vibration level (*V*0, *V*1, *V*2, *V*3). The frequency counts were evaluated with *χ*2-tests. The *χ*2-statistics and *p*-values are given for each category; *p* < 0.05 is highlighted in bold.

Although we did not alter the sound throughout the experiment, Table 1 shows that the vibration levels (*V*0, *V*1, *V*2, *V*3) have an influence on the pianists' sound perception. Pairwise testing with Bonferroni correction showed that the significance of the *χ*2-test in the sound category arises from the difference between *V*0 and *V*2.

The phenomenon of vibrotactile feedback causing a difference in sound perception is known as integration of auditory and tactile information [14–16] or weak synesthesia [17]. The preference of *V*2 over *V*3, confirming the evaluation of the preference ratings in Figure 2, was surprising because vibration level *V*3 is closer to the levels of acoustic instruments (see Figure 1). An explanation lies in five statements about vibration level *V*3, which were classified as negative. These statements criticized the balance of the perceived sound as having "too much bass" or being "unbalanced".

The results for control and feel are not significant according to the *χ*2-tests. For vibration level *V*1 there were twice as many "no difference" statements than for all other levels, which indicates that *V*1 is most difficult to differentiate.

Only two participants consciously noticed a change in vibration levels during the experiment, when vibration level *V*3 was compared to *V*0 and vice versa. Both recognized the vibrations during the last part, when the keyboard was the focal point.

To find possible explanations for the above presented differences, we analyzed the verbal self-reports of the participants in more specific categories. We observed that the key vibrations influence the timbre and the perceived loudness of the bass keys. Also, the timbre of treble notes was judged more pleasant when playing with *V*2 or *V*3. Some participants also noted a sensation of space when playing with higher vibration levels (*V*2, *V*3) and described it as room or reverb effect of the sound. In contrast when pianists played with vibration levels *V*0 or *V*1 the sound was sometimes described as dry. Comparisons to acoustic instruments and e-pianos also align with this observation.

A critical aspect for discussion is that the sequence of parts (A, B, and C) was the same for all participants, which might have influenced the results. We designed the experiment protocol to steer the players' attention to different multi-modal aspects to discover unexpected connections between vibrotactile feedback and the players' judgment of the instrument quality.

In part A the participants could freely describe their impressions and we did not suggest any quality criteria for the comparisons. Unbiased comparisons are only possible within the first part of the experiment. Ten out of eleven participants naturally made statements about sound and control in the first part, which justifies the suggestion of these criteria in the following parts.

We decided to put part C at the end of the experiment session, because we did not want to risk that the participants are already consciously aware of the vibrations, when comparing the levels with regard to sound. In part C the participants focused on the keyboard. Therefore we expected that it is most likely that the participants recognize the vibrations in this part (which happened in two cases).

Finally, also the difficulty of the evaluation task might have altered the judgments of the pianists. In each trial, the participant played the music piece excerpt for a duration of around 30 s, communicated his or her impression with regard the previous setting and sometimes also answered clarifying questions from the experimenter. Thereafter he or she performed the music piece for the next comparison.

#### *3.2. Influence of Piano Key Vibration Levels on Musical Performance*

To analyze the MIDI-based performance data, a custom data structure was used. The structure groups notes played at the same time (±40 ms) into clusters, removes accidentally played wrong notes, and assigns if the note was played by the left or right hand.

Key velocity *υ*, a measure of a keystroke's excitation strength, was directly extracted from the MIDI messages. For the calculation of the inter-onset interval *τ*—the time interval between two subsequent note onsets—only notes played by the left hand were considered. The tempo of each trial was normalized.

To compare trials by the distribution of key velocity *υ* (analogously for inter-onset interval *τ*) histogram intersection was used. Histogram intersection was introduced by Swain and Ballard [18] to identify objects by color similarity in computer vision. Histogram intersection is defined as [18]

$$H\_1(\upsilon) \cap H\_2(\upsilon) = \sum\_{i=1}^n \min\left(h\_{1i}(\upsilon), h\_{2i}(\upsilon)\right),\tag{1}$$

where *H*1 and *H*2 represent two trials by normalized discrete distributions of key velocity *υ* with *n* bins *h*1*i*, *h*2*i*. Equation (1) measures the overlap of two histograms in the range [0, 1]. The number '1' corresponds with perfect overlap; '0' means no overlap. In contrast to an evaluation based on mean values only, histogram intersection also identifies differences in the distributions' shape or offset.

To judge significance, two tests were demanded to reject the null hypothesis: non-parametric Friedman analysis of variance with a 95% confidence level in combination with pairwise Wilcoxon signed-rank tests with Bonferroni correction. We used non-parametric tests, because the evaluated quantities do not necessarily follow a normal distribution and because of the sample size.

Two approaches were used to compare the distributions of key velocity *υ* and inter-onset interval *τ* by histogram intersection. The distribution of both parameters was calculated for each trial and was analyzed for each participant separately.

The first approach considered if the pianists adapted their playing to the vibration levels (e.g., if a pianist perceived an overemphasis of bass notes and therefore played the bass notes with less finger force than before). This force adaption manifests in the shape of the distribution of key velocity *υ*. An increase in the "amount of adaption" was expected with increasing vibration levels. The "amount of adaption" for key velocity *υ* was measured as follows.

Let *HPk*,*V*-,*i*(*υ*) denote the normalized histograms describing the distribution of key velocity *υ* for trial *i* ∈ {1, ... , 9}, vibration level *V* with - ∈ {0, 1, 2, 3}, and pianist *Pk* with *k* ∈ {1, ... , 11}. Then the "amount of adaption" *APk*,*Vn* of pianist *Pk* to vibration levels *Vn* with *n* ∈ {1, 2, 3} was estimated as the histogram difference relative to *V*0 condition *HPk*,*V*0,*i*(*υ*) ∩ *HPk*,*Vn*,*j*(*υ*) for all combinations of trials *i*, *j* ∈ {1, . . . , 9} and key velocity *υ*. For the inter-onset interval *τ* the same procedure was conducted.

Figure 3 shows the "amount of adaption" of the pianist's playing to the feedback levels for both performance parameters. The differences in Figure 3 are not significant. There is no general tendency that the participants adapt their playing to the key vibration level. Nonetheless, by analyzing the

influence of the vibration levels per participant individually, three participants showed significant differences for key velocity *υ* and three for inter-onset interval *τ*.

(**a**) "Amount of adaption" *APk* ,*Vn* (*υ*) of key velocity *υ*

(**b**) "Amount of adaption" *APk* ,*Vn* (*τ*) of inter-onset interval *τ*

**Figure 3.** "Amount of adaption" to the feedback levels *Vn* (*n* = {1, 2, 3}) relative to the no-vibration condition *V*0 across all pianists *Pk* and for both performance parameters. The differences in the amount of adaption for all vibration levels are not significant for both parameters. The line in the center of the box-plot marks the median, the box extends from the first to the third quartile, and the whiskers mark the value range.

Possible explanations for the majority of pianists not adapting their playing to the feedback levels include that the combined task of playing, judging the impression, and adapting their playing was too difficult, that the levels were too small to cause a reaction, or that the method was not accurate enough to unveil such differences.

The second approach investigated how accurately the pianists could repeat the music piece excerpt when playing with different vibration levels. If key vibrations support the precise control of finger forces, a lower variance in the distribution of key velocity *υ* could be expected. In consequence the shape of the distribution of key velocity *υ* would be altered and hence a difference in repeatability could be detected. Likewise, if the pianist's tempo was more stable, a different shape of the distribution of the inter-onset interval *τ* would occur. Indirect causes are also possible (e.g., the pianist feels more comfortable to play and therefore plays with higher repeatability). The time-point of the trials during the experiment was not taken into account. We decided to analyze and present the data for each pianist individually hereafter, because we observed a strong dependency of the repeatability on the player.

The repeatability *RPk*,*V*0 (*υ*) for participant *Pk* playing with vibration level *V*0 was computed by comparing the distributions of key velocity *υ* by *HPk*,*V*0,*i*(*υ*) ∩ *HPk*,*V*0,*j*(*υ*) for all combinations of trials *i*, *j* ∈ {1, ... , 9}, where *i* = *j*. For vibration levels *V*1, *V*2, and *V*3, and for inter-onset interval *τ* similar procedures were conducted.

The resulting repeatability per vibration level and per participant is presented in Figure 4 for key velocity *υ* and inter-onset interval *τ*. The differences in repeatability were significant for a majority of the participants. However, no consistent tendency or pattern in repeatability occurred among the pianists in Figure 4.

Therefore, the vibration levels of our experiments do not have a conclusive influence on the pianists' repeatability, and the measured MIDI data do not support the hypothesis that key vibrations assist the precise control of finger force. Consequently the observations do not confirm the results of Ahmaniemi [2] or Galica et al. [19], for the piano playing case. Galica et al. [19] showed that even unconscious vibratory stimulation applied to the soles of the feet can cause lower variance in kinematic interactions.

**Figure 4.** Estimated repeatabilities per participant *Pk* and vibration level *V* for both performance parameters. No consistent tendency occurred among the pianists, but the vibration levels had a significant influence (marked with \*) on repeatability for a majority of the participants.

In summation, the repeatability estimates *RPk*,*V*- (*υ*), *RPk*,*V*- (*τ*) for vibration levels (*V*0, *V*1, *V*2, *V*3) depended on the player. The pianists in this study were more accurate in repeating key velocity *υ* (median of *RPk*,*V*- (*υ*) ∈ [0.85, 0.88], ∀*k*, -) than in repeating inter-onset interval *τ* (median of *RPk*,*V*- (*τ*) ∈ [0.66, 0.77], ∀*k*, -). For the inter-onset interval *τ* the intra-individual variance was also considerably larger (see Figure 4), although we normalized the tempo of each trial before the analysis.

As a concluding aspect of interest, no categorical differences (in repeatability or playing adaption to feedback levels) were found between the group of students (*P*1 to *P*7 in Figure 4) and the group of professional pianists (*P*8 to *P*11 in Figure 4).

#### **4. Conclusions**

By systematically investigating the players' personal judgment of the instrument quality of the vibration level in the keys, we observed strong integration effects of auditory and tactile information. The results give an illustration of the strong multi-modal effects in piano playing. The subjects perceived differences in the sound of the instrument when the vibration level in the keys was changed. The preference of vibration level *V*2 over *V*3 indicates an optimum or a "sweet spot" of piano key vibration levels, which depends on the instrument's sound and sound balance.

In line with the results of Keane and Dodd [8], the vibration levels in this experiment significantly affected the instrument's judged sound quality but not its control and feel. However, in contrast to the design of the present experiment, Keane and Dodd [8] reduced the level of the broadband part of piano key vibrations of an acoustic instrument.

An interesting direction for future research is to determine the vibration level differences that pianists can differentiate. This would help to understand if an instrument can be identified based on its vibrotactile feedback only, while the instrument's auditory and force feedback are kept constant. Some participants in this experiment perceived a spatial sensation and described it as room or reverb effect on the sound when playing with higher vibration levels (*V*2, *V*3). For several applications it could be interesting to understand the conditions that can cause such an illusion.

We did not find systematic differences by analyzing the influence of the vibration levels on the timing and dynamics accuracy of the pianists' musical performances. We can not exclude that such an influence exists but with the proposed measures, "amount of adaption" and repeatability, we could not measure such a relation. Furthermore, the basic results of Ahmaniemi [2], that vibrotactile feedback assists the precise control of finger forces could not be confirmed in our case. For future studies, we suggest to include a larger number of participants. This could help to identify groups reacting similarly to key vibrations. Future studies might also include multiple experiments over a certain range of time to exclude influences of physical and mental state on the day of testing. Finally, an analysis on a note-by-note basis could clarify if, for example, the vibrotactile feedback of a long-lasting bass note helps the precise control of the dynamics in subsequent keystrokes. We could not generalize such a relation with the data of the presented experiment.

If our results can be confirmed on acoustic instruments, our findings of the perception part of the experiment suggest that piano manufacturers should design the vibrations in the piano keys in balance with the sound of the lower notes of the instrument. Furthermore, it would be interesting to investigate the just-noticeable difference of piano key vibration levels, which might possibly be around 6 dB. Further research in this area could help to answer the question, if the tonal parts of the Steinway and Sons and the Yamaha concert grand pianos (the tonal parts for notes A2, A3, and A4 differ by more than 6 dB [7]) can be differentiated based on their vibrotactile feedback only by the player. Of course in such an experiment the vibrotactile feedback should be rendered on the same instrument, otherwise cues from the auditory or kinematic sensations might dominate the perceived impression.

**Author Contributions:** M.F. designed, accomplished, and evaluated this study. T.G. and G.T. contributed as consultants. All discussed the results.

**Funding:** This research has been pursued as part of the "Musician's behavior based on multi-modal real-time feedback" project, Grant No. 166588, funded by the Swiss National Science Foundation (SNSF).

**Acknowledgments:** The authors are indebted to Anders Askenfelt for contributing to the experiment's design, for inspiring discussions, and for offering advice about the topics of the presented study. They would also like to thank Yamaha for generously providing the AvantGrand N3X.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Appendix A**

**Figure A1.** Music sheet of the study. The excerpt was taken from Klage by composer Gretchaninov [11]. The excerpt was edited to cover a broad dynamic range and also to include accents.

#### **References**

c 2018 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

### *Article* **MRI Compatible Planar Material Acoustic Lenses**

### **Daniel Tarrazó-Serrano 1 , Sergio Castiñeira-Ibáñez 1 , Eugenio Sánchez-Aparisi 2, Antonio Uris 1 and Constanza Rubio 1,\***

Received: 10 October 2018; Accepted: 13 December 2018; Published: 15 December 2018

**Abstract:** Zone plate lenses are used in many areas of physics where planar geometry is advantageous in comparison with conventional curved lenses. There are several types of zone plate lenses, such as the well-known Fresnel zone plates (FZPs) or the more recent fractal and Fibonacci zone plates. The selection of the lens material plays a very important role in beam modulation control. This work presents a comparison between FZPs made from different materials in the ultrasonic range in order to use them as magnetic resonance imaging (MRI) compatible materials. Three different MRI compatible polymers are considered: Acrylonitrile butadiene styrene (ABS), polymethyl methacrylate (PMMA) and polylactic acid (PLA). Numerical simulations based on finite elements method (FEM) and experimental results are shown. The focusing capabilities of brass lenses and polymer zone plate lenses are compared.

**Keywords:** MRI; Zone Plates; ultrasonic lenses

#### **1. Introduction**

The development of modulating and focusing energy systems has been a field of study of great interest for scientist and engineers. The lens is a devices that is able to perform this energy modulation. Lenses allow beam forming, control propagation and focusing the energy that impinges on them. These effects are produced by refractive and diffractive phenomena. Transmission efficiency is one of the most important aspects, particularly when low impedance contrast is presented between the lens and the host medium. Due to the wide versatility of the lenses, they have been used in different areas. For example, they have been applied in sonochemistry [1], construction [2] and the pharmaceutical industry [3].

The acoustic lenses, depending on the physics involved in the beam formation, can be divided into different groups, including refractive lenses and diffractive lenses. One example of lenses based on the refraction phenomenon are sonic crystal lenses made of periodic distributions of rigid cylinders [4]. Due to the subsonic sound speed inside the crystal, these lenses act similar to those in optical systems. Another example of this typology of acoustic lenses are those which modify the refractive index using labyrinths. These type of lenses are the so-called Gradient-Index lenses [5–7].

The other subtype of lenses, based on the diffractive phenomenon, conducts its behavior on the constructive interferences of the pressure field. An example of these types of lenses is the fractal lenses, which are able to generate different foci depending on their fractal geometrical properties [8]. Fresnel Zone Plates (FZP) have an improved focusing capacity. Among the different ways to implement FZPs, one of the most common and easiest is to alternate transparent and blocking zones, which results in a Soret type FZP [9]. To obtain these blocking areas, materials that are opaque to sound are required. This fact is accomplished by selecting materials that have a high impedance contrast with the host

medium. There are studies that have implemented Soret FZP (SZP) by ultrasounds based on these type of lenses [10].

A material that has a high impedance with respect to water and that allows for the creation of opaque zones to achieve a Soret type lens is brass. However, this type of material has limitations, especially when used in fields such as bioengineering. The use of acoustic lenses in medicine for high intensity focused ultrasounds (HIFU) treatment is one of the current lines of research. magnetic resonance imaging (MRI) is the technique that is most used for guiding HIFU treatment [11].

MRI is a technique used for soft tissue structure imaging in a non-invasive way. The image is obtained by aligning and relaxing the magnetic moments of the atoms of the introduced elements in the MRI. Tissues are exposed to a strong external time-independent magnetic field. Thus, metallic elements cannot be introduced in the resonance zone due to their interference in the image and because they could damage MRI-systems. To avoid interaction with the electromagnetic field, non-metallic materials should be used in the construction of lenses. One of these materials is polylactic acid (PLA) [12]. The MRI environment requires materials such as PLA for medical instruments and patient supports. Recently, PLA has been used for this purpose and its reliability has been shown [13]. The HIFU transducer is embedded within a specially designed table that fits into the MRI device. This integrated system, has a degassed water bath where the transducer is located. The patient lies over this system on [14]. Although the transducer and the lens must be immersed in this water bath, degradation of the PLA will occur over long-term immersion. PLA degrades in water after a period ranging from months to a year [15]. Therefore it must be taken into account that, in MRI systems, the lens and the transducer are not permanently submerged. After 20 to 25 min, the system is extracted from the water bath, and for this reason, the time of degradation due to being immersed in water can be prolonged considerably.

In this work, three lenses with three types of compatible materials with MRI environments are compared. In this sense, acrylonitrile butadiene styrene (ABS), polymethyl methacrylate (PMMA) and polylactic acid (PLA) materials are used. Furthermore, a SZP built in brass is compared. Although, this material is not MRI compatible, it is the nearest to the ideal SZP that can be implemented in real projects. In the comparison, a not compatible with MRI lens built in brass and an ideal Soret lens will be added. Results are obtained and compared both numerically and experimentally. Numerical results have been obtained using the commercial software COMSOL Multiphysics 4.3a by COMSOL Inc. (Sweden) [16]. In this work, it has been verified that the ratio of the transmission capacity that is related to the ratio of impedances of the medium and the lens, directly influences the focusing capacity.

#### **2. Methodology and Theoretical Analysis**

Fresnel zone plates are circular concentric structures, which are known as Fresnel regions. Every consecutive region has a *π* phase shift between them. This fact makes a coherent contribution to obtain high intensity levels at focal length (*FL*), which is the location in the axial coordinate where the focus is placed. The number of Fresnel regions is defined as *N*, this includes both opaque and transparent acoustic sections. The working frequency is defined as *f*0 and radial distances (*rn*) of each Fresnel zone can be obtained using Equation (1) valid for plane wave incidence.

$$r\_n = \sqrt{n\lambda F\_L + \left(\frac{n\lambda}{2}\right)^2} \qquad n = 1, 2, \ldots, N \tag{1}$$

In this work, underwater transmission is considered and lenses are designed for ultrasound applications. Therefore, piston sources have to be considered when FZPs are implemented. Due to spherical wave incidence consideration, Equation (2) has been used where *d* is the separation between the point source and the lens.

$$d + F\_L + \frac{n\lambda}{2} = \sqrt{d^2 + r\_n^2} + \sqrt{F\_L^2 + r\_n^2} \tag{2}$$

*Appl. Sci.* **2018**, *8*, 2634

The acoustic wave has to propagate through the host medium, then cross Fresnel regions and afterwards continue through the host medium. A three-layer configuration has to be considered (Figure 1). Acoustic impedance (*Z*) is defined as the product of the medium density (*ρ*) and the sound propagation velocity (*c*) in it. Therefore, it is necessary to consider the input (*Zin*) and output (*Zout*) acoustic impedance and the transmission pressure coefficients must be calculated (*t*). This coefficient is a clear indicator of the blocking capacity of the elements of the FZP. Hence, *t* is defined as the relation between the transmitted field and the incident field. Density (*ρ*), sound propagation velocity (*c*) and acoustic impedance (*Zmat*) values have been shown in Table 1. Using these values in Equation (3), *Zin* could be obtained [17].

$$Z\_{in} = Z\_{\text{mat}} \frac{Z\_{\text{out}} + jZ\_{\text{mat}} \tan(k\_{\text{mat}}d)}{Z\_{\text{mat}} + jZ\_{\text{out}} \tan(k\_{\text{mat}}d)} \tag{3}$$

where *km* is the wave number, defined as *km* = *ω*/*c*. Considering *ω* = 2*π f*0. Once *Zin* is obtained, reflection coefficient is defined in Equation (4).

$$r\_{in} = \frac{Z\_{in} - Z\_{water}}{Z\_{in} + Z\_{water}} \tag{4}$$

The equation that relates the field balance as a function of the impedance and reflection coefficient of the system is defined in Equation (5) and gives *t* values depending on the material.

$$|t| = \frac{|p\_t^+|}{|p\_{in}^+|} = \sqrt{(1 - |r\_{in}|^2)}\tag{5}$$

**Figure 1.** Transmission diagram of the implemented lenses.

**Table 1.** Density and sound speed values. Acrylonitrile butadiene styrene (ABS), polylactic acid (PLA), and polymethyl methacrylate (PMMA).)

Considering the transmission coefficient values obtained (0.23 for brass, 0.51 for PLA-Air-PLA and more than 0.95 for ABS, PLA and PMMA), it can be affirmed that full implemented MRI compatible material lenses will focus less energy at the *FL* if it is compared to brass FZP or ideal SZP. Therefore, one solution is proposed to obtain the desired impedance contrast. A FZP that includes an air chamber inside the structure has been implemented by using a 3D-printer. Thus, both lenses, full-PLA and air-chamber, have been compared.

#### *2.1. Numerical Model*

The finite elements method (FEM) has been used to obtain a numerical solution of the physical problem. The finite elements method allows us to study the physical phenomena involved in the interaction of waves with FZPs. Therefore, a mathematical model that replicates the conditions of the problem has been implemented. This method also allows us to determine the pressure distribution of the diffracted fields generated by the FZP when there is a piston emitter, causing interference phenomena. From the mesh generated by FEM, a partial differential equation solution is obtained for each node [18]. In this case, acoustic Helmholtz equation is considered (Equation (6)). To solve the Helmholtz equation, standard values of water such as density of the medium (*ρ* =1000 kg/m3) and sound propagation velocity (*c* = 1500 m/s) have been considered. The working frequency of the FZPs is 250 kHz and it can be found by its relation with the angular velocity (*ω*). Finally, *p* corresponds to the acoustic pressure.

$$\nabla \cdot \left( -\frac{1}{\rho\_0} (\nabla p) \right) = \frac{\omega^2 p}{\rho\_0 c^2} \tag{6}$$

If a 3D model is considered, this will require high computational resources. To simplify the model and reduce this computational cost, as shown in previous works [19,20], the geometrical properties of the model are used taking advantage of its axisymmetry. Therefore, the model is simplified by implementing a semi-lens only. A complete solution is obtained by rotating it from its symmetry axis. This procedure achieves a reduction of the degrees of freedom necessary to obtain the results of the numerical simulation and thus significantly diminishing the calculation time.

The boundary conditions defined in the numerical models are explained below as seen in Figure 2. The contours of the model are defined as wave radiation condition boundary to emulate an infinitely large medium and therefore the Sommerfeld condition is satisfied. Acoustic impedance domain definition has been used for all opaque Fresnel regions for each lens. In the case of the SZP lens, the contours are considered infinitely rigid, applying the Neumann condition (the sound velocity in the contour is zero).

**Figure 2.** Scheme of the finite element method (FEM) conditions.

#### **3. Experimental Set-Up**

It is required to validate the results obtained from the theoretical models with other solutions such as numerical models and experimental measurements. In this sense, obtaining experimental results is fundamental to validate the numerical models. A complex measurement and acquisition system is needed to perform the experiments given the technical difficulties to control the underwater devices. The Center for Physics Technologies: Acoustics, Materials and Astrophysics of the Universitat Politècnica de València has a robotized and automated system for high precision ultrasound measurements. The robot is built based on the size of the immersion tank where the tests and experiments are carried out, which contains distilled and degassed water, with dimensions of 0.5 m wide by 0.5 m high by 1 m long. These dimensions suppose that the immersion tank must contain around 200 L of distilled water to be functional, and allow both transducers and devices to be completely submerged, and avoid reflections due to the impedance changes produced by the change medium.

The measurement system is composed by a fixed emitter and a hydrophone coupled to the robotic system. This system obtains reliable and precise results that allow for the evaluation of the acoustic phenomena involved in these types of lenses. A plane immersion piston transducer built by Imasonic with 250 kHz of central working frequency and an active diameter of 32 mm has been used as the emitter. Also, a Precision Acoustics hydrophone, model 1.0 mm Needle Hydrophone is used as the receiver. This hydrophone is capable of measuring high frequencies, even if they have a very weak signal level. The sensitivity of the hydrophone is 850 nV/Pa (−241.4 dB 1V/μPa) with a tolerance of ±3 dB. The frequency response is flat ±2 dB between 3 and 12 MHz and ±4 dB between 200 kHz and 15 MHz. The bandwidth ranges from 5 kHz to 15 MHz. Figure 3 shows an experimental set-up in a measurement. Two types of different configurations are used to generate and amplify the signals. The first one is to use an external function generator connected to a high power amplifier. The second configuration is to use a pulse generator (5077PR of Panametrics) with integrated amplifier. This generator and amplifier allows generating pulses with frequencies between 100 kHz and 20 MHz, a pulse repetition frequency (PRF) from 100 Hz to 5 kHz and a pulse amplitude between 100 and 400 V.

**Figure 3.** Experimental set-up.

All the results shown below are obtained for a working frequency of 250 kHz. For the experimental comparison, three lenses have been implemented, two made of PLA and one of brass. Every lens considered in this work was designed with 11 Fresnel zones and an outer radius of 88.8 mm. The thickness of the brass lens was 1 mm. For manufacturing reasons, the rest of the lenses had a total thickness of 5 mm. Figure 4 shows both PLA and brass lenses. Both PLA lenses are identical, the only difference being an inner air chamber to achieve a higher impedance contrast. As described in the previous section, it is not possible to differentiate them by the naked eye.

**Figure 4.** Implemented lenses, (**a**) PLA and (**b**) brass.

#### **4. Results**

Intensity gain for longitudinal axis cuts and maps have been calculated to compare all the lenses coherently. One parameter, which can be used to evaluate the focusing capacity of a lens is the intensity gain (*G*). The intensity gain is related with the intensity with both the intensity with lens (*I*) and intensity without lens (*I*0), as shown in Equation (7).

$$G(dB) = 10 \cdot \log\_{10}(I/I\_0) \tag{7}$$

Intensity gain values have been calculated from Equation (7). Figure 5, shows the intensity gain for longitudinal cuts on the Z axis for both, numerical and experimental results. It can be seen from Figure 5a, that higher impedance contrast, as in the case of the ideal SZP or brass FZP lens, gives rise to a higher gain levels. As expected, the lower gains are obtained with those materials with impedance contrast values between 1 and 2 and for impedance contrast values lower than 1, the intensity gain increases. ABS, PLA, PLA-Air-PLA, and PMMA polymers, according with the values showed in Table 1, are not able to achieve enough intensity gain as brass FZP. By comparing the experimental results obtained for brass and PLA (see Figure 5b) with numerical ones (see Figure 5a) it can be seen that there is a good agreement. From Figure 5b, it is observed that the air chamber PLA lens has higher intensity gain than full PLA lens. This fact can be explained by the introduction of an air layer. This layer, due to its low acoustic impedance, can block the transmission of the ultrasonic waves approaching its behavior to an ideal SZP. Nevertheless, a focal length displacement of 1.66*λ* is observed in the FZP lens built with air chamber and PLA. In this case, the displacement is due to the new three-layer configuration (PLA-Air-PLA). The resolution of the 3D-printer and the wall width needed to avoid porosities means that there is an interface between the host medium and the air chamber.

**Figure 5.** Intensity gain longitudinal cuts for (**a**) FEM results and (**b**) experimental results.

Figure 6 shows four intensity gain maps, the first three obtained experimentally and the fourth numerically. The experimental ones correspond to PLA, PLA-Air-PLA and brass, while the numerical one has been obtained using an ideal SZP. It has been verified how the results obtained with brass resembles the ideal SZP. This is due to the rigidity of the material. On the other hand, in PLA lens results, a diminishing intensity gain is observed. This intensity gain level can be improved using a PLA-Air-PLA lens. All the lenses are designed with a focal length located at 8.33*λ* for a working frequency of 250 kHz. When the lens is able to block destructive interference, it is possible to locate the focus in *FL*. This occurs in brass and the ideal SZP case. Resulting from the lack of blocking capability, the full PLA FZP could not impede the incident pressure wave transferal generating aberrations in the *FL*.

**Figure 6.** Intensity gain maps for experimental measurements and ideal SZP numerically obtained (FEM).

#### **5. Conclusions**

Non-metallic materials can be used for the construction of acoustic lenses. Three alternative materials, compatible with magnetic resonance, have been proposed instead of brass lenses. It has been possible to verify that the higher the impedance contrast of the materials, the higher the intensity gain levels. The PMMA lens has higher intensity level than ABS and PLA ones, because it has a slightly higer impedance contrast value than ABS or PLA. However, the use of an air chamber inside the PLA lens increases the intensity gain levels, due to the fact that values of impedance contrast less than one means blocking of the waves. PLA is a biocompatible material and is cheaper than PMMA. 3D printers give open field of new lens design MRI compatible. Moreover, since PLA is a biodegradable material, it is a environmental friendly material. This point is important in procedures that generate waste. Nevertheless, the manufacture of PLA lenses require great care because of microporosities that could appear. The appearance of pores can cause water to enter into the lens, drastically reducing the

blocking capacity. In addition, polymers such as PLA, ABS or PMMA are more affordable than metal plates. This will lower the costs in the production of HIFU treatment devices based on acoustic lenses. For this reason, PLA is proposed as an MRI compatible material with great potential for therapeutic applications of ultrasound focusing.

**Author Contributions:** A.U. and C.R. coordinated the theoretical development, participating in the establishment of the theory principles used in this work, as well as in the drafting of the manuscript. D.T.-S. coordinated experimental development. S.C.-I. developed part of the theory used and designed some characterization. E.S.-A. participated in the analysis of the state of art.

**Funding:** This research was funded by spanish Ministerio de Economía y Competitividad (MINECO) TEC2015-70939-R.

**Acknowledgments:** This work has been supported by Spanish MINECO (TEC2015-70939-R).

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

c 2018 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

#### *Article*

### **Theoretical and Numerical Estimation of Vibroacoustic Behavior of Clamped Free Parabolic Tapered Annular Circular Plate with Different Arrangement of Stiffener Patches**

#### **Abhijeet Chatterjee 1,\*, Vinayak Ranjan 2, Mohammad Sikandar Azam 1 and Mohan Rao 3**

Received: 13 September 2018; Accepted: 2 November 2018; Published: 8 December 2018

**Abstract:** This paper compares the vibroacoustic behavior of a tapered annular circular plate having different parabolic varying thickness with different combinations of rectangular and concentric stiffener patches keeping the mass of the plate and the patch constant for a clamped-free boundary condition. Both numerical and analytical methods are used to solve the plate. The finite element method (FEM) is used to determine the vibration characteristic and both Rayleigh integral and FEM is used to determine the acoustic behavior of the plate. It is observed that a Case II plate with parabolic decreasing–increasing thickness variation for a plate with different stiffener patches shows reduction in frequency parameter in comparison to other cases. For acoustic response, the variation of peak sound power level for different combinations of stiffener patches is investigated with different taper ratios. It is investigated that Case II plate with parabolic decreasing–increasing thickness variation for an unloaded tapered plate as well as case II plate with 2 rectangular and 4 concentric stiffeners patches shows the maximum sound power level among all variations. However, it is shown that the Case III plate with parabolically increasing–decreasing thickness variation with different combinations of rectangular and concentric stiffeners patches is least prone to acoustic radiation. Furthermore, it is shown that at low forcing frequency, average radiation efficiency with different combinations of stiffeners patches remains the same, but at higher forcing frequency a higher taper ratio causes higher radiation efficiency, and the radiation peak shifts towards the lower frequency and alters its stiffness as the taper ratio increases. Finally, the design options for peak sound power actuation and reduction for different combinations of stiffener patches with different taper ratios are suggested.

**Keywords:** thick annular circular plate; Rayleigh integral; finite element modeling; rectangular and concentric stiffener patches; taper ratio; thickness variation

#### **1. Introduction**

Tapered annular circular plates with different combinations of rectangular and concentric patches has many engineering applications. They are used in many structural components i.e., building, design, diaphragms and deck plates in launch vehicles, diaphragms of turbines, aircraft and missiles, naval structures, nuclear reactors, optical systems, construction of ships, automobiles and other vehicles, the space shuttle etc. These tapering plates with different combinations of rectangular and concentric patches are found to have greater resistance to bending, buckling and vibration in

comparison to plates of uniform thickness. It is interesting to know that tapered plates with different thickness variation have drawn the attention of most of the researchers in this field. However, tapered plates with different combination of rectangular and concentric patches can alter the dynamic characteristic of structures with a change in stiffness. Hence, for practical design purposes, the vibration and acoustic characteristics of such tapered plates are equally important. In comparison to the present study, several existing works are presented where the researchers have investigated the vibration response [1–9] of circular or annular plates of tapered or uniform thickness. But in terms of acoustic behavior, many researchers have contributed most. Lee and Singh [10] used the thin and thick plate theories to determine the sound radiation from out-of-plane modes of a uniform thickness annular circular plate. Thompson [11] used the Bouwkamp integral to determine the mutual and self-radiation impedances both for annular and elliptical pistons. Levine and Leppington [12] analyzed the sound power generation of a circular plate of uniform thickness using exact integral representation. Rdzanek and Engel [13] determined the acoustic power output of a clamped annular plate using an asymptotic formula. Wodtke and Lamancusa [14] minimized the acoustic power of circular plates of uniform thickness using the damping layer placement. Wanyama [15] studied the acoustic radiation from linearly-varying circular plates. Lee and Singh [16] used the flexural and radial modes of a thick annular plate to determine the self and mutual radiation. Cote et al. [17] studied the vibro acoustic behavior of an unbaffled rotating disk. Jeyraj [18] used an isotropic plate with arbitrarily varying thickness to determine its vibro-acoustic behavior using the finite element method (FEM). Ranjan and Ghosh [19] studied the forced response of a thin plate of uniform thickness with attached discrete dynamic absorbers. Bipin et al. [20] analyzed an isotropic plate with attached discrete patches and point masses with different thickness variation with different taper ratios to determine its vibro acoustic response. Lee and Singh [21] investigated the annular disk acoustic radiation using structural modes through analytical formulations. Rdzanek et al. [22] investigated the sound radiation and sound power of a planar annular membrane for axially-symmetric free vibrations. Doganli [23] determined the sound power radiation from clamped annular plates of uniform thickness. Nakayama et al. [24] investigated the acoustic radiation of a circular plate for a single sound pulse. Hasegawa and Yosioka [25] determined the acoustic radiation force used on the solid elastic sphere. Lee and Singh [26] used a simplified disk brake rotor to investigate the acoustic radiation through a semi-analytical method. Thompson et al. [27,28] analyzed the modal approach for different boundary conditions to calculate the average radiation efficiency of a rectangular plate. Rayleigh [29] determined the sound radiation from flat finite structures. Maidanik [30] analyzed the total radiation resistance for ribbed and simple plates using a simplified asymptotic formulation. Heckl [31] used the wave number domain and Fourier transform to analyses the acoustic power. Williams [32] determined the wave number as a series in ascending power to estimate the sound radiation from a planar source. Keltie and Peng [33] analyzed the sound radiation using the cross- modal coupling from a plane. Snyder and Tanaka [34] demonstrated the importance of cross-modal contributions for a pair of modes through total sound power output using modal radiation efficiency. Martini et al. [35] investigated the structural and elastodynamic analysis of rotary transfer machines by a finite element model. Croccolo et al. [36] determined the lightweight design of modern transfer machine tools using the finite element model. Martini and Troncossi [37] determined the upgrade of an automated line for plastic cap manufacture based on experimental vibration analysis. Pavlovic et al. [38] investigated the modal analysis and stiffness optimization: the case of a tool machine for ceramic tile surface finishing using FEM.

While reviewing the literature, this comes to a discussion at a common point that has inspired the present paper based on a comparison of vibroacoustic behavior of a parabolic tapered annular circular plate with attached rectangular and concentric patches at different positions. The paper is significant for the analysis of the comparison of vibroacoustic behavior of such clamped free tapered plate, which is done by keeping the mass of the plate and patch constant. Therefore, this paper is based on the vibroacoustic analysis of a clamped free parabolic tapered annular circular plate with different attachments of rectangular and concentric stiffener patches for different positions with different taper ratios under time-varying harmonic excitations.

#### **2. Mathematical Modeling and Analysis**

#### *2.1. Plate Free Vibration*

Let us considered a plate with outer radius 'a' and inner radius 'b' as shown in Figure 1. In this paper, the modal analysis is performed to estimate the natural frequency and modes shape of the plate is given by the following equation:

$$\left( [K] - \omega^2 [M] \right) \psi\_{\text{mu}} = 0 \tag{1}$$

where [*M*] is the mass matrix and [*K*] is the stiffness matrix where as *ψmn* is the mode shape and *ω* is the respective natural frequency of the plate in rad/sec. The non-dimensional frequency parameter *λ*2 is given by the following equation:

$$
\lambda^2 = \omega a^2 \sqrt{\frac{\rho h}{D}} \tag{2}
$$

where D, the flexure rigidity = *Eh*3 12(1−*υ*2) , a = outer radius, E = Young's modulus of elasticity, *υ* = Poisson's ratio, h = thickness of the plate and *ρ* = density of plate.

**Figure 1.** Sound radiation investigated for thick annular circular plate in Z direction enclosed in a sphere.

#### *2.2. Analytical and Numerical Formulation for Acoustic Radiation from Tapered Annular Circular Plate*

It is considered that an annular circular plate of inner radius 'b' and outer radius 'a' in flexural vibration is set on flat rigid baffle having infinite extent as reported in Figure 1. Acoustic scattering of the edges of a vibrating structure is neglected in this study. Let P be the sound pressure amplitude, Ss be the surface of the sound source, *q* be the Green methods function in free field, *ls* and *lp* be the position vectors of source and receiver and the surface normal vector at *ls* be *f*; then structure sound radiation can be obtained by the Rayleigh integral [10] as given by Equation (3):

$$P(l\_{\mathbb{P}}) = \int\_{S\_s} \left( \frac{\partial q}{\partial f}(l\_{p\prime}, l\_s) P(l\_{\mathbb{P}}) - \frac{\partial P}{\partial f}(l\_s) q\left(l\_{p\prime}, l\_s\right) \right) ds(l\_s) \tag{3}$$

*Appl. Sci.* **2018**, *8*, 2542

The sound pressure, radiated from non-planar source in far and free field environment based on plane wave approximation can be expressed by Equation (4):

$$P(l\_p) = \frac{\rho\_0 c\_0 B}{4\pi} \int\_{S\_s} \frac{e^{iB|l\_p - l\_s| \cdot \mathcal{U}} (l\_s)}{|l\_p - l\_s|} (1 + \cos \eta) dS \tag{4}$$

Let *ρ*0 be the mass density of air, *c*0 be the speed of sound in air, B be the corresponding acoustic wave number, and . U and . u be the corresponding vibratory velocity amplitude and spatial dependent vibratory velocity amplitude in the *z* direction at *ls*, then from a normal plane [10], the modal sound pressure Pmn for an annular plate with (m, n)th mode is obtained from simplifying Equation (4) with Hankel transform and is expressed by Equations (5) and (6):

$$P\_{\rm mm}(R, a, \beta) = \frac{\rho\_0 c\_0 B e^{i R\_{\rm mm} R\_d}}{2 R\_d} \cos n\beta (-i)^{n+1} A\_n \left[ u(l) \right] (1 + \cos \eta) \tag{5}$$

$$A\_f\left[u(l)\right] = \int\_0^\infty u(l)l\_n(B\_l l)l dl\_l \text{ } lb = B\sin\theta ; \ R\_d = \left|l\_p - l\_s\right|\tag{6}$$

where Xn is Bessel function of order n, (α, β) are the cone and azimuthal angles of the observation positions, respectively, η is the angle between the surface normal vector and the vector from source position to receiver position, and A is the Hankel transform. According to the far field condition, Rd in the denominator is approximated by R whereR=|*lp*| is considered to be radius of the sphere. Consider that on a sphere Sv the observation positions are represented by some points having equal angular increments (Δ*ϕ*, Δ*α*). If 'Δ*ϕ*' represents the small increment in the circumferential direction of the plate, at all of the observation positions, the sound pressures is given by Equations (4)–(6). The modal sound power Smn for the (m, n)th mode [10,16] from the far-field is given by Equation (7):