Shape-preserving average frequency response curves using rational polynomials: A case study on human stapes vibration measurements

.


Introduction
To investigate the behavior of the middle ear under incident sound, the velocity response of the ossicles is often measured.Using methods such as Laser Doppler Vibrometry (LDV), the vibration amplitude and phase due to incident sound can be recorded.As responses differ between individuals, average response curves must be computed from the LDV data to obtain a generic middle ear response (ME).Ideally, these average response curves can then be used to detect pathological ears.Moreover, computer simulations of the ME rely on these data to validate the accuracy of the models [3].However, there is no well-defined method for considering the variations between individuals to compute such average response curves.
Response curves of the ME of individuals show similar characteristics.For example, we show in Fig. 1 human stapes footplate response curves measured using LDV by Niklasson et al. [8].Fig. 1A shows the velocity amplitude of each of the six footplates measured in Ref. [8].Fig. 1B shows the corresponding phase of the motion.It can be observed that each curve, shown in gray, responds similarly to sound.At low frequencies, the velocity amplitude increases gradually and reaches a resonance around 1 kHz.The velocity remains constant up to frequencies of 3 kHz, but then sharply increases again.Above 5 kHz, all six samples show a rapid decrease in velocity amplitude.The velocity phase shows a similar behavior (Fig. 1B).In all samples, the input sound pressure wave leads by 0.25 cycles at low frequencies.A slight drop in the phase is seen when the first resonance around 1 kHz is reached.When the second resonance between 3 and 4 kHz is reached, the individual curves show a sharp change in the phase of about − 0.5 cycles.The experimental procedure used in Ref. [8] had a partially intact ear canal.Therefore, two resonances are present in the data.The first resonance around 1 kHz is that of the hinging motion of the middle ear ossicles [2,5,8], and the second resonance around 3 kHz is that of the ear canal [6,8,9].Fig. 1 illustrates why the traditional average of the velocity curves is a poor predictor of the generic response of the footplate.The solid black line represents the average curve of the experimental data in Fig. 1.Due to the differences in the location of the resonance peaks on the frequency axes, the amplitude response of the individual curves is 'washed-out' in the average.Similarly, the sharp phase jump between 3 and 4 kHz of the individual curves is lost when taking the average on each frequency.
Gladiné and Dirckx [4] developed a semi-automatic landmark assignment (SLA) method to compute shape-preserving average response curves.SLA attempts to align the LDV curves before averaging.However, shifting the curves does not allow aligning all resonance peaks, since the distance between the peaks in one curve on the frequency axis can differ between measurements.Therefore, a warping function was defined to rescale the data along the frequency axis.First, the authors defined so-called landmarks on the individual curves (e.g., at resonance peaks).Then, the average location of the landmarks on the frequency axis for all measurements was determined.Finally, the individual landmarks were shifted to the mean location, and the response curve in-between the landmarks was deformed using the warping function.While SLA gave good results, the authors of [4] discussed some drawbacks.First, the number of features needs to be chosen, which requires knowledge of the expected behavior of the response curves.Secondly, the data is used piecewise over the frequency axis, and the intervals and the amount of datapoints between measurements needs to be the same.A final remark is that all curves contribute equally to the average.The authors of [4] manually removed one of the curves to account for a measurement that deviated significantly from the other five measurements.Therefore, it would be valuable to develop a method which does not require a priori knowledge of the behavior of the response curve (the locations of the landmarks), and achieves curve alignment with fewer degrees of freedom to minimize the effect of measurement noise on the results.
This paper presents a method to determine the average response curve of frequency response data, which preserves the shape of the individual measurement curves.The procedure is based on frequency response fitting (FR-fit), and the theory of transfer functions is used.Transfer functions describe the behavior of oscillating systems using rational polynomials (see Materials and Methods).An automated way to choose the best order of the polynomials during the data fitting is presented.It is shown that no a priori knowledge of the system's behavior is required, and that the method allows identification of the location of the resonance peaks on the frequency axis.The stapes velocity data of Niklasson et al. [8] is used as a real-life illustrative example.Still, the FR-fit method is general, and any vibrating system can be analyzed using the method presented here.Additionally, we compare our FR-fit method with the results of the SLA approach of Gladiné and Dirckx [4] (see Discussion).

