Some Theoretical and Experimental Extensions Based on the Properties of the Intrinsic Transfer Matrix

The work approaches new theoretical and experimental studies in the elastic characterization of materials, based on the properties of the intrinsic transfer matrix. The term ‘intrinsic transfer matrix’ was firstly introduced by us in order to characterize the system in standing wave case, when the stationary wave is confined inside the sample. An important property of the intrinsic transfer matrix is that at resonance, and in absence of attenuation, the eigenvalues are real. This property underlies a numerical method which permits to find the phase velocity for the longitudinal wave in a sample. This modal approach is a numerical method which takes into account the eigenvalues, which are analytically estimated for simple elastic systems. Such elastic systems are characterized by a simple distribution of eigenmodes, which may be easily highlighted by experiment. The paper generalizes the intrinsic transfer matrix method by including the attenuation and a study of the influence of inhomogeneity. The condition for real eigenvalues in that case shows that the frequencies of eigenmodes are not affected by attenuation. For the influence of inhomogeneity, we consider a case when the sound speed is varying along the layer’s length in the medium of interest, with an accompanying dispersion. The paper also studies the accuracy of the method in estimating the wave velocity and determines an optimal experimental setup in order to reduce the influence of frequency errors.


Introduction
A transfer matrix is an important tool in the study of wave and pulse propagation in finite and infinite homogeneous and inhomogeneous elastic media. The method is widely used in computer simulation or in elastic characterization of different kinds of elastic media [1][2][3].
The main characteristic of matrix methods with respect to other approaches is their simple and compact analytical form and the ease of application in obtaining theoretical results [4]. This has further advantages in the development of numerical methods, where a compact encoding of wave propagation and scattering allows building clear and efficient simulations, which are also easy to test and modify. Computational requirements are also relatively low, and the approach is easy to extend. The downside of matrix methods is that, while they are well adapted to 1D wave propagation, for larger dimensions they require complicated transfer matrices which are harder to work with and interpret in terms of involved wave phenomena. Nevertheless, in many applications and experimental setups, considering a 1D-wave propagation is sufficient, which allows a refined and straightforward approach by a suitable matrix method in order to simulate, explain and characterize the considered system [5].
The elements of the transfer matrix are obtained as a consequence of the boundary conditions and propagation mechanism [6,7]. For the 1D propagation and solid elastic media, boundary conditions imply continuity of the stress and of the wave function at the interface [8][9][10].
Especially in the study of multilayered media, the transfer matrix method is used in sound insulation or transmission loss factor in order to evaluate sound attenuation [11][12][13], in automotive design of the interior of the car [1,4] or the design of multiple connected mufflers [14,15]. A large volume of research refers to sonic crystals' behavior and optimization using the transfer matrix method [16][17][18]. By using the transfer matrix, formalism is possible to model the behavior of elastic media with inhomogeneities. In that case, the transfer matrix approach combines with the finite element method in order to describe homogeneous and non-homogeneous acoustic absorbing materials. The characterized acoustic materials are mainly metamaterials made of multiple layers, where at least one layer consists of a non-homogeneous material. The equivalent transfer matrix of the non-homogeneous material is determined and, by using the equivalent transfer matrix of the nonhomogeneous material coupled analytically in a series with other transfer matrices, complex multilayer systems can be modeled easily and quickly in configurations wherein the use of finite element calculation alone will be more expensive and time consuming [19]. Other methods have been developed to numerically simulate waves in complex materials, some of which take advantage of matrix formulations [20][21][22][23]. Especially in the case of polycrystalline materials, velocity surfaces [24] or slowness surfaces [25,26] are used to characterize local anisotropy. To take into account inhomogeneity and anisotropy is important in sound wave propagation [27,28] as it gives a clearer view on wave propagation and scattering and also on the emerging sound dispersion in composite or multilayered structures.
Intrinsic transfer matrix represents a special kind of transfer matrix written for amplitudes of the Fourier components of the waves confined in an embedded elastic system. Combining the properties of the intrinsic transfer matrix with a corresponding modal analysis, we can determine some elastic constants of the system. This application of the method is presented in Section 2. By imposing the condition that the eigenvalues of the intrinsic transfer matrix are real at resonance, we can generate another form of the resonance condition for the considered elastic system. This extension is presented in Section 3.1. The resonance condition may be generalized if we consider the attenuation of the amplitudes of the involved waves. The obtained results in case of attenuation were compared with those obtained by using a simulation based on the model of the coupled oscillators (this method of simulation is very often used for the study of multilayered materials). This is presented in Section 3.2. By using computer simulation, the intrinsic transfer matrix was also applied to study the behavior of a multilayered medium with inhomogeneities, by considering the case in which one layer consists of a non-homogeneous material. This is presented in Section 3.3. One possible effect of inhomogeneity is the appearance of dispersion, i.e., frequency-dependent sound velocity. This is applied in the case of multilayered media, and also for polycrystalline materials composed of anisotropic grains such as metals, or for metamaterials in a larger sense. The simulations are presented in Section 3.3. In Section 3.4, we approach an optimization procedure for the design of a multilayer medium, where the purpose is to place the sample of interest such as to minimize the sound velocity errors due to frequency determination. The case of a ternary elastic system is being analyzed, but the study may be generalized to other complex multilayered media.

