function [ experimentTag ] = waymarkModelMain( experimentName, srcData, environmentName, trajectoryFunction, surfaceParams, surfaceFunction, simulationTime, waymarkPeriod, misc)
%WAYMARKMODELMAIN This is the top level for the Waymark channel model
% This is the main loop for the Waymark channel model as described in [1] and [2].

% References
% ==========
% [1] C. Liu, Y. V. Zakharov, and T. Chen, ?Doubly Selective Underwater Acoustic Channel Model for a Moving Transmitter/Receiver,? IEEE Transactions on Vehicular Technology, vol. 61, no. 3, pp. 938?950, 2012.
% [2] B. Henson, J. Li, Y. V. Zakharov, and C. Liu, ?Waymark baseband underwater acoustic propagation model,? in Underwater Communications and Networking (UComms), 2014. IEEE, 2014, pp. 1?5.
% [3] M. B. Porter, ?The bellhop manual and user?s guide: Preliminary draft,? Heat, Light, and Sound Research, Inc., La Jolla, CA, USA, Tech. Rep, 2011.
% [4] (2010, Accessed 3, December, 2016) Heat, Light & Sound Research, Inc. website. [Online]. Available: http://oalib.hlsresearch.com/Modes/AcousticsToolbox/

% This code is based code written by Chunshan Liu.
% The Bellhop ray tracing acoustic simulation is from the acoustic toolkit written by Michael Porter [3][4].

% Start stopwatch timer
timerWaymark = tic;

% Setup Parameters
% ================
% Experiment setup
experimentTag = [experimentName, '_', srcData.waveformName,'_', datestr(now, 'yyyymmddTHHMMSS')];
mkdir(experimentTag);
experimentTagDir = ['./', experimentTag, '/'];

% Input waveform
% ==============
% Read the information from the input audio file, check the sample rate.  We will need the TotalSamples later.
disp(' ')
disp('=== Input waveform ===')
disp(['Loading ' srcData.waveformName '.' srcData.waveformExt ' source waveform']);
srcInfo			= audioinfo([srcData.waveformName '.' srcData.waveformExt]);
[inputData, ~] = audioread([srcData.waveformName '.' srcData.waveformExt]);

disp(['Sample rate:   ' num2str(srcInfo.SampleRate) 'Hz'])
disp(['Sample period: ' num2str(1/srcInfo.SampleRate) 's'])
disp(['Length:        ' num2str(srcInfo.TotalSamples(1))])

if srcData.fSample ~= srcInfo.SampleRate
	warning('WarnSetup:fSample', 'The Fs of the input file and the Fs specified do not agree, using Fs of the input file');
end
F_SAMPLE		= srcInfo.SampleRate;									% Specify the sample rate for the system
T_SAMPLE		= 1/F_SAMPLE;

FREQ_SHIFT	= srcData.fCarrier;

% Misc. environment parameters
% ============================
fSrcBaseBW			= srcData.bandwidth/2;
maxDoppler			= (misc.RXTX_SPEED_MAX/misc.c0)*(srcData.fCarrier + fSrcBaseBW);	% Calculate the maximum Doppler shift expected at the highest passband frequency

% Decimation parameter
% ====================
% Calculate how much the signal may be decimated by.
disp(' ')
disp('=== Down shift and Decimation parameters ===');
disp(['Source carrier frequency:   '	num2str(srcData.fCarrier) 'Hz']);
disp(['Source bandwidth:           '	num2str(srcData.bandwidth) 'Hz']);
% Divide the sampling frequency by the baseband bandwidth and add a 100% margin to make sure not to attenuate the signal with the filter.
F_DECI_MULT		= 4;														% The margin between the baseband bandwidth and the new decimated sampling frequency 
DECI_STEP		= floor(F_SAMPLE/(F_DECI_MULT*fSrcBaseBW));	% Calculate the decimation step
F_DECI_SAMPLE	= F_SAMPLE/DECI_STEP;								% Decimated sampling frequency
T_DECI_SAMPLE	= 1 / F_DECI_SAMPLE;									% Decimated sampling period
disp(['Decimation step:            ' num2str(DECI_STEP)]);
disp(['Decimated sample frequency: ' num2str(F_DECI_SAMPLE) 'Hz']);
disp(['Decimated sample period:    ' num2str(T_DECI_SAMPLE) 's']);

