ECE 403  HWd – Version 1.14 April 10, 2017  Spring 2017
Univ. of Illinois  Lec 16, Disc: Lec 17, Due: Lec 20  Prof. Allen

Topic of this homework: Webster Horns: 1D, 3D & Exp Horns; Reflectance; Thermal noise; Hilbert Transform

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 doc files.

Some credit is better than NO credit.

1 Wave equation

1.1 History of the wave equation

1. What year did d’Alembert derive his solution to the wave equation? Solution: d’Alembert first proved this in 1747.

2. What is the form of D’Alembert’s solution? Solution: f(t − x/c) + g(t + x/c)

3. Who was the first person to calculate the speed of sound, and what was the result? Solution: Newton did this in 1648. His formula was in error due to the dynamic stiffness of air, which is γP_{o}. His result was too small by the factor of \sqrt{1.4}.

1.2 The Webster wave equation:

The Webster Horn equation may be written in the time domain as 1D transmission line equation:

\frac{∂}{∂x} [\begin{matrix}p(x, t) \\ ν(x, t) \end{matrix}] = − [\begin{matrix}0 & \frac{ρ_{o}}{A(x)} \\ \frac{A(x)}{γP_{o}} & 0 \end{matrix}] \frac{∂}{∂t} [\begin{matrix}p(x, t) \\ ν(x, t) \end{matrix}] ,  (1)

where ν(x, t) = A(x)u(x, t) is the volume velocity, more generally defined as the integral over the normal component of the particle velocity u(x, t), over the cross-sectional area A(x) of the tube. Transforming to the frequency domain we have

\frac{d}{dx} [\begin{matrix}P (x, ω) \\ V (x, ω) \end{matrix}] = − [\begin{matrix}0 & Z_{s}(s, x) \\ Y_{s}(s, x) & 0 \end{matrix}] [\begin{matrix}P (x, ω) \\ V (x, ω). \end{matrix}]  (2)

Here we use the complex Laplace frequency s when referring to the per-unit impedance

Z_{s}(s, x) ≡ s\frac{ρ_{o}}{A(x)} = sM(x)  (3)

and per-unit admittance

Y_{s}(s, x) ≡ s\frac{A(x)}{γP_{o}} = sC(x),  (4)

where M(x) = ρ_{o}/A(x) is the horn’s per-unit-length mass, and C(x) = A(x)/γP_{o} per-unit-length compliance, to remind ourselves that these functions must be causal, and except at their poles, analytic in s.

The horn of a loudspeaker cone is either conical A(x) = A_{0}(x/x_{0})^{2}, or in the shape of an exponential.

1.3 To Do:

1. Assuming a conical horn, having area A(x) = A_{o}(x/x_{o})^{2} with A_{o} ≤ 4π, rewrite these equations as a second order equation solely in terms of the pressure P (remove U), and thereby find the frequency domain solutions P_{±}(x, s) for the conical horn equation. Solution: If we let P_{x} ≡ ∂P/∂x (i.e., the partial with respect to space) then

P_{x} + ZU = 0  (5)

and

U′ + YP = 0.  (6)

Taking the partial wrt x of the first equation, and then using the second, gives

P_{xx} + Z_{x}U + ZU_{x} = P_{xx} − \frac{Z_{x}}{Z} P_{x} − ZYP = 0.  (7)

Using the relation

Z_{x}/Z = \frac{d}{dx} ln Z = −\frac{d}{dx} ln A(x),  (8)

with A = A_{o}(x/x_{o})^{2}, we find

P_{xx} + \frac{2}{x}P = \frac{s^{2}}{c^{2}} P.  (9)

Just as it was important to replace real frequency ω with the Laplace frequency s, since the roots are typically complex, using the same reasoning, we replace the “real wave number” k = ω/c with a complex wave number as κ(s). Only in the case of non-dispersive waves (e.g., plane waves) is κ(s) = s/c.

Since the wave equation is

\frac{∂^{2}P}{∂x^{2}} = \frac{1}{c^{2}} \frac{∂^{2}P}{∂x^{2}}  (10)

by inspection we see that ω^{2}/c^{2} = ZY, which results in the final formula for the speed of sound.

2. Starting from the transmission line equations given above, and assuming an exponential area

function

A(x) = A_{0}e^{2mx}

(m is a positive constant, called the horn flair parameter), derive the exponential horn equation for the pressure