The Idea of an Intrinsic Transfer Matrix
If we consider a single elastic homogeneous layer with width l placed between two semi-infinite elastic media, the transfer matrix in the Fourier space that describes the propagation of a progressive and a regressive wave with respective amplitudes A, B is obtained as: In Formula (1) Z in and Z out are the elastic impedances of the two semi-infinite media, Z is the elastic impedance of the medium of interest, ω is the angular frequency and c Materials 2022, 15, 519 3 of 15 the phase velocity of the wave. The transfer matrix in this simple case is a product of a propagation matrix P and two discontinuity matrices D, with general expressions: with k = ω/c as the wavenumber and z = Z 1 /Z 2 as the relative impedance of media 1 and 2.
In the case of a standing wave, when the wave is confined only within the sample, the Fourier components of the waves inside the sample are determined only by the propagation matrix P(ω, l) in (2), which has two eigenvalues given by: As it is known, if the standing wave consists of a superposition of its own eigenmodes, for which we have ω n = nπc l , n = 1, 2 . . . . . ., it is obvious from (3) that in this special case the eigenvalues become real. Experiments confirm that we can extend the reasoning to complex embedded elastic systems and consider that for eigenfrequencies the eigenvalues become real. This important behavior of the eigenvalues can be used experimentally to find elastic constants of materials. The implied procedure is, to calculate the intrinsic transfer matrix and its eigenvalues, to experimentally find the frequencies of eigenmodes and to evaluate by a numerical analysis for which values of the wave velocity the eigenvalues are real.
The intrinsic transfer matrix method can become very useful in finding the longitudinal wave velocity in solid elastic media, especially for small samples that are not suitable to be measured by classical resonance methods. Such a sample is embedded in an elastic system with free ends containing the sample of interest and gauge materials. The entire system is excited using an impact hammer that generates an approximate Gaussian elastic pulse so that a standing wave is generated in the embedded system. The standing wave in the system is highlighted by a noncontact technique using a Doppler vibrometer, in our case an Ometron VQ-400A vibrometer. The analog signal proportional to the vibration velocity of the surface at the end of the sample is acquired with an acquisition board and spectrally processed using an analysis software, in our case LabView (or directly using an FFT analyzer). A simple FFT power spectrum allows for the determination of the frequencies of the system's eigenmodes. Knowing the dimensions of the embedded components in the elastic system determined by length measurements, the mass density determined by weighting, and asserting the condition that the eigenvalues are real, we can estimate by a numerical analysis the wave velocity in the sample of interest.
As an example of the application of the intrinsic matrix method, in Figure 1 we illustrated the Fourier spectrum of the eigenmodes for a system consisting of two aluminum rods as gauge materials and a sample of interest of wood. The two spectra (red and black) refer to the same wood species but the two samples were cut differently, one along the fiber and the other one radially.
In Figure 2, we plotted the theoretical values of eigenfrequencies obtained from the condition that the eigenvalues are real for the ternary system illustrated in Figure 1, which consists of two aluminum rods as gauge materials (l 1 = 150 mm, l 3 = 300 mm) and the radially cut wood sample (l 2 = 40.4 mm). In Figure 2, we plotted the theoretical values of eigenfrequencies obtained from the condition that the eigenvalues are real for the ternary system illustrated in Figure 1, which consists of two aluminum rods as gauge materials ( = 150 , = 300 ) and the radially cut wood sample ( = 40.4 ). There is a good concordance between theory and experiment, especially for the first eigenmode of the embedded system, and as a result this recommends it in practical applications for determining the velocity of longitudinal elastic waves in small samples unsuitable for classical resonance measurements. For applications, it is important for the design of the experimental setup to be correlated to the possibility of the measuring system and the desired accuracy.

