% FD scheme for the Webster equation
% Technologies for Singing workshop, October 2017, Detmold, Germany
% Charlotte Desvages

%% INITIALISATION

% Initial commands: make sure the workspace is clear.
clear all; % close all;
input       = 'pulse'; % pulse, audio
boreProfile = 'cyl';  % cyl, cone, exp, neck

% Flags: turn on/off certain features.
flag.Plots = 1;	% live plots on/off
flag.Sound = 0; % sound output
%    global StopButton;
global ButtonStatus
ButtonStatus = 0;
%    StopButton =0; % interrupts calculation

% Input parameters: simulation parameters, physical parameters.
	Fs = 22050;		% sample rate (Hz)
	Tf = 2.0;			% simulation duration (s)
	lambda = 1.00;		% Courant number
	
	L = 2;			% tube length (m)
	c = 340;			% speed of sound in the air (m/s)
	r0 = 0.01;		% tube radius at the mouthpiece/left end (m)

% Derived parameters: calculate useful coefficients from the input parameters.
	k = 1/Fs;						% time step (s)
	Nf = floor(Fs*Tf);	% total number of samples

	gamma = c/L;
	h = (gamma*k)/lambda;	% grid spacing (m)
	N = floor(L/h);			% number of grid points
	L = h*N;				% adjusted L, to fit an integer number of grid segments
	x = [0:h:L]';			% make a spatial vector with N grid points