\frac{∂^{2}p(x, t)}{∂x^{2}} + 2m\frac{∂p(x, t)}{∂x} = \frac{1}{c^{2}} \frac{∂^{2}p(x, t)}{∂^{2}t}  (11)

Solution: Starting from the basic definitions with A(x) = A_{0}e^{2mx} along with the basic equation for a horn, explicitly write out the two equations

\frac{dP}{dx} + s\frac{ρ_{o}}{A_{0}} e^{−2mx}V = 0  (12)

\frac{dV}{dx} + s\frac{A_{0}e^{2mx}}{γP_{o}} P = 0  (13)

Next solve for the pressure (remove V):

\frac{d^{2}P}{dx^{2}} + s\frac{ρ_{o}}{A_{0}} (e^{−2mx} \frac{dV}{dx} + V \frac{d}{dx}e^{−2mx}) = 0  (14)

or

\frac{d^{2}P}{dx^{2}} + s\frac{ρ_{o}}{A(x)} (\frac{dV}{dx} − 2mV ) = 0.  (15)

Going back to the basic equations, again removing V

\frac{d^{2}P}{dx^{2}} + s\frac{ρ_{o}}{A(x)} (−s\frac{A(x)}{γP_{o}}P + 2m\frac{A(x)}{sρ_{o}} \frac{dP}{dx} ) = 0  (16)

which simplifies to the requested result

\frac{d^{2}P}{dx^{2}} + 2m\frac{dP}{dx} − s^{2} \frac{ρ_{o}}{γP_{o}} P = 0.  (17)

since −s^{2} = ω^{2} and γP_{o} = ρ_{o}c^{2}.

3. In general the soluton to a wave equation is of the form

p(x, t) = P^{+}(κ, s)e^{κx}e^{st} + P^{−}(κ, s)e^{−κx}e^{st}

where s = σ + ȷω is the Laplace frequency and κ(s) is the complex “wave number.” What is κ(s) for the

(a) conical horn? Solution: κ(s) = ±s/c.

(b) exponential horn? Solution: κ(s) = −m ± \sqrt{m^{2} + (s/c)^{2}}.

(c) What is the significance of κ(s)? Solution: It is called the dispersion relation of the differential equation, that relates the wavelength to the frequency.

(d) Why is it a function of s? Solution: In general the wavelength is a function of frequency. This is a generalization of the plane wave relation λf = c, or κ(s) = s/c.

(e) What is the role of P^{±}(κ, s)? Solution: These are aplitudes of the forward and reverse traveling waves, that are to be determined by applying boundary conditions to the solution.

4. Show that the solution to Eq. 11 is of the form

P^{±}(x, s) = e^{−mx}e^{∓\sqrt{m^{2}+(s/c)^{2}} x}.

Solution: This equation is an ordinary constant coefficient equation. Accordingly, substitution of p(x, t) = P (κ, s)e^{−κx}e^{st}, one may find κ = −m ± \sqrt{m^{2} + (s/c)^{2}}, giving the solution, in the freq. domain, of the requested form.

Discussion: This solution is loss-less under all conditions. For high frequencies, when ω > mc, the pressure changes by e^{−mx} as it propagates with a frequency dependent delay, making it an acoustic transformer. At very high frequencies, when m can be ignored with respect to ω/c, the wave propagates without dispersion, but still with a decay. It is still loss-less. The decay in amplitude is due to the change in area. At very low frequencies, the wave solution is still causal, but no longer decays spacial [the exp() term is canceled out]. To see this requires taking the inverse Laplace transform of P (x, s). Note the interesting and seemingly related relations (s is the Laplace frequency and ↔ represents the Laplace transform)

sinh^{−1}(s) = log_{e}[s + \sqrt{s^{2} + 1}]  (18)

(see Matlab’s doc asinh) and the low-pass sinc()-like relation

\frac{J_{1}(t)}{t}U(t) ↔ \sqrt{s^{2} + 1} − s  (19)

1.4 Reflectance:

1. Find (derive) the formula for the “input” impedance of a transmission line, having characteristic impedance z_{0}(x, s), in terms of the reflectance. Define all the terms. Hint, I did this in class several times. Solution: Define z_{0} ≡ \frac{P_{+}}{U_{+}} and R ≡ P_{−}/P_{+} = U_{−}/U_{+}, then

