% ZinWaveguide.m
%
% This matlab script calculates the impedance of a cylindrical duct as
% modeled by a digital waveguide system.  The open end of the duct is
% terminated with a digital filter approximation to the analytic
% solution of Levine and Schwinger.  Thermo-viscous losses along the
% pipe walls are ignored.  No delay line interpolation is implemented.
%
% by Gary P. Scavone

% Physical Constants and evaluation frequencies

N = 20000;
fs = 44100;
c = 347.23;    % speed of sound (m/sec)
rho = 1.1769;  % density of air (kg/m^3) %

% Pipe radius (m)

ra = 0.016 / 2;

% Length of pipe (m)

L = 0.35;

% Delay Line Length

delay = round( L * fs / c);

% Pipe characteristic impedance

Zc = rho * c / (pi * ra^2);

% Open-end Reflection Filter - based on Levine & Schwinger results
% and designed using 'openpipe' matlab function from scavone
% thesis.  The following values were found for a radius of 8 mm.
% As me for the "openpipe" script if you want to experiment with
% other sizes.

b = 0.98 * [-0.221 -0.4108];
a = [1.0 -0.3801 0.0119];
z_l = [0 0];

% Initialize delay line and output vector

dp = zeros(1, delay);
dm = zeros(1, delay);
pointer = 1;
p = zeros(1,N);
u = [Zc, zeros(1, N-1)];

% *************************************** %
%                                         %
%             Delay Line(d)               %
%    |-------------------------------|    %
%    ^                                    %
%  pointer                                %
%                                         %
%     >>--- pointer increments --->>      %
%                                         %
% *************************************** %
%
% The pointer initially points to the delay line output.
% We can take the output and calculate a new input value
% which is placed where the output was taken from.  The
% pointer is then incremented and the process repeated.

% ******* Run Loop Start ******* %

for i = 1:N,

  p(i) = dm(pointer);          % p- traveling into input
  
  % Open-end filtering
  [dm(pointer), z_l] = filter(b, a, dp(pointer), z_l);
  
  dp(pointer) = p(i) + u(i);   % new p+ @ input
  p(i) = p(i) + dp(pointer);   % physical pressure @ input

  pointer = pointer + 1;       % increment "pointer"
	
  if pointer > delay           % check delay-line limits
    pointer = 1;
  end
  
end

% Clear Figure Window and Plot
clf
subplot(2,1,1)
T = 1/fs;
t = (0: T : (N-1)*T) * 1000;
plot( t, p / Zc )
grid
xlabel('Time (ms)')
ylabel('Gain / Zc')
title('Impulse Response of Pipe')
xlim([0 20]);

subplot(2,1,2)
f = 0: fs/(N-2) : fs/2;
P = fft( p );
plot( f, abs( P(1:N/2) / Zc ) )
grid
xlabel('Frequency (Hz)')
ylabel('Magnitude / Zc')
title('Input Impedance of Pipe')
v = axis;
axis( [0 10000 v(3) v(4) ] );

%sound(y, fs)
