Horn Loudspeaker Simulation part 1: Radiation and T-Matrix

The first part of the series will cover basic radiation and horn throat impedance.

Radiation Impedance

We will work our way from the horn mouth, which means we start with the radiation impedance. The simplest model, which is used in Hornresp for many of the horn flares, is the circular piston in an infinite baffle. This is a reference model that is easy to calculate using built-in functions. The equation for the acoustic impedance for a piston of radius a and area S is

It uses two so-called special fuctions, the Bessel function of order 1, and the Struve function of order 1. The Bessel function is a built-in function in Octave/Matlab, but the Struve function can fortunately be approximated using other standard functions:

Our first functions will therefore be

function Znorm = circularPistonIB(ka)
R = 1 - besselj(1,2*ka)./(ka);
X = struveH1(2*ka)./(ka);
Znorm = R+1i*X;

function H1 = struveH1(x)
% Using a Struve function approximation from http://mathworld.wolfram.com/StruveFunction.html
H1 = 2/pi - besselj(0,x) + (16/pi - 5)*sin(x)./x + (12-36/pi)*(1-cos(x))./x.^2;

These can be placed in the same file since H1 is't used for other things. Note that the function returns the normalized impedance, and the argument is ka, not frequency and area. We will see why later.

We can now test our radiation impedance function easily like this:

ka = logspace(-1,2,500);

Z = circularPistonIB(ka);

figure(1);
loglog(ka, real(Z), ka, imag(Z));

This will produce the following plot:

CircPistonZaPlot

The blue line is the real part of the impedance, the orange line is the reactive part.

T-Matrix

Relations between the pressure and volume velocity at one end of the horn to pressure and volume velocity at the other end can be describes by transfer matrices (T-matrices):

The numbers of the matrices can be read like this: T_12 is the matrix producing p_1 and U_1 given p_2 and U_2. We can find the values of the matrix elements for a number of horns: conical, exponential, hypex and bessel. One advantage of the T-matrices is that they can be multiplied together to form composite matrices for multisegment horns, or a horn with no analytical solution can be represented by a large number of small (typically conical) segments approximating the horn profile.

Exponential Horn

Here we use the following definition for the exponential horn:

where

The reason for using 2mx rather than mx in the exponent is that it makes m equal to the cutoff wavenumber of the horn. It's my personal preference, but it makes the equations cleaner; some derivations define a second variable equal to m/2 and so on to simplify notation.

Apart from the cylindrical tube, the expontial horn has one of the simplest expressions for the T-matrix:

ExpoTmatrix

where

which becomes imaginary below cutoff. Fortunately, Octave/Matlab handles complex numbers directly without us having to think about it, so this poses no problem. We can implement this as a function:

function [a,b,c,d] = expoHornMatrix(k,Zrc,S1,S2,L)

m = log(S2/S1)/(2*L);
gamma = sqrt(k.^2-m^2);
gL = gamma*L;

singl = sin(gL);
cosgl = cos(gL);

emL = exp(m*L);
a = emL*(cosgl-m./gamma.*singl);
b = emL*1j*Zrc/S2*k./gamma.*singl;
c = emL*1j*S1/Zrc*k./gamma.*singl;
d = emL*S1/S2*(cosgl+m./gamma.*singl);

Octave/Matlab handles matrix multiplication directly, so we could generate 2x2 matrices, and multiply them together as needed. Instead I have decided to implement the horn matrices as vectors of matrix entries, and do the matrix multiplication "by hand". It is possible to implement "vectorized matrix multiplication", where you have a 3D matrix consisting of an array of matrices, and make a functions that multiply them together and other operations, but using vectors like here utilises the vectorized operations in Octave/Matlab in, I think, a more efficient way.

Now that we have our T-matrix and radiation impedance, we can find the throat imedance as

We are now ready to calculate the throat impedance of an exponential horn.

% Define medium properties (Hornresp default values)

rho = 1.205;
c = 344;

% Define frequency range
fmin = 10;
fmax = 20e3;
freq = logspace(log10(fmin), log10(fmax), 533);
k = 2*pi*freq/c;

% Define horn dimensions
S1 = 80e-4;
S2 = 5000e-4;
L12 = 150e-2;

% Calculate radiation impedance
a = sqrt(S2/pi);
Z2 = rho*c/S2 * circularPistonIB(k*a);

% Calculate horn matrix
Zrc = rho*c;
[a12,b12,c12,d12] = expoHornMatrix(k,Zrc,S1,S2,L12);

% Calculate and normalize throat impedance
Z1 = (a12.*Z2 + b12) ./ (c12.*Z2 + d12);
Z1norm = Z1*S1/(rho*c);

figure(1);
semilogx(freq, real(Z1norm), 'k', freq, imag(Z1norm), 'r');
xlim([fmin, fmax]);
ylim([-0.5, 2.5]);
xlabel('Frequency (hertz)');
title('Acoustical impedance');
grid

The result is shown below, together with the same horn simulated in Hornresp:

OctaveZth1

HornrespReferenceZth1

I have stretched the Octave figure to match the aspect ratio of the Hornresp plot. The impedance has been normalised before plotting, since this is what is done in Hornresp. Normalized impedance always approaches 1 at high frequencies, and it makes it easier to se how the ripple changes for different profiles, even if they have different throat areas.

The .m files for this post can be found here.