Z = P/U = \frac{P_{+} + P_{−}}{U_{+} − U_{−}} = z_{0} \frac{1 + R}{1 − R}.  (20)

2. Find the formula for the reflectance R(s) in terms of the load impedance Z_{L}(s) and the characteristic impedance z_{0} if:^{1} Solution: Solve above for R(s) in terms of Z_{L} gives:

R = \frac{Z_{L} − z_{0}}{Z_{L} + z_{0}}  (21)

(a) Z_{L}(x, s) = r [Nt-s/m^{5}]

(b) Z_{L}(x, s) = 1/sC [Nt-s/m^{5}]

(c) Z_{L}(x, s) = r||sM [Nt-s/m^{5}] Solution: Let Z_{L} = \frac{rsM}{r+sM} then

R = \frac{rsM − z_{0}(r + sM)}{rsM + z_{0}(r + sM)} = \frac{(r − z_{0})Ms − z_{0}r}{(r + z_{0})Ms + z_{0}r} .  (22)

(d) Two transmission lines are in cascade, the first one having an area of 1 [cm^{2}] and a second having an area of 2 [cm^{2}], with lengths L_{1} and L_{2} respectively, terminated with a resistor r = ρc/A, where A = 2 × 10^{−4} [m^{2}]. Find R(x = 0, s). Solution: Since the second line is terminated in its own impedance, it is just a resistor at its input, which makes the problem very simple. As a result

R(s) = \frac{1/1 − 1/2}{1/1 + 1/2}e^{−s2L_{1}/c} = 1/3e^{−s2L_{1}/c},  (23)

where L_{1} is the length of the first TL. Note that if the line were not matched at the end, the story would be very different.

(e) What is the inverse Laplace transform of

i. H(s) = 1/(s + 1)? Find h(t). Solution: h(t) = e^{−t}U(t)

ii. R(s) = \frac{Z−z_{0}}{Z+z_{0}} where L = 1, Z = 1 and z_{0} = 2? Find r(t) at the input. Solution:

r(t) = (1 − 2)/(1 + 2)δ(t − 2L/c) = −δ(t − 2)/3 iii. H(s) = s/(s + 1)? Solution: H(s) = 1-1/(s+1) ↔ h(t) = \frac{d}{dt}e^{−t}U(t) = δ(t) − e^{−t}U(t)

2 Nyquist Thm on Thermal noise

The purpose of this problem is to do a simulation of Harry Nyquist’s famous result on the noise of a resistor. The experiment is to use a Thevenin equivalent model of a resistor as a resistance R in series with a voltage source. A stub of transmission line having characteristic impedance z_{0} is terminated in each end with this Thevenin model, with R = z_{0}. Then at t = 0, the resistance is short or open circuited. If open circuited you may watch the voltage at one end, and if short circuited, you may watch the current. Let’s monitor the voltage with an open circuit. This setup is shown in the figure.

1Note that r, C and M represent an acoustic resistance, compliance and mass. Namely they are positive constants.

R  R V_{m}(x = 0, t < 0)  V_{m}(L, t < 0) z_{0}  z_{0}  V_{m}(t > 0)

The transmission line stores the voltage at t = 0 once the switch is opened removing the resistors (and also the Thevenin source). At that point voltage V_{m}(L, t) becomes periodic with a period of 2L/c, where L is the length of the line and c is the speed of the wave.

Define an array that has a duration of the period, and load it with thermal random samples [x=randn(1,N) with N the number of noise samples]. Set the RMS to 1. The RMS may be computed using std(x,1). Use a sampling period T = 1/f_{s} with f_{s} = 10 × 10^{3} Hz for this experiment, and let c = 345 [m/s] be the speed of sound (an acoustic transmission line), with L = 10 [m].

1. What is the fundamental period of the noise? Solution: The period is the round trip delay T = 2L/c ≈ 58 [ms].

2. Every periodic signal has a Fourier series. If the period is T [s], what are the Fourier series frequecies? You know the answer to the question, but you may need to think about it. This is not difficult. Solution: If a function is periodic with period T, then this may be indicated as f((t))_{T}. This notation with the double parentheses is meant to indicate that the function is periodic, and the subscript T is the period. Any periodic funciton has harmonics at the Fourier series frequencies given by f_{k} = k/T [Hz], where k is an integer. When k = 0 it is called DC (zero frequency). All the other frequencies are called AC. Such terms are jargon. The fundamental frequency is corresponds to k = 1, namely f_{1} = 1/T, also known as the the first harmonic. This may be confusing, until you think about it a minute. The second harmonic corresponds to k = 2 with f_{2} = 2f_{0}. The third harmonic is k = 3 with f_{3} = 3f_{0}, etc. You first learned this in Physics 101, or maybe in high school. The normal modes of an organ pipe and of a guitar string obey this relationship (approximately).