A New Form of Resonance Condition
A multilayer medium with layers 1, 2, ..., n has the intrinsic transfer matrix = , … , . The propagation of a progressive and a regressive wave with amplitudes = ( ), = ( ) through the medium is generally described by: where , are complex-valued coefficients depending on and the layer properties and * , * are their complex conjugates. If the medium has free ends, then at resonance = , = and (4) becomes an eigenvalue equation: * * In Figure 2, we plotted the theoretical values of eigenfrequencies obtained from the condition that the eigenvalues are real for the ternary system illustrated in Figure 1, which consists of two aluminum rods as gauge materials ( = 150 , = 300 ) and the radially cut wood sample ( = 40.4 ). There is a good concordance between theory and experiment, especially for the first eigenmode of the embedded system, and as a result this recommends it in practical applications for determining the velocity of longitudinal elastic waves in small samples unsuitable for classical resonance measurements. For applications, it is important for the design of the experimental setup to be correlated to the possibility of the measuring system and the desired accuracy.

A New Form of Resonance Condition
A multilayer medium with layers 1, 2, ..., n has the intrinsic transfer matrix = , … , . The propagation of a progressive and a regressive wave with amplitudes = ( ), = ( ) through the medium is generally described by: where , are complex-valued coefficients depending on and the layer properties and * , * are their complex conjugates. If the medium has free ends, then at resonance = , = and (4) becomes an eigenvalue equation: There is a good concordance between theory and experiment, especially for the first eigenmode of the embedded system, and as a result this recommends it in practical applications for determining the velocity of longitudinal elastic waves in small samples unsuitable for classical resonance measurements. For applications, it is important for the design of the experimental setup to be correlated to the possibility of the measuring system and the desired accuracy.

A New Form of Resonance Condition
A multilayer medium with layers 1, 2, . . . , n has the intrinsic transfer matrix T = P N D N,N−1 . . . P 2 D 1,2 P 1 . The propagation of a progressive and a regressive wave with amplitudes A = A(ω), B = B(ω) through the medium is generally described by: where a, b are complex-valued coefficients depending on ω and the layer properties and a * , b * are their complex conjugates. If the medium has free ends, then at resonance A in = B in , A out = B out and (4) becomes an eigenvalue equation: Then, at resonance λ is real-valued and from (5) the resonance condition is: Generally, (6) is a nonlinear equation in ω and it must be solved numerically in order to obtain resonance frequencies ω res , or the layers' parameters if ω res is known.

Intrinsic Transfer Matrix in Case of Attenuation
To study the influence of attenuation it is necessary to introduce attenuation coefficients of the implied materials. By considering a ternary system, where β 1 , β 2 , β 3 are the attenuation coefficients of amplitude and c 1 , c 2 , c 3 are the corresponding wave velocities, the intrinsic transfer matrix in the presence of attenuation will be: The condition for real eigenvalues shows that the frequencies of eigenmodes are not affected by attenuation. The equations show that the most probable condition to have real eigenvalues is: and consequently: where n and m are integers. A numerical simulation by the coupled oscillators method confirms the independence of eigenfrequencies of the attenuation coefficients. For this, a ternary brass-aluminum-brass system with free ends was modeled as an ensemble of 5000 elements with mass connected by springs and a short impulse was applied to one end. A Fourier transform was applied to the signal at the ends. Figure 3 shows the eigenfrequencies without and with attenuation for the ternary system.
Then, at resonance is real-valued and from (5) the resonance condition is: Generally, (6) is a nonlinear equation in and it must be solved numerically in order to obtain resonance frequencies , or the layers' parameters if is known.

Intrinsic Transfer Matrix in Case of Attenuation
To study the influence of attenuation it is necessary to introduce attenuation coefficients of the implied materials. By considering a ternary system, where , , are the attenuation coefficients of amplitude and , , are the corresponding wave velocities, the intrinsic transfer matrix in the presence of attenuation will be: The condition for real eigenvalues shows that the frequencies of eigenmodes are not affected by attenuation. The equations show that the most probable condition to have real eigenvalues is: and consequently: where and are integers. A numerical simulation by the coupled oscillators method confirms the independence of eigenfrequencies of the attenuation coefficients. For this, a ternary brass-aluminum-brass system with free ends was modeled as an ensemble of 5000 elements with mass connected by springs and a short impulse was applied to one end. A Fourier transform was applied to the signal at the ends. Figure 3 shows the eigenfrequencies without and with attenuation for the ternary system.

Influence of the Inhomogeneity Studied by Intrinsic Transfer Matrix
A special application refers to an inhomogeneous medium. Based on the intrinsic transfer matrix for a ternary system, we considered the sample of interest with a randomly varying sound speed c 2 along the middle layer, which was divided into 100 uniform slices. We numerically generated profiles for c 2 by two methods: as random values with a Gaussian distribution, or as a random walk along the layer ( Figure 4). Figure 3. The eigenfrequencies of a brass-aluminum-brass system without (above) and with (below) attenuation.

