function [stim, unfilt] = genstim(varargin)
%GENSTIM Generate stimulus
% 	[S, U] = GENSTIM(stimlen, envlen, separation)
% 	Generate S, a series of 4 white noise stimuli, each of 
% 	length <stimlen> (msec) with a sin-squared envelope of 
% 	length <envlen> (msec).  U will be a vector with 2
%	identical white noise stimuli, to be used as reference.  Each
%	stimulus will be separated by a <separation> (msec) length
%	silence.  Stimuli are monaural, 16-bit, with a sampling
%	rate of 44.1kHz.  Stimuli will be filtered by the FIR
%	filter saved in the .mat file noisefiltB.mat, which should be
%	located in the current directory.
% 	Currently, this is a bandpass roughly from 200Hz-17000Hz;
% 	All arguments are optional, and they will default to 
% 	stimlen=750 (msec), envlen=20 (msec) and separation=1250 (msec)
% 	S is a column vector.

% Look for input arguments or use defaults
if (nargin < 3)
	filler_msec = 1250;
else
	filler_msec = varargin{3};
end

if (nargin < 2)
	envlen_msec = 20;
else
	envlen_msec = varargin{2};
end

if (nargin < 1)
	stimlen_msec = 750;
else
	stimlen_msec = varargin{1};
end


Fs = 44100;	% Sampling rate
plot_stim = 0;	% debugging flag that plots stimulus

stimlen_samp = (stimlen_msec/1000) * Fs;

% Generate Gaussian white noise
stim = randn(stimlen_samp, 1);

load noisefiltB;

% Bandpass filter noise
stim = filter(noisefiltB, 1, stim);

%Generate envelope
envlen_samp = (envlen_msec/1000) * Fs;
env = (sin((1:envlen_samp) / envlen_samp).^2)';

% Apply envelope
index = stimlen_samp-envlen_samp+1;

stim(1:envlen_samp) = stim(1:envlen_samp) .* env;
stim(index:end) = stim(index:end) .* flipud(env);

% Generate "filler" (silence) vector
filler = zeros(Fs*(filler_msec/1000),1);

% Generate total stimulus by alternating noise and silence
totalstim = [stim; filler; stim; filler; stim; filler; stim; filler];

unfilt = [stim; filler; stim; filler; zeros(44100,1)];

% Normalize
totalstim = totalstim / max(totalstim);
unfilt = unfilt / max(unfilt);

stim = totalstim;

if (plot_stim == 1)
	figure(2)
	subplot(221)     
	plot(env)
	title('envelope');

	subplot(222)
	plot(totalstim);
	title('total stimulus')

	subplot(223)
	plot(totalstim(1:Fs))
	title('one second of stimulus')

	subplot(224)
	plot(real(fft(totalstim)));
	title('magnitude DFT of total stimulus');

	figure(1)
	plot(unfilt)
	title('unfilt');

end;