% Waymark interval
% ================
disp(' ')
disp('=== Waymark interval ===')
% Round the Waymark period down to the nearest number of decimated sample periods
NUM_DECI_SAMPLES_WAYMARK	= floor(waymarkPeriod/T_DECI_SAMPLE);
T_WAYMARK						= NUM_DECI_SAMPLES_WAYMARK*T_DECI_SAMPLE;
NUM_SAMPLES_WAYMARK			= floor(T_WAYMARK/T_SAMPLE);
disp(['Input Waymark interval:                      ' num2str(waymarkPeriod) 's']);
disp(['Decimated samples in rounded waymark period: ' num2str(NUM_DECI_SAMPLES_WAYMARK)]);
disp(['Rounded waymark period:                      ' num2str(T_WAYMARK) 's']);

% Raised cosine parameters
% ========================
disp(' ')
disp('=== Low pass filter parameters ===')
F_FILT_MULT				= 2+(1+misc.ROLL_OFF);													% The margin between the baseband bandwidth and the corner frequency of the filter
SAMPLES_IN_SYMBOL		= floor(F_SAMPLE/(F_FILT_MULT*fSrcBaseBW));
F_SYMBOL					= F_SAMPLE/SAMPLES_IN_SYMBOL;
FILT_LEN_MULT			= 5;																	% Make the total filter length a multiple of the SAMPLES_IN_SYMBOL

NUM_POINTS_DECI_SMP	= ceil(SAMPLES_IN_SYMBOL*FILT_LEN_MULT/DECI_STEP);		% Calculate the number of decimated samples needed for the filter, this is for the sampling at the output
NUM_POINTS				= (NUM_POINTS_DECI_SMP*DECI_STEP)+1;						% The number of points to be generated for the filter.
if rem(NUM_POINTS,2)==0																		% If the number of points is even then add one so that there is a centre sample.
    NUM_POINTS = NUM_POINTS + 1;
end
% Generate the raised-cosine filter
% =================================
timeRaisedCos	=  -(T_SAMPLE*(floor(NUM_POINTS/2))) : T_SAMPLE : (T_SAMPLE*(floor(NUM_POINTS/2)));	% Generate the timeRaisedCos steps to generate the raised-cosine.  Starts at T_SAMPLE because t=0 element is added below
raisedCos		=	((sin(pi * F_SYMBOL * timeRaisedCos) ./ (pi * F_SYMBOL * timeRaisedCos))) .* (cos(pi() * F_SYMBOL * misc.ROLL_OFF * timeRaisedCos)) ./ (1 - ((2 * F_SYMBOL * misc.ROLL_OFF * timeRaisedCos).^2) + eps);			% +eps to avoid divide by zero error
raisedCos(timeRaisedCos == 0) = 1;
raisedCos		= raisedCos/sum(raisedCos, 2);

disp(['Symbol rate:            ' num2str(F_SYMBOL) 'Hz']);
disp(['Num samples per symbol: ' num2str(SAMPLES_IN_SYMBOL)]);
disp(['Filter length:          ' num2str(NUM_POINTS)]);

T_LOWPASS_NEG_START				= ((NUM_POINTS-1)/2)*T_SAMPLE;					% Calculate the time added onto the signal as it passes through the decimating low pass filter
WAYMARKS_LOWPASS_NEG_START		= T_LOWPASS_NEG_START/T_WAYMARK;					% Calculate the length of half the lowpass filter in units of waymark time, this is the number of extra Waymarks needed to accomodate the filter

disp(['Filter negative start (Waymarks): ' num2str(WAYMARKS_LOWPASS_NEG_START)]);

% Channel impulse response parameters
% ===================================
disp(' ')
disp('=== Channel impulse response ===')
T_IMPULSE_MIN					= 3;															% Desired impulse response time in seconds
NUM_DECI_SMP_IMP				= ceil(T_IMPULSE_MIN/T_DECI_SAMPLE);				% Record the number of decimated samples in the impulse response, rounding up to the nearest sample
T_IMPULSE						= NUM_DECI_SMP_IMP*T_DECI_SAMPLE;					% Find the duration of the rounded impulse response
NUM_DECI_SMP_IMP_NEG_START	= floor((T_IMPULSE*0.1)/T_DECI_SAMPLE);			% Move the time of interest back to around 10% of T_IMPULSE but make it a round number of decimatied samples so that zero time is aligned.
T_IMPULSE_NEG_START			= NUM_DECI_SMP_IMP_NEG_START*T_DECI_SAMPLE;		% Calculate the time that the T_IMPULSE was moved back
NUM_WAYMARK_IMPULSE_LAG		= 1+1+floor(NUM_DECI_SMP_IMP_NEG_START/NUM_DECI_SAMPLES_WAYMARK);	% Num. of Waymarks that the impulse response calculation needs to lag the input sample generation so that the T_IMPULSE_NEG_START
																															% can be accommodated.  
																															% +1, the calculation of the impulse response convolution is done in Waymark blocks, therefore there needs to be an extra one for the spline support.
																															% +1+floor(), this is for consistancy in the case where there is an exact number of decimated samples in the impulse response.