Influence of the Inhomogeneity Studied by Intrinsic Transfer Matrix
A special application refers to an inhomogeneous medium. Based on the intrinsic transfer matrix for a ternary system, we considered the sample of interest with a randomly varying sound speed along the middle layer, which was divided into 100 uniform slices. We numerically generated profiles for by two methods: as random values with a Gaussian distribution, or as a random walk along the layer (Figure 4). The first method may model fractures and inclusions into the medium, while the second method generates continuous profiles and can model the influence of time-varying manufacturing processes on the material properties. Each generated profile ( ) was normalized such that its mean is 〈 〉 = 5018 / . Profiles with different value spreads, as estimated by the mean square root error ( ) = ( ), were obtained and the first eigenfrequency was computed from the intrinsic transfer matrix condition ( ) = 0. Figure 5 shows the influence of the inhomogeneity on for different ( ).  The first method may model fractures and inclusions into the medium, while the second method generates continuous profiles and can model the influence of time-varying manufacturing processes on the material properties. Each generated profile c 2 (x) was normalized such that its mean is c 2 = 5018 m/s. Profiles with different value spreads, as estimated by the mean square root error MSE(c 2 ) = σ(c 2 ), were obtained and the first eigenfrequency f 1 was computed from the intrinsic transfer matrix condition Im(λ) = 0. Figure 5 shows the influence of the inhomogeneity on f 1 for different MSE(c 2 ).  Figure 3. The eigenfrequencies of a brass-aluminum-brass system without (above) and with (below) attenuation.

Influence of the Inhomogeneity Studied by Intrinsic Transfer Matrix
A special application refers to an inhomogeneous medium. Based on the intrinsic transfer matrix for a ternary system, we considered the sample of interest with a randomly varying sound speed along the middle layer, which was divided into 100 uniform slices. We numerically generated profiles for by two methods: as random values with a Gaussian distribution, or as a random walk along the layer (Figure 4). The first method may model fractures and inclusions into the medium, while the second method generates continuous profiles and can model the influence of time-varying manufacturing processes on the material properties. Each generated profile ( ) was normalized such that its mean is 〈 〉 = 5018 / . Profiles with different value spreads, as estimated by the mean square root error ( ) = ( ), were obtained and the first eigenfrequency was computed from the intrinsic transfer matrix condition ( ) = 0. Figure 5 shows the influence of the inhomogeneity on for different ( ).  It can be seen that there is a similar variation of the first resonance frequency for the two methods for MSE(c 2 ) smaller than about 1000 m/s, with f 1 decreasing by about 1 Hz for MSE(c 2 ) = 500 m/s, or about 10% of c 2 . Thus, inhomogeneities in the sound speed have a small effect on the resonance frequency of the medium.
One possible effect of inhomogeneity is the appearance of dispersion, i.e., frequencydependent sound velocity. This is the case for multilayer media, and also for polycrystalline materials composed of anisotropic grains, such as metals [29][30][31], or for metamaterials in a larger sense. In order to follow such an effect, we completed a series of simulations where layer 2 of the ternary brass-aluminum-brass medium is composed of such anisotropic grains and acquires a frequency-dependent sound velocity c 2 = c 2 ( f ). Layers 1 and 3 were considered homogeneous and isotropic. To induce graininess, layer 2 was decomposed in N F = 100 longitudinal "fibers" (propagation directions along its length), each being composed of an equal number N G of grains with random lengths l ij having an exponential distribution [31,32]. We considered sound propagation only along the fibers, independently on each fiber. The anisotropy of grains is determined by a velocity surface v = v 0 ·F v (θ, ϕ) giving the sound velocity along direction (θ, ϕ) in spherical coordinates. Here, v 0 = 5018 m/s is the normal sound velocity in layer 2 (when considered homogneous) and F v (θ, ϕ) is a direction function; in the simplest case, this is a Fourier series in θ, ϕ of the form: but it can have other expressions too. For example, C v is a constant determined by the condition: Different other conditions can be imposed to the function F v (θ, ϕ), e.g., F v (π − θ, π + ϕ) = F v (θ, ϕ), which ensures equal sound velocities in opposite directions. However, for metamaterials, this condition may not hold. In the simulation, for a given material in layer 2, a specific expression for F v (θ, ϕ) was applied to all the grains and fibers; for each grain, the direction of propagation was chosen at random, as (θ, ϕ) = (acos(1 − 2r 1 ), 2πr 2 ), where r 1 , r 2 are uniformly distributed random numbers between 0 and 1. Four velocity surfaces were used in (Figure 6), given by: F v (θ, ϕ) = 1 + (0.05 cos(4θ) − 0.025 cos(8θ))· cos(4ϕ); (9c) F v (θ, ϕ) = 1 − (π − θ)e θ−π ·(1 + (0.05 cos(4θ) − 0.025 cos(8θ))· cos(4ϕ)).
Surface (a) is suitable for modeling metals with cubic crystals [30]. Surface (d) in particular is highly asymmetric and could be achieved in a metamaterial.
At the input of layer 2 a pair of waves was applied with A in = 1, B in = 0 for a range of frequencies. The outputs for all fibers were averaged. An effective sound velocity for layer 2 c 2e f f was computed from the average output: Examples of effective velocities for multigrain materials can be seen in Figure 7. Additionally, resonance frequencies for the whole ternary medium were obtained and their variation with the grain number N G was studied (Figure 8). Examples of effective velocities for multigrain materials can be seen in Figure 7. Additionally, resonance frequencies for the whole ternary medium were obtained and their variation with the grain number was studied (Figure 8).  Figure 6, for different grain number along fibers. The case = 1 is for the homogeneous (monocrystalline) layer 2. The legend applies to all graphs. (a-d) correspond to Equations (9a)-(9d) respectively. Figure 6. The four velocity surfaces F v (θ, ϕ) in 3D that were used to describe anisotropy in sound wave propagation. On the left, the surface colors represent the distance from the center, from small (blue) to large (yellow). On the right, the meshes represent the anisotropic (dark red) and isotropic (green) surface velocities. (a-d) correspond to Equations (9a)-(9d) respectively. Examples of effective velocities for multigrain materials can be seen in Figure 7. Additionally, resonance frequencies for the whole ternary medium were obtained and their variation with the grain number was studied (Figure 8).   Figure 6, for different grain number N G along fibers. The case N G = 1 is for the homogeneous (monocrystalline) layer 2. The legend applies to all graphs. (a-d) correspond to Equations (9a)-(9d) respectively.

