ECE 403  HWe – Version 1.11 April 10, 2017  Spring 2017
Univ. of Illinois  Disc: Lect 23, Due: April 13 (Lect 24)  Prof. Allen

Topic of this homework: Model of the human middle ear Deliverable: Show your work. If you hand it in late, you will get zero credit. I would like a paper copy, with your name on it.

No files.doc, please.

Some credit is better than NO credit.

Note: Due: April 13 (Lect 24). Each person is to do their own final writeup, but obviously you can discuss it as much as you like between yourselves. However, you’re crossing the line if you share computer files. The general rule is you need to process all the words you write through your own eyes and fingers. When ever you use material from someone else, you must give them credit, as I do here.

This homework is important because the middle ear is the acoustic load on an in the ear earphone (earbud). The middle ear model has many thing in common with a loudspeaker, but it runs in reverse.

1 Transducer Thevenin Parameters

The Hunt model of a electromechanical motor was first introduced by Kennelly and Nukiyama (1919), and fully developed by Hunt (1954). It contains a gyrator, proposed by Tellegen (1948). The history is carefully discussed in Chapter 1 of Hunt’s book. The equations are also presented in our lab manual, and HW-C.

B_{o}l  ^{f(t)}^{u(t)}^{i(t)}^{ϕ(t)} Z_{e}(s)  Z_{m}(s)

Figure 1: The Hunt model is composed of an electrical impedance Z_{e}(s), a Gyrator having parameter T = lB_{o}, and a mechanical impedance Z_{m}(s). The electrical side is driven by a voltage ϕ(t) and a current i(t), while the mechanical side produces a force f(t) and a velocity v(t).

1.1 To do:

1. Define and then find the ABCD (Transmission) T(s) matrix.

2. Define and then find the impedance matrix Z(s).

3. Compute Δ_{Z} and Δ_{T}. What is the significance of each of these determinants?

4. Define, and then find the Thevenin equivalent force F_{thev}(s) [N], given an input voltage Φ_{o}.

5. Define, and then find the Thevenin equivalent mechanical impedance in mechanical Ohms, Z_{thev}(s) [Ω].

2 Simulation of the middle ear

Acoustic constants: ρ_{o}c = 407 [Rayls], c = 345 [m/s].

The middle ear may be treated as an L = 2.5 cm long stub of transmission line, terminated in a series combination of a stiffness, shown as a capacitor C_{al} and damping, shown as a resistor R_{c}. The experimental details of the cat ear have been extensively explored by Guinan and Peake (1967), and first modeled by Zwislocki (1962, 1957) and much later by Lynch et al. (1982). The model shown here is a significant simplification, but adequate for our purpose. This homework was the basis for the masters thesis of Pierre Parent (Parent and Allen, 2007).

The diameter of the ear canal is d_{c} = 0.75 = 2 ∗ r_{c} cm (area A_{c} = πr_{c}^{2}). As shown in the figure, the free field sound pressure, defined as P_{0}(ω) acts as a source in series with the radiation resistance R_{rad}. The total radiation impedance Z_{rad}(s, A_{rad}) is a combination of the resistance and a reactive component L_{rad}, which represents the local stored field. The two impedances are in parallel

Z_{rad}(s) = sL_{rad}R_{rad}/(sL_{rad} + R_{rad}) = 1/Y_{rad}(s).

where s = σ + jω is the Laplace complex frequency variable. The radiation admittance for a sphere is

Y_{rad} = 1/Z_{rad} = \frac{A_{rad}}{sr_{c}ρ_{o}} + \frac{A_{rad}}{ρ_{o}c} = \frac{1}{sL_{rad}} + \frac{1}{R_{rad}} ,

where r_{c} is the radius of the sphere and A_{rad} is the effective area of the radiation.

If we assume that the pinna (the flap of skin we call our ear) changes the radiation area, then we can account for this pinna-horn effect (a transformer), changing the effective area. The effective area is equivalent to resonant scattering characterized by kr_{rad} ≈ 1, where k = ω/c = 2π/λ. Thus the larger the effective area the smaller the impact of the mass reactance L_{rad} on the radiation efficiency. Namely, the size of L_{rad} decreases as the effective area is made larger via the action of the horn transformer. We can chose the area to match the frequency where the radiation load switches from a mass to a resistance. This frequency was found to be 3 kHz in the cat ear (Rosowski et al., 1988, Fig. 10).

z_{0}^{→}^{U_{ec}(ω)}^{P_{tragus}(ω)}^{P_{0}(ω)}^{R_{rad}}  L  C_{al} R_{c} P_{c}(ω)^{L_{rad}}