3. Plot two periods of the time domain signal. Using fft(), find the spectrum of the periodic noise, for 1, 4 and 10 periods. Be sure to properly label all axes (with units)! For each of the FFT plots, show one figure of the full FFT (from 0 Hz up to the half-sample-rate), and one figure zoomed in on the range 0 to ∼300 Hz (linear axis, not log axis). Comment/explain any observations.

Note: Let’s say your array length is N. When creating 10 periods of random noise, don’t use randn(1,10*N). Use randn(1,N) to randomly populate 1 period, and then repeat that noise 10 times, so that it is identical. Remember, its in a delay line. Solution: There will be harmonics at f_{n} = n/T = 17[Hz], and the bandwidth of the noise lines become narrow as more periods are included.

4. Why won’t the values of your spectral peaks be the same as your fellow students? Average the values of your spectral peaks over many noise samples (that is, many initializations of the randn(1,N) array), and plot the resulting spectral average. Just do this for the case of 10 periods.

Solution: The average harmonic amplitude may be computed using Parseval’s formula, since the power in the noise is independent of how it is calculated, either in the time domain, or in the frequency domain. Parseval’s Thm is

\frac{1}{T} ∫_{0}^{T} |w(t)|^{2}dt = ∑_{k=−∞}^{∞} |W (2πk/T)|^{2}.

On the average the spectral level should be |W (ω)|, but because of the finite sample of the noise, for any given frozen sample, the spectrum will not be constant. This is very important point that many students (and graduates) seem to miss. It is the ensemble that is constant, not the FT of a fixed sample. No matter now long the noise sample is, the spectral samples |W (ω_{k})| will never become the same, except in the average.