Figure 8.
Graphs of relative variations of resonance frequencies for the ternary medium as a function of for the velocity surfaces in Figure 6 and velocity profiles in Figure 7. Here, and are resonance frequencies with homogeneous ( = 1) and, respectively, grainy layer 2. The first resonance for = 1 is at 2068.36 Hz, while the next resonances are approximate multiples of it. The resonances No. 1,2,5,10,100 and 1000 were considered. (a-d) correspond to Equations (9a)-(9d) respectively.
Anisotropy induces a frequency-dependent sound velocity (Figure 7), and this effect is stronger for a small number of grains (large grain size), while for small grains the effective velocity is almost constant. The average effective velocity is close to that of the homogeneous layer, except for the velocity surface (d), which induces a strong anisotropy.
The considered dispersion has a minor effect on the resonances of the ternary medium ( Figure 8). This is partly due to the fact that layer 2 is much shorter than layers 1 and 3, but also to the fact that the considered surface velocities are close to the isotropic one. The relative variation of the resonance frequency is quite small (around 10 , or 0.2 Hz for the first resonance), and it is smaller for large resonance frequencies. A relatively large change in the resonance frequency occurs for the velocity surface (d) (close to 10 , or 2 Hz for the first resonance). The case where layers 1 and 3 are also anisotropic was not considered.

A Proposed Optimization Alghorithm Based on the Intrinsic Transfer Matrix
An optimization procedure for the design of a multilayer medium is described, where the purpose is to place the measured sample such as to minimize the sound velocity errors due to frequency determination. Typically, this is applied to a ternary (or more complex) multilayer medium consisting of a known material 1 and the sample material 2. The density is known for both; the sound velocity is known for material 1 only. The first N resonance frequencies obey a resonance condition (6) of the general form: ( , ) = 0, = 1,2, … , For example, for a binary medium with material order 1-2 and a ternary medium 1-2-1, with layer lengths , = / and = / , (10) we have the formulas:  Figure 6 and velocity profiles in Figure 7. Here, f i0 and f i are resonance frequencies with homogeneous (N G = 1) and, respectively, grainy layer 2. Anisotropy induces a frequency-dependent sound velocity (Figure 7), and this effect is stronger for a small number of grains (large grain size), while for small grains the effective velocity is almost constant. The average effective velocity is close to that of the homogeneous layer, except for the velocity surface (d), which induces a strong anisotropy.
The considered dispersion has a minor effect on the resonances of the ternary medium ( Figure 8). This is partly due to the fact that layer 2 is much shorter than layers 1 and 3, but also to the fact that the considered surface velocities are close to the isotropic one. The relative variation of the resonance frequency is quite small (around 10 −4 , or 0.2 Hz for the first resonance), and it is smaller for large resonance frequencies. A relatively large change in the resonance frequency occurs for the velocity surface (d) (close to 10 −3 , or 2 Hz for the first resonance). The case where layers 1 and 3 are also anisotropic was not considered.