Figure 2: Model of the ear canal, terminated by the radiation impedance Z_{rad}(s) at the tragus (x = 0), and by the eardrum and cochlea at x = L. The horn transformation is the transformer that represents the concha, that converts the area of the canal to the area of the pinna/concha opening. This transformer ratio seems necessary to reduce the effective radiation mass seen by the ear canal, to improve the matching of energy between the canal and free-space. It does this by making the area defining L_{rad} larger by the transformer turns ratio. The effective turns ratio needs to be estimated from measured data (Rosowski et al., 1988, Fig. 10).

The model: To model the middle ear we first need to simulate the transmission line by using d’Alembert solution of the wave equation. This simulation will operate in the time domain.

Create two arrays that represent the forward and backward velocity traveling waves u_{+}(x − ct) and u_{−}(x + ct). Let c = 345 m/s.

The boundary conditions at each end of the line will be realized via reflection coefficients R(L, s) at the cochlear end and R(0, s) at the input. Initially assume that P_{0} is zero, and drive the middle ear with a velocity source U_{ec} at the tragus.

On the left we terminate the ear in a radiation impedance

Z_{rad}(s) = \frac{sL_{rad}R_{rad}}{sL_{rad} + R_{rad}} .

At the cochlear end we terminate the line with an impedance

Z_{c} = R_{c} + 1/(sC_{al}),

where R_{c} is the cochlear impedance and C_{al} is the stiffness of the annular ligament, which is the ligament that holds the stapes in the oval window. The cochlear resistance (R_{c}) is assumed to be twice the characteristic impedance of the ear canal.

Each impedance may be converted into a reflection coefficient, defined as

R(s) ≡ \frac{U_{−}}{U_{+}} = \frac{Z(s) − z_{0}}{Z(s) + z_{0}} ,

where z_{0} = ρ_{o}c/A_{c} is the characteristic impedance of the ear canal, having area A_{c}. For example, the reflection coefficient at x = L is the transfer function between the wave reflected U_{−}(L, s) over the incident wave U_{+}(L, s), as defined above. The two reflection coefficients are a function of complex frequency s.

Use the bilinear z transform^{1} to convert R(s) into a discrete time IIR filter having numerator polynomial N(z) (zeros) and denominator polynomial D(z) (poles), so that the reflectance filters (e.g., R(0, s) and R(L, s)) are implemented in the time domain. I recommend that you implement the filter manually rather than use the filter() command. [If you chose to use the filter() command, you must save state between calls. See help filter for the details on how to do this, and carefully verify that it is working.]

Finally is the question of the boundary conditions. This is where things are a little tricky. These are dealt with by defining the end reflection coefficients

R(L, s) ≡ \frac{U_{−}(L, s)}{U_{+}(L, s)}

and

R(0, s) ≡ \frac{U_{+}(0, s)}{U_{−}(0, s)}.

In the frequency domain we can find

U_{−}(L, s) = R(L, s)U_{+}(L, s)

and

U_{+}(0, s) = R(0, s)U_{−}(0, s).

These relations give the junction reflected wave in terms of the incident wave. For example, at x = L the incident wave is U_{+} and the reflected wave is U_{−}. At x = 0 the signs must change. Since this relation is in the frequency domain, it represents a convolution in the time domain. Since we wish to run the middle ear model in the time domain, we must design the filter represented by these reflection coefficients. To do this we may use the bilinear z transform.

Each reflection coefficient R is given by

R(s) = \frac{Z(s) − z_{0}}{Z(s) + z_{0}} .

1substitute s = \frac{2}{T_{s}} \frac{1−z^{−1}}{1+z^{−1}} where T_{s} is the sampling frequency.

I derived this formula in class and showed that Z(s) is the load impedance and z_{0} = ρ_{o}c/A is the characteristic impedance of the transmission line. Taking the case of x = L we find

R(L, s) = \frac{R_{c} + 1/(sC_{al}) − z_{0}}{R_{c} + 1/(sC_{al}) + z_{0}} .

When using the bilinear z transform one substitutes 2F_{s} \frac{z−1}{z+1} for s. Try doc bilinear from

Matlab^{\begin{matrix}c \\ ◯ \end{matrix}}.

Show your work. Check all the filters and all the components carefully. Show your work. For example, show the impulse response for each of these filters r(0, t) ↔

R(0, s) and r(L, t) ↔ R(L, s). Also test your time domain transmission line impulse response by setting the reflection coefficients to 1, and then compute the impulse response of p_{c}(L, t) with an impulse at the input (e.g., find p_{ec}(0, t) = δ(t)). Show your work. Write this up with a full description of what you did, so that I can easily understand how you got your results. I will grade this homework partially on the readability of the presentation.

To do:

1. Find z_{0}, the characteristic impedance of the ear canal. (1 min)

2. Assume that the cochlear load resistor R_{c} is equal to twice the characteristic impedance of the ear canal. Find C_{al}, by assuming that the impedance of this element (Z_{al} ≡ 1/sC_{al}) is equal to the cochlear impedance R_{c} at a frequency of 0.800 kHz. Show all your work. (1 min)

