% Pipe.m; resonance of a tube
clear
T = input('Enter air temperature in degrees Fahrenheit: ');
T = (5/9)*(T - 32) + 273; % convert to Kelvin
v = 20.1*sqrt(T); % calculate sound speed in air
L = input('Enter length of tube in inches: ');
L = L*.0254; % convert to m
N = input('Enter number of harmonics to sum: ');
rate = 1000; % sample rate (22254 for higher fidelity)
t = linspace(0,1,rate); % use one second of sound
accum = zeros(1,length(t));
OC = input('Enter O for open-open tube or C for closed-open tube: ','s');
i = 1;
while i<=N
if OC == 'O' % open-open tube
n = i;
nstring = int2str(n);
f = n*v/(2*L); % calculate each frequency
fstring = num2str(f); % convert to string
ampl=input(['Enter relative amplitude of harmonic ',...
nstring,' at ',fstring,' Hz: ']);
y = ampl*sin(2*pi*f*t); % generate the appropriate sinusoid
accum = accum + y; % accumulate sum
sound(accum,rate); % play accumulated sound
pause(0.5);
i = i + 1; % next harmonic
elseif OC == 'C' % open-closed tube
n = 2*i - 1; % odd harmonics only
nstring = int2str(n);
f = n*v/(4*L); % calculate each frequency
fstring = num2str(f); % convert to string
ampl=input(['Enter relative amplitude of harmonic ',...
nstring,' at ',fstring,' Hz: ']);
y = ampl*sin(2*pi*f*t); % generate the appropriate sinusoid
accum = accum + y; % accumulate sum
sound(accum,rate); % play accumulated sound
pause(0.5);
i = i + 1; % next harmonic
else disp('Not a valid choice for a tube!');
i = N; % end loop
end
end