A Proposed Optimization Alghorithm Based on the Intrinsic Transfer Matrix
An optimization procedure for the design of a multilayer medium is described, where the purpose is to place the measured sample such as to minimize the sound velocity errors due to frequency determination. Typically, this is applied to a ternary (or more complex) multilayer medium consisting of a known material 1 and the sample material 2. The density is known for both; the sound velocity is known for material 1 only. The first N resonance frequencies ω i obey a resonance condition (6) of the general form: For example, for a binary medium with material order 1-2 and a ternary medium 1-2-1, with layer lengths l i , k i = ω/c i and z = Z 1 /Z 2 , (10) we have the formulas: C 2 : (z + 1) sin(k 1 l 1 + k 2 l 2 ) + (z − 1) sin(k 1 l 1 − k 2 l 2 ) = 0 (11) Typically, the experimental resonances ω exp, i = 2π f exp,i , i = 1, . . . , N will have an error due to the discrete Fourier transform used in their determination. To further estimate c 2 , theoretical resonances, ω th, i (c 2 ) are obtained from (10) which then are matched to ω exp, i by minimizing the convex error: Thus, the estimated c 2 will have an error due to resonance frequency errors. The error of c 2 also depends on the position of sample material 2 in the multilayer medium, which allows one to minimize this error.
For simplicity, let us consider a ternary medium where the sample 2 must be placed at optimal position l 1 . Furthermore, we assume ω exp, i are uniformly distributed in intervals ω exp, i ∈ [ω th, i (c 2 ) − ∆ω; ω th, i (c 2 ) + ∆ω] with a fixed ∆ω. Then, the standard deviation of ω exp, i is σ ω exp, i = ∆ω/ √ 3 and assuming this is small compared to ω exp, i , the standard deviation of c 2 = c 2 ω exp,1 , ω exp,2 , . . . , ω exp,N is estimated by error propagation: where v i is the Euclidean norm of a vector v. Thus l 1 is given by the minimum of (15) with c 2 from (14): Steps (14) and (16) can be repeated to further improve l 1 ; usually one pass and going through (14) again is enough to minimize (15).
Generally, estimating numerical derivatives and performing, e.g., a least squares minimization in (14) and (16) may be time and computation intensive. Obtaining an analytical formula of ∂c 2 /∂ω i from (10) and (16) speeds up computation. Differentiating (10) with respect to c 2 twice at a fixed i, with ω i = ω i (c 2 ) one obtains: From (14), after labeling ω i (c 2 ) = ω th, i (c 2 ), it follows: To determine ∂c 2 /∂ω i in (15), we apply a small perturbation δω exp, i in ω exp, i with a corresponding perturbation δc 2 of c 2 . Then (19) becomes: By retaining terms up to first-order perturbations, we obtain: For fixed k, take in the right-hand side δω exp, i = δω exp, k if i = k; else δω exp, i = 0. Then: where c 2 is determined from (14) and ω i = ω th,i (c 2 ) from (10). The derivatives with respect to c 2 are given by (17) and (18). In particular, if one solves (16) with just the theoretical resonance spectrum, i.e., ω exp, i = ω th,i then: A numerical analysis was performed to illustrate the approach. A ternary medium made of brass (layers 1 and 3) and aluminum (layer 2) was considered, with layer lengths l 1 = 544 mm, l 2 = 18.16 mm, l 3 = 251.5 mm and densities ρ 1 = ρ 3 = 8315 kg/m 3 , ρ 2 = 2713 kg/m 3 . The sound speed in brass is considered known c 1 = c 3 = 3373 m/s while c 2 is to be determined from the resonance spectrum; typically, a value c 2 = 5018 m/s was obtained and is used in the analysis below. This medium has the first 4 resonances at 2104.03 Hz, 4205.93 Hz, 6303.47 Hz and 8394.18 Hz. Next, the frequency error was considered fixed ∆ f = 10 Hz, and the optimum distance l 1 that minimizes (15) was obtained from (16). When modifying l 1 in the analysis below the total length of the medium was kept constant.
The norm (22) as a function of l 1 was of interest (Figure 9). If the sample is placed at the ends of the medium, c 2 has very large errors, while inside it lies close to the measured value c 2 = 5018 m/s. The positions of minimum error for c 2 lie at about l 1 = 10 cm and l 1 = 70 cm, but also positions at about l 1 = 30 cm and l 1 = 50 cm trigger similar errors. The average of c 2 obtained for different resonance spectra approximates well the known value ( Figure 10), except at the ends, where σ(c 2 ) is very large.
where is determined from (14) and = , ( ) from (10). The derivatives with respect to are given by (17) and (18). In particular, if one solves (16) with just the theoretical resonance spectrum, i.e., , = , then: A numerical analysis was performed to illustrate the approach. A ternary medium made of brass (layers 1 and 3) and aluminum (layer 2) was considered, with layer lengths = 544 mm , = 18.16 mm , = 251.5 mm and densities = = 8315 kg/m 3 , = 2713 kg/m 3 . The sound speed in brass is considered known = = 3373 m/s while is to be determined from the resonance spectrum; typically, a value = 5018 m/s was obtained and is used in the analysis below. This medium has the first 4 resonances at 2104.03 Hz, 4205.93 Hz, 6303.47 Hz and 8394.18 Hz. Next, the frequency error was considered fixed ∆ = 10 Hz, and the optimum distance that minimizes (15) was obtained from (16). When modifying in the analysis below the total length of the medium was kept constant.
The norm (22) as a function of was of interest (Figure 9). If the sample is placed at the ends of the medium, has very large errors, while inside it lies close to the measured value = 5018 m/s. The positions of minimum error for lie at about = 10 cm and = 70 cm , but also positions at about = 30 cm and = 50 cm trigger similar errors. The average of obtained for different resonance spectra approximates well the known value ( Figure 10), except at the ends, where ( ) is very large.   (22), not to scale. The total length of the medium is constant. Figure 10. At different positions , 20 sets of resonance spectra were randomly generated with the same ∆ and for each set the optimum was obtained from (14) (the average of the 20 values of is plotted to the left); the standard deviation of (right) was computed numerically (blue dots) and theoretically from (15) (red curve).
We obtained the dependence of (22) with respect to and (Figure 11), and the number of resonances ( Figure 12). The norm (22) decreases with decreasing velocities , increasing sample lengths , and with increasing . For a fixed , the norm has local minima, most of which have about the same minimum value. The absolute minimum is reached at the outermost local minima, and this minimum moves closer to the ends when increases, which could make it difficult to place the sample there. On the other hand, placing the sample at the very end would give the greatest errors for . So, a binary medium is not very appropriate for sound velocity measurements.   We obtained the dependence of (22) with respect to c 2 and l 2 (Figure 11), and the number of resonances N ( Figure 12). The norm (22) decreases with decreasing velocities c 2 , increasing sample lengths l 2 , and with increasing N. For a fixed N, the norm has N local minima, most of which have about the same minimum value. The absolute minimum is reached at the outermost local minima, and this minimum moves closer to the ends when N increases, which could make it difficult to place the sample there. On the other hand, placing the sample at the very end would give the greatest errors for c 2 . So, a binary medium is not very appropriate for sound velocity measurements.  Figure 9. Optimum values of (blue dots) obtained from (14) for random sets of N = 4 resonances with frequency error ∆ = 10 Hz, for different positions of the sample. The red line is from (22), not to scale. The total length of the medium is constant. Figure 10. At different positions , 20 sets of resonance spectra were randomly generated with the same ∆ and for each set the optimum was obtained from (14) (the average of the 20 values of is plotted to the left); the standard deviation of (right) was computed numerically (blue dots) and theoretically from (15) (red curve).
We obtained the dependence of (22) with respect to and (Figure 11), and the number of resonances ( Figure 12). The norm (22) decreases with decreasing velocities , increasing sample lengths , and with increasing . For a fixed , the norm has local minima, most of which have about the same minimum value. The absolute minimum is reached at the outermost local minima, and this minimum moves closer to the ends when increases, which could make it difficult to place the sample there. On the other hand, placing the sample at the very end would give the greatest errors for . So, a binary medium is not very appropriate for sound velocity measurements.     (22), not to scale. The total length of the medium is constant. Figure 10. At different positions , 20 sets of resonance spectra were randomly generated with the same ∆ and for each set the optimum was obtained from (14) (the average of the 20 values of is plotted to the left); the standard deviation of (right) was computed numerically (blue dots) and theoretically from (15) (red curve).
We obtained the dependence of (22) with respect to and (Figure 11), and the number of resonances ( Figure 12). The norm (22) decreases with decreasing velocities , increasing sample lengths , and with increasing . For a fixed , the norm has local minima, most of which have about the same minimum value. The absolute minimum is reached at the outermost local minima, and this minimum moves closer to the ends when increases, which could make it difficult to place the sample there. On the other hand, placing the sample at the very end would give the greatest errors for . So, a binary medium is not very appropriate for sound velocity measurements.   A 5-layer medium was also investigated, with layer ordering 1-2-1-2-1, with l 2 = l 4 = 18.16 mm and l 1 + l 3 + l 5 = 795.5 mm fixed. The norm (22) was minimized and optimal values for l 1 and l 3 were obtained ( Figure 13). Thus, the optimization algorithm can be applied even for complex experimental multilayer media.
A 5-layer medium was also investigated, with layer ordering 1-2-1-2-1, with = = 18.16 mm and + + = 795.5 mm fixed. The norm (22) was minimized and optimal values for l1 and l3 were obtained ( Figure 13). Thus, the optimization algorithm can be applied even for complex experimental multilayer media. Figure 13. Map of the norm (22) in ( ) as function of and for a five-layer medium with two identical samples, with = 4 resonances. The total length was fixed. The high values are trimmed at the edges. The absolute minimum is marked with a red "plus" sign.

