% ZinMatrix.m
%
% This matlab script calculates both the theoretical input impedance
% and impulse response of an open cylindrical pipe. The calculation
% ignores viscous and thermal losses. It uses the Levine & Schwinger
% unflanged tube radiation approximation. The input impedance is
% calculated using a transmission matrix. The impulse response is
% then calculated as the IFFT of the input impedance.
%
% by Gary P. Scavone
% Physical Constants and evaluation frequencies
fs = 44100;
fmax = 22050; % maximum evaluation frequency (Hz)
N = 10000; % number of frequencies for evaluation (even)
finc = fmax / (N-1);
f = 0:finc:fmax; % evaluation frequencies
omega = 2 * pi * f; % radian frequencies
c = 347.23; % speed of sound (m/sec) %
rho = 1.1769; % density of air (kg/m^3) %
k = omega/c;
% Bore dimensions
ra = 0.016 / 2; % radius of pipe (m) %
% Length of pipe (m)
L = 0.35;
delays = round( 2 * L * fs / c); % compute equivalent length is unit delays and then round
L = c * delays / fs / 2; % compute corresponding length from rounded delay length
% Pipe characteristic impedance
Zc = rho * c / (pi * ra^2);
% Lossless Cylindrical Transmission Matrix Elements
sinL = sin(L*k);
cosL = cos(L*k);
A = cosL;
B = 1i * Zc * sinL ;
C = 1i * sinL / Zc ;
D = cosL;
% Radiation impedance (added 0.98 factor for extra losses)
%Z_L = 0; % Ideal open-end approximation
z = k * ra;
[r,l] = reflpoly(z); % Levine & Schwinger approximation
R = -0.98 * r .* exp(-2*1i*z.*l); % Reflectance with small extra loss factor
Z_L = Zc * (1+R) ./ (1-R);
% Matrix Multiplications
Zin = (B + A.*Z_L) ./ (D + C.*Z_L);
% Clear figure window, bring it to front, and plot
clf
% Plot Impedance
subplot(2,1,1)
plot( f, abs(Zin / Zc) );
grid
xlabel('Frequency (Hz)')
ylabel('| Impedance | / Zc')
v = axis;
axis( [0 10000 v(3) v(4) ] );
% Calculate the impulse response
% In order to obtain a real valued IFFT of the input
% impedance, it is necessary to make H conjugate symmetric.
% Since h is technically not time-limited, there ends up
% being some high frequency aliasing which appears as a
% ripple in the beginning of the impulse response.
H = Zin;
H = [ H, conj( fliplr( H(2:end-1) ) ) ]; % make conjugate symmetric
h = real( ifft( H ) );
T = 1 / (2 * fmax );
t = (0 : T : (2*N-3)*T) * 1000; % time in milliseconds %
subplot(2,1,2)
plot( t, h );
grid
xlabel('Time (ms)');
ylabel('Impulse Response');
xlim([ -0.5 20 ])