Describing frequency-response data using rational polynomials
As seen in Fig. 1, the velocity response of the stapes is complicated and highly frequency dependent.When systems such as the linear harmonic oscillator are studied analytically, the solutions to the equation of motion can often be written as a rational polynomial [1].The motion's amplitude and phase are then computed for each frequency via a ratio of two polynomials.This function is often called the transfer function of the system.Mathematically, the system's transfer function H(a, b, s) can be expressed as In Eq. ( 1), the vector a holds all the coefficients a i of the numerator with order n z .The order n z determines the number of zeros in the polynomial of the numerator.Similarly, the order n p determines the number of zeros in the denominator.The coefficients b j are grouped in the vector b.Since zeros in the denominator increase the response amplitude, these zeros are called poles.For the present paper, the poles' locations provide insight into the locations of increased stapes velocity and will define the peaks in Fig. 1.The value of s in Eq. ( 1) is the input frequency.However, transfer functions are computed in analytic problems via a Fourier (or Laplace) transform of the original differential equation [1].Therefore, we do not use the original frequency f as the function input, but the frequency expressed in radians per second ω = 2πf.Moreover, the data are transformed to the complex domain, so the final relation becomes s = iω, with i the imaginairy unit.Since s is complex, the value of H(a, b, s) is also complex.The amplitude of H(a, b, s) provides the velocity amplitude, and the phase angle provides the phase of the motion.

Curve fitting procedure and determining average frequency response
Each of the six stapes velocity curves of Niklasson et al. (2018) given in Fig. 1 will be denoted by the complex value h k (s) in what follows.In brief, our FR-fit method consists out of the following steps: by minimizing the squared difference of h k (s) and H(a k , b k , s).This solution depends on the initial choice of n z and n p .
2. Next, compute the average response curve H(ã, b, s) using the individual curve coefficients a k and b k .
3.Then, compute the deviation of the average curve H(ã, b, s) from the individual curves.In this way, the goodness of fit of the average curve with the entire dataset is known.4. Repeat steps 1-3 for different choices of n z and n p .Higher n z and n p allow to model more complicated behavior, but may lead to overfitting and will be more susceptible to noise in the data.5. Finally, choose the best result of step 4. The optimal choice of n p determines the number of poles and thus the number of resonance peaks.From the locations of the poles H(ã, b, s) on the frequency axes, the system's resonant frequencies can be determined.
In step 1, the MatLab command 'tfest' was used to determine the  1) is non-linear in the terms a k and b k , an iterative algorithm was used to determine the best possible solution.We refer to the MatLab documentation of 'tfest' [7] for more details on the exact fitting algorithm.
In step 2, the average response H(ã, b, s) needs to be computed based on the H(a k , b k , s) curves.To achieve this goal, the following computations were performed Eq. ( 2) determines the average coefficients ã and b by weighted averaging over the individual polynomial coefficients a k and b k .Then, Eq. ( 3) defines the new average response curve as the rational polynomial with coefficients ã and b.The weights w i are computed as follows.First, we compute the roots (poles) of the b k polynomial for each dataset.Next, we compute the distance of the poles to the median pole location of all curves.Using the median allows for a robust estimation of the true pole location, since the number of curves is often limited.Finally, we take the inverse of the distances as our weights w i , so curves with large distances to the median poles contribute less to the average curve.
Step 3 determines the goodness of the fit with the current choice of n z and n p .Hence, a measure of the error φ k (s) between the average curve and the individual curves is required.While it is possible to use the squared difference of the H (a k , b k , s) and h k (s) as in Step 1, we found better results by using The asterisk in Eq. ( 4) denotes the complex conjugate.By comparing the logarithm of the squared amplitudes in Eq. ( 4), the high-and lowfrequency ranges also contribute meaningfully to the error measure φ k (s).When a linear scale is used, the large amplitude of the resonance peaks causes the contribution of the low and high frequencies to be negligible.When a range of n z and n p values are computed, the lowest value of φ k (s) gives the optimal combination of poles and zeros to describe the current dataset.

Results
Using Eq. ( 4) and steps 1-5 outlined previously, it was found that the best combination of n z and n p was {4,4} for the current dataset.Thus, the order of the numerator and denominator were the same and of order 4. In general, the number of poles would then be four.However, the polynomial coefficients of H(a k ,b k ,s), and of the average curve H(ã, b,s), are real-valued and not complex-valued.Therefore, the poles of the individual curves H(a k , b k , s) will be either real or appear in complexvalued pairs.It was found that for all of the six H(a k , b k , s) response curves, the poles appeared in two complex conjugate pairs.Therefore, there were only two unique poles on the frequency axis.
Fig. 2 shows a boxplot to inspect the distribution of the two pole frequencies for all of the six curves.Since the two poles differ in frequency, we standardized the frequency values.Standardization was performed by subtracting the median pole frequency from each of the two sets of poles and dividing by the median absolute deviation.
Standardization accounts for differences in the spread of the pole frequencies due to the larger frequency value of the second pole.
In Fig. 2, each measurement is shown via a gray circle.It can be seen that the variation on the first pole is more prominent than on the second pole.For the first pole, measurement 5 produces the largest estimation of the pole frequency, while measurement 3 results in the smallest frequency estimation.For the second pole, measurement 6 results in the largest frequency estimation.Measurement 5 results in the smallest estimation of the second pole frequency.Measurement 5 produced the most spread for both poles on the results in Fig. 2. Gladiné and Dirckx [4] argued that measurement 5 may be an outlier based on similar observations and manually removed it from it their computations.In the method presented here, such manual intervention is not required.The weights w i of data 5 will be small due to the large distance to the median pole value.
In Fig. 3, the results of our FR-fit approach are visualized.Fig. 3A shows the amplitude, and Fig. 3B shows the phase.Two curves are shown using our FR-fit approach.The first curve (in red with rightpointing triangles) was computed by setting all w i weights to a value of one, i.e., with all data contributing equally.The second curve (in green with left-pointing triangles) used the automatically computed weighted average as outlined above.Additionally, we show the experimental data of Niklasson et al. [8] (solid gray lines) and the results of the SLA method [4] (blue curve with square markers).The black curve with circular markers denotes the regular average computed per frequency value.The two vertical dashed lines show the median pole locations.The median pole values found were 1.14 ± 0.13 kHz and 3.61 ± 0.43 kHz, respectively.The coefficients of the weighted curve are given in the appendix in Table .1 for researchers wishing to replicate our results.

Discussion
Fig. 3 shows the proposed steps of Sec.3.2 lead to an accurate description of the average stapes response.It was found that the combination n z = n p = 4 resulted in the lowest error via Eq.( 4).This led to the identification of two unique poles.The presence of two poles is expected from the literature, since both the middle ear resonance [2,5,8] and the ear canal resonance [6,8,9] should be present in the data.The median pole values found were 1.14 ± 0.13 kHz and 3.61 ± 0.43 kHz, respectively.When the weighting w i of Eq. ( 2) was omitted, the resulting FR-fit (with right-pointing triangles) gave poles at slightly higher frequencies than with the weighting (green with left-pointing triangles).The clearest effect is that the weighted curve has a larger amplitude at the second pole location, which almost coincides with the result of the SLA method at that frequency.
Compared to SLA, our FR-fit method has some important advantages.In Ref. [4], it was argued that curve 5 was an outlier.However, SLA needs to define an average landmark location to determine the shift of each of the measurement curves.Therefore, one needs to select the curves used in the average before an additional curve can be compared.This approach assumes that curve 5 is an outlier to omit it from the average.Our method does not have this drawback.Since we first determine the poles of each curve, see Fig. 2, we can look at the individual pole locations on the frequency axis.By using the distance of the poles of the data to the median pole frequency value, curves with large deviations contribute little to the average response curve.A second related advantage is that our method immediately identifies the poles of the response curve.When computer models of the middle ear are built, the resonance location often changes via changing the model parameters, such as the stiffness.The pole locations are immediately returned using FR-fit, while SLA would still require some additional steps to estimate the pole locations.
One consequence of our method is that it has a smoothing effect.Compared to SLA, the degrees of freedom of our method are much lower.We used ten polynomial coefficients in total, and primarily the five denominator coefficients determine the resonance behavior.SLA used 100 equally spaced frequency intervals with a cubic spline interpolation, which resulted in 102 degrees of freedom.Note that this is the major difference between the method presented here and SLA.Our method uses all velocity data over the frequencies measured to determine the transfer function.I.e., no resampling and subdivision over the frequency axis are required.The result is expressed on a polynomial basis, requiring the ten polynomial coefficients (see appendix).SLA requires the data to be resampled and subdivided over the frequency range, which may introduce interpolation artifacts.The result of SLA is expressed on a cubic piecewise basis, which is more complicated to share and implement for other users.
SLA may be more advantageous than the present method depending on the application.If large amplitude or phase variations are expected on a small portion of the frequency axis, the order of the polynomials in the FR-fit would need to be high, which lowers the numerical stability.SLA would not suffer from a poorer fit, since it uses the original data curves to align it to form the average.However, SLA may have difficulty assigning the landmarks before the curves are aligned.If the response changes slowly along the frequency axis (as in Fig. 3 up to 6 kHz) FR-fit can better smooth out minor local variations and provide a clear overview of the average system response.Compared with the classical average (Fig. 3, black-dotted curve), both SLA and FR-fit provide a better description of the average system response.

Conclusion
This paper presented a novel method to determine the average frequency response curve with preservation of response characteristics of individual experimentally acquired curves.We used rational polynomials to perform a weighted frequency response curve fit (FR-fit) on human stapes velocity data.An order of four in both the numerator and denominator polynomials gave the best FR-fit results, which was expected from the literature.The two unique poles (zeros of the denominator polynomial) directly identified the two resonance locations for Fig. 3. A) Stapes velocity amplitude in dB SPL.B) Velocity phase in cycles.The experimental data of Niklasson et al. [8] is shown in solid gray.Our frequency response fit (FR-fit) with fixed weights of one is shown in red with right-pointing triangles.The weighted FR-fit curve is shown in green with left-pointing triangles.The blue line with square markers denotes the SLA curve of Gladiné and Dirckx [3].The black curve with circular markers denotes the arithmetic average per frequency value.The two vertical dashed lines show the median pole locations of 1.14 ± 0.13 kHz and 3.61 ± 0.43 kHz, respectively.(For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)each of the six individual curves.
The average curve was determined from the individual measured curves by weighted averaging individual FR-fit polynomial coefficients.It was found that the stapes vibration showed two resonance peaks at 1.14 ± 0.13 kHz and 3.61 ± 0.43 kHz.Compared to the classical average, our method provides a more realistic description of the general population's generic average frequency response of the footplate.Compared to other methods, FR-fit does not require a priori knowledge of the locations of the peaks in the data and is fully automatic.However, a smoothing effect is present and rapid variations in velocity may be missed using FR-fit.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Fig. 1 .
Fig. 1.A) Stapes velocity amplitude at different frequencies of sound measured by Niklasson et al. [8].Gray lines denote individual measurements.The black curve shows the average stapes response computed by averaging the velocities for each frequency.A "washed-out" view of the stapes response is seen, which does not represent the behavior of the individual curves.B) The corresponding phase of the stapes velocity.Similarly, the average smooths out the sharp phase gradients and does not accurately represent the stapes motion.

Fig. 2 .
Fig. 2. Boxplot of the two unique poles found by fitting the six H(a k , b k , s) curves.The frequency values of the two poles were standardized to account for scaling effects due to the higher value of the second pole on the frequency axis.The (blue) box contains the 25th to 50th percentile of the data.The (red) line inside the box denotes the median value of the pole, which is zero due to the standardization.Each individual measurement is shown via a gray circle and numbered accordingly.(For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)