Conclusions
This work presents a method of determining the longitudinal wave velocity at a normal incidence in a sample which cannot be measured by the classical resonance method. The method suggests embedding the sample in an elastic system and then determining the frequencies of the eigenmodes of the elastic system. Taking into account that for a standing wave in the system, for such a wave, only the intrinsic matrix can be considered in a wave transfer; for simple systems it can be calculated, and we can find an analytical expression for it. Generalizing for the elastic systems the behavior of the eigenvalues of the intrinsic transfer matrix, which for eigenfrequencies become real, we can determine by a numerical analysis the wave velocity in the sample of interest and then the corresponding elastic constants. From the experimental point of view, a ternary system is preferred because such a system preserves much better the longitudinal plan wave, a special case for which the transfer matrix has a simple analytical form. An example of the application of the method is illustrated for a ternary aluminum-woodaluminum system, which shows the sensitivity of the method.
The paper presents also some experimental and computational studies that refer to the possibilities and accuracy of the method. The method can be generalized so as to take into account the attenuation. In that case, the condition for real eigenvalues shows that the frequencies of eigenmodes are not affected by attenuation.
The case of an inhomogeneous medium is also analyzed by considering two kinds of inhomogeneities. The first can model fractures and inclusions into the medium, while the second generates continuous profiles and it can model the influence of time-varying manufacturing processes on the acoustical properties. Inhomogeneity may also induce dispersion in the medium, and this was analyzed too; a small influence on resonance frequencies was detected.
The paper also proposes an optimization procedure based on the properties of the intrinsic transfer matrix that makes it possible to place the measured sample in a multilayer medium in such a way as to minimize the errors for the measured sound velocity due to frequency errors. It is very clear that, from an experimental point of view, good estimations of the wave velocity can be obtained only by a very precise Fourier analysis; hence, a good frequency resolution will lead to a good estimation of the phase velocity of Figure 13. Map of the norm (22) in σ(c 2 ) as function of l 1 and l 3 for a five-layer medium with two identical samples, with N = 4 resonances. The total length was fixed. The high values are trimmed at the edges. The absolute minimum is marked with a red "plus" sign.