When you use the FFT, it is necessarily over a finite bandwidth, determined by the sampling rate you chose. Thus the above formula is not the one to use for the case of a band-limited signal (its not wrong, just not the most relevant, given the more detailed assumption of a band-limited and sampled noise. Thus the formula to use is Parseval’s Thm for the DFT

∑_{n=0}^{N−1} x[n]^{2} = \frac{1}{N} ∑_{k=0}^{N−1} |X(e^{j2πk/(N−1)})|^{2}.

Now the problem with this formula is that the normalization is wrong, as it will not give you the true ”spectral level” you would actually measure with a volt meter placed at the output of a filter around one of the spectral lines (i.e., what you would see coming out of a spectral analyzer, calibrated in volts/\sqrt{Hz}.

In other words, the solution to this problem is a nightmare of details, and the books (at least that I read) don’t tell you any of this. I prefer to proceed with a numerical solution, for speed and accuracy. Namely, make a narrow band filter around one of the spectral lines, and measure the level of the line. Do this for many noise samples, and take an average. That would be a good problem. Maybe next time.

5. Given T = 300 degrees Kelvin (300-273 = 27 degrees C) and k = 1.38 × 10^{−23} [J/degree K], what would the RMS value of the voltage be? Justify your answer. Hint: this is also know as “Johnson/Nyquist Noise”. By sampling we inherently band-limit the signal. Solution: Use the formula from the paper, and you get that the RMS level is V_{rms} = \sqrt{4kT RB}, which is the same as V_{rms}^{2}/R = 4kT B.

3 Hilbert transform

In all parts of this problem h(t) ↔ H(s) and H(ω) = H(s)|_{s=jω} Analyze the real impulse response

h(t) = e^{−t/τ_{0}}u(t),

with τ_{0} = 10 [ms], in terms of its Hilbert transform (integral) relations.

1. Find H(s), the Laplace transform of h(t). Solution: Let a = 1/τ_{0}. Then e^{−at}u(t) ↔ \frac{1}{s+a}.

2. Where are the poles of H(s) ↔ e^{−t/τ_{0}}u(t)? Solution: Let a = 1/τ_{0}. Since e^{−at}u(t) ↔ \frac{1}{s+a} H(s) has a simple pole at s = −1/τ_{0}.

3. Evaluate the following:

(a) Determine g(t) ≡ h(t) ∗ δ(t). (∗ represents convolution.) Solution: Convolution with a delta function is an identity, thus g(t) = h(t) ∗ δ(t) = h(t)

(b) Determine h(t)δ(t−1). Solution: Multiplying by a delta function evaluates the function at the time of the delta function. Thus h(t)δ(t − 1) = e^{−1/τ_{0}}δ(t − 1).

(c) Determine h(0). Solution: h(0^{−}) = 0 and h(0^{+}) = 1. I would conclude that at t = 0, h(t) is not defined. What does this say about delta functions?

4. Find the real and imaginary parts of H(ω) ≡ H(s)|_{s=jω}. Solution: First rationalize the

denominator:

H(s) = \frac{s^{∗} + a}{(s + a)(s^{∗} + a)} = \frac{σ − jω + a}{σ^{2} + ω^{2} + a^{2}} .

Next take the real ℜ and imaginary ℑ parts, and evaluate s on the ω axis (set σ = 0):

ℜH(ω) = a/(a^{2} + ω^{2}), ℑH(ω) = −ω/(a^{2} + ω^{2}).

The real part is constant below the cutoff (resonant) frequency ω = a, and goes at -12 dB/oct above the cutoff. The imaginary part is bandpass with ±6 dB/Oct above and below the resonance frequency.

5. Write out the symmetric h_{e}(t) and antisymmetric h_{o}(t) functions.

Solution: 2h_{e}(t) = h(t) + h(−t) ≡ e^{−at}u(t) + e^{at}u(−t), while 2h_{o}(t) = h(t) − h(−t) ≡ e^{−at}u(t) − e^{at}u(−t). It trivially follows that h(t) = h_{e}(t) + h_{o}(t).

6. Find the Fourier transforms of h_{e}(t) ↔ H_{e}(ω) and h_{o}(t) ↔ H_{o}(ω).

Solution: Since h(−t) ↔ H^{∗}(ω), a symmetric time function is real in the frequency domain,

2h_{e}(t) = h(t) + h(−t) ↔ H(ω) + H^{∗}(ω) = 2ℜH(ω),

thus h_{e}(t) ↔ ℜH(ω). In a similar fashion, an antisymmetric time function is pure imaginary

h_{o}(t) ↔ jℑH(ω).

Again with a ≡ 1/τ_{0}:

H_{e}(ω) = ℜH(ω) = \frac{a}{ω^{2} + a^{2}} ,

H_{o}(ω) = jℑH(ω) = \frac{−jω}{ω^{2} + a^{2}} ,

thus

H(ω) = H_{e}(ω) + H_{o}(ω) ↔ h(t) = h_{e}(t) + h_{o}(t).

The inverse Fourier transform of H_{o}(ω) is zero at t = 0, which makes it very different from the inverse Laplace transform, which is not defined at t = 0. What is the inverse FT of H_{e}(ω)? Be sure to discuss what happens at t = 0.

7. Find the Hilbert (integral) relations between H_{r} ≡ ℜH(ω) (real part) and H_{i} ≡ ℑH(ω) (imag part) of H(ω).

Solution: It follows from the above results that

jH_{i}(ω) = \frac{1}{jπ} ∫ \frac{H_{r}(ω′)}{ω − ω′}dω′  (24)

which for the case at hand is

\frac{ω}{ω^{2} + a^{2}} = \frac{1}{π} ∫ \frac{a dω}{(ω′ − ω)(ω^{′2} + a^{2})}.  (25)

A second derivation of the requested integrals may be found from

h(t) = h(t)u(t),  (26)

(note this is not exactly true at t = 0) which after a FT, results in

H(ω) = \frac{1}{2π} H(ω) ⋆ (πδ(ω) + \frac{1}{jω} ) ,  (27)

which may be rewritten as

H_{r}(ω) = \frac{1}{2}H_{r}(ω) + H_{i}(ω) ⋆ \frac{1}{2πω} .  (28)

The final relations are [Papoulis (1977), Signal Analysis, McGraw Hill, page 251]

H_{r}(ω) = \frac{1}{π} ∫ \frac{H_{i}(ω′)}{ω − ω′}dω′  and  H_{i}(ω) = −\frac{1}{π} ∫ \frac{H_{r}(ω′)}{ω − ω′}dω′  (29)