NUM_WAYMARK_IMPULSE			= 1+floor((NUM_DECI_SMP_IMP-NUM_DECI_SMP_IMP_NEG_START)/NUM_DECI_SAMPLES_WAYMARK) + NUM_WAYMARK_IMPULSE_LAG;	% Find the number of decimated samples needed for the impulse response,
																																														% this inclused the rounding for the negative start and then the rounding for the rest
																																														% +1+floor(), this is for consistancy in the case where there is an exact number of decimated samples in the impulse response.
disp(['Impulse response length:                        ' num2str(T_IMPULSE) 's']);
disp(['Decimated samples in rounded Impulse response:  ' num2str(T_IMPULSE/T_DECI_SAMPLE)]);
disp(['Impulse response length before first arrival:   ' num2str(T_IMPULSE_NEG_START) 's']);
disp(['Decimated samples before first arrival:         ' num2str(NUM_DECI_SMP_IMP_NEG_START)]);
disp(['Waymark periods before first arrival:           ' num2str(NUM_DECI_SMP_IMP_NEG_START/NUM_DECI_SAMPLES_WAYMARK) ' (' num2str(NUM_WAYMARK_IMPULSE_LAG) ')']);

% Start time adjustment
% =====================
disp(' ')
disp('=== Processing length and start adjustment ===')
% This section calculates the start time for the trajectory and the amount that would be added to the signal to make it up to a round number of Waymarks in length

T_SIGNAL						= simulationTime+(T_SAMPLE*(NUM_POINTS-1)+(T_SAMPLE*T_IMPULSE_NEG_START));			% Time to run the experiment for.  This is the time for the simulation, the time through the rasied cosine filter and the delay in the channel filter
SINGLE_SPLINE_SUPPORT	= 3;																			% The number of data points required to build a single spline
SPLINE_SUPPORT				= 2*SINGLE_SPLINE_SUPPORT;												% The number of data points required to build enough splines to interpolate a point
WAYMARK_SPLINE_PAD.pre	= SINGLE_SPLINE_SUPPORT + ceil(WAYMARKS_LOWPASS_NEG_START);	% Number of waymarks added before the signal to allow for spline interpolation and the filter start
WAYMARK_SPLINE_PAD.post	= SINGLE_SPLINE_SUPPORT + 1;											% Number of waymarks added after the signal to allow for spline interpolation
NUM_WAYMARKS				= ceil(T_SIGNAL/T_WAYMARK) + WAYMARK_SPLINE_PAD.pre + WAYMARK_SPLINE_PAD.post +1;	% Number of waymarks to generate, +1 for the zeroth element
PROP_DELAY_WM_NUM			= 1;																											% Specify which waymark is to be used initially for the common propagation delay
disp(['Total num. Waymarks:             ' num2str(NUM_WAYMARKS)]);
disp(['Num. Waymarks added at begining: ' num2str(WAYMARK_SPLINE_PAD.pre)]);
disp(' ')

if (WAYMARK_SPLINE_PAD.pre > SPLINE_SUPPORT)
	warning('WarnSetup:delayCalc', 't=0 is not reached in the Waymark-spline initalisation.  The delays will relative to the first Waymark not the Waymark at t=0');
end

% B-Spline generation
% ===================
% B-Spline generation, calculate the four quarters of the spline into a matrix for later use
t = [-2 -1 0 1 2];
pp = bspline(t);
bSplineMat = zeros(NUM_SAMPLES_WAYMARK, 4);
for i = 1:4
	bSplineMat(:,i) = (fnval(pp,((2-i) : 1/NUM_SAMPLES_WAYMARK : (3-i)-1/NUM_SAMPLES_WAYMARK))).';
end

