% 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 ])
