function[sweep,inverseSweep] = logSweepGen(start,stop,duration,Fs);
% function[sweep,inverseSweep] = logSweepGen(start,stop,duration,Fs);
% function to produce a logarithmic frequency sweep rising from the start
% frequency to stop frequency (both in Hz) for the duration T (in seconds).
% The sample rate (Fs) must also be specified
% See page 6 of "Simultaneous measurement of impulse response and distortion
% with a swept sine technique" - Angelo Farina for the equation used.
% The inverse sweep is generated by time-reversing the sweep and applying a
% -6dB per octave envelope
% start: start frequency of sweep
% stop: stop frequency of sweep
% duration: duration of sweep in seconds
% Fs: sample rate
% convolution of the sweep with its inverse will produce an impulse
% Note that a 100 ms ramp in and out is applied to the sweep to prevent start
% and stop transients, therefore the sweep should be as long as possible (e.g.
% many seconds)
% This implementation by Jez Wells, University of York, August 2006

% convert from Hz to radians per second
start = start *2*pi;
stop = stop*2*pi;

% produce swept sinusoid
temp1 = log(stop/start);
temp2 = exp(((linspace(0,duration,Fs*duration)./duration)*temp1))-1;
sweep = sin(((start*duration)/temp1)*temp2);
arg =((start*duration)/temp1)*temp2;

% apply ramp in and ramp out over 100 ms
sweep(1:round(0.1*Fs)) = sweep(1:round(0.1*Fs)).*linspace(0,1,round(0.1*Fs));
sweep(end-round(0.1*Fs)+1:end) = sweep(end-round(0.1*Fs)+1:end).*linspace(1,0,round(0.1*Fs));

% to generate inverse time-reverse the excitation and apply a -6dB octave
% envelope
inverseSweep = rot90(sweep,-2);
envelope = logspace(log10(1),log10(start/stop),Fs*duration);
inverseSweep = envelope.*inverseSweep;