% Input Data weighting matrix
SPLINE_WEIGHTS_MAT	 = [	-1/6  4/3 -1/6    0	  0    0 ;
										0 -1/6  4/3 -1/6    0    0 ;
										0    0 -1/6  4/3 -1/6    0 ;
										0    0    0 -1/6  4/3 -1/6 ].';


% Frequency response generation
% =============================
% Frequency sweep generation, this is used to find the delays for each waymark at the baseband frequency.  This should accommodate the anticipated Doppler shift at the carrier frequency.
FREQ_SWEEP_SEARCH.freqLow	= -1024;
FREQ_SWEEP_SEARCH.freqHigh	= 1023;
FREQ_SWEEP_SEARCH.res		= 8;
FREQ_SWEEP_SEARCH.freqStep	= 1/(FREQ_SWEEP_SEARCH.res*T_IMPULSE);													% The frequency step needs to be small enough to capture the lowest frequency that can be represented by within the impulse response time.
freqSearch = (FREQ_SWEEP_SEARCH.freqLow : FREQ_SWEEP_SEARCH.freqStep : FREQ_SWEEP_SEARCH.freqHigh);	% Generate the frequency sweep values

% Field Calculation
% =================
% Store the results from the field calculation long enough to be used for the impulse response interpolation
SIGNAL_CALC_DELAY		= 3;									% This is the number of Waymarks that it takes to calculate the delay compensated signal
FIELD_STORE_LEN		= NUM_WAYMARK_IMPULSE_LAG;
fieldStore.delay		= zeros(FIELD_STORE_LEN,[]);	% Delays of the multipath
fieldStore.amp			= zeros(FIELD_STORE_LEN,[]);	% Complex amplitude of the multipath
DELAY_STORE_LEN		= max([(NUM_WAYMARK_IMPULSE_LAG) SPLINE_SUPPORT]);
delayComp				= zeros(1, DELAY_STORE_LEN);	% Store enough of the waymark delays for the spline calculation
delayCompT0				= 0;									% Store the delay compensation to Waymark t=0, this is so that it can be applied as an offset so that tfirst sample is aligned with t=0
fieldPropDelayFirst	= 0;									% Store the group propagation delay to the first Waymark, this is for the calculation of the delays
fieldPropDelayT0		= 0;									% The group propagation delay associated with the first arrival at Waymark 0, this is for building the impulse responses so that they align with t=0

% Input signal delay compensation
% ===============================
INPUT_BUFF_LEN				= 5;
SAMPLE_SPLINE_WEIGHTS	= [0 -1/6 4/3 -1/6 0];
SAMPLE_SPLINE_SUPPORT	= size(SAMPLE_SPLINE_WEIGHTS, 2);			% The number of data points required to build a single spline for the sample delay interpolation 

% Impulse response generation
% ===========================
% Frequency sweep generation, this is used to construct the impulse response.  This should accomodate the anticipated Doppler shift at the carrier frequency.
FREQ_SWEEP.freqLow	= -(fSrcBaseBW + maxDoppler);
FREQ_SWEEP.freqHigh	= (fSrcBaseBW + maxDoppler);
FREQ_SWEEP.res			= 8;
FREQ_SWEEP.freqStep	= 1/(FREQ_SWEEP.res*T_IMPULSE);								% The frequency step needs to be small enough to capture the lowest frequency that can be represented by within the impulse response time.
freq = (FREQ_SWEEP.freqLow : FREQ_SWEEP.freqStep : FREQ_SWEEP.freqHigh);	% Generate the frequency sweep values