3. Determine the radiation impedance looking out the ear canal, and compute the resonant frequency, defined as the frequency f_{c} ≡ s/2πj, for which the real and imaginary parts are equal. For an example of the frequency response of this impedance look at Rosowski et al. (1988, Fig. 10) provided on the class web site. Hint: A_{rad} should be approximately the surface area of a half-sphere.

4. The simulation is to be done in the time domain. Define two arrays in Matlab^{\begin{matrix}c \\ ◯ \end{matrix}} called u_{+}(x−ct) and u_{−}(x+ct), which represent the velocity of the waves in the ear canal, traveling in the two directions. At each time step, the data in the first array is shifted to the right, while the second is shifted to the left. The distance between two points corresponds in the array is dx = cT where sampling rate T = 1/F_{s}.Choose the sampling rate F_{s} such that there are at least 5 points along the length of the ear canal (you can modify the length of the canal slightly to make the samples come out at the ends).

5. Find the formula for the input impedance Z(0, s) of the middle ear at the entrance of the ear canal, when the cochlea is “blocked” (Z_{c} = ∞ or R(L, s) = 1)? (Hint, this has a simple answer that you can easily derive).

6. Find the two reflection coefficients and plot their magnitude frequency response.

7. What are the coefficients of the time domain filters that implement the reflection coefficients?

(Hint 1: Use Matlab’s bilinear() command. Hint 2: here is the answer to 1 decimal place, to help you debug up to this point. I’ve used a sampling rate of ≈70 kHz to compute this: R(0, z^{−1}) = (.5 − .08z^{−1})/(1 − .4z^{−1}) and R(L, z^{−1}) = (0.3 − 0.3z^{−1})/(1 − 0.95z^{−1}).)

8. What are the poles and zeros in s and in z^{−1} planes? List them, or plot them using ‘x‘ to represent poles, and ‘o’ to represent zeros.

9. Compute and plot the input impedance (numerically derived at x = 0), using your Matlab model at the point labeled P_{ec}(ω) (i.e., at x = 0). Use log-frequency [Hz] and dB for the magnitude. Use z_{0} as the dB reference (e.g., plot 20 log_{10}(|Z_{c}(0, f)|/z_{0})). Don’t forget that Z ≡ P/U where p = (u^{+} − u^{−})/z_{0} and p ≡ u^{+} − u^{−}. Compute the time sequence for p and u and then FFT these and take the ratio. Be sure that each function has settled to zero or else you will alias the answer. I am not asking for the transmission line model, I’m looking for the numerical solution.

10. Using the time domain model, compute the tragus transfer function of the middle ear, defined as the ratio of the cochlear pressure to the ear canal pressure P_{c}(ω)/P_{ec}(ω). Pressure P_{ec} is more precisely P_{tragus}.

11. Compute the free field transfer function of the middle ear, defined as the ratio of the cochlear pressure P_{c}(ω)/P_{0}(ω). I suggest you work backward, and find the pressure P_{0}, the assumed source, required to give the given source U_{ec}.

Note that this question may be just too hard, and take to long, so do this problem last (or just skip it if you don’t have the time to deal with it. (30-60 min)

12. A 16 electrical ohm earphone delivers 120 dB SPL into the ear canal with 1 volt RMS input at 1 kHz. Using your middle ear model, calculate the efficiency (Acoustic/Electrical power ratio, in %) of the earphone, at 1 kHz.

References

Guinan, J. J. and Peake, W. T. (1967), “Middle-ear characteristics of anethetized cats,” J. Acoust. Soc. Am. 41, 1237–1261.

Hunt, F. V. (1954), Electroacoustics: The Analysis of Transduction and Its Historical Background (Harvard University Press, Cambridge, MA).

Kennelly, A. E. and Nukiyama, H. (1919), “Electromagnetic theory of the telephone receiver with special reference to motional impedance,” The 348th Meeting of the American Institute of Electrical Engineers .

Lynch, J. T., Nedzelnitsky, V., and Peake, W. T. (1982), “Input impedance of the cochlea in cat,” J. Acoust. Soc. Am. 72, 108–130.

Parent, P. and Allen, J. B. (2007), “Wave model of the cat tympanic membrane,” J. Acoust. Soc. Am. 122(2), 918–931.

Rosowski, J. J., Carney, L. H., and Peake, W. T. (1988), “The radiation impedance of the external ear of cat: Measurements and applications,” J. Acoust. Soc. Am. 84, 1695–1708.

Tellegen, B. D. H. (1948), “The Gyrator: A new electric network element.” Philips Research Reports 3, 81–101.

Zwislocki, J. (1957), “Some impedance measurements on normal and pathological ears,” J. Acoust. Soc. Am. 29, 1312–1317.

Zwislocki, J. (1962), “Analysis of the middle-ear function. Part I: Input impedance,” J. Acoust. Soc. Am. 34, 1514–1523.