% Building a bore profile: choose between 'cyl','cone','exp','neck'

	switch boreProfile

	case 'cyl'		% cylindrical tube
		r = r0*ones(N+1,1);

	case 'cone'		% conical tube
		rL = 2.0*r0;		% tube radius at the bell/right end (m)
		r = r0 + (rL - r0)*x;

	case 'exp'		% exponential tube
		flare = 10;			% how flared is the tube. try different coefficients!
		rL = 20*r0;			% tube radius at the bell/right end (m)
		A = (rL - r0)./(exp(flare*L) - 1);
		r = A*(exp(flare*x) - 1) + r0;

	case 'neck'		% cylindrical tube with a narrow neck
		neckDepth = 0.8;		% prominence of the neck, as a fraction of r0. Try different values between 0 and 1.
		neckPos = 0.7;			% location of the neck, as a fraction of L
		neckWidth = 0.2;		%	width of the neck, as a fraction of L

		neckSamp = round(neckWidth*(N+1));								% number of samples in the initial pulse
		neckStart = floor((neckPos-0.5*neckWidth)*(N+1));	% sample number of start of pulse
		neckEnd = neckStart + neckSamp - 1;								% sample number of end of pulse

		r = r0*ones(N+1,1);
		r(neckStart:neckEnd) = 0.5*r0*(2 - neckDepth + neckDepth*cos(2*pi*(0:neckSamp-1)'/neckSamp));

	otherwise
		error('Please select one of the suggested bore profiles.');
	end
	
	S0 = pi*r.^2;			% cross-section area along the tube (m2)
	S = [S0(1); (S0(3:N+1) + 2*S0(2:N) + S0(1:N-1))./4; S0(N+1)];	% discrete approximation to the bore profile (mu_tt S)
	
	alpha1 = 2.0881*L*sqrt(1/(S(1)*S(N+1)));	% radiation coefficients
	alpha2 = 0.7407*lambda;

% Compute scheme coefficients
	B = 2*(1 - lambda^2);
	C = 0.5*lambda^2*((S0(2:N) + S0(1:N-1))./S(2:N));
	D = 0.5*lambda^2*((S0(2:N) + S0(3:N+1))./S(2:N));
	E = -(gamma^2*k^2/(h*S0(1)))*(3*S0(1) - S0(2));
	F1 = alpha1*gamma^2*k^2*(1.5*S0(N+1) - 0.5*S0(N))/(h*S(N+1));
	F2 = alpha2*gamma^2*k*(1.5*S0(N+1) - 0.5*S0(N))/(h*S(N+1));
	G1 = 2*lambda^2/(1 + F1 + F2);
	G2 = -(1 + F1 - F2)/(1 + F1 + F2);

% Initialise arrays. We only need the full state for 3 time steps at a time.
	phi  = zeros(N+1,1);	% future time step			phi^n+1
	phi1 = zeros(N+1,1);	% current time step			phi^n
	phi2 = zeros(N+1,1);	% previous time step		phi^n-1

% Initialise our output vector, which we can save as a .wav file at the end.
	outPos = N+1;			% sample position of the readout location
	out = zeros(Nf,1);
	out(1) = phi2(outPos);	% fill out first 2 values with our initial conditions
	out(2) = phi1(outPos);

% Initialise the input waveform: choose between 'pulse', 'audio'
	
	switch input
	
	case 'pulse'	% single pulse
		pulseDepth = 0.1;		% amplitude of the pulse, as a fraction of r0
		pulseWidth = 0.002;	% duration of the pulse (s)

		pulseSamp = round(pulseWidth*Fs);									% number of samples in the initial pulse

		vIn = zeros(Nf,1);
		vIn(1:pulseSamp) = 0.5*pulseDepth*r0*(1 - cos(2*pi*(0:pulseSamp-1)'/pulseSamp));

	case 'audio'	% audio signal, from a .wav file
		%[vIn,FsIn] = audioread('TestGuitarPhraseMono.wav');
		[vIn,FsIn] = audioread('astrud.wav');
		if FsIn~=Fs
			error('Pick an audio file with the right sample rate, or change your simulation sample rate.');
		elseif size(vIn,2)~=1
			vIn = vIn(:,1);		% convert to mono by picking the left channel
		end

		vIn = 0.0001*r0*vIn;

	otherwise
		error('Please choose a valid input.');
	end

% Initialise figure
if flag.Plots
	fig1 = figure(1);
	ax1 = axes(fig1);
    clf
    maxrax=max(r);
    yyaxis right
    tubeDisp = plot(x,r,'k-',x,-r,'k-', 'LineWidth',2);
    ylim([-1.1*maxrax 1.1*maxrax])
    titlehdl=title(['Time: 0 ms'],'Color','black','FontSize',14);
end

%uibutton('Text', {'Stop','Simulation'}, 'Position',[100 100 100 40]);

%% MAIN LOOP

for n=1:Nf
	% Recursion equation
		phi(2:N) = B*phi1(2:N) + C.*phi1(1:N-1) + D.*phi1(3:N+1) - phi2(2:N);

	% Left boundary: excitation
		phi(1) = B*phi1(1) + 2*lambda^2*phi1(2) - phi2(1) + E*vIn(n);

	% Right boundary: loss
		phi(N+1) = G1*phi1(N) + G2*phi2(N+1);

	% Live plotting
	if flag.Plots
        yyaxis left
%		hold on
        maxlax=max(abs((phi-phi1)/k));
		pressureDisp = plot(x,(phi-phi1)/k,'Color','b','LineWidth',1.5);
%		hold off
		ylim([-maxlax-1e-10 maxlax+1e-10])
        xlabel('Distance (m)')
        ylabel('Pressure, Distance (Pa, m)')

        drawnow
        StopButton(fig1);
%        StopButton=cla(b);
    switch ButtonStatus
        case 1
%            disp('Running...');            
        case 2
            disp('Exiting loop ...');            
            break;
        case 3
            disp('Press key to continue');
            pause;
            ButtonStatus = 1;
    end
    titlehdl.String=['Time: ', num2str(n/Fs*1000,'%.1f'),' ms'];
		%pause
	end

	% Readout pressure
	out(n) = (phi(outPos) - phi1(outPos))/k;

	% Update the state of the potential
	phi2 = phi1; phi1 = phi;
end

%% PLOT RESULTS

% Compute and plot spectrum of the output
figure(1);
clf
plot([0:length(out)-1]/Fs,out)
xlabel('Time (s)')
ylabel('Amplitude (Pa)')

figure(2);
outSpectrum = fft(out);
spectrum = plot([0:Nf-1]./(Nf*k),20*log10(abs(outSpectrum)),'k');
axis([0 Fs/2 -80 0]);
xlabel('Frequeny (Hz)')
ylabel('Level (dB)')

if flag.Sound
    soundsc(out,Fs)
end
function StopButton(FigureHandle)
global ButtonStatus
b = uicontrol(FigureHandle,'Style','popupmenu');
b.Position = [60 50 50 20];
b.String = {'Run','Stop','Pause'};
b.Callback = @selection;

    function selection(src,event)
        val = b.Value;
        str = b.String;
        str{val};
        disp(['Selection: ' str{val}]);
        ButtonStatus = val;
    end
%ButtonStatus = c.Callback;
end

%function StopButton=cla(b)
%global StopButton
% hObject    handle to pushbutton1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
%    display('Goodbye');
%    StopButton=get(b,'Value');
%end