timeChannel	= (0 : T_DECI_SAMPLE : T_IMPULSE-T_DECI_SAMPLE) - T_IMPULSE_NEG_START;					% Defining the length of channel impulse response (after delay adjustment).  Move the period back to before 0 time.
impulseResp	= zeros(length(timeChannel), SPLINE_SUPPORT);
invFftMatix	= 2*(1/FREQ_SWEEP.res)*(1/size(timeChannel, 2))*exp(1i*2*pi*(freq.'*timeChannel));	% Matrix used for inverse Fourier transform, specifically for the frequencies of interest

% Initialise control variables and constants
% ==========================================
freqRespWaymarkA				= [];
signalDeci						= zeros(1, NUM_DECI_SAMPLES_WAYMARK);
NUM_SIG_DECI_BLOCKS			= 3;															% The number of decimated signal waymark blocks to store, this is three because the delay may be +ve or -ve and so the 
																									% interpolation of the central block may require samples from either block
signalDelayT0Block			= zeros(1, (1+NUM_WAYMARK_IMPULSE_LAG)*NUM_SAMPLES_WAYMARK);
signalDComp						= zeros(1, NUM_DECI_SAMPLES_WAYMARK);
signalDCompStore				= zeros(1, NUM_WAYMARK_IMPULSE*NUM_DECI_SAMPLES_WAYMARK);
NUM_WAYMARK_HALF_POINTS		= (ceil(ceil(NUM_POINTS_DECI_SMP/2)/NUM_DECI_SAMPLES_WAYMARK));
NUM_SIG_CHNL_BLOCKS			= 1 + 2*NUM_WAYMARK_HALF_POINTS;												% Calculate the number of decimated sample Waymark blocks that need to be kept to fill half the low pass filter
																															% either side of the block of interest.  The low pass filter is in output samples so we need to makes sure there
																															% is enough decimated samples.  So 1 for the centre block and two rounded up halfs
signalChnl						= zeros(1, NUM_DECI_SAMPLES_WAYMARK);
signalChnlBlock				= zeros(1, NUM_SIG_CHNL_BLOCKS*NUM_DECI_SAMPLES_WAYMARK);

% Upsample and interpolate
% ========================
interpSignalUp			= zeros(1, NUM_SAMPLES_WAYMARK);
INTERP_STEP				= DECI_STEP;
WAYMARK_UP_LAG			= SIGNAL_CALC_DELAY + NUM_WAYMARK_IMPULSE_LAG;								% The lag for the waymarks that is being calculated as the output
interpSignalUShift	= [];

% Write out the results
% =====================
N_BITS        = 32;		% The bit depth for each wave file sample



% =========
% Main Loop
% =========
% This is the main loop running the simulation in steps of a Waymark period
for waymarkIndex = -WAYMARK_SPLINE_PAD.pre : NUM_WAYMARKS
	

	% ==========================
	% Pressure field calculation
	% ==========================
	% Generate the time for the current waymark, this is relative the the Waymark at t=0
	rxTxTime	= waymarkIndex * T_WAYMARK;
	
	% Waymark trajectory
	% ------------------
	% Calculate the Waymark trajectory point

	[rxTxTraj] = trajectoryFunction(rxTxTime);
		% rxTxTraj is a structure to hold the trajectory data.  Each run of Bellhop accepts a rx depth, a tx depth, a range (this is a much simplified way to use Bellhop)
		% rxTxTraj.rxDepth
		% rxTxTraj.txDepth
		% rxTxTraj.rxTxRange

	% Waymark surface
	% ---------------
	% Generate the surface for each waymark
	[surface] = surfaceFunction(rxTxTime, surfaceParams);


	% Pressure field calculation
	% --------------------------
	% Run external program to build pressure field
	[field] = waymarkFieldGenBellhop(rxTxTraj, environmentName, surface);
	

	
	% Store the field calculation in a buffer ready to be used by the impulse response interpolation
	% This buffer is required because the impulse response calculation needs to happened after the input signal has been generated
	fieldStore.delay(1:end-1,:)							= fieldStore.delay(2:end, :);
	fieldStore.delay(end, 1:numel(field.delay))		= field.delay(1,:);
	
	fieldStore.amp(1:end-1,:)						= fieldStore.amp(2:end, :);
	fieldStore.amp(end, 1:numel(field.amp))	= field.amp(1,:);


	% Calculate the common propagation delay, firstly from the first Waymark, then from the Waymark at t=0
	if (waymarkIndex == -WAYMARK_SPLINE_PAD.pre)
		fieldDelayNoZeros			= field.delay(PROP_DELAY_WM_NUM, :);
		fieldDelayNoZeros(fieldDelayNoZeros == 0)	= [];						% Remove any zeros delays, this can happen if there is a varying number of multipath components
		fieldMinDelay				= min(fieldDelayNoZeros,[],2).';			% For each waymark find the minimum delay, only delays greater than zero should be included from the Field generation
		fieldPropDelayFirst		= floor(fieldMinDelay/T_SAMPLE)*T_SAMPLE;		% Record the common propagation time (the time to the first waymark), this will be subtracted from arrival delays.  Make the delay a round number of input samples.
		fieldPropDelayT0			= floor(fieldMinDelay/T_SAMPLE)*T_SAMPLE;		% Initialise fieldPropDelayT0 to the first delay in order to create the initial impulse responses.  Make the delay a round number of input samples
		disp(['Common propagation delay for Waymark ' num2str(waymarkIndex) ': ' num2str(fieldPropDelayFirst) 's']);		
	end
	
	[freqRespWaymarkA, waymarkDelay] = waymarkDelayCalc(	misc.c0, ...
																			misc.RXTX_SPEED_MAX, ...
																			(waymarkIndex == -WAYMARK_SPLINE_PAD.pre), ...
																			freqRespWaymarkA, ...
																			field, ...
																			fieldPropDelayFirst, ...
																			T_DECI_SAMPLE, ...
																			T_WAYMARK, ...
																			freqSearch);


	delayComp	= [delayComp(2:end) (waymarkDelay+delayComp(end))];	% Store enough of the waymark delays for the spline calculation
																							% Store the composite delay (i.e. the  running delay total), this being the relative delay from the first Waymark

	if (waymarkIndex == 0)		% The next Waymark result after the pre-Waymarks have been calculated
		fieldDelayNoZeros			= field.delay(PROP_DELAY_WM_NUM, :);
		fieldDelayNoZeros(fieldDelayNoZeros == 0)	= [];						% Remove any zeros delays, this can happen if there is a varying number of multipath components
		fieldMinDelay				= min(fieldDelayNoZeros,[],2).';			% For each waymark find the minimum delay, only delays greater than zero should be included from the Field generation
		fieldPropDelayT0			= floor(fieldMinDelay/T_SAMPLE)*T_SAMPLE;		% Record the common propagation time (the time to the Waymark 0), this will be subtracted from arrival delays when calculatin the impulse response.  Make the delay a round number of input samples
		disp(['Common propagation delay for Waymark ' num2str(waymarkIndex) ': ' num2str(fieldPropDelayT0) 's']);	

		delayCompT0	= delayComp(end);												% Record the composite delay at the Waymark corresponding to t=0
		disp(['Delay compensation relative to Waymark ', num2str(waymarkIndex), ' at ', num2str(delayCompT0) 's']);	
	
	end



	% ===============================
	% Input signal delay compensation
	% ===============================
	% Interpolate delays at decimate sample rate
	% ------------------------------------------
	% Use the local B-Spline to interpolate the delay for each sample
	signalDelay				= bSplineMat*(delayComp(end-(SPLINE_SUPPORT-1):end)*SPLINE_WEIGHTS_MAT)';

	signalDelayT0			= signalDelay - delayCompT0;																			% Readjust delays to be relative to the Waymark associated with t=0
	signalDelayT0Block	= [signalDelayT0Block(NUM_SAMPLES_WAYMARK+1 : end) signalDelayT0'];
	
	signalDelayDeciSamp	= downsample(signalDelayT0, DECI_STEP)/T_DECI_SAMPLE;											% Get the delays in units of decimated samples

	[~, sigDelayWaymarkOffsetMax]	= max(abs(fix((-1*signalDelayDeciSamp)/NUM_DECI_SAMPLES_WAYMARK)));						% Work out if the delay will be over the number of decimated samples in the Waymark, first look at the maximum movement,
	sigDelayWaymarkOffset			= fix(-1*signalDelayDeciSamp(sigDelayWaymarkOffsetMax)/NUM_DECI_SAMPLES_WAYMARK);	% then get the direction and round towards zero
	signalDelayDeciSampOff			= signalDelayDeciSamp+(sigDelayWaymarkOffset*NUM_DECI_SAMPLES_WAYMARK);					% Find the offset for the Waymarks to select for the signal generation



	% =====================================
	% Input signal downshift and decimation
	% =====================================
	% Input sample window
	signalWinStart	= ((waymarkIndex-INPUT_BUFF_LEN+sigDelayWaymarkOffset)*NUM_SAMPLES_WAYMARK) - ((NUM_POINTS-1)/2);			% From the Waymark count calculate the start of the window for the input data, (NOTE: The Waymark number is reduced so that
																																									% the delay in calculating the sample delays is taken into account
	signalWinEnd	= ((waymarkIndex+sigDelayWaymarkOffset)*NUM_SAMPLES_WAYMARK) + ((NUM_POINTS-1)/2) - 1;							% Calculate the end of the window for the input, -1 because it is up to and not including the next Waymark


	if	((signalWinEnd < 0) || (signalWinStart >= srcInfo.TotalSamples(1)))						% If the index is not in the file then zero the input
		inputSignalWin = zeros(1, (INPUT_BUFF_LEN*NUM_SAMPLES_WAYMARK) + (NUM_POINTS-1));	% Zero the input data buffer
	else																												% else, calculate the file index.  The signal point counts from zero, so add one to make it into an index into the file.
		if (signalWinStart < 0)
			fileIndexStart	= 1;
			inputSignalWinStart = -signalWinStart;										% Sort out where to place the loaded data in the input buffer
		else
			fileIndexStart	= signalWinStart+1;
			inputSignalWinStart = 0;
		end
		
		if (signalWinEnd >= srcInfo.TotalSamples(1))
			fileIndexEnd	= srcInfo.TotalSamples(1);
		else
			fileIndexEnd	= signalWinEnd+1;
		end

%		% Read the input .wav here if the input is streamed
%		[inputData, ~] = audioread([srcData.waveformName '.' srcData.waveformExt], [fileIndexStart, fileIndexEnd]);	% Load the window of data from the file
		inputSignalWin = [zeros(1, inputSignalWinStart) inputData(fileIndexStart : fileIndexEnd, 1)' zeros(1, ((INPUT_BUFF_LEN*T_WAYMARK)/T_SAMPLE) + (NUM_POINTS-1) - (fileIndexEnd-fileIndexStart+1) - inputSignalWinStart)];	% Arrange the data in the input buffer,
																																																																										% accommodating a short file and the begining and end 
	end

	timeIpSignal = signalWinStart*T_SAMPLE : T_SAMPLE : signalWinEnd*T_SAMPLE;
	signalDShift = inputSignalWin.*exp( -1i*2*pi*FREQ_SHIFT*timeIpSignal );		% Shift the spectrum of the test signal down by the carrier frequency

	% Pass the input signal through the filter and decimate
	for i = 1 : INPUT_BUFF_LEN*NUM_DECI_SAMPLES_WAYMARK
		signalDeci(i) = signalDShift(	((i-1)*DECI_STEP)+1 : (((i-1)*DECI_STEP)+1 + size(raisedCos,2)-1)) * raisedCos.';
	end

	signalDeciBlock = signalDeci;							% Add new decimated input samples to buffer

	% This section could be made more efficient by storing the previously calculated delayed samples
	% signalDeciBlock = [signalDeciBlock(NUM_DECI_SAMPLES_WAYMARK+1 : NUM_SIG_DECI_BLOCKS*NUM_DECI_SAMPLES_WAYMARK) signalDeci];							% Add new decimated input samples to buffer



	% ==========================
	% Delay spline interpolation
	% ==========================
	for i = 1 : NUM_DECI_SAMPLES_WAYMARK
		signalDComp(i) = signalBSpline(	signalDeciBlock(		((floor(INPUT_BUFF_LEN/2)*NUM_DECI_SAMPLES_WAYMARK)+i - (SAMPLE_SPLINE_SUPPORT-1) + floor(-1*signalDelayDeciSampOff(i))) : ((floor(INPUT_BUFF_LEN/2)*NUM_DECI_SAMPLES_WAYMARK)+i - (SAMPLE_SPLINE_SUPPORT-1) + floor(-1*signalDelayDeciSampOff(i)) + (2*SAMPLE_SPLINE_SUPPORT)-1)), ...
													SAMPLE_SPLINE_WEIGHTS, ...
													mod(-1*signalDelayDeciSampOff(i),1) ...
												);
	end

	signalDCompStore = [signalDComp(end : -1 :1) signalDCompStore(1 : (NUM_WAYMARK_IMPULSE-1)*NUM_DECI_SAMPLES_WAYMARK)];		% Store the generated delay compensated input signals until we have got enough for the negative start on the impulse response,
																																										% reverse the order for convolution


	% ===========================
	% Impulse response generation
	% ===========================
	freqRespCompDelay			= fieldStore.amp(1, :) * exp(-1i*2*pi*(fieldStore.delay(1, :) - fieldPropDelayT0 - (delayComp(end-(NUM_WAYMARK_IMPULSE_LAG-1))-delayCompT0)).'*(freq+FREQ_SHIFT));	% Calculate the impulse response for this Waymark,
	impulseResp(:,1:end-1)	= impulseResp(:,2:end);
	impulseResp(:,end)		= (freqRespCompDelay*invFftMatix).';																																											% and store it for the interpolation

	impulseRespInterpWeight = impulseResp*SPLINE_WEIGHTS_MAT;

	for i = 1 : NUM_DECI_SAMPLES_WAYMARK
		impulseRespInterp = bSplineMat(((i-1)*DECI_STEP)+1,:)*impulseRespInterpWeight.';
		% Convolve the individually delayed input samples with the interpolated impulse response, move to the start point which is part way through a Waymark and an extra Waymark as the index move back through the signal
		% The (NUM_DECI_SAMPLES_WAYMARK*2) is to move on an extra waymark, the first is for the spline support and therefore the data position is moved on one.
		signalChnl(i) = signalDCompStore((NUM_DECI_SAMPLES_WAYMARK*2) - mod(NUM_DECI_SMP_IMP_NEG_START, NUM_DECI_SAMPLES_WAYMARK) - (i-1): (NUM_DECI_SAMPLES_WAYMARK*2) - mod(NUM_DECI_SMP_IMP_NEG_START, NUM_DECI_SAMPLES_WAYMARK) + NUM_DECI_SMP_IMP-1  - (i-1))*impulseRespInterp.';
	end
	
	signalChnlBlock = [signalChnlBlock(NUM_DECI_SAMPLES_WAYMARK+1 : NUM_SIG_CHNL_BLOCKS*NUM_DECI_SAMPLES_WAYMARK) signalChnl];							% Add new decimated input samples to buffer







	% ========================
	% Upsample and interpolate
	% ========================
	signalChnlUp		= upsample(signalChnlBlock, INTERP_STEP);

	for i = 1 : NUM_SAMPLES_WAYMARK
		interpSignalUp(i) = signalChnlUp((NUM_SAMPLES_WAYMARK-((NUM_POINTS-1)/2)) + i : ((NUM_SAMPLES_WAYMARK-((NUM_POINTS-1)/2)) + i + NUM_POINTS-1)) * (raisedCos*INTERP_STEP).';
	end

	% Calculate the output sample time delay
	timeOpSignal =  ((waymarkIndex-WAYMARK_UP_LAG)*T_WAYMARK : T_SAMPLE : ((waymarkIndex-WAYMARK_UP_LAG+1)*T_WAYMARK) - T_SAMPLE)	- signalDelayT0Block(1 : NUM_SAMPLES_WAYMARK);

	if (waymarkIndex >= WAYMARK_UP_LAG)		% Wait that appropriate delay in Waymarks and start recording the output
														% Shift the spectrum of the test signal up by the carrier frequency, taking into account the appropriate delays
		interpSignalUShift = [interpSignalUShift real(interpSignalUp.*exp( 1i*2*pi*FREQ_SHIFT*timeOpSignal ))];
	end

	if (mod(waymarkIndex, 10) == 0)
		disp(['Saving output waveform, Waymark ' num2str(waymarkIndex) ' of ' num2str(NUM_WAYMARKS) ', ' experimentTag '.wav']);
		fieldPropDelayT0InSamples = round(fieldPropDelayT0/T_SAMPLE);		% Calculate the number of zero samples to add at the beginning of the file to make it the correct length with the common delay added.
		audiowrite([experimentTagDir experimentTag '.wav'], [zeros(1,fieldPropDelayT0InSamples) interpSignalUShift], F_SAMPLE, 'BitsPerSample', N_BITS);
	end	

end


% ==================
% Record the results
% ==================
outputSignalTime	= size(interpSignalUShift,2)*T_SAMPLE;
fieldPropDelayT0InSamples = round(fieldPropDelayT0/T_SAMPLE);		% Calculate the number of zero samples to add at the beginning of the file to make it the correct length with the common delay added.
% wavwrite can't append so write the whole file here
audiowrite([experimentTagDir experimentTag '_' num2str(outputSignalTime) 's.wav'], [zeros(1,fieldPropDelayT0InSamples) interpSignalUShift], F_SAMPLE, 'BitsPerSample', N_BITS);

% Stop stopwatch timer
calcRunTime = toc(timerWaymark);
disp(['Run time, using simple tick-toc counter: ' num2str(calcRunTime/60) ' mins'])

% Save the whole variable set
save ([experimentTagDir, experimentTag, '.mat']);
waymarkPlot;



end

