## Introduction

_{ea}(s) and T

_{ae}(s), show anti-reciprocal characteristics of the system. For example, under the DC status (ignoring acoustic components), the system’s electro-magnetic coupling coefficient is a constant, -T

_{ea}= T

_{ae}= -T

_{a}= B

_{0}l where B

_{0}and l are the DC magnetic field and length of wire, respectively. The T

_{ea}(s) and T

_{ae}(s) have the same magnitudes but opposite signs.

_{Z}=Z

_{e}Z

_{a}-T

_{ea}T

_{ae}. Note that if C=0, Z does not exist.

## Theoretical Background

#### Two-port network theory

#### Transfer function and equivalent forms

#### Differential equation form of H(s)

#### Rational polynomial fraction (pole-zero) form of H(s)

_{i}) and infinity (s=d

_{k}), are called the system’s zero and pole, respectively.

#### Partial fractions with residual form of H(s)

_{m}is pole and r

_{m}is residual; d and e are real values, whereas p

_{m}and r

_{m}can be either real or complex values. If e is zero, then the number of poles and zeros in the system is the same.

#### Stability of LTI system

## Pole-Zero Fitting Algorithm

#### Gustavsen’s method: rational approximation

_{i}, Q

_{i}, D, and E in a residual expansion of the form

_{s=jω}. A two-step process is performed in this fitting procedure, pole identification and residual identification. This procedure is accomplished by converting a nonlinear least squares problem to a linear least squares problem by introducing an unknown scaling function σ(s) and known poles (q

_{i}in Eq. 17).

_{i}and p

_{i}, d, and e terms must be found. Multiplying data F(ω) with the scaling function σ(s=jω) and evaluating over the many available frequency data points in F(ω), we can overdetermine this linear problem to estimate the unknown values of p

_{i}, b

_{i}, e, and d. The mathematical expression yields,

_{i}, b

_{i}, e, and d in Eq. 17, the equation can be redefined as a product form

_{i}in Eq. 19) become the poles of the model function

_{i}=Q

_{i}. Gustavsen and Semlyen [8] report that the residuals for

_{i}, D, and E) can be calculated from p

_{i}, d, e, and bi by using zi as the starting poles for

_{i}). Once this procedure is completed, we can approximate our data F(s) from the fitting function

_{i}=p

_{i}, D=d, and E=e (Eq. 16 and Eq. 17).

_{i}(Eq. 19) order changes, they converge to the true poles. Then the scaling function σ(s) will become 1, meaning that the calculated zeros z

_{i}of σ(s) have values arbitrarily close to the values of the starting poles, q

_{i}. This single process can iterate over the data’s frequency range until the fitting reaches the best fit. When it returns

_{i}>0 referring to an unstable pole in

_{z}. For example, N

_{z}is the same as the number of poles (N) when e is zero but d is not zero. If e and d are zero, we can have one fewer zeros than the poles because the numerator order is one less than the denominator (see Eq. 17). Similarly, when d and e are both nonzero, N

_{z}is N+1. In this study, we use the fitting condition with e zero and d nonzero, allowing us the same number of poles and zeros. We reach a slightly different fitting result by varying the initial value of e and d, but this does not critically affect the overall performance of the fitting result. Eventually, when this augmented equation is fitted to the given data F(ω) iterating toward σ(s)=1, the order of the poles for the fitting function meets the given error tolerance. The error is calculated using a mean squared error (MSE), in decibels, relative to the L2 norm of the signal:

#### Prony’s method

_{i}, σ

_{i}, φ

_{i}, and f

_{i}are amplitude, damping coefficient, phase, and frequency of component i, respectively. Also, L is the total number of the damped exponential components. This function estimates given data y(t), which consists of N samples y(t

_{k})=y[k], where k=0, 1, 2, ..., N-1 are evenly spaced.

_{i}t+φ

_{i}) term can be rewritten as a sum of exponentials:

_{i}and μ

_{i}in three steps: First, it solves the linear prediction model, which is constructed by the observed data set. Second, it finds roots of characteristic polynomials formed from the linear prediction coefficients. Third, it solves the original set of linear equations to obtain estimates of the exponential amplitude and sinusoidal phase.

^{2}) via the inverse fast Fourier transform. The impulse response shows that our data are periodic in time. Therefore, Gustavsen’s method is appropriate in our case. As a result, we achieve a better data fit result from Gustavsen’s method. Details of the result are discussed later in this document.

## Application with Experiments

#### Electrical input impedance measurements

_{in}has been computed based on Ohm’s law. The signal-to-noise ratio is enhanced by averaging the response over more than 10 consecutive frames, and the first response block is dropped to consider the steady-state response.

_{in}using Hunt parameters (Eqs. 2-5) is

_{L}is

_{L}≈Z

_{0}coth (a·tube length) where Z

_{0}and a are the characteristic impedance of a tube and the complex propagation function, respectively. They are computed considering the viscous and thermal losses by Keefe [12] for calculation c=334.8 m/s, assuming 20°C room temperature.

_{in}is sensitive to the acoustic load conditions; thus, each condition of Z

_{in}can be denoted as Z

_{in}|

_{A}, Z

_{in}|

_{B}, ..., Z

_{in}|

_{F}excluding the no-tube condition. There are three unknowns (Hunt parameters) in Eq. 25; to compute one set of unknowns, three conditions of Z

_{in}must be selected. The specific computation of the Hunt parameters using three input impedances is explained in Weece and Allen [11] and Kim and Allen [3].

_{in}in the BAR are randomly chosen to calculate the Hunt parameters (please see a Fig. 11 of the published paper by Kim and Allen [3]). Hunt parameters explain the intrinsic characteristics of the system; thus, calculating Hunt parameters is one way to calibrate the system. By manipulating these data, the system’s transfer function is computed as described in the next section.

#### Transfer function from Hunt’s parameters

_{in}) across the two electrical terminals of ED7045 (Knowles, Itasca, IL, USA) (U=0). Other than this direct pressure measurement, all responses are derived from the Hunt parameter calculation, using the electrical input impedance measurements data with various acoustical loads.

## Results

#### The pole-zero fitting

^{2}are much better than those of H(ω). Thus, H(ω)

^{2}is used as our fitting data. The π jump happens during the Hunt’s parameter calculation procedure in MATLAB programming when computing T

_{a}from T

^{2}. As a result, we can expect double poles and zeros for H(ω)

^{2}compared to H(ω), but this does not affect the system analysis as the double poles and zeros are on top of each other. Fig. 5 shows the pole-zero plot of H(ω)

^{2}.

#### All-pass/minimum-phase decomposition of H(ω)

_{AP}(ω) and H

_{MP}(ω) are an all-pass filter and a minimum-phase filter, respectively. Without changing the amplitude of the original filter, H(ω), the maximum-phase zeros in RHP of H(ω), are reflected in LHP to compose H

_{MP}(ω), whereas H

_{AP}(ω) has the original maximum-phase zero (Fig. 6).

_{AP}(ω), has a frequency response where the magnitude is always 1,

_{AP}(ω)=φ(ω). Based on Fig. 7, we can approximate φ(ω)=ωT, yielding

_{g}(ω) of H

_{AP}(ω)

^{2}can be calculated as follows:

_{g}(ω) to an acoustic delay of the sound. Considering the squaring factor, the acoustic group delay is approximately 50 μs below 3.5 kHz (Fig. 8).