Conclusions
This work presents a method of determining the longitudinal wave velocity at a normal incidence in a sample which cannot be measured by the classical resonance method. The method suggests embedding the sample in an elastic system and then determining the frequencies of the eigenmodes of the elastic system. Taking into account that for a standing wave in the system, for such a wave, only the intrinsic matrix can be considered in a wave transfer; for simple systems it can be calculated, and we can find an analytical expression for it. Generalizing for the elastic systems the behavior of the eigenvalues of the intrinsic transfer matrix, which for eigenfrequencies become real, we can determine by a numerical analysis the wave velocity in the sample of interest and then the corresponding elastic constants. From the experimental point of view, a ternary system is preferred because such a system preserves much better the longitudinal plan wave, a special case for which the transfer matrix has a simple analytical form. An example of the application of the method is illustrated for a ternary aluminum-wood-aluminum system, which shows the sensitivity of the method.
The paper presents also some experimental and computational studies that refer to the possibilities and accuracy of the method. The method can be generalized so as to take into account the attenuation. In that case, the condition for real eigenvalues shows that the frequencies of eigenmodes are not affected by attenuation.
The case of an inhomogeneous medium is also analyzed by considering two kinds of inhomogeneities. The first can model fractures and inclusions into the medium, while the second generates continuous profiles and it can model the influence of time-varying manufacturing processes on the acoustical properties. Inhomogeneity may also induce dispersion in the medium, and this was analyzed too; a small influence on resonance frequencies was detected.
The paper also proposes an optimization procedure based on the properties of the intrinsic transfer matrix that makes it possible to place the measured sample in a multilayer medium in such a way as to minimize the errors for the measured sound velocity due to frequency errors. It is very clear that, from an experimental point of view, good estimations of the wave velocity can be obtained only by a very precise Fourier analysis; hence, a good frequency resolution will lead to a good estimation of the phase velocity of the wave. Nevertheless, precision can be increased by choosing an optimal position of the sample inside the experimental setup with respect to frequency error propagation to other measures of interest that are determined. Author Contributions: N.C. did conceptualization and methodology, experimental measurements, part of theoretical description and applications, and part of writing; M.-I.P. did simulations, part of theoretical description and applications, and part of writing; H.S.A.P. did some simulations and experiments refeering to data acquisition and Doppler interferometry. All authors have read and agreed to the published version of the manuscript.