WO2020257188A1 - System and method for measuring causality-constrained acoustic reflectance and reflection function of the ear - Google Patents

Abstract

Methods and systems are used to determine a causal response function of the ear. A probe assembly is inserted into each of one or more acoustic calibration waveguides. The probe assembly includes an acoustic source and a microphone. Calibration signals are received from the microphone when the probe assembly is inserted into each of the acoustic calibration waveguides and the acoustic source generates sound responsive to an electrical input signal. Time domain and/or frequency domain calibration waveguide model response functions are used that satisfy causality and are indicative of acoustic transfer characteristic with viscothermal effects for acoustic transmission in each acoustic calibration waveguide. The calibration waveguide model response functions are used to determine measurement system parameters that satisfy causality. The measurement system parameters are used with an acoustical response measured by the microphone in the ear canal to determine the causal response function of the ear.

Description

SYSTEM AND METHOD

FOR MEASURING CAUSALITY-CONSTRAINED ACOUSTIC REFLECTANCE AND

REFLECTION FUNCTION OF THE EAR

This application claims the benefit of U.S. Provisional Application Serial No.

62/862,448, filed June 17, 2019 and U.S. Provisional Application Serial No. 63/019,113, filed May 1, 2020, both of which are incorporated by reference herein.

SUMMARY

Embodiments described below include systems and methods for measuring causality- constrained acoustic reflectance and reflection function of the ear. In one embodiment, a probe assembly is inserted into each of one or more acoustic calibration waveguides. The probe assembly includes an acoustic source and a microphone. Calibration signals are received from the microphone when the probe assembly is inserted into each of the acoustic calibration waveguides and the acoustic source generates sound responsive to an electrical input signal. One or more calibration waveguide model response functions are determined in the time domain that each satisfies causality for a corresponding acoustic calibration waveguide. Each calibration waveguide model response function is indicative of an acoustic transfer characteristic with viscothermal effects that is a function of one or more waveguide model parameters, with each waveguide model parameter value initialized to a plausible estimate.

An outer loop calculation is performed to determine values of final measurement system parameter signals that are causal. Each iteration within the outer loop involves: determining values of causal measurement system parameter signals in the time domain by solving a first over-determined matrix equation, which includes functions of the calibration signals and the calibration waveguide model response functions, with each function calculated in terms of the one or more waveguide model parameter values; and determining an update of the waveguide model parameter values by minimizing a residual norm of a second over-determined matrix equation, which includes functions of the calibration signals, the measurement system parameter values, and the waveguide model response functions. The outer loop calculation concludes when the residual norm is minimized. Values of causal final measurement system parameters signals in the time domain are determined as values of the measurement system parameters determined at a conclusion of the outer loop calculation.

A measurement signal is received from the microphone when the probe assembly is inserted into the ear canal and the acoustic source generates sound responsive to the electrical input signal. The causal response function of the ear is determined in the time domain based on the measurement signal in the ear and the final measurement system parameters.

In another embodiment, a probe assembly is inserted into each of one or more acoustic calibration waveguides. The probe assembly includes an acoustic source and a microphone. Calibration signals are received from the microphone when the probe assembly is inserted into each of the acoustic calibration waveguides and the acoustic source generates sound responsive to an electrical input signal. Calibration waveguide model response functions in the frequency domain are determined for each of the acoustic calibration waveguides that are indicative of an acoustic transfer characteristic with viscothermal effects for acoustic transmission in each acoustic calibration waveguide. Each calibration waveguide model response function in the frequency domain is a function of one or more waveguide model parameters. Final measurement system parameters are determined in the frequency domain such that each satisfies causality by solving system equations that contain the calibration signals and the calibration waveguide model response function for each of the calibration signals.

A measurement signal in the ear is received from the microphone after the probe assembly is inserted into the ear canal and the acoustic source generates sound with input from the electrical input signal. A causal response function of the ear is determined in the frequency domain based on the measurement signal in the ear and the final measurement system

parameters.

BRIEF DESCRIPTION OF THE DRAWINGS

FIG. l is a block diagram of a system according to an example embodiment;

FIG. 2 is a block diagram of a probe assembly according to an example embodiment;

FIG. 3 is a set of graphs showing reflectance functions of acoustic calibration waveguides according to an example embodiment; and FIGS 4-5 are flowcharts of procedures according to example embodiments.

DETAILED DESCRIPTION

U.S. Patent No. 8,840,565 B2 to Keefe and entitled TEST BATTERY SYSTEM AND METHOD FOR ASSESSMENT OF AUDITORY FUNCTION, U.S. Patent No. 5,792,072 to Keefe and entitled SYSTEM AND METHOD FOR MEASURING ACOUSTIC

REFLECTANCE, U.S. Patent Application No. 62/862,448 to Keefe entitled SYSTEM AND METHOD FOR MEASURING CAUSALITY-CONSTRAINED ACOUSTIC REFLECTANCE AND REFLECTION FUNCTION OF THE EAR filed June 17, 2019, and U.S. Patent

Application No. 63/019,113 to Keefe entitled SYSTEM AND METHOD FOR MEASURING CAUSALITY-CONSTRAINED ACOUSTIC REFLECTANCE AND REFLECTION

FUNCTION OF THE EAR filed May 1, 2020, are herein incorporated by reference.

The various embodiments of techniques described herein can be utilized with any suitable system or assembly, e.g., the systems described in U.S. Patent Nos. 8,840,565 B2 and 5,792,072. For example, FIGS. 1-2 are various views of one embodiment of a system 200 that can utilize or incorporate the various techniques described herein. The system 200 includes a probe assembly 225, one or more acoustic calibration waveguides 227 (FIG. 2), a stimulus signal generator 213, a main signal processor 210, and processing circuitry (e.g., as a part of the main signal processor).

The probe assembly 225 can be adapted to be inserted into an ear canal of a patient. The assembly 225 can include an acoustic source 216 (DRIVER) adapted to generate an acoustic stimulus in response to an electrical input signal, and a microphone 223, which can also be an acoustic transducer, to convert acoustic pressure to a detected electrical signal.

The main signal processor 210 can generate a digital electrical input, or stimulus, signal 214 that can be converted to an analog signal 215, e.g., by a digital to analog converter no. 1 (DAC1) 213. In one or more embodiments, the system 200 may lack the DAC1 213, and the analog signal 215 is a copy of the stimulus signal 214. The analog signal 215 is input to a driver 216 (e.g., a miniature loudspeaker).

In one or more embodiments, the acoustic source 216 and the microphone 223 can be disposed in the probe assembly 225. In one or more embodiments, the acoustic source 216 and microphone 223 can be separate from the probe assembly 225 and connected to the probe assembly using any suitable technique or techniques. In one or more embodiments, each calibration waveguide 227 has an internal air volume closed at its far end that is a cylindrical volume with waveguide model parameters of cylindrical length and radius, and is thus termed a calibration tube. The probe assembly 225 is inserted into the open end of each calibration tube in a leak-free manner during the calibration procedure. For example, the acoustic source 216 can be the source of acoustic stimulus signal that is delivered to the ear through a small tube or set of tubes 217 to a first part of a probe assembly 225 and then through a second tube or set of tubes 218 in a second part of the probe assembly 225, from which the acoustic stimulus 219 is delivered into a calibration waveguide or the ear canal of a human patient. Prior to the measurement of a response to the acoustic stimulus 219, the probe assembly 225 can be inserted in a substantially leak-free manner in the calibration waveguide or ear canal. In one or more embodiments, the first and second parts of the probe assembly 225 are identical and the tube or set of tubes 217 is identical to the tube or set of tubes 218.

In one or more embodiments, the signal processor 210 is adapted to receive a set of detected calibration electrical signals from the microphone when the probe assembly 225 is inserted into each of the one or more acoustic calibration waveguides or tubes 227 and receive a detected measurement electrical signal when the probe assembly is inserted into the ear canal. In one or more embodiments, an acoustic response 220 of the external and middle-ear, or of a calibration waveguide, to the acoustic stimulus 219 is transmitted into an inlet tube or set of tubes 221 of the probe assembly 225 and then through a tube or set of tubes 222 to the input of the microphone 223. The electrical output 224 of the microphone 223 is a detected electrical analog signal that is converted into a digital response signal 212 by an analog to digital converter (ADC) 211. The response signal 212 is recorded by the main signal processor 210, which can be synchronized in time to the digital electrical input signal 214. In one or more embodiments, the system 200 lacks the ADC 211, and the response signal 212 is a copy of the electrical output 224. In one or more embodiments, the inlet tube or set of tubes 221 is identical to the tube or set of tubes 222.

In one or more embodiments, the signal processor 210 is further adapted to receive the set of detected calibration electrical signals recorded in the set of calibration waveguides 227 and detect in at least one of the set of detected calibration electrical signals an incident signal, the incident signal being separable in time from a first reflected signal in the same detected calibration electrical signal.

Substantially simultaneously with the presentation of the stimulus signal 214 and the recording of the response signal 212, the main signal processor 210 can optionally transmit an intended air-pressure signal 205 to an optional pump control unit 204, and the pump control unit 204 can transmit a driver signal 207 to the electrical input of a pressure pump 208. The pump output generates a volume displacement of air into the tube or set of tubes 209 that is coupled to the probe assembly 225 sealed in the ear canal. This volume displacement of air generates an air pressure within the ear canal, which is sensed by a pressure transducer 201 via inlet tubes 202 and 209 through the probe assembly 225 that is coupled to the volume of air in the ear canal. The term“air pressure within the ear canal” refers to the equilibrium air pressure within the ear canal over time scales longer than the time scale of the acoustical disturbance of the air, which is conveyed by the term“acoustic pressure within the ear canal.” The air-pressure signal can vary with time to produce a desired time-varying air pressure during any measurement, e.g., in a pressure sweep from positive to negative air pressures (relative to the atmospheric pressure of air), or from negative to positive air pressures.

Multiple drivers may be utilized when the probe assembly 225 is used to measure an acoustic reflectance response with an otoacoustic emission (OAE) response or with middle-ear muscle reflex (MEMR) threshold and supra-threshold responses. Reflectance, OAE and MEMR responses of the ear are useful diagnostic acoustic test responses to detect ear disorders. The system 200 can measure an acoustic response function of the ear for purposes of calibrating the sound field in the ear canal at high frequencies up to approximately 20 kHz to specify the acoustic input in other electrophysiological or behavioral tests, e.g., an auditory brainstem response test, a behavioral audiogram, and a behavioral test of spatial hearing abilities using measurements in both ears of a subject.

The operation of the system 200 includes a set of measured sound pressure signals in one or more waveguides and in the test ear, a set of calibration waveguide model response functions, and a causal response function of the ear. Each such response function is a signal. Any signal can be specified in the time or frequency domain. Time-domain signals may be represented in continuous-time or discrete-time descriptions, and frequency-domain signals may be represented in continuous-frequency or discrete-frequency descriptions. Those skilled in the art are familiar with discrete-time sampling of a continuous-time signal at a fixed sample rate to determine a discrete-time signal, from which the corresponding frequency-domain signals can be calculated in a continuous-frequency representation using the Fourier transform (FT) or in a discrete- frequency representation using the discrete Fourier transform (DFT), In this specification, time- domain signals are named in lower case and frequency-domain signals in upper case, except for matrix and vector notations in several places. The sample period T s the inverse of the sample rate.

The probe assembly in system 200 is inserted into an open end of one or more acoustic calibration waveguides. The inputs in system 200 include a microphone measurement of sound pressure in one or more waveguides of known dimensions in the form of a detected calibration electrical signal that is output by the microphone pre-amplifier, in which the probe assembly is inserted into one open end of the waveguide and its opposite end is closed. The system 200 includes for each waveguide a calibration waveguide model response function in the time domain or in the frequency domain for each of these waveguides. Each such model response function is specified in terms of one or more waveguide parameters whose values are determined by system 200. The calibration waveguide model response function in the time domain satisfies causality, as its response is zero at times earlier than time zero. The calibration waveguide model response function in the frequency domain can be constrained to satisfy causality as explained below. It has been used in previous systems without this constraint, thereby adding a source of error to the measurement system.

As shown in FIG. 4, the operation of the system 200 can include a calibration procedure 400 followed by an evaluation procedure in 401 to calculate a causal response function of the ear in the frequency domain. Corresponding calibration and evaluation procedures in FIG. 5 are used to calculate a causal response function of the ear in the time domain. The inputs to the calibration procedure in FIG. 5 include the set of detected calibration electrical signals 600 and the set of waveguide model parameters 612 and 606 used to calculate the set of calibration waveguide model response functions 602. The outputs from the calibration procedure are the final measurement system parameters 608, which then serve as inputs to the evaluation procedure, which is shown below the dashed line in FIG. 5, along with the detected measurement signal in the ear 609. The output of the evaluation procedure 610 is a causal response function of the ear in the time domain, which can be converted to the frequency domain using a DFT. The calibration waveguide 227 can have an open first end to permit the insertion of the probe assembly 225 and a closed second end that reduces the effects of ambient noise and can be modeled. However, the calibration waveguide 225 may be virtually any shape or dimension and have a closed second end or an open second end so long as the acoustical transfer characteristics of the calibration waveguide can be modeled.

A specific embodiment of a calibration waveguide model causal response function applies for a calibration waveguide 227 in the form of a straight hard-walled cylindrical tube closed at its end opposite to the end in which the probe assembly is inserted. A calibration waveguide model response function is constructed to satisfy causality, whether expressed analytically in the time domain or frequency domain. If the calibration model causal response function is expressed analytically in the frequency domain, then it may also be expressed numerically in the time domain using a waveguide digital filter based on an input invariance design method. The preferred embodiment is a calibration waveguide model causal response function representing the sound field in a closed cylindrical tube expressed analytically in the time domain, which may optionally include an analytical model component to describe evanescent-mode effects in the tube near the probe assembly. In one or more embodiments, each calibration waveguide model causal response function includes an evanescent-mode contribution near a tip of the probe assembly 225 that is a function of the geometry of the probe assembly and the one or more waveguide model parameters. This tip is the location in the ear canal into which the acoustic stimulus 219 is delivered and the acoustic response is transmitted to the microphone 223.

The processing circuitry can be adapted to determine a causal response function of the ear in the time domain using any suitable technique or techniques. In one or more embodiments, the processing circuitry can be adapted to determine a causal response function of the ear in the time domain based on the set of detected calibration electrical signals, the detected measurement electrical signal in the ear, and a calibration waveguide model causal response function for each of the one or more acoustic calibration waveguides that is indicative of a causal acoustic transfer characteristic with viscothermal effects for acoustic transmission, and such that the calibration waveguide model causal response function in the time domain is a function of one or more waveguide model parameters The processing circuitry can also be adapted to determine one or more measurement system parameters or calibration outputs using any suitable causally-constrained technique or techniques, e.g., one or more of the techniques described in Keefe, D.H (2020),“Causality- constrained measurements of aural acoustic reflectance and reflection functions,” J. Acoust. Soc. Am 147(1), 300-324, with supplementary material published by the journal on line (hereinafter “Keefe, 2020”), which is included in this application and incorporated by reference. In one or more embodiments, the final measurement system parameters 608 of a calibration in the time domain in FIG. 5 can be determined by calculations in an iterative outer loop designed to calculate a final set of waveguide model parameters. The set of calibration waveguide model response functions 602 for each calibration waveguide is calculated as functions of the waveguide model parameters 606, whose values are varied with the outer loop. The set of calibration waveguide functions and set of calibration signals compose a first over-determined matrix equation, which can be solved for the measurement system parameters 604 using a direct or iterative regularized deconvolution approach.

The calculations within the outer loop vary the values of the waveguide model parameters 606 in each iteration so as to calculate different values of the calibration waveguide model response functions 602 and different values of the measurement system parameters 604 in each iteration. The outer loop procedures minimizes the residual norm of a second over determined matrix equation, which is a function of the calibration signals (that are the same on each iteration), the measurement system parameters 604 obtained by solving the first over determined matrix equation in each iteration, and the calibration waveguide model response functions 602 in each iteration. If the residual norm is not yet minimized in 607, then the waveguide model parameters 606 are updated and are used in the next iteration of the outer loop. If the residual norm is minimized in 607 according to a calculated criterion, the final

measurement system parameters are determined as those calculated in the final iteration of the outer loop.

The causality constraints imposed in the calibration procedure are that: (1) each calibration waveguide model response is causal in the time domain and the causality-constrained in the frequency domain, which means that the frequency domain response is the FT or DFT of a causal response, and (2) each wave measurement system parameter, which is a time-domain signal, is causal. The first and second over-determined matrix equations can be solved using techniques of regularized deconvolution, including a criterion to choose a regularization parameter from a range of values within the regularization procedure that achieves a balance of minimizing the residual norm and avoiding too large of a solution norm. Further, in one or more embodiments, the processing circuitry is further adapted to determine a causal response of the ear in the frequency domain using the determined causal response of the ear in the time domain.

The processing circuitry can determine final measurement system parameters or final calibration outputs in the frequency domain in 401 in FIG. 4 such that each final measurement system parameter is constrained to be a causal signal by solving system equations that contain the set of detected calibration electrical signals and each calibration waveguide model response 400 in the frequency domain for each detected calibration electrical signal of the set of detected calibration signals. The processing circuitry can also determine a causal response of the ear in the frequency domain by calculating a set of equations that includes the electrical input signal and the causally-constrained final measurement system parameters or final calibration outputs.

Further, in one or more embodiments, the processing circuitry is adapted to determine a causal response of the ear in the time domain using the determined causal response of the ear in the frequency domain.

Similarly, the processing circuitry can determine a tympanometric causal response function of the ear in at least one of the time domain or frequency domain. A causal response function of the ear in the frequency domain can be calculated as the discrete Fourier transform of the causal response function of the ear in the time domain. A tympanometric causal response function of the ear in the frequency domain can be calculated as the discrete Fourier transform of the tympanometric causal response function of the ear in the time domain.

The acoustic characteristic impedance at the probe tip is the input impedance of the waveguide for a case in which there is a forward or outgoing sound wave from the probe tip with no reverse or incoming sound wave from elsewhere within the waveguide. The characteristic impedance Z_c of the waveguide is Z_c = ci S, in which the equilibrium density of air is p, the free-space phase velocity of air is c, and the circular cross-sectional area of the waveguide (e.g., a cylindrical tube) is S. Viscothermal, or lossy, corrections to this expression for Z_c are much smaller than those associated with the acoustic propagation wavenumber and can be neglected in one or more embodiments. Various examples of measurement system parameters in system 200 include the incident pressure and source reflectance of the probe at its probe tip specified in terms of the

characteristic impedance ¾of the waveguide, from which can be calculated in the frequency domain the Thevenin pressure P_T(f) and Thevenin impedance P_T{f) of the probe, and the

Norton volume velocity U_N(f) and Norton admittance Y_N(f) of the probe. These Thevenin and

Norton variables are generally specified in the frequency domain, but their corresponding time- domain variables can be formally calculated using the inverse Laplace or inverse Fourier transform, or using a digital waveguide filter. Alternatively, they can also be directly calculated in terms of the incident pressure, source RF and characteristic impedance. Thus, the

measurement system parameters are interchangeable in the frequency domain between the reflectance, Thevenin and Norton pairs of variables. The corresponding transformations are obtained in the time domain using operations of convolution and deconvolution. The

deconvolution calculation can be numerically unstable, and direct or iterative regularization procedures are available to obtain numerical results from a deconvolution that are more stable. Use of an inverse Laplace or inverse Fourier transform procedure, or of a digital waveguide filter, to calculate measurement system parameters in the time domain can result in non-causal responses, thereby introducing a source of error into a calculation of a response function of the ear in the time domain.

Reflection Functions and the Reflectance Waveguide Equation

An example is described for the case that each acoustic calibration waveguide is a cylindrical tube closed at the end opposite to the end in which the probe is inserted. A cylindrical tube includes the case in which a cylindrical volume of air of finite length with an open end is contained in some material, for example, as a cavity of air.

The reflectance waveguide equation (RWE),

describes sound transmission in the continuous-time domain in a single propagating mode with convolution operator * in terms of the Dirac delta function 5(t) . These waveforms, which are referenced to the probe-tip location within the waveguide, include waveguide pressure p(t) , incident pressure /¾(¾ source RF T₀(i), and waveguide RF r(t). The measurement system parameters associated with the RWE are incident pressure and source RF. Once these are determined, the RF of an ear r(t) may be calculated using the RWE by measuring the pressure p(t) at the probe tip, once inserted in a leak-free manner into the ear canal.

The number of unknown measurement system parameter signals is reduced from two to source reflectance alone if a sufficiently long waveguide or tube is used such that the incident pressure is separable in time from the first-reflected sound. This is the“known-po” type of RWE, as contrasted with the“unknown-po” type when incident pressure cannot be directly estimated from a single tube recording due to temporal overlap of multiple internal reflections. An analogous discrete-time version of the RWE is used for practical calculations.

The calibration cylindrical tube length L is sufficiently short that its round-trip delay time t = 2 L i e is shorter than the measurement duration. Viscous and thermal boundary layers at the wall are each small relative to the tube radius a. In this“small boundary layer” limit, the model RF r{t) of a calibration tube closed at its far end is,

in which o (t) is the Heaviside step function defined as zero for t <0 and one for t > 0.

Thermodynamic constants of air include the shear viscosity h, thermal conductivity K, specific heat at constant pressure c _, and ratio g of c _p to the specific heat at constant volume. This r(i ) in this example is the calibration waveguide model response function in the time domain and satisfies causality through the action of the Heaviside step function.

The RWE measurement system parameters in the frequency domain are the incident pressure spectrum P_Q(f) and source reflectance -¾(/), which are the Fourier transforms (FTs) of ¾( ^ar|d X respectively. Any FT can be written either as a function of frequency or radian frequency w = 2p/ for frequency /, e.g., as -¾(/ or l^i /J). With a known

characteristic impedance, the RWE measurement system parameters may be used to calculate the Thevenin and Norton measurement system parameters. For a waveguide or test ear, the pressure reflectance R (co ) is the FT of r(t ) .

In the frequency domain, the acoustic transmission line for a calibration cylindrical tube is characterized by its characteristic impedance and complex propagation constant G = G_ί, + jT_i (unit of radian per length). In a small boundary-layer limit applicable to ear response

measurements, the real part G_;. of the propagation constant in this same approximation varies with the square root of the radian frequency and the imaginary part G is,

For any waveguide with plane wave sound transmission, the input impedance Z (co ) is,

1+K(w)

Z(w) = Z_c (4) l-R{a>)

The pressure reflectance K(w ) at the entrance to a cylindrical tube of constant cross-section area S and length L that is closed at its opposite end is,

R(c ) = e ^2l (5) and its input impedance

In the low-amplitude limit of acoustical theory applicable to ear response measurements, a response function of the ear is substantially independent of stimulus level except for any small nonlinear effects due to OAEs or other nonlinear acoustic responses of the ear. The response functions of the ear are acoustic response functions that can be measured using the pressure measured by the probe microphone after insertion into the ear canal of the ear. Examples of response functions of the ear include the following transfer functions defined within the ear canal at the probe tip: pressure reflectance, reflection function, admittance and impedance. The pressure reflectance of the ear is termed the reflectance.

The acoustic response function of the ear can be measured by inserting the probe assembly 225 in the ear canal in a leak-free manner and conducting measurements using the system 200. The driver 216 can produce a short-duration sound field (i.e., an acoustic transient signal 219). The microphone 223 measures the sound pressure, which is the acoustic response 200, and the system 200 can calculate a reflectance or reflection function of the ear from the microphone response. The reflectance or reflection function assesses functional properties of the ear, and especially of the ear canal and middle ear. In addition to the acoustic wave that propagates down the acoustic waveguide, the acoustic response includes a response from non propagating modes. Each non-propagating mode, also called an evanescent mode, refers to an acoustic signals that is non-propagating at sufficiently low frequencies such that the acoustic wavelength is small relative to the circumference of the flexible tubes 218, 221 or ear canal. A sum of responses over multiple evanescent modes can describe the localized acoustic field near the probe assembly 225. Any localized, non-propagating acoustic field caused by evanescent modes in the vicinity of the probe assembly 225 can be attenuated by restricting the frequency content of the external stimulus or by other well-known methods such as drawing the

microphone probe 223 slightly beyond the plane of the driver probe 216, as shown in one or more embodiments in FIG. 2. One or more other embodiments can draw the tube 218 slightly beyond the plane of the driver probe 216. Effects of evanescent modes on the sound field near the probe assembly can also be described by an analytical model, which is known to those skilled in the art, with parameters set based on the geometrical dimensions of the probe and the estimated dimensions of the ear canal near the probe assembly.

Causality-Constrained Calibration in the Frequency Domain

FIG. 2 illustrates an acoustic calibration waveguide 227 (or ear canal) whose cross- section area is similar to the ear-canal area in a patient to be tested, or the entryway area of any other unknown system. The calibration waveguide 227 has an open first end in which the probe assembly 225 is inserted and a second end opposite the first end. The acoustical characteristics of the calibration waveguide 227 are derived from a model of a cylindrical tube. An example of the calibration waveguide 227 is a straight hard-walled cylindrical tube with a circular cross-section and specified value of length and radius. Other examples use a flexible-walled tube or coiled cylindrical tube. Those skilled in the art can appreciate that the calibration waveguide 227 can be an acoustic waveguide of virtually any shape, such as square tube, oval tube, conical tube or the like, whose cross-section as a function of position along the length of the calibration waveguide is known and whose acoustical properties, including viscothermal effects, can also be described by a calibration waveguide model response function. The calibration waveguide model response function can be described in the time domain or in the frequency domain. This disclosure describes the use of determining causal measurement system parameters by use of a calibration waveguide model causal response function. Examples of such causal signals are the incident pressure and source reflection function or source reflectance of the probe, the Thevenin pressure and Thevenin impedance of the probe, and the Norton volume velocity and Norton admittance of the probe.

Measurements are performed in M_t ³ 1 calibration waveguides using any calibration procedure equivalent to the RWE. Such equivalent procedures include use of the Thevenin or Norton acoustical equivalent circuit. The calibration waveguide model response function is calculated in terms of the set of waveguide model parameters for each waveguide, with an initial plausible estimate of each parameter value. The RWE is used as an example with a cylindrical tube for each calibration waveguide. The model pressure reflectance

to tube m with

1 < _m < M _t is calculated using Eq. (5) based on the estimated length L_m of tube w, which is a waveguide model parameter.

In the known-po approach, the incident pressure spectrum P₀[k\ in the k-th frequency bin is inferred by separability of the incident response in time after a first reflected response in the mean pressure measured in a tube 1 (or generally waveguide 1) . The frequency-domain RWE, and the Thevenin and Norton circuits, take the following form for the k-th frequency bin,

A Wm⁼bM (6) in which the matrix components for the RWE in this example are,

The pressure spectrum measured in tube m is

which is calculated from the discrete DFT of the discrete-time pressure waveform P_m[n measured in the tube at time sample n with 0 £ n £ N - 1 . The A(xl column vectors A\k\ and b\k | are defined that contain the data and models for all M_t tubes. Their transposes are row vectors: ^[k\ = A k ...,A_N \k\ j, b^r[k\— (¾[&], .. b_N^ [k]j . The resulting over-determined matrix equation for the unknown x| k \ in general form is,

A k x k = b[k ]. (8)

The matrix A^H\k\ is the Hermitian of A\k |. which is the complex conjugate of the transpose of A[k], and is a lxA^ row vector. The least-squares solution x[k^~\ , which minimizes the residual norm The term in square

brackets in the previous equation is the pseudoinverse Af\fc\ of matrix A\k\ . The residual norm II A[k]x[k] - b\k | || quantifies the errors associated in the least-squares solution in the k-th frequency bin. An overall real non-negative error function E² of the waveguide model parameters is defined by summing the square of the residual norms over all analysis frequencies,

The waveguide model parameters include the M_t tube lengths L_m, and can include other parameters such as those related to viscothermal effects or evanescent-mode effects.

A constrained or unconstrained nonlinear multivariable optimization can be performed to minimize the square of the residual norm E² by varying the waveguide parameter values. The constrained optimization is preferred in that the range of allowed parameter values can be specified in the initialization. The optimization procedure can satisfy inequality constraint equations on each iteration. A sequential quadratic programming algorithm can be used. This algorithm is analogous for constrained optimization to how Newton’s method works for unconstrained optimization. A numerical approximation to the Hessian of the Lagrange function including constraints is calculated in a major iteration using a quasi-Newton method. Each major iteration includes multiple minor iterations using a line-search optimizer. A preferred algorithm uses an active-set algorithm can result in a smaller residual norm. This algorithm uses the Kuhn- Tucker equations that cancel at the solution point the gradients of the objective function E with the gradients of the inequality constraints represented by positive Lagrange multipliers. The active-set algorithm appeared more effective at forcing an iterative solution out of a local minimum of the objective function so as to arrive at a more accurate solution. In the unknown-po RWE approach, the unknowns to be determined at the -th frequency are P₀[k and -¾[ ], for which 2x1 column vector x[k] is defined:

The frequency-domain RWE is re-arranged for tube Wl at the k-th frequency into the general form of Eq.(6), except that AJ c for tube Wl is now a 1x2 matrix and bjji\ is a lxl column vector defined by:

At each frequency, a M_t x 2 matrix A\k \ and A-^xl column vector b\k | are defined that contain the data and models for all N_t tubes. They are defined by inserting this new ^] into row Wl of A[k], and this new b_m\k\ into row Wl of h \ k \ . The 2 M_t Hermitian matrix A^H\k\ is defined as the complex conjugate of Ap\k\. The least-squares solution x[k] of this new form of Eq. (8) is calculated. The new form of the overall real non-negative error function E² is calculated using Eq. (9). The waveguide model parameters remain the same in the unknown-po and known-po approaches. The solution procedure using a constrained nonlinear multivariable optimization procedure to minimize the square of the residual norm E² is as in the above.

This solution to the unknown-po approach is algebraically identical to the calibration obtained by substituting the Thevenin pressure P_T\_k\ and Thevenin impedance A_j\k\ (or the

Norton volume velocity and Norton admittance) in place of the variables P₀[k and k\ . Thus, the RWE in the unknown-po approach is equivalent to using the Thevenin circuit. Numerical results can differ since optimization procedures differ across the choice of approach due to effects of the degree of numerical precision in the calculations.

A signal such as acoustic pressure is causal if its response does not occur earlier in time than the stimulus eliciting the response. If the transfer function is specified in the frequency domain, then a transfer function is causal if the real and imaginary parts of the transfer function are related to one another via a Hilbert transform, which is well understood to one skilled in the art. Such a transfer function is a causal response function of the ear in the frequency domain. The latter means that the calibration phase is causally constrained, so that the measurement system parameters or calibration outputs are a set of causal signals.

Each detected electrical signal from the microphone has real values and is of finite duration /), i.e., the signal is transient. Suppose a real continuous-time signal is w(i) . In one example of a discrete-time description, this signal is sampled at a fixed sample rate such that the sample period T is the inverse of the sample rate. The resulting real, discrete-time, transient signal is w(nT) , in which integer n varies from 0, 1, 2, N-l over the duration D = NT . A real sequence w[n ] = w(nT) satisfies a periodic causality condition if w[n ] = 0 for n < 0 and n³N . This transient signal is embedded in a N-point, periodic sequence [«], which has a discrete Fourier transform (DFT) W\k\ = W_R\k\ + jW_j[k\ for frequency index 0 < k < N. Such a periodic causal sequence is termed herein a causal signal. The DFT is expressed in terms of its real part ¾[^]and imaginary part W_I\k\ with unit imaginary number j. A N-point periodic sequence is periodic causal if w[n ] = 0 for all n modulo Ain N / 2 < n < N , i.e., w[n ] = 0 for -A / 2 < n < 0. Even and odd parts of w[n\ for -A / 2 < n < A / 2 are w_e [n] = ( w[n\ + w[-«]) / 2 and w_o[n] = (w[«] -w[-«]) / 2, respectively. For real u|«|,

is the inverse discrete Fourier transform (IDFT) of ¾[ ] and W₀[ /?] is the IDFT of jW_j\k\ . Moreover, ¾[ ] is even in & and W_j\k\ is odd. It follows that a period-N causal w[«] is fully specified by its even part, namely, n] = _nw_e[nl (12) in which z₀ =1 and z_h = 2 for 0 < n £ A / 2. The mathematical formalism in this paragraph is described in A.V. Oppenheim and R.W. Schafer,“Discrete-Time Signal Processing,” (Prentice Hall, New York, 1989).

A least-squares reflectance calibration has N_c unknowns. Because each unknown variable is physically known to be real and causal, then the imaginary part of each spectral variable is known once the real part of each spectral variable is calculated. It suffices to construct an unknown column vector x[k] as a concatenation of the real parts of the spectra of each unknown variable. These unknown variables are estimated by solving the following over determined system of equations:

A [k ] x[k ] = b [k ] , (13) in which Am and b [k ] are defined using the set of calibration signals from the microphone and the calibration waveguide model response function. The solution in terms of the pseudoinverse matrix A ⁺ [k] of Am is,

x [k ] = A ⁺ [k ]b [k ] (14)

Equation (13) is an over-determined set of equations that is solved in Eq. (14) by a least- squares procedure for x| k | using the pseudoinverse. Suppose that x| k | contains N_c variables

W_jR[k for 1 £ i £ N_c such that W_iR\fc\ is the real part of the spectrum W^k\ of the i-th periodic causal variable to be determined along with its unknown imaginary part W_a\_k]. The estimated values of the complex spectra W^k\ and waveform w [/?] for the i-th output variable are calculated from x[k ] by a causality-constrained procedure having these steps:

1. Calculate Am and /,^' _| A· _| from measured data and model with current waveguide

model parameter values,

2. Calculate x[k] as the least-squares solution of Eq. (13),

3. Minimize E² = ^|| A[k]x[k] -b[k] ||² over waveguide model parameters as in Eq.(9), k

4. For each spectral variable ¾[ ] in x[k]

• Calculate even part _{w , e [« ]} as the IDFT of W_jR\fc\ ,

• Calculate _t\ri\ from

using Eq. (12),

• Calculate W k\ as the DFT of W [/?] .

Iterate on steps 1-3 until the stopping criterion of the minimization is attained. Once the least-squares solution for x[k ] is obtained for all k , step 4 calculates each periodic causal waveform w- [/?] and its complex spectrum

The spectrum W_a\k\ is obtained as the imaginary part of W \ . A major benefit of the approach is that the solution is always causal and its form in the frequency is said to be causally-constrained. This benefit is not necessarily so for a standard calibration approach. Errors can remain due to any inaccuracies in the least-squares approximate solution.

In FIG. 4, a flowchart summarizes the steps in the prodecure above according to an example embodiment. Calibration waveguide model response functions are determined 400 in the frequency domain for each of the acoustic calibration waveguides (step 1 above). The calibration waveguide model response functions are indicative of an acoustic transfer characteristic with viscothermal effects for acoustic transmission in each acoustic calibration waveguide, and such that the calibration waveguide model response function in the frequency domain is a function of one or more waveguide model parameters. Final measurement system parameters in the frequency domain are determined 401 such that each satisfies causality by solving system equations that contain the calibration signals and the calibration waveguide model response function for each of the calibration signals (steps 2-4 above).

In real measurements, the length N of the sequence Ή [/^z] can be chosen sufficiently large so that the resulting transient w[/?] satisfies the periodic causal condition for each i. A major benefit of the causally-constrained approach is that the solution always causal, in contrast to the least-squares methods in standard use based on the RWE, Thevenin or Norton circuits. Errors can remain due to inaccuracies in any least-squares approximate solution and

measurement errors in the recordings.

For the causality-constrained known-po calibration, the calibration outputs to be determined are the set of ·¾ ] = ^[ ] over the analysis band of frequencies. The hatted variables using the above procedure with Eq. (7) for tube m is,

such that the real and imaginary parts of -¾[ ] are -¾[ ] and -¾[ ] , respectively, the real and imaginary parts of A J\ are

and ^[^] , respectively, and the real and imaginary parts of b_m\_k\ are b_mR[k\ and _/[^], respectively. The solution is given by Eq. (14) for use in the causality-constrained procedure. For the causality-constrained unknown-po calibration, the solution for tube m in the causality-constrained procedure is expressed in terms of,

such that the real and imaginary parts of P₀[k are P_0R\k\ and P_QI\k\ . For the causality- constrained unknown-po calibration, the unknown vector x[k ] in Eq. (16) contains _R \k and ^_M^k it is convenient to normalize all acoustic pressures by the peak amplitude /7_llK of the measured pressure waveforms in all waveguides at all sample times. The variable k’_t/[k] can be eliminated from the system equations, and the remaining three variables are related by for waveguide m by,

with coefficients,

¾W = (1 + R_mR[k])² + R_m ² _I[k],

r_mm = ( R_mR[k]P_mI[k] - R_mI[k]P_mR[k] - (R_m ² _R[k] + R_m ² _I[k])P_mI[k]) / P_max,

L\k\ = {( ⁺ A^)^rA^⁺ l^^rA^ ΐ p_max·

(The variable can alternatively be eliminated from the system equations to obtain an equivalent set of relations to the above.) Each row fit of a M _x 2 matrix im-_| is denoted by

A_m[k], which is defined below along with row ffl of column vector _{h | k} \ :

Given data recorded in M_t tubes, there are M_t =M_t \/[(M_t— 2)!2!] pair-wise equations of the above form, or six pairs in the example with

=4. The matrix A_[k] and vector t\_{k ]} are constructed as in the above in terms of the set of A k] and b_m[k] across all tubes. This is the desired relation in terms of two unknown real parts and -¾*[ ]. Although this im-_| has more than M_t rows with additional dependencies of terms across rows, the least-squares approach can still be used.

The acoustical Thevenin circuit includes a Thevenin pressure P_T\k\ ⁼ P_TM\ ⁺ jP_jM\ in series with a Thevenin impedance Z_T\k\ =R_r[k + jX_T\k to which is connected in series an input impedance Z_m[k\ = R_m[k\+ jX_m[k\ for waveguide m, in which the pressure spectrum

P_m[P is measured during calibration. The input impedance waveguide is the calibration waveguide model response function for the Thevenin circuit. Note in this section only that upper-case R is used to denote frequency-domain acoustic resistance and lower-case r to denote time-domain resistance, rather than their uses elsewhere for pressure reflectance and RF, respectively. Upper- and lower-case X denote acoustic reactance in frequency and time domains, respectively. The pressure divider equation expresses the relationship between the Thevenin circuit parameters and the pressure and impedance in waveguide m.

PM _ z_m[k]

P Z +Z

The standard acoustical Thevenin calibration procedure uses data acquired in more than two tubes (C > 2) to construct an overdetermined matrix equation to solve for the unknown transposed vector _{c ]} ^G _{= [}p_{T [}k _] z _{T [}k _{] ]} in the least-squares approximation. The model parameters are adjusted by minimizing a residual. Once P_T and Z_r are calculated, the pressure divider equation is again used to evaluate the impedance Z of an unknown ear or other waveguide in terms of its measured pressure and Thevenin parameters.

A causally-constrained Thevenin calibration uses an over-determined matrix equation of the general form of Eq. (13), in which the unknown vector x[k ] at the k-th frequency contains the real parts of the Thevenin pressure PJJ\ and Thevenin impedance Z_IR\k\. The pressure divider equation for tube M is expressed in terms of the real and imaginary parts of all quantities. The following coefficients are defined in terms of the waveguide data and model impedance and takes the form (the index [ k ] is omitted in intermediate calculations), a_, =(^ + )/(¾ ),

The unknown column vector at the k-th frequency is defined by,

Row ftt of the _{M ' x} 2 matrix im- 1 is denoted by A_m[k], which is defined below along with row fit of column vector b by data from tubes i and /:

K ⁼ 4 ^_f -4_f ^y _e ·

Each of the M _t rows of A_[k] and _{h |} /- 1 contain one of the pair-wise equations for i from 1 to

M_t, and over all j < i. This completes construction of the over-determined matrix equation in

Eq. (13) for the causal Thevenin calibration with the final measurement system parameters given in the least-squares solution in Eq. (14). A causally-constrained Norton circuit calibration is obtained in a similar manner to the RWE and Thevenin circuit examples based on the well- known relation between Norton and Thevenin circuit parameters.

Causal Calibration in the Time Domain

Measurements can be performed using a discrete-time implementation with sample period T and sample rate HT. Discrete-time waveforms are constructed by sampling the waveform r(ί_h) at times t_n =nT with taking on values between 0 and N- 1. The discrete time waveform is p[n\ at sample n, and similarly for other signals and transfer functions in time.

A discrete-time convolution b = x * y of transient causal signals x[n] and y[n\ can be expressed in matrix form:

b = Ax. (17) The components in row n of column vectors b and x are b[n ] and x[n ] , respectively, at sample n and sample period T. The matrix A has coefficients

in row m and column n£m given by,

A_m =y{m-n\ . (18)

The remaining terms of A_mn = 0 for ΪUKUI. The resulting A is a lower triangular Toeplitz matrix. All convolutions of causal variables are represented in matrix form by a lower triangular Toeplitz matrix. The transpose w ^T of any column vector W of length A is a row vector, and the product of w ^Tw is

which is a non-negative scalar. The Euclidean norm || w | | of W is w⁷w. Other norms can be used here and elsewhere.

The simplest discretization procedure directly samples the model RF at the audio sample rate so that [/?] =r(t_n). This procedure is inaccurate near the peak of the model RF just after time ¾ in Eq. (2). A more accurate discretization procedure is to sample the model RF at a rate Q > 1 times faster than the audio sample rate \/T. Upsampling by Q > 1 is computationally expensive as deconvolutions would be calculated at an increased sample rate relative to the audio sample rate. This would require up-sampling of the measured data by a factor of Q , and down- sampling the output of each deconvolution back to the audio sample rate.

A third procedure, the decimation approach, initially discretizes the RF at an interpolation factor Q > 1 to better approximate the continuous-time RF. Next, the RF is filtered with a linear- phase digital lowpass filter to the maximum analysis frequency, and then decimated, or down- sampled, back to the audio sample rate. The resulting model RF is band-limited to the audio bandwidth so that the deconvolution was calculated at the audio sample rate. A Kaiser window can be used with a linear-phase FIR filter using a least-squares error minimization, because it has a more narrow frequency range between pass and stop bands than do Blackman, Hamming, Hanning and rectangular windows, when the filters are compared in the condition of having the same stop-band level. The impulse response length of the FIR filter is relatively short, e.g., 41 samples (0.85 ms) in Keefe, 2020 to limit the duration of the audio-rate model RF. This RF occurs in the deconvolution described below only as a convolution with the tube pressure p, i.e., as r*p, in which the measured tube pressure is sampled at the audio sample rate. A short- duration FIR filter window ensures that the convolution of the RF with the incident-pressure part of the click remains well-separated from time zero or the convolution of the RF with the first- reflected click. A fourth procedure is to represent the RF using a waveguide filter.

A lossy waveguide (LW) in a fifth approach is the preferred embodiment to calculate a discrete-time model RF for a closed cylindrical tube of any length such that the resulting model RF is causal. All calculations are performed in the time domain. This solves a problem described in Keefe, 2020 that calculated a discrete-time model RF with the decimation approach for tubes having a length of at least 90 cm. The modeling approach was deemed inadequate for much shorter tubes, for which the calculated model did not satisfy causality. This limited the ability to measure a causal RF of the ear in the time domain based on a calibration using only long tubes. The LW approach calculated a discretized model RF of a closed tube that is based on a smoothed model RF in continuous time. The resulting discretized model RF can be used for a closed tube of any length.

A finite-duration, quasi-impulsive, transient d_Q(t) is defined as equal to 1/7" for the audio-rate sample period T centered at time zero, and zero for all other times. Its integral over a single-period sample centered at time zero is one. This transient is,

in which the Heaviside step function s(t) is one for positives times and zero for negative times. As explained following A. Papoulis,“The Fourier Integral and Its Applications” (McGraw-Hill, New York, 1962), its FT D₀(a> ) = sin( / 2) / (col· / 2) is attenuated a higher frequencies although it is not band-limited according to the Nyquist sampling theorem for the audio sample rate 1/7". A smoothed continuous-time RF r(t) = ( d₀(t ) * r(t))(t) is defined as the convolution of i ₀ (0 with r(t ) in Eq. (2), and is evaluated at sample times t_n =nT as r(t_n) = r\n\ = r[n\/ T.

This discrete-time RF is calculated in terms of the error function erf (x) as,

with x_T = A / jT and e[h ] equal to one for n > 0, zero for n < 0, and 1/2 for n = 0. A practical approximation for n > /^ + ^Tl^^L is,

The system equations involving r[n ] always occur with a convolution with a band-limited signal, as in the discrete-time form of Eq. (2). The output of the convolution of r[n ] with a band- limited signal is band-limited, even though the smoothed continuous-time f (/) is not band- limited. The choice of smoothing function d_Q(t) is preferred, but other smoothing functions can be used.

For example, model RF responses r[n ] are shown in the left panel of FIG. 3 for a cylindrical tube of radius a= 4 mm and a length of 90.4 cm. This corresponds to the tube model for the shortest tube 4 used in Keefe, 2020. Calculation procedures from this article were used except for a more accurate interpolation factor of Q = 15 to define the high sample rate. The LW model was more accurate than the decimation model.

The LW model RF in Eq. (20) has no restriction on tube length. An example is shown in the right panel of FIG. 3 for a short tube of length 2 cm with radius of 4 mm on a logarithmic amplitude scale (i.e., 201 Qg₁₀ |/?]), because of the large peak amplitude (the early components with zero values are not displayed). The peak RF value from the continuous-time solution was 50 dB larger than the peak of the LW RF at the audio rate, and the time of the peak for the continuous-time solution was slightly larger than the round-trip time T_2L = 0.116 ms and intermediate between adjacent sample times.

In a known-po calibration procedure in the time domain, the incident pressure p₀ is assumed known so that only the source RF ri, is calculated. Using tube 1 as an example, the RWE is expressed as a convolution product with unknown input ·¾[/?] = r₀[/?] convolved with a known <¾[«] to output a known b^ri\ ^'.

a_l *x=b_l , with

The N x 1 column vector _x holds the N waveform values of x[n ] with x[0] in row 1 down to x[N - 1] in row N. The Mx 1 column vector b, holds the M waveform values of bpri] with b ] in row 1 down to b^M—X in row M . The M x N matrix A is defined with coefficients

_{Ai mn = a}^_{m - n}\ in row M and column fl. The resulting matrix is lower triangular Toeplitz due to causality. The above relations are expressed as the matrix equation,

x ⁼b r (22)

A similar relation holds for each of the N_t tubes. A (M_f -M) 1 column vector b is defined by vertically stacking the individual column vectors b into a single column vector. That is, rows ( m - \)M ₊ 1 to rrM of b hold the components of b_m for all M from tubes 1 to M_t. The (M_t · M) x N matrix A is defined by vertically stacking the individual matrices A_m into a single matrix having the same number of columns. That is, rows ( m - 1 )M ₊ 1 to mM of A hold the components of matrix A_m for all M from tubes 1 to N_t. The resulting matrix equation for the over-determined system to measure x is,

Ax = b. (23)

Equation (23) can be solved for 2C = r₀ using regularized deconvolution.

In the unknown-po calibration procedure in the time domain, the discrete-time RWE for tube 1 is expressed as a matrix convolution product as follows:

By analogy with the known-po approach, the unknown-po discrete-time approach in Eq. (24) also takes the matrix forms of Eqs. (22) and (23). The resulting Eq. (23) can be solved for the unknown RWE outputs p₀ and contained in vector fusing a regularized deconvolution.

Using standard transformations, the known Cl_x and b , and the unknown L, in Eq. (24) in terms of the RWE can be expressed in terms of the Thevenin impulse response z_T , whose FT is the Thevenin impedance Z^’, _T, and Thevenin pressure p_T, whose FT is the Thevenin pressure

P_T , in place of p₀ and T₀. The RWE terms can similarly be replaced by time-domain forms of the Norton variables. The resulting equations that are obtained can be written in the matrix form of Eq. (23) and can be solved using regularized deconvolution.

The (M_t ·AR)^cN matrix A in the general-form Eq. (23) in the time domain has no inverse because (M_t’M) > N, so that an appropriate goal is to calculate the least-squares solution vector _x that minimizes the norm of the residual vector || Ax - b || . This corresponds to deconvolution in the time domain. An overall real non-negative error function E) over all tube variables and times is defined as the squared norm of the residual vector

This is the error function that is minimized. Numerically stable deconvolutions can be obtained using iterative or direct regularization techniques. Regularized deconvolution procedures

One deconvolution procedure is based on the singular value deconvolution (SVD) of matrix A. The SVD of a M x N matrix A with M > N is,

A=UåV (26) in which U and V are orthogonal matrices of size MxM and N x N, respectively, and superscript T represents the transpose operation. Columns of U are eigenvectors of AA^T and columns of V are eigenvectors of A' A The å is a M x N matrix that is diagonal in its first N rows such that ¾ =cr in row i and zero in the remaining rows. The r <N non-zero singular values cr (with 1 < i < r ) are square roots of the eigenvalues of both AA^T and A^TA . Singular values are ordered in a non-increasing manner, i.e., <J_l £ <J₂ £...£cr . In numerical deconvolutions, typically r = N with very small values of cr _at large i on the order of the floating-point precision or other types of noise.

For singular value cr the eigenvector of U is U and the eigenvector of V is V . The SVD takes the form,

Terms with large singular values (e.g., <J_X) make a large contribution to the SVD sum of A, whereas terms with small singular values (e.g., cr ) make a small contribution, i.e., the presence of system noise has only a small perturbative effect on calculating A.

The least-squares solution minimizes the norm of the error of the residual ||z4x -b|| in terms of the pseudoinverse of A. Assuming that all singular values are non-zero, this solution is,

This relation is unduly sensitive to small singular values because cr occurs in the denominator.

Small levels of system noise produce large effects in the solution to the inverse problem; see, e.g., P. C. Hansen,“Rank-Deficient and Discrete Ill-Posed Problems: Numerical Aspects of Linear Inversion.” Society for Industrial and Applied Mathematics, Philadelphia, 1998.

(hereinafter“Hansen, 1998”).

The matrix A may be contaminated by noise. Regularization refers to procedures used to calculate solutions with a reduced sensitivity to such errors. SVD regularization procedures omit or attenuate terms in Eq. (28) that are associated with smaller singular values. The presence of smaller singular values dominated by noise on the solution otherwise leads to large-amplitude fluctuations in its solution norm ||x|| .

Tikhonov regularization minimizes a weighted sum of the residual and solution norms with a real, non-negative weighting factor l termed the regularization parameter. The Tikhonov procedure minimizes the following expression to calculate an approximate solution ||x|| based on a pre-selected value of l :

The resulting solution vector with lowpass filter function f\ⁱ i is,

The filter function attenuates terms with smaller singular values (J <

The generalized SVD (GSVD) is a refinement of the Tikhonov procedure to minimize,

A general preconditioner matrix L has ( rows and N columns with N_L £N. It is selected to smooth the solution vector such that ||AX|| is better behaved numerically than ||x||. Two commonly used preconditioners can be used:

with N_L = N— 1 is a finite-difference approximation to the first time derivative, and

with N_L = N— 2 approximates the second derivative. Each is a Toeplitz matrix defined by values in the first column and first row. The first column vector of both preconditioners has 1 as the first component and 0 otherwise. Row 1 of matrix is expressed in component form as L^{ j ) for its value in row 1 and column j. The vector holds all zeros except:

A(U)=l· A( )=-i^. (32)

The row 1 vector in matrix

holds all zeros except: 4(1,1)=1, 4( l2)=-2. 4(1,3) =1. (33)

Singular values are expected to decay faster for the GSVD solution than for the Tikhonov SVD solution, thus achieving faster convergence towards the best solution.

The best choice of L differs according to the structure of A and b in an individual problem, so that one tries several preconditioners, or none, to determine which gives the best solution, e.g., one in which the unknown(s) in ||x|| are smooth functions of time with the smallest residual error. A trade-off exists in regularization between minimizing the residual norm

|z4x-b| while also controlling for a solution norm ||x|| or || X|| that is not too large, which would indicate the presence of large-amplitude fluctuations in the solution.

This trade-off is represented by a log-log plot of solution norm versus the residual norm. Each pair of norms varies with l in the case of Tikhonov regularization. Such a plot is called a L-curve; see. e.g., P. C. Hansen,“The use of the L-curve in the regularization of discrete ill- posed problems.” SIAM J. Sci Comput., 14:1487-1503, 1993. A generic L-curve has two sub curves. One part initially grows from zero to become approximately horizontal for a range of larger l , while the other part is approximately vertical for a range of smaller l such that more terms are included in the sum of Eq. (30) with noise-dominated singular values. The corner or “elbow”' of the L-curve, which corresponds to the point on the L-curve with maximum intrinsic curvature, can be selected to be the best solution as it has achieved both a relatively small residual norm without a substantial increase in its solution norm. In the case that the solution is exact, there is no unusual growth in the solution norm so no elbow is observed, but only a curve that remains approximately horizontal.

The calculation of a SVD is slow for large matrices, especially for large N. The conjugate gradient (CG) technique is an iterative deconvolution method that can provide high- quality results with substantially fewer floating point operations than SVD in acoustic waveguide problems to calculate RFs. A CG least squares (CGLS) algorithm is well-suited to solve deconvolution problems by iteratively solving the normal equation A^TAx = A^Tb for unknown _x (Hansen, 1998). Here, A^TA is a square N x N matrix. Performing A iterations of CGLS results in the SVD solution (except for floating point errors) and is the maximum number of iterations in any CG algorithm. Key properties are that successive iterations of CGLS regularize, or smooth, the solution, and relatively few iterations may be needed to calculate the solution _x compared to SVD.

In a preconditioned CGLS (PCGLS) method, a preconditioner matrix L such as defined in Eqs. (32)-(33) can be added to further smooth the solution. The PCGLS method outputs a new estimate of a solution vector _x, from which the residual norm | d x - b | and preconditioned solution norm ||LX|| are calculated. The L-curve of solution norm versus residual norm can be calculated as a function of iteration number, and the regularized solution _x is obtained at the elbow.

A regularized deconvolution procedure refers to any of the above procedures along with other well-known direct and iterative deconvolution procedures with optional regularization and/or preconditioner stages. A regularized deconvolution procedure also includes use of a criterion to the choose the best solution among the trial solutions labelled in a direct procedure by a parameter such as l for the SVD and GSVD procedures or labelled in an iterative procedure by the iteration number such as used in the CG, CGLS and PCGLS procedures. The choice of this criterion can be based on the value on a L-curve having the maximum intrinsic curvature. Other criteria can be used, e.g., those described in Hansen (1998).

An procedure for determining measurement system parameters in the time domain (e.g., p_Q and f₀ ) that satisfy causality is shown according to an embodiment in the flowchart of FIG.

5. Calibration signals 600 are received from the microphone of a probe when the probe assembly is inserted into each of one or more acoustic calibration waveguides and the acoustic source generates sound responsive to the electrical input signal. Each calibration signal is causal. At the start of the procedure, the waveguide model parameters 612 used by the calibration waveguide model response functions are initialized to plausible estimates. As in the frequency-domain optimization, initial acoustical estimates of tube lengths in the set of calibration waveguide model parameters 612 can be calculated from the measured tube data using a short-time autocorrelation procedure.

The procedure enters an outer optimization loop which reduces the overall error in estimating the measurement system parameters by varying each waveguide model parameter in the parameter set such that the residual norm of a second over-determined matrix equation 605 is minimized over all waveguide model parameter values. This is an iterative loop calculation, as indicated by blocks 602-607 in FIG. 5. The calibration waveguide model response function 602 is calculated for each calibration waveguide as a function of the set of waveguide model parameters 612 or 606. The first over-determined matrix equation has as known inputs the set of calibration signals 600 and the set of calibration waveguide model response functions 602. An approximate solution _x of the first over-determined matrix equation 603 contains the

measurement system parameters 604.

An example of a first over-determined matrix equation 603 takes the form of Eq. (23), which can be solved using a regularized deconvolution. This regularization can adopt a stopping criterion based on a L-curve or other procedure to output a best solution _x with the measurement system parameters 604. The notation of Eq. (23) is extended by adding a superscript to each matrix variable for the iteration number in the outer loop, so that A^(J) represents the system matrix A in iteration j, vector b ^(;) represents the vector b in iteration j, and vector x ⁾ represents the least-squares approximate solution _x of the first over-determined matrix equation 603 in iteration j. This least-squares approximation solution does not satisfy ri⁰⁾x^(;) = b^(;) because of measurement errors and modeling errors, which include errors in the set of calibration waveguide model response functions 602. The purpose of the outer loop is to reduce the overall least-squares approximation error by improving the accuracy of the set of calibration waveguide model response functions 602. This accuracy is improved by the choice and value of each waveguide model parameter used in the respective model response function.

A set of waveguide model parameters is represented as an abstract vector variable Q, in which the m-th parameter q_hi of this vector includes the m-th tube length for 1 £ m <M_t and can include additional parameters 0_i for i > M_t holding any additional viscothermal parameter or evanescent-mode parameter that might further improve model accuracy. These are the waveguide model parameters that are updated in box 606 after completing each iteration in the outer loop. The residual error vector in iteration j of the outer loop is ri⁰⁾x⁰⁾ - b^(;) . This residual error comprises a second over-determined matrix equation that is a function of the set of waveguide model parameters q^w. These are varied in each iteration of the loop so that A¹ ' "¹ in iteration j+1 differs from A ¹ ' ¹ in iteration j, and b°⁺¹⁾ in iteration j+1 differs from b ^(;) in iteration j. It follows that the least-squares solution x°⁺¹⁾ from solving the first over-determined matrix equation in iteration j+1 differs from the solution x ⁾ in iteration j. The residual norm II ^0⁾ _cϋ⁾— 1}^{( /)} II

_seconci over-determined matrix equation is minimized in the outer loop by varying the waveguide model parameters so as to generate a more accurate solution as the iteration number j increases from 1 to its final value.

The outer optimization loop in boxes 602 through 607 can use a constrained nonlinear multivariable optimization based on the active-set algorithm. Inequality constraints can be adopted to bracket the admissible range for each waveguide model parameter value in each iteration. The optimization includes a criterion for stopping at the numerically optimum solution based on calculated results in all iterations. After the residual norm is minimized (block 607 returns‘yes’), the final measurement system parameters 608 are returned as the values of the updated measurement system parameters 604 calculated on the final iteration of the calculations in boxes 604-605. Causality of all signals is maintained throughout.

Causality lends an additional constraint to limit the space of possible errors in a least- squares approximation. Limiting the space of errors is expected to increase the accuracy of the reflectance and RF calibrations which produce causal outputs for the reflectance and RF measured in an ear canal.

After calibration, the system is used to receive a measurement signal 609 from the microphone when the probe assembly is inserted in a leak-free manner into the ear canal and the acoustic source generates sound responsive to the electrical input signal 215. The causal response function 611 of the ear is determined in the time domain based on the measurement signal in the ear and the final one or more measurement system parameters 608. This determination is based on solving a matrix equation of the general form of Eq. (23) with modifications in the system matrix and known and unknown vectors. The method of solution can use a direct or iterative regularized deconvolution with a L-curve stopping criterion.

Evanescent-Mode Effects in Causally-Constrained Calibrations

The front surface of a typical ear-canal probe has a narrow orifice through which an acoustic source delivers sound, and a small duct through which sound pressure acting over its surface is transduced by a microphone. This section describes additional effects of the localized sound field generated by the effects of the sound-source and microphone apertures on the surface of the probe. The waveguide into which the probe is inserted has area S = ncr . Sound exiting the source spreads out to form a one-dimensional uniparametric wave that propagates along the longitudinal waveguide axis. The localized sound field can be expressed as acting at low frequencies as an additional inertive impedance in series with the input impedance of the tube.

The volume velocity delivered by the probe equals the volume velocity through the inertance and into the waveguide. For a short, narrow tube of air of length and area S , the inertance I e

with the latter normalized by the characteristic impedance. The frequency-domain evanescent mode impedance Z_e{cd) =jCol_e =2jCOT_IZ_{c ,} with characteristic time 7} = /, / 2c. For example, in the publication by K. R. Norgaard, E. Fernandez-Grande, and S. Laugesen,“ Incorporating evanescent modes and flow losses into reference impedances in acoustic Thevenin calibration,”

J. Acoust. Soc. Am., 142:3013-3024, 2017 (hereinafter Norgaard 2017), an inertance of I_e of 37 kg/m⁴ was calculated for a probe inserted into a cylindrical tube of radius a = 4 mm. Their resulting l_e = 1.6 mm with l_e / <2 = 0.40. The impedance Z_e(⁽D) is in series with the input impedance Z{w) of the tube into which the probe is inserted. The total input impedance of the tube including evanescent modes is Z(w) = Z(w) +2ZJCOT_i. The magnitude of the evanescent- mode impedance is much less than that of the tube impedance except in the vicinity of the minima in the tube impedance magnitude. The corresponding pressure reflectance k _{(w )} is then calculated in terms of z (_{® )} as,

This is an example of a calibration waveguide model response function 400 in the frequency domain that includes an evanescent-mode contribution. Adding one or more evanescent-mode parameters to a calibration waveguide model response function may lead to more accurate determinations of the final measurement system parameters 401, and especially so for a relatively short cylindrical tube as the calibration waveguide (Keefe, 2020). The^ _{(® )} includes any evanescent-mode contribution that can be used to determine to determine the final measurement system parameters in the frequency domain that satisfy causality. For an anechoic tube with R(co) = 0, this relation simplifies to,

With _{t >} o_, the corresponding causal, time-domain solution r(t) to R _{(® )} in Eq. (36) is obtained by integrating the acoustical differential equations as,

That this is the solution is confirmed by the fact that it satisfies causality and its FT is R _{(® )} in Eq. (36). This time-domain solution avoids use of the inverse Laplace transform, although its use leads to the same result. The anechoic-tube Eq. (37) describes the model RF seen by the outgoing incident pressure signal prior to any reflection received from the opposite closed end of the tube.

In Norgaard 2017, a negative inertance corresponding to 7} <0 was measured for an ear- canal probe inserted into tubes of radius a = 2 (the corresponding I_e =— 101.5kg/m⁴), and confirmed its negative value using a computational model. The resulting l_e was about -1.1 mm so that IJ Cl was about -0.53. However, a negative inertance causes the RF in Eq. (37) to diverge exponentially. Such behavior is unphysical. This conclusion is not ameliorated by using other correction terms added to describe evanescent-mode effects.

In K. R. Norgaard, S. T. Neely, and D. M. Rasetshwane,“Quantifying undesired parallel components in Thevenin-equivalent acoustic source parameters,” J. Acoust. Soc. Am.,

143: 1491-1503, 2018, an alternative circuit model was described for the negative junction inertance I _e in series with the impedance of the closed tube of length L. In the absence of viscothermal losses, a short length l_e of the transmission line has a series impedance I_e = pl_e l S and shunt compliance c_e = si_{e / pc} ² For l_e >0 and with a reasonably small value (i.e., l_e ! a < \ ), the combination of the series and shunt impedance is an excellent approximation to a short segment of the acoustic transmission line of length l_e. For the case of a negative junction inertance l_e =-|/,|, the equivalent circuit with the tube impedance in series with the negative inertance is augmented by a pair of shunt compliances in parallel that leave the overall circuit response unchanged, one with a negative compliance si _{p c} ¹ and the other with a positive compliance C_e of,

The combination of negative series inertance and negative shunt compliance terminated by the impedance of the tube of length L is replaced by an equivalent transmission line for a slightly shorter tube of length L— |/ . The remaining positive shunt compliance acts in parallel with this modified transmission line. The tube length is much larger than the magnitude of the negative junction length. With reference to Eq. (2), the value of A in the transmission line of shorter length is adjusted to a_vt ( L - 1/,|) / a , and ¾ is replaced by T_2L = I(L - |/,|) / 2. Otherwise, this model continues to be denoted by r(t) with R( ) in the subsequent development.

It is convenient to define a characteristic time T_C =|/,| / 2c. In the frequency domain, the reflectance R _{(® )} corresponding to this parallel combination of the junction compliance and the reflectance R(co) of the tube is (Keefe, 2020),

This is an example of a calibration waveguide model response function 400 in the frequency domain that includes an evanescent-mode contribution. This relation simplifies for an anechoic tube with R(co) = 0 to,

With _{t >} o_, the corresponding causal, time-domain solution f (t) to

obtained by integrating the acoustical differential equations as, r(t) = -5(t) +— e^~t/Tcs(t).

T_r That this is the solution is confirmed by the fact that it satisfies causality and its FT is _{(w )} in Eq. (40). In general, T_c =—T_j , so the strategy is to adopt the RF model having the non-negative time delay parameter. The anechoic-tube Eq. (37) describes the model RF seen by the outgoing incident pressure signal prior to any reflection received from the opposite closed end of the tube.

This example uses the reflectance k _{(w )} in Eq. (35) or (39) as the calibration waveguide model response function in the frequency domain that includes an evanescent-mode contribution depending on whether either T or T_c, respectively, is positive for a particular choices of the probe and waveguide dimensions. This example uses the RF f(t) in Eq. (37) or (41) as the calibration waveguide model response function in the time domain that includes an evanescent mode contribution in the incident pressure, depending on whether either T or T_c, respectively, is positive for a particular choices of the probe and waveguide dimensions. The corresponding RF r(t) at times after the first reflection returns to the probe tip can be obtained by numerical integration of the appropriate choice of Eq. (37) or (41) in the acoustical field equations with the corresponding reflection characteristics embodied in r(t) for the calibration waveguide model response function in the time domain in the absence of evanescent-mode effects.

Ear Evaluation Procedures

After calibration is completed so that P₀[k and

are known, the frequency-domain pressure reflectance R[k] of an ear or other waveguide with measured pressure spectrum P[k] in response to the same DAC stimulus used in calibration is evaluated using the RWE as an example as,

The area of the ear (or other waveguide) generally differs from that of the calibration tubes, while the reflectance is initially evaluated as having the same area as the calibration tubes. An ear-adjusted reflectance

of the ear can be calculated in terms of R[k] and a knowledge of the ear-canal area S_e at the probe tip. The measured ear pressure P[k ] is causal because it is the DFT of a causal waveform. Suppose a causality-constrained calibration is used to calculate P₀\_k\ and k\ at all analysis frequencies. Then, all variables on the right-hand side of Eq. (42)

, so that the ear response functions R[k], and P_L\_k\ derived from R\k |. are also causal. This procedure determines a causal response function of the ear in the frequency domain, which in this example is either of R[k] or R R\ based on the detected electrical signal from the microphone in the ear that measures P[k] and the final measurement system parameters P₀[k\

. Other causal response functions of the ear in the frequency domain are the acoustic impedance or acoustic admittance, which can be calculated from the reflectance R[k] of the ear and the area of the calibration waveguide at the probe tip, or equivalently from the ear-adjusted reflectance R R\ of the ear and the ear-canal area. For example, the relation between the acoustic impedance and R[k] is given in Eq. (4), and acoustic admittance is the inverse of acoustic impedance. These other response functions of the ear are causal because either of R[k] or causal. In summary, the use of a causally-constrained calibration results in evaluations of reflectance, impedance and admittance in an ear that are also causal.

In the time domain, the RF of the ear r[n ] defined below is an example of a causal response function of the ear in the time domain. This r[n ] can be calculated in the time domain based on the measured pressure p[n\ at the probe tip in the ear, which is in the form of a detected measurement electrical signal in the ear at the output of the microphone contained in the probe assembly, and the final measurement system parameters

^ar|d Ά in the time domain that are causal signals.

The RWE in Eq. (1) is re-written with suppressed time-dependence as,

b = a * r,

a = Po + r₀ * P (43)

P max

b P - Po

F’max

Equation (43) may be expressed in discrete time by the following over-determined matrix equation,

Ax = b (44) The lower-triangular Toeplitz matrix A has coefficients A_tm =a[m—n\. Because r[n] is zero for negative times, it is evaluated at discrete time steps from 0 £n £ (N_e—1) in the column vector r . The N_e is sufficiently large that r[n] is of negligible amplitude at times later than L{7( with the exception that any small otoacoustic emission components generated within the cochlea are neglected. Deconvolutions of Eq. (44) to solve for _r can be performed based on techniques such as SVD or CG with and without regularization and with or without pre-conditioners. The preferred embodiment to calculate the RF of the ear uses the PCGLS algorithm with

pre conditioner and L-curve criterion. This contrasts with the preferred embodiment to calculate the final measurement system parameters, in which p₀[n is calculated using the separability of incident pressure in a long cylindrical tube, and ft is calculated using the PCGLS algorithm with I pre-conditioner and L-curve criterion. Additional details are in Keefe 2020.

This invention to determine a causal response function of the ear in the time domain differs from the previous art by the causality constraint. This is facilitated by using a calibration waveguide model response function in the time domain that satisfies causality.

The causal RF r[n ] in the rows of the solution vector r in Eq. (44) is based on an evaluation waveguide with known area S₀ equal to that of the calibration waveguides at the probe tip. Further processing is needed to calculate the causal RF t\ of the ear (or other waveguide with unknown area) with known area S_e differing from

It can be appreciated that the solution technique for ft\ is not unique. An example of a RWE solution technique is described based on the RWE and the linearized equations of acoustics in an ear canal, namely, the RF r_e[ri of the ear satisfies,

These equations are re-written in matrix form as,

^Ae^re = K (46)

The lower-triangular Toeplitz matrix A_e has coefficients A _{e m n} = a _e [m - « ] . Because f,[/7] is zero for negative times, it is evaluated at discrete time steps from 0 £n £ (N_e—1) in the column vector r . The N_e is sufficiently large that

is of negligible amplitude at times later than N_cj with the exception that any small otoacoustic emission components generated within the cochlea are neglected. Deconvolutions of Eq. (46) to solve for r_e can be performed based on techniques such as SVD or CG with and without regularization and with or without pre-conditioners.

Tympanometric Causal Response Functions of the Ear

When the main signal processor 210 in FIG. 1 transmits an air-pressure signal 205 to a pump control unit 204, and the pump control unit 204 transmits a driver signal 207 to the electrical input of the pressure pump 208, then the resulting acoustic response function of the ear is measured after inserting the probe assembly 225 in the ear canal in a leak-free manner. The ear response 220 is measured in response to the sound stimuli 219 while the pump either maintains a fixed air pressure differing from atmospheric pressure or varies the air pressure as a function of time with respect to time. The resulting measured by system 200 provides an acoustic

tympanometric response function of the ear that varies with the air pressure produced by the action of the pressure pump and varies with either the time or frequency of the acoustic stimulus. The air pressure within the ear canal is measured by the pressure transducer 201, which transmits the pressure value to the main signal processor. The main signal processor 210 adjusts air pressure within the ear canal via the pump controller 204, pressure pump 208, and pressure transducer 201 obtain the desired air pressure during the test, whether that desired air pressure is constant with time, a sweep between initial and final air pressure values, or some other time variation. An acoustic tympanometric response function of the ear is an acoustic tympanometric causal response function of the ear if the transfer function, which includes reflectance, reflection function, impedance and admittance, is causal. The resulting acoustic tympanometric causal response function of the ear can be analyzed in the time domain 611 or the frequency domain 401 to provide clinically useful information to screen for, or diagnose, middle-ear dysfunction in patients across the age range from newborns to older adults.

All headings provided herein are for the convenience of the reader and should not be used to limit the meaning of any text that follows the heading, unless so specified.

The terms“comprises” and variations thereof do not have a limiting meaning where these terms appear in the description and claims. Such terms will be understood to imply the inclusion of a stated step or element or group of steps or elements but not the exclusion of any other step or element or group of steps or elements.

The words“preferred” and“preferably” refer to embodiments of the disclosure that may afford certain benefits, under certain circumstances; however, other embodiments may also be preferred, under the same or other circumstances. Furthermore, the recitation of one or more preferred embodiments does not imply that other embodiments are not useful, and is not intended to exclude other embodiments from the scope of the disclosure.

In this application, terms such as“a,”“an,” and“the” are not intended to refer to only a singular entity, but include the general class of which a specific example may be used for illustration. The terms“a,”“an,” and“the” are used interchangeably with the term“at least one.” The phrases“at least one of’ and“comprises at least one of’ followed by a list refers to any one of the items in the list and any combination of two or more items in the list.

The phrases“at least one of’ and“comprises at least one of’ followed by a list refers to any one of the items in the list and any combination of two or more items in the list.

As used herein, the term“or” is generally employed in its usual sense including“and/or” unless the content clearly dictates otherwise.

The term“and/or” means one or all of the listed elements or a combination of any two or more of the listed elements.

As used herein in connection with a measured quantity, the term“about” refers to that variation in the measured quantity as would be expected by the skilled artisan making the measurement and exercising a level of care commensurate with the objective of the measurement and the precision of the measuring equipment used. Herein,“up to” a number (e.g., up to 50) includes the number (e.g., 50).

The complete disclosure of the patents, patent documents, and publications cited herein are incorporated by reference in their entirety as if each were individually incorporated. In the event that any inconsistency exists between the disclosure of the present application and the disclosure(s) of any document incorporated herein by reference, the disclosure of the present application shall govern.

Claims

What is claimed is:

1. A system comprising:

a probe assembly adapted to be inserted into an ear canal of a patient, the probe assembly comprising an acoustic source and a microphone;

one or more acoustic calibration waveguides having known dimensions;

a signal generator coupled to deliver an electrical input signal to the acoustic source; a processor coupled to the probe assembly and the signal generator and adapted to:

receive calibration signals from the microphone when the probe assembly is inserted into each of the acoustic calibration waveguides and the acoustic source generates sound responsive to the electrical input signal; and

determine one or more calibration waveguide model response functions in the time domain that each satisfies causality for a corresponding acoustic calibration waveguide, each calibration waveguide model response function being indicative of an acoustic transfer characteristic with viscothermal effects that is a function of one or more waveguide model parameters, with each waveguide model parameter value initialized to a plausible estimate;

perform an outer loop calculation to determine values of final measurement system parameter signals that are causal with the processing on each iteration within the outer loop comprising:

determining values of causal measurement system parameter signals in the time domain by solving a first over-determined matrix equation, which includes functions of the calibration signals and the calibration waveguide model response functions, with each function calculated in terms of the one or more waveguide model parameter values; and

determining an update of the waveguide model parameter values by minimizing a residual norm of a second over-determined matrix equation, which includes functions of the calibration signals, the measurement system parameter values, and the waveguide model response functions, the outer loop calculation concluding when the residual norm is minimized; determine values of causal final measurement system parameters signals in the time domain as values of the measurement system parameters determined at a conclusion of the outer loop calculation;

receive a measurement signal from the microphone when the probe assembly is inserted into the ear canal and the acoustic source generates sound responsive to the electrical input signal; and

determine the causal response function of the ear in the time domain based on the measurement signal in the ear and the final measurement system parameters.

2. The system of claim 1, wherein the processor is further adapted to detect in at least one of the calibration signals an incident signal, wherein the incident signal is separable in time from a first reflected signal in the at least one calibration signal.

3. The system of any one of claims 1-2, wherein each of the acoustic calibration

waveguides comprises an internal air volume closed at its far end that is a cylindrical volume with waveguide model parameters of cylindrical length and radius.

4. The system of any one of claims 1-3, wherein each calibration waveguide model response function includes one or more waveguide model parameters to describe an evanescent mode contribution in the sound field adjacent to a tip of the probe assembly.

5. The system of any one of claims 1-4 further comprising:

an air pressure pump controlled by the processor; and

an air pressure transducer that transmits the value of a measured air pressure within the ear canal to the processor, wherein the processor is further adapted to:

adjust the air pressure within the ear canal based on a desired air pressure and the measured air pressure within the ear canal; and

at the desired air pressure, determine a tympanometric causal response function of the ear in the time domain based on the measurement electrical in the ear and the final measurement system parameters.

6. The system of any one of claims 1-5, wherein the measurement system parameters comprise the incident pressure and source reflectance as functions of time at a tip of the probe.

7. The system of any one of claims 1-6, wherein the second over-determined matrix equation further includes functions of at least one of viscothermal parameters and evanescent mode parameters of the acoustic calibration waveguides.

8. A system comprising:

a probe assembly adapted to be inserted into an ear canal of a patient, the assembly comprising: an acoustic source and a microphone;

one or more acoustic calibration waveguides having known dimensions;

a signal generator coupled to deliver an electrical input signal to the acoustic source; a processor coupled to the probe assembly and the signal generator and adapted to:

receive calibration signals from the microphone when the probe assembly is inserted into each of the acoustic calibration waveguides and the acoustic source generates sound responsive to the electrical input signal;

determine a calibration waveguide model response functions in the frequency domain for each of the acoustic calibration waveguides that are indicative of an acoustic transfer characteristic with viscothermal effects for acoustic transmission in each acoustic calibration waveguide, and such that each calibration waveguide model response function in the frequency domain is a function of one or more waveguide model parameters; determine final measurement system parameters in the frequency domain such that each satisfies causality by solving system equations that contain the calibration signals and the calibration waveguide model response function for each of the calibration signals;

receive a measurement signal in the ear from the microphone after the probe assembly is inserted into the ear canal and the acoustic source generates sound with input from the electrical input signal; and

determine a causal response function of the ear in the frequency domain based on the measurement signal in the ear and the final measurement system parameters.

9. The system of claim 8, wherein the signal processor is further adapted to receive the calibration signals and detect in at least one of the calibration signals an incident signal, wherein the incident signal is separable in time from a first reflected signal in the at least one calibration signal.

10. The system of any one of claims 7-8, wherein the processing circuitry is further adapted to determine a causal response function of the ear in the time domain using the determined causal response function of the ear in the frequency domain.

11. The system of any one of claims 7-9, wherein each calibration waveguide model response function in the frequency domain includes an evanescent-mode contribution adjacent to a tip of the probe assembly that is a function of the geometry of the probe assembly and the one or more waveguide model parameters.

12. The system of any one of claims 7-10 further comprising an air pressure pump and an air pressure transducer, which transmits the value of a measured air pressure within the ear canal to the processor, wherein the processor is further adapted to determine a tympanometric causal response function of the ear in the frequency domain based on the measurement electrical signal and the final measurement system parameters.

13. The system of claim 12, wherein a tympanometric causal response function of the ear in the time domain is determined from the tympanometric causal response function of the ear in the frequency domain.

14. A method comprising:

inserting a probe assembly into each of one or more acoustic calibration waveguides, wherein the probe assembly comprises an acoustic source and a microphone;

receiving calibration signals from the microphone when the probe assembly is inserted into each of the acoustic calibration waveguides and the acoustic source generates sound responsive to an electrical input signal; determining one or more calibration waveguide model response functions in the time domain that each satisfies causality for a corresponding acoustic calibration waveguide, each calibration waveguide model response function being indicative of an acoustic transfer characteristic with viscothermal effects that is a function of one or more waveguide model parameters, with each waveguide model parameter value initialized to a plausible estimate; performing an outer loop calculation to determine values of final measurement system parameter signals that are causal with the processing on each iteration within the outer loop comprising:

determining values of causal measurement system parameter signals in the time domain by solving a first over-determined matrix equation, which includes functions of the calibration signals and the calibration waveguide model response functions, with each function calculated in terms of the one or more waveguide model parameter values; and determining an update of the waveguide model parameter values by minimizing a residual norm of a second over-determined matrix equation, which includes functions of the calibration signals, the measurement system parameter values, and the waveguide model response functions, the outer loop calculation concluding when the residual norm is minimized;

determining values of causal final measurement system parameters signals in the time domain as values of the measurement system parameters determined at a conclusion of the outer loop calculation;

receiving a measurement signal from the microphone when the probe assembly is inserted into the ear canal and the acoustic source generates sound responsive to the electrical input signal; and

determining the causal response function of the ear in the time domain based on the measurement signal in the ear and the final measurement system parameters.

15. The method of claim 14, further comprising detecting in at least one of the calibration signals an incident signal, wherein the incident signal is separable in time from a first reflected signal in the at least one calibration signal.

16. The method of claims 13 or 14, wherein each of the acoustic calibration waveguides comprises an internal air volume closed at its far end that is a cylindrical volume with waveguide model parameters of cylindrical length and radius.

17. The method of any one of claims 13-15, wherein each calibration waveguide model response function includes an evanescent-mode contribution adjacent to a tip of the probe assembly that is a function of the geometry of the probe assembly and the one or more waveguide model parameters.

18. The method of any one of claims 13-16 further comprising:

adjusting an air pressure within the ear canal based on a desired air pressure and the measured air pressure within the ear canal; and

at the desired air pressure, determining a tympanometric causal response function of the ear in the time domain based on the measurement electrical in the ear and the final measurement system parameters.

19. A method, comprising:

inserting a probe assembly into each of one or more acoustic calibration waveguides, wherein the probe assembly comprises an acoustic source and a microphone;

receiving calibration signals from the microphone when the probe assembly is inserted into each of the acoustic calibration waveguides and the acoustic source generates sound responsive to an electrical input signal;

determining a calibration waveguide model response functions in the frequency domain for each of the acoustic calibration waveguides that are indicative of an acoustic transfer characteristic with viscothermal effects for acoustic transmission in each acoustic calibration waveguide, and such that each calibration waveguide model response function in the frequency domain is a function of one or more waveguide model parameters;

determining final measurement system parameters in the frequency domain such that each satisfies causality by solving system equations that contain the calibration signals and the calibration waveguide model response function for each of the calibration signals; receiving a measurement signal in the ear from the microphone after the probe assembly is inserted into the ear canal and the acoustic source generates sound with input from the electrical input signal; and

determining a causal response function of the ear in the frequency domain based on the measurement signal in the ear and the final measurement system parameters.