rmeddis@0: function [audio, msg, time]=stimulusCreate(globalStimParams,...
rmeddis@0:     stimComponents, doPlot)
rmeddis@0: % updated June 2007
rmeddis@0: % the only authorised version of stimulus create is the version to be found
rmeddis@0: % in MAP1_6.  Other versions are copies!!
rmeddis@0: %
rmeddis@0: % for a simple tone you need
rmeddis@0: %
rmeddis@0: % % Mandatory structure fields
rmeddis@0: %  globalStimParams.FS=100000;
rmeddis@0: %  globalStimParams.overallDuration=.1;  % s
rmeddis@0: % doPlot=1;
rmeddis@0: %
rmeddis@0: % stim.type='tone';
rmeddis@0: % stim.phases='sin';
rmeddis@0: % stim.toneDuration=.05;;
rmeddis@0: % stim.frequencies=500;
rmeddis@0: % stim.amplitudesdB=50;
rmeddis@0: % stim.beginSilence=.01;
rmeddis@0: % stim.endSilence=-1;
rmeddis@0: % stim.rampOnDur=.002;
rmeddis@0: % stim.rampOffDur=-1;
rmeddis@0: % [audio, msg]=stimulusCreate(globalStimParams, stim, doPlot);
rmeddis@0: %
rmeddis@0: % % or for multi-component stimuli
rmeddis@0: % % Mandatory structure fields
rmeddis@0: %  globalStimParams.FS=100000;
rmeddis@0: %  globalStimParams.overallDuration=.1;  % s
rmeddis@0: % ear=1;
rmeddis@0: % componentNo=1;
rmeddis@0: %
rmeddis@0: % stimComponents(ear, componentNo).type='tone';
rmeddis@0: % stimComponents(ear, componentNo).phases='sin';
rmeddis@0: % stimComponents(ear, componentNo).toneDuration=.05;;
rmeddis@0: % stimComponents(ear, componentNo).frequencies=500;
rmeddis@0: % stimComponents(ear, componentNo).amplitudesdB=50;
rmeddis@0: % stimComponents(ear, componentNo).beginSilence=.01;
rmeddis@0: % stimComponents(ear, componentNo).endSilence=-1;
rmeddis@0: % stimComponents(ear, componentNo).rampOnDur=.002;
rmeddis@0: % stimComponents(ear, componentNo).rampOffDur=-1;
rmeddis@0: %
rmeddis@0: % % All components are forced to have the same overall duration and sample rate
rmeddis@0: %
rmeddis@0: %
rmeddis@0: % [audio, msg]=stimulusCreate(globalStimParams, stimComponents);
rmeddis@0: %
rmeddis@0: %  Optional fields
rmeddis@0: %  .ears overides ear setting by component
rmeddis@0: %  globalStimParams.ears='monoticL'; % 'monoticL', 'monoticR', 'diotic', 'dichotic'
rmeddis@0: %
rmeddis@0: %  globalStimParams.filter = [leftfl leftfu rightfl right fu]
rmeddis@0: %    filter is applied separately to left and right combined sounds
rmeddis@0: %
rmeddis@0: %  correction factor is applied to all signals to compensate for differences in output devices.
rmeddis@0: % audioOutCorrection is a scalar
rmeddis@0: %  globalStimParams.audioOutCorrection=2;
rmeddis@0: %
rmeddis@0: %  globalStimParams.FFT= 1; % {0, 1}
rmeddis@0: %  globalStimParams.doPlot=1; % {0, 1}
rmeddis@0: %  globalStimParams.doPlay=1; % {0, 1}
rmeddis@0: %
rmeddis@0: %  stimComponents(ear, componentNo).filter=[100 10000 2] % [lower, upper, order] applies to an
rmeddis@0: % individual component
rmeddis@0: %
rmeddis@0: % Output arguments:
rmeddis@0: %  audio is a stereo signal, a 2-column vector
rmeddis@0: %
rmeddis@0: %
rmeddis@0: % stim.type='noise';    % {'IRN', 'irn', 'noise', 'pinkNoise'}
rmeddis@0: % stim.toneDuration=.05;;
rmeddis@0: % stim.amplitudesdB=50;
rmeddis@0: % stim.beginSilence=.01;
rmeddis@0: % stim.endSilence=-1;
rmeddis@0: % stim.rampOnDur=.002;
rmeddis@0: % stim.rampOffDur=-1;
rmeddis@0: % [audio, msg]=stimulusCreate(globalStimParams, stim);
rmeddis@0: %
rmeddis@0: % % for IRN only
rmeddis@0: % stim.type='IRN';
rmeddis@0: % stim.toneDuration=.05;;
rmeddis@0: % stim.amplitudesdB=50;
rmeddis@0: % stim.beginSilence=.01;
rmeddis@0: % stim.endSilence=-1;
rmeddis@0: % stim.rampOnDur=.002;
rmeddis@0: % stim.rampOffDur=-1;
rmeddis@0: % stim.niterations = 8;   %0 for white noise
rmeddis@0: % stim.delay = 1/150;
rmeddis@0: % stim.irnGain = 1;
rmeddis@0: % [audio, msg]=stimulusCreate(globalStimParams, stim);stimComponents.clickWidth;
rmeddis@0: %
rmeddis@0: % stim.type='clicks';    % uses frequencies for duty cycle
rmeddis@0: % stim.toneDuration=.05;;
rmeddis@0: % stim.frequencies=500;
rmeddis@0: % stim.amplitudesdB=50;
rmeddis@0: % stim.beginSilence=.01;
rmeddis@0: % stim.endSilence=-1;
rmeddis@0: % stim.rampOnDur=.002;
rmeddis@0: % stim.rampOffDur=-1;
rmeddis@0: % [audio, msg]=stimulusCreate(globalStimParams, stim);
rmeddis@0: % stimComponents.clickWidth=.0001;    % default= dt
rmeddis@0: % stimComponents.clickHeight= 20e-6;  % default= 28e-6 * 10^(stimComponents.amplitudesdB/20);
rmeddis@0: % stim.clickTimes=[.4 .6];            % times in cylce when clicks occur, default = 1
rmeddis@0: %
rmeddis@0: %
rmeddis@0: 
rmeddis@0: msg=''; %error messages; no message is a good message
rmeddis@0: 
rmeddis@0: % plotting can be set as a separate argument or as a globalstimParams
rmeddis@0: % variable. this is for backwards compatibility only
rmeddis@0: if nargin>2,
rmeddis@0:     globalStimParams.doPlot=doPlot;
rmeddis@0: end
rmeddis@0: 
rmeddis@0: % stimComponents(1,1).endSilence=-1;  % end silence is always computed
rmeddis@0: 
rmeddis@0: % perform checks and set defaults
rmeddis@0: [globalStimParams, stimComponents]=checkDescriptors(globalStimParams, stimComponents);
rmeddis@0: 
rmeddis@0: 
rmeddis@0: % create empty stereo audio of appropriate length
rmeddis@0: audio=zeros(globalStimParams.nSignalPoints, 2);
rmeddis@0: % signal=zeros(globalStimParams.nSignalPoints, 1);
rmeddis@0: dt=globalStimParams.dt;
rmeddis@0: 
rmeddis@0: [Nears nComponentSounds]=size(stimComponents);
rmeddis@0: for ear=1:Nears % left/ right
rmeddis@0:     % combinedSignal is the sum of all sounds in one ear
rmeddis@0:     % it is a column vector
rmeddis@0:     combinedSignal=zeros(globalStimParams.nSignalPoints,1);
rmeddis@0: 
rmeddis@0:     % find valid components
rmeddis@0:     % if a component is empty, it is not a validComponent and is ignored
rmeddis@0:     validStimComponents=[];
rmeddis@0:     for i=1:nComponentSounds
rmeddis@0:         if ~isempty(stimComponents(ear,i).type)
rmeddis@0:             validStimComponents=[validStimComponents i];
rmeddis@0:         end
rmeddis@0:     end
rmeddis@0: 
rmeddis@0:     for componentNo=validStimComponents
rmeddis@0:         % compute individual components before adding
rmeddis@0:         stim=stimComponents(ear,componentNo);
rmeddis@0:         switch stim.type
rmeddis@0:             case 'tone'
rmeddis@0:                 stimulus=maketone(globalStimParams, stim);
rmeddis@0: 
rmeddis@0:             case 'fmTone'
rmeddis@0:                 stimulus=makeFMtone(globalStimParams, stim);
rmeddis@0: 
rmeddis@0:             case 'OHIO'
rmeddis@0:                 stim.beginSilence=0;
rmeddis@0:                 stimulus=makeOHIOtones(globalStimParams, stim);
rmeddis@0: 
rmeddis@0:             case 'transposedStimulus'
rmeddis@0:                 stim.beginSilence=0;    % necessary because of recursion
rmeddis@0:                 stimulus=makeTransposedStimulus(globalStimParams, stim);
rmeddis@0: 
rmeddis@0:             case { 'noise', 'pinkNoise'}
rmeddis@0:                 stimulus=makeNoise(globalStimParams, stim);
rmeddis@0: 
rmeddis@0:             case { 'whiteNoise'}
rmeddis@0:                 stimulus=makeWhiteNoise(globalStimParams, stim);
rmeddis@0: 
rmeddis@0:             case {'IRN', 'irn'}
rmeddis@0:                 stimulus=makeIRN(globalStimParams, stim);
rmeddis@0: 
rmeddis@0:             case {'RPN'}
rmeddis@0:                 stimulus=makeRPN(globalStimParams, stim);
rmeddis@0: 
rmeddis@0:             case 'clicks'
rmeddis@0:                 stimulus=makeClicks(globalStimParams, stim);
rmeddis@0: 
rmeddis@0:             case 'PressnitzerClicks'
rmeddis@0:                 % kxx clicks
rmeddis@0:                 % k is 1/clickRepeatFrequency
rmeddis@0:                 stimulus=makePressnitzerClicks(globalStimParams, stimComponents);
rmeddis@0: 
rmeddis@0:             case 'PressnitzerABxClicks'
rmeddis@0:                 % kxx clicks
rmeddis@0:                 % k is 1/clickRepeatFrequency
rmeddis@0:                 stimulus=makePressnitzerABxClicks(globalStimParams, stimComponents);
rmeddis@0: 
rmeddis@0:             case 'ABxClicks'
rmeddis@0:                 % A=rand*k,  B=k-A, x=rand*k.
rmeddis@0:                 stimulus=makeABxClicks(globalStimParams, stimComponents);
rmeddis@0: 
rmeddis@0:             case 'YostClicks'
rmeddis@0:                 % kxx clicks
rmeddis@0:                 % k is 1/clickRepeatFrequency
rmeddis@0:                 stimulus=makeYostClicks(globalStimParams, stimComponents);
rmeddis@0: 
rmeddis@0:             case 'kxxClicks'
rmeddis@0:                 % kxx clicks
rmeddis@0:                 % k is 1/clickRepeatFrequency
rmeddis@0:                 stimulus=makeKxxClicks(globalStimParams, stimComponents);
rmeddis@0: 
rmeddis@0: 
rmeddis@0: %             case 'babble'
rmeddis@0: %                 % NB random start in a long file
rmeddis@0: %                 [stimulus sampleRate]= wavread('babble');
rmeddis@0: %                 nPoints=round(sampleRate*...
rmeddis@0: %                     stimComponents(ear,componentNo).toneDuration);
rmeddis@0: %                 start=round(rand(1,1)*(length(stimulus)-nPoints));
rmeddis@0: %                 stimulus=stimulus(start:start+nPoints-1);
rmeddis@0: %                 rms= 20*log10((mean(stimulus.^2).^0.5)/20e-6);
rmeddis@0: %                 dBSPLrms=stimComponents(ear,componentNo).amplitudesdB;
rmeddis@0: %                 gain=10.^((dBSPLrms-rms)/20);
rmeddis@0: %                 stimulus=stimulus'*gain;
rmeddis@0: 
rmeddis@0:             case 'speech'
rmeddis@0:                 [stimulus sampleRate]= wavread('speech');
rmeddis@0:                 stimulus=stimulus';
rmeddis@0:                 nPoints=sampleRate*stimComponents(ear,componentNo).toneDuration;
rmeddis@0:                 if nPoints > length(stimulus)
rmeddis@0:                     initialSilence=zeros(1,nPoints-length(stimulus));
rmeddis@0:                 else
rmeddis@0:                     initialSilence=[];
rmeddis@0:                     start=round(rand(1,1)*(length(stimulus)-nPoints));
rmeddis@0:                     stimulus=stimulus(start:start+nPoints-1);
rmeddis@0:                 end
rmeddis@0:                 rms= 20*log10((mean(stimulus.^2).^0.5)/20e-6);
rmeddis@0:                 dBSPLrms=stimComponents(ear,componentNo).amplitudesdB;
rmeddis@0:                 gain=10.^((dBSPLrms-rms)/20);
rmeddis@0:                 stimulus=stimulus*gain;
rmeddis@0:                 stimulus=[stimulus initialSilence ];
rmeddis@0: 
rmeddis@0: 
rmeddis@0:             case 'file'
rmeddis@0:                 % component already read from file and stored in stimulus. Insert it here
rmeddis@0:                 % additional code for establishing signal rms level
rmeddis@0:                 % NB signal is always mono at this stage
rmeddis@0: 
rmeddis@0:                 stimulus=stim.stimulus;
rmeddis@0:                 dBSPL=stim.amplitudesdB;
rmeddis@0: 
rmeddis@0:                 nPoints=round(stim.toneDuration/dt);
rmeddis@0:                 [r c]=size(stimulus);
rmeddis@0:                 if r>c, stimulus=stimulus'; end    % secure horizontal vector
rmeddis@0:                 stimulus=stimulus(1,1:nPoints); % only mono taken from file
rmeddis@0: 
rmeddis@0:                 try
rmeddis@0:                     % dB rms
rmeddis@0:                     rms= 20*log10((mean(stimulus.^2).^0.5)/20e-6);
rmeddis@0:                     % special request to set fixed rms for stimulus
rmeddis@0:                     dBSPLrms=stimComponents(ear,componentNo).amplitudesdBrms;
rmeddis@0:                     if ~(dBSPLrms==-1)
rmeddis@0:                         gain=10.^((dBSPLrms-rms)/20);
rmeddis@0:                         stimulus=stimulus*gain;
rmeddis@0:                     end
rmeddis@0:                 catch
rmeddis@0:                     % If no request for rms level is made
rmeddis@0:                     % set dB as peak amp
rmeddis@0:                     [stimulus gain]= normalize(stimulus);
rmeddis@0:                     dBSPL=stimComponents(ear,componentNo).amplitudesdB;
rmeddis@0:                     if ~(dBSPL==-1)
rmeddis@0:                         amplitude=28e-6*10.^(dBSPL/20);
rmeddis@0:                         stimulus=stimulus*amplitude;
rmeddis@0:                     end
rmeddis@0:                 end
rmeddis@0: 
rmeddis@0:             case 'none'
rmeddis@0:                 % no stimulus
rmeddis@0:                 stimulus=zeros(1,round(stim.toneDuration/dt));
rmeddis@0: 
rmeddis@0:             case 'digitStrings'
rmeddis@0:                 % select a digit string at random anduse as target
rmeddis@0:                 files=dir(['..' filesep '..' filesep 'multithresholdResources\digitStrings']);
rmeddis@0:                 files=files(3:end);
rmeddis@0:                 nFiles=length(files);
rmeddis@0:                 fileNo=ceil(nFiles*rand);
rmeddis@0:                 fileData=files(fileNo);
rmeddis@0:                 fileName=['..\..\multithresholdResources\digitStrings\' fileData.name];
rmeddis@0:                 [stimulus sampleRate]=wavread(fileName);
rmeddis@0:                 stimulus=stimulus(:,1)';  % make into a row vector
rmeddis@0:                 % estimate the extend of endsilence padding
rmeddis@0:                 nPoints=sampleRate*...
rmeddis@0:                     stimComponents(ear,componentNo).toneDuration;
rmeddis@0:                 if nPoints > length(stimulus)
rmeddis@0:                     endSilence=zeros(1,nPoints-length(stimulus));
rmeddis@0:                 else
rmeddis@0:                     % or truncate the file
rmeddis@0:                     endSilence=[];
rmeddis@0:                     stimulus=stimulus(1:nPoints);
rmeddis@0:                 end
rmeddis@0:                 % compute rms before adding silence
rmeddis@0:                 rms= 20*log10((mean(stimulus.^2).^0.5)/20e-6);
rmeddis@0:                 dBSPLrms=stimComponents(ear,componentNo).amplitudesdB;
rmeddis@0:                 gain=10.^((dBSPLrms-rms)/20);
rmeddis@0:                 stimulus=stimulus*gain;
rmeddis@0:                 stimulus=[stimulus endSilence];
rmeddis@0:                 global stimulusParameters
rmeddis@0:                 stimulusParameters.digitString=fileName(end-7:end-5);
rmeddis@0: 
rmeddis@0:             otherwise
rmeddis@0:                 switch stim.type(end-5:end)
rmeddis@0:                     % any file name with 'Babble' at the end is a
rmeddis@0:                     % multiThreshold file
rmeddis@0:                     case 'Babble'
rmeddis@0:                         % one of the many babbles is selected.
rmeddis@0:                         % NB random start in a long file
rmeddis@0:                         %  stim.type should contain the name of the babble file
rmeddis@0:                         fileName=['..' filesep '..' filesep ...
rmeddis@0:                             'multithresholdResources' filesep ...
rmeddis@0:                             'backgrounds and maskers'...
rmeddis@0:                              filesep stim.type];
rmeddis@0:                         
rmeddis@0:                         [stimulus sampleRate]= wavread(fileName);
rmeddis@0:                         if ~isequal(sampleRate, globalStimParams.FS)
rmeddis@0:                             % NB the file will read but will disagree with
rmeddis@0:                             % tone stimuli or other files read
rmeddis@0:                             msg= ['error: file sample rate disagrees ' ...
rmeddis@0:                                 'with sample rate requested in paradigm'...
rmeddis@0:                                 ' file (' ...
rmeddis@0:                             num2str(globalStimParams.FS) ').'];
rmeddis@0:                             error(msg);                            
rmeddis@0:                         end
rmeddis@0:                         nPoints=round(sampleRate*...
rmeddis@0:                             stimComponents(ear,componentNo).toneDuration);
rmeddis@0:                         start=round(rand(1,1)*(length(stimulus)-nPoints));
rmeddis@0:                         stimulus=stimulus(start:start+nPoints-1);
rmeddis@0:                         rms= 20*log10((mean(stimulus.^2).^0.5)/20e-6);
rmeddis@0:                         dBSPLrms=stimComponents(ear,componentNo).amplitudesdB;
rmeddis@0:                         gain=10.^((dBSPLrms-rms)/20);
rmeddis@0:                         stimulus=stimulus'*gain;
rmeddis@0:                         
rmeddis@0:                     otherwise
rmeddis@0:                         % stim.type may be the name of a file to be read
rmeddis@0:                         % play from beginning for stimulus duration
rmeddis@0:                         try
rmeddis@0:                             [stimulus sampleRate]= wavread([stim.type '.wav']);
rmeddis@0:                         catch
rmeddis@0:                             error(['stimulusCreate: unrecognised stimulus type -> '...
rmeddis@0:                                 stim.type])
rmeddis@0:                         end
rmeddis@0:                         if ~isequal(sampleRate, globalStimParams.FS)
rmeddis@0:                             % NB the file will read but will disagree with
rmeddis@0:                             % tone stimuli or other files read
rmeddis@0:                             msg= ['error: file sample rate disagrees ' ...
rmeddis@0:                                 'with sample rate requested in paradigm'...
rmeddis@0:                                 ' file (' ...
rmeddis@0:                             num2str(globalStimParams.FS) ').'];
rmeddis@0:                             error(msg);
rmeddis@0:                         end
rmeddis@0:                         stimulus=stimulus';  % make into a row vector
rmeddis@0:                         % estimate the extend of endsilence padding
rmeddis@0:                         nPoints=sampleRate*...
rmeddis@0:                             stimComponents(ear,componentNo).toneDuration;
rmeddis@0:                         if nPoints > length(stimulus)
rmeddis@0:                             endSilence=zeros(1,nPoints-length(stimulus));
rmeddis@0:                         else
rmeddis@0:                             % or truncate the file
rmeddis@0:                             endSilence=[];
rmeddis@0:                             stimulus=stimulus(1:nPoints);
rmeddis@0:                         end
rmeddis@0:                         % compute rms before adding silence
rmeddis@0:                         rms= 20*log10((mean(stimulus.^2).^0.5)/20e-6);
rmeddis@0:                         dBSPLrms=stimComponents(ear,componentNo).amplitudesdB;
rmeddis@0:                         gain=10.^((dBSPLrms-rms)/20);
rmeddis@0:                         stimulus=stimulus*gain;
rmeddis@0:                         stimulus=[stimulus endSilence];
rmeddis@0:                 end
rmeddis@0:         end
rmeddis@0: 
rmeddis@0:         % NB stimulus is a row vector now!
rmeddis@0:         % audio and combinedSignal were column vectors
rmeddis@0:         % signal will be a row vector
rmeddis@0: 
rmeddis@0:         % filter stimulus
rmeddis@0:         try
rmeddis@0:             % if filter field is present, [lower upper order]
rmeddis@0:             if stim.filter(1)>0	% 0 means don't filter
rmeddis@0:                 stimulus=Butterworth (stimulus, dt, stim.filter(1), ...
rmeddis@0:                     stim.filter(2), stim.filter(3));
rmeddis@0: 
rmeddis@0:             end
rmeddis@0:         catch
rmeddis@0:         end
rmeddis@0: 
rmeddis@0: 
rmeddis@0:         % apply amplitude modulation
rmeddis@0:         if isfield(stim,'AMfrequency') & isfield(stim,'AMdepth')
rmeddis@0:             if stim.AMfrequency>0 & stim.AMdepth>0
rmeddis@0:                 time=dt:dt:dt*length(stimulus);
rmeddis@0:                 modulator=sin(2*pi*stim.AMfrequency*time);
rmeddis@0:                 modulator=modulator*stim.AMdepth/100 + 1; % depth is percent
rmeddis@0:                 stimulus=stimulus.*modulator/2;
rmeddis@0:             end
rmeddis@0:         end
rmeddis@0: 
rmeddis@0:         % Basic stimulus is now created.
rmeddis@0:         % Add trappings, ramps, silences to main stimulus
rmeddis@0:         rampOnTime=0;		 %ramp starts at the beginning of the stimulus
rmeddis@0:         rampOffTime=stim.toneDuration-stim.rampOffDur;
rmeddis@0:         if stim.rampOnDur>0.0001
rmeddis@0:             stimulus=applyRampOn(stimulus, stim.rampOnDur, rampOnTime, 1/dt);
rmeddis@0:             stimulus=applyRampOff(stimulus, stim.rampOffDur, rampOffTime, 1/dt);
rmeddis@0:         end
rmeddis@0:         if stim.rampOnDur<-0.0001   % apply Gaussian ramp
rmeddis@0:             stimulus=applyGaussianRamps(stimulus, -stim.rampOnDur, 1/dt);
rmeddis@0:         end
rmeddis@0: 
rmeddis@38:         % begin silence
rmeddis@0:         % start with a signal of the right length consisting of zeros
rmeddis@0:         signal=zeros(1, globalStimParams.nSignalPoints);
rmeddis@0:         % identify start of stimulus
rmeddis@0:         insertStimulusAt=round(stim.beginSilence/dt)+1;
rmeddis@0:         % add stimulus
rmeddis@0:         endOfStimulus=insertStimulusAt+length(stimulus)-1;
rmeddis@0:         if endOfStimulus<=globalStimParams.nSignalPoints
rmeddis@0:             signal(1, insertStimulusAt: endOfStimulus)=stimulus;
rmeddis@0:         else
rmeddis@0:             error('stimulusCreate: tone too long to fit into the overall duration')
rmeddis@0:         end
rmeddis@0: 
rmeddis@0:         % time signature
rmeddis@0:         time=dt:dt:dt*length(signal);
rmeddis@38:         % figure(22), plot(signal), title([num2str(ear) ' - ' num2str(componentNo)]),pause (1)
rmeddis@0: 
rmeddis@0:         try
rmeddis@0:             % create a column vector and trap if no signal has been created
rmeddis@0:             signal=signal';
rmeddis@0:             % also traps if signals are not the same length
rmeddis@0:             combinedSignal=combinedSignal+signal;
rmeddis@0:             %             figure(21), plot(combinedSignal), title([num2str(ear) ' - ' num2str(componentNo)]),pause (1)
rmeddis@0:         catch
rmeddis@0:             % dump everything because there is a problem
rmeddis@0:             globalStimParams
rmeddis@0:             [ear  componentNo]
rmeddis@0:             stim
rmeddis@0:             size(combinedSignal)
rmeddis@0:             size(signal)
rmeddis@0:             [   size(initialSilence)            size(signal)            size(endSilence)]
rmeddis@0:             [ ear  componentNo...
rmeddis@0:                 round(globalStimParams.overallDuration*globalStimParams.FS)...
rmeddis@0:                 round(stim.beginSilence*globalStimParams.FS)...
rmeddis@0:                 round(stim.toneDuration*globalStimParams.FS) ...
rmeddis@0:                 round(stim.endSilence*globalStimParams.FS)...
rmeddis@0:                 (                       round(globalStimParams.overallDuration*globalStimParams.FS)...
rmeddis@0:                 -round(stim.beginSilence*globalStimParams.FS)...
rmeddis@0:                 -round(stim.toneDuration*globalStimParams.FS) ...
rmeddis@0:                 -round(stim.endSilence*globalStimParams.FS))...
rmeddis@0:                 ]
rmeddis@0:             error(' trouble in stimulusCreate: signals are the wrong length ')
rmeddis@0:         end
rmeddis@0: 
rmeddis@0:     end % component no
rmeddis@0: 
rmeddis@0:     audio(:,ear)= combinedSignal;
rmeddis@0: 
rmeddis@0:     % FFT
rmeddis@0:     try
rmeddis@0:         if globalStimParams.FFT
rmeddis@0:             showFFT (audio, dt)
rmeddis@0:         end
rmeddis@0:     catch
rmeddis@0:     end
rmeddis@0: 
rmeddis@0: 
rmeddis@0: end % ear
rmeddis@0: 
rmeddis@0: switch globalStimParams.ears
rmeddis@0:     % normally the signals are created in appropriate ears but .ears can
rmeddis@0:     % overide this to produce a mono signal.
rmeddis@0:     case 'monoticL';
rmeddis@0:         % combine left and right ears to make a mono signal in the left ear
rmeddis@0:         audio(:,1)=audio(:,1)+audio(:,2);
rmeddis@0:         audio(:,2)=zeros(globalStimParams.nSignalPoints, 1);
rmeddis@0: 
rmeddis@0:     case 'monoticR';
rmeddis@0:         % combine left and right ears to make a mono signal in the right ear
rmeddis@0:         audio(:,2)=audio(:,1)+audio(:,2);
rmeddis@0:         audio(:,1)=zeros(globalStimParams.nSignalPoints, 1);
rmeddis@0: 
rmeddis@0:     case 'diotic';
rmeddis@0:         % combine left and right ears to make a mono signal in both ears
rmeddis@0:         bothEarsCombined=audio(:,1)+audio(:,2);
rmeddis@0:         audio(:,2)=bothEarsCombined;
rmeddis@0:         audio(:,1)=bothEarsCombined;
rmeddis@0: 
rmeddis@0:     otherwise
rmeddis@0:         % Any other denomination produces no effect here
rmeddis@0: end
rmeddis@0: 
rmeddis@0: % Plotting as required
rmeddis@0: if globalStimParams.doPlot
rmeddis@0:     figure(9), clf
rmeddis@0:     plot(time,audio(:,1)'), hold on,
rmeddis@0:     if Nears>1
rmeddis@0:         offSet=(max(audio(:,1))+max(audio(:,2)))/10;
rmeddis@0:         offSet=2*offSet+max(audio(:,1))+max(audio(:,2));
rmeddis@0:         plot(time,audio(:,2)'+offSet,'r'), hold off
rmeddis@0:     end
rmeddis@0:     ylabel('left                             right')
rmeddis@0: end
rmeddis@0: 
rmeddis@0: % Calibration
rmeddis@0: % all signals are divided by this correction factor
rmeddis@0: % peakAmp=globalStimParams.audioOutCorrection; % microPascals = 100 dB SPL
rmeddis@0: 
rmeddis@0: audio=audio/globalStimParams.audioOutCorrection;
rmeddis@0: 
rmeddis@0: if globalStimParams.doPlay
rmeddis@35:         if ispc
rmeddis@35:             wavplay(audio,globalStimParams.FS)
rmeddis@35:         else
rmeddis@35:             sound(audio,globalStimParams.FS)
rmeddis@35:         end
rmeddis@35: 
rmeddis@0:     wavplay(audio,globalStimParams.FS)
rmeddis@0: end
rmeddis@0: % all Done
rmeddis@0: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
rmeddis@0: 
rmeddis@0: 
rmeddis@0: 
rmeddis@0: %---------------------------------------------------- maketone
rmeddis@0: function tone=maketone(globalStimParams, stimComponents)
rmeddis@0: % maketone generates a stimComponents tone
rmeddis@0: % tone=maketone(dt, frequencies, toneDuration, dBSPL, phases)
rmeddis@0: % tone is returned in Pascals
rmeddis@0: % frequencies is a list of frequencies
rmeddis@0: % phase is list of phases or 'sin', 'cos', 'alt'
rmeddis@0: %
rmeddis@0: % defaults:
rmeddis@0: %  phase = sin
rmeddis@0: %  dBSPL=56 dB SPL
rmeddis@0: 
rmeddis@0: dt=globalStimParams.dt;
rmeddis@0: frequencies=stimComponents.frequencies;
rmeddis@0: toneDuration=stimComponents.toneDuration;
rmeddis@0: dBSPL=stimComponents.amplitudesdB;
rmeddis@0: phases=stimComponents.phases;
rmeddis@0: 
rmeddis@0: if ischar(phases)
rmeddis@0:     switch phases
rmeddis@0:         case 'sin'
rmeddis@0:             phases= zeros(1,length(frequencies));
rmeddis@0:         case 'cos'
rmeddis@0:             phases= pi/2*ones(1,length(frequencies));
rmeddis@0:         case 'alt'
rmeddis@0:             phases= repmat([0 pi/2], 1, floor(length(frequencies)/2));
rmeddis@0:             if length(phases)<length(frequencies)
rmeddis@0:                 phases=[phases 0];
rmeddis@0:             end
rmeddis@0:         case {'ran', 'rand'}
rmeddis@0:             phases= 2*pi*rand(1,length(frequencies));
rmeddis@0:     end
rmeddis@0: end
rmeddis@0: 
rmeddis@0: if length(phases)==1, phases=repmat(phases(1), 1, length(frequencies)); end
rmeddis@0: if length(phases)<length(frequencies)
rmeddis@0:     error('makeTone:phase specification must have the same length as frequencies')
rmeddis@0: end
rmeddis@0: 
rmeddis@0: if length(dBSPL)==1, dBSPL=repmat(dBSPL(1), 1, length(frequencies)); end
rmeddis@0: if length(dBSPL)<length(dBSPL)
rmeddis@0:     error('makeTone:dBSPL specification must have the same length as frequencies')
rmeddis@0: end
rmeddis@0: 
rmeddis@0: time=dt:dt:toneDuration;
rmeddis@0: amplitudes=28e-6* 10.^(dBSPL/20);
rmeddis@0: 
rmeddis@0: tone=zeros(size(time));
rmeddis@0: for i=1:length(frequencies)
rmeddis@0:     frequency=frequencies(i);
rmeddis@0:     phase=phases(i);
rmeddis@0:     tone=tone+amplitudes(i)*sin(2*pi*frequency*time+phase);
rmeddis@0: end
rmeddis@0: 
rmeddis@0: % figure(1), plot(time, signal)
rmeddis@0: 
rmeddis@0: 
rmeddis@0: %---------------------------------------------------- makeOHIOtones
rmeddis@0: function stimulus=makeOHIOtones(globalStimParams, stimComponents)
rmeddis@0: 
rmeddis@0: % Generates a stimulus consisting of one or more 10-ms tones
rmeddis@38: %  The length of the list of frequencies determines the numberof tones
rmeddis@38: % Tones are either presented at 10-ms intervals or simultaneously
rmeddis@38: % all tones are individually ramped
rmeddis@0: % Each tone has its own amplitude and its own ramp.
rmeddis@0: 
rmeddis@0: frequencies=stimComponents.frequencies;
rmeddis@0: amplitudesdB=stimComponents.amplitudesdB;
rmeddis@0: nFrequencies=length(frequencies);
rmeddis@0: 
rmeddis@38: if amplitudesdB==-100
rmeddis@38:     % catch trial
rmeddis@38:     amplitudesdB=repmat(-100,1,nFrequencies);
rmeddis@38: end
rmeddis@38: 
rmeddis@0: dt=globalStimParams.dt;
rmeddis@38: 
rmeddis@0: toneDuration=.010;
rmeddis@0: time=dt:dt:toneDuration;
rmeddis@0: 
rmeddis@38: % fixed ramp, silenceDuration, toneDuration
rmeddis@38: rampDuration=0.005;
rmeddis@0: rampTime=dt:dt:rampDuration;
rmeddis@38: ramp=[0.5*(1+cos(2*pi*rampTime/(2*rampDuration)+pi)) ...
rmeddis@38:     ones(1,length(time)-length(rampTime))];
rmeddis@0: ramp=ramp.*fliplr(ramp);
rmeddis@0: 
rmeddis@38: silenceDuration=0.010;
rmeddis@38: silenceDurationLength=round(silenceDuration/dt);
rmeddis@38: initialSilence=zeros(1,silenceDurationLength);
rmeddis@38: 
rmeddis@38: silenceToneDuration=toneDuration + silenceDuration;
rmeddis@38: silenceToneDurationLength=round(silenceToneDuration/dt);
rmeddis@38: 
rmeddis@38: % OHIO spect produces simultaneous tones
rmeddis@38: if strcmp(stimComponents.OHIOtype,'OHIOspect')
rmeddis@38:     totalDuration=silenceToneDuration;
rmeddis@38: else
rmeddis@38:     totalDuration=silenceToneDuration*nFrequencies;
rmeddis@38: end
rmeddis@38: 
rmeddis@38: totalDurationLength=round(totalDuration/dt);
rmeddis@38: stimulus=zeros(1,totalDurationLength);
rmeddis@38: toneBeginPTR=1;
rmeddis@0: 
rmeddis@0: for i=1:nFrequencies
rmeddis@0:     frequency=frequencies(i);
rmeddis@0:     dBSPL=amplitudesdB(i);
rmeddis@0:     amplitude=28e-6* 10.^(dBSPL/20);
rmeddis@0:     tone=amplitude*sin(2*pi*frequency*time);
rmeddis@0:     tone=tone.*ramp;
rmeddis@38:     % stimulus is normally zeros except for OHIOspect
rmeddis@38:     stimulus(toneBeginPTR:toneBeginPTR+silenceToneDurationLength-1)=...
rmeddis@38:         [initialSilence tone]+...
rmeddis@38:         stimulus(toneBeginPTR:toneBeginPTR+silenceToneDurationLength-1);    
rmeddis@38:     if ~strcmp(stimComponents.OHIOtype,'OHIOspect')
rmeddis@38:         toneBeginPTR=toneBeginPTR+silenceToneDurationLength;
rmeddis@38:     end
rmeddis@0: end
rmeddis@0: % figure(2), plot( stimulus')
rmeddis@0: 
rmeddis@0: 
rmeddis@0: %---------------------------------------------------- makeFMtone
rmeddis@0: function tone=makeFMtone(globalStimParams, stimComponents)
rmeddis@0: % maketone generates a stimComponents tone
rmeddis@0: % tone=maketone(dt, frequencies, toneDuration, dBSPL, phases)
rmeddis@0: % tone is returned in Pascals
rmeddis@0: % frequencies is a list of frequencies
rmeddis@0: % phase is list of phases or 'sin', 'cos', 'alt'
rmeddis@0: %
rmeddis@0: % defaults:
rmeddis@0: %  phase = sin
rmeddis@0: %  dBSPL=56 dB SPL
rmeddis@0: 
rmeddis@0: dt=globalStimParams.dt;
rmeddis@0: frequencies=stimComponents.frequencies;
rmeddis@0: toneDuration=stimComponents.toneDuration;
rmeddis@0: dBSPL=stimComponents.amplitudesdB;
rmeddis@0: phases=stimComponents.phases;
rmeddis@0: fmDepth=stimComponents.fmDepth;
rmeddis@0: fmFrequency=stimComponents.fmFrequency;
rmeddis@0: 
rmeddis@0: if ischar(phases)
rmeddis@0:     switch phases
rmeddis@0:         case 'sin'
rmeddis@0:             phases= zeros(1,length(frequencies));
rmeddis@0:         case 'cos'
rmeddis@0:             phases= pi/2*ones(1,length(frequencies));
rmeddis@0:         case 'alt'
rmeddis@0:             phases= repmat([0 pi/2], 1, floor(length(frequencies)/2));
rmeddis@0:             if length(phases)<length(frequencies)
rmeddis@0:                 phases=[phases 0];
rmeddis@0:             end
rmeddis@0:     end
rmeddis@0: end
rmeddis@0: 
rmeddis@0: if length(phases)==1, phases=repmat(phases(1), 1, length(frequencies)); end
rmeddis@0: if length(phases)<length(frequencies)
rmeddis@0:     error('makeTone:phase specification must have the same length as frequencies')
rmeddis@0: end
rmeddis@0: 
rmeddis@0: if length(dBSPL)==1, dBSPL=repmat(dBSPL(1), 1, length(frequencies)); end
rmeddis@0: if length(dBSPL)<length(dBSPL)
rmeddis@0:     error('makeTone:dBSPL specification must have the same length as frequencies')
rmeddis@0: end
rmeddis@0: 
rmeddis@0: time=dt:dt:toneDuration;
rmeddis@0: amplitudes=28e-6* 10.^(dBSPL/20);
rmeddis@0: 
rmeddis@0: tone=zeros(size(time));
rmeddis@0: for i=1:length(frequencies)
rmeddis@0:     frequency=frequencies(i);
rmeddis@0:     phase=phases(i);
rmeddis@0:     tone=tone+amplitudes(i)*sin(2*pi*frequency*time+phase + fmDepth*sin(2*pi*fmFrequency*time));
rmeddis@0: end
rmeddis@0: 
rmeddis@0: % figure(1), plot(time, signal)
rmeddis@0: 
rmeddis@0: %----------------------------------------------------makeTransposedStimulus
rmeddis@0: function [transposedStimulus]=makeTransposedStimulus(globalStimParams, stim)
rmeddis@0: % globalStimParams.FS=100000;
rmeddis@0: % globalStimParams.overallDuration=.1;  % s
rmeddis@0: % globalStimParams.doPlot=1;
rmeddis@0: %
rmeddis@0: % stim.type='transposedStimulus';
rmeddis@0: % stim.phases='sin';
rmeddis@0: % stim.toneDuration=.05;;
rmeddis@0: % stim.frequencies=500;
rmeddis@0: % stim.amplitudesdB=50;
rmeddis@0: % stim.beginSilence=.01;
rmeddis@0: % stim.endSilence=-1;
rmeddis@0: % stim.rampOnDur=.002;
rmeddis@0: % stim.rampOffDur=-1;
rmeddis@0: %
rmeddis@0: % stim.transCarrierFreq=4000;
rmeddis@0: % stim.transModFreq=100;
rmeddis@0: 
rmeddis@0: transposedStimulus=[];
rmeddis@0: % make envelope of transposed tone
rmeddis@0: for i=1:length(stim.transModFreq)
rmeddis@0:     stim.type='tone';
rmeddis@0:     stim.frequencies=stim.transModFreq(i);
rmeddis@0:     stim.endsilence=-1;     stim.beginSilence=0;
rmeddis@0:     [envelope, msg]=stimulusCreate(globalStimParams, stim);    % NB recursive
rmeddis@0:     envelope=envelope(:,1);   % mono
rmeddis@0:     % HW rectify
rmeddis@0:     envelope(find(envelope<0))=0;
rmeddis@0:     % LP filter
rmeddis@0:     maxEnvelope=max(envelope);
rmeddis@0:     envelope=UTIL_Butterworth (envelope, globalStimParams.dt, 10, .2*stim.transModFreq(i), 2);
rmeddis@0:     envelope=envelope*maxEnvelope/max(envelope);
rmeddis@0: 
rmeddis@0:     % make the carrier
rmeddis@0:     stim.frequencies=stim.transCarrierFreq(i);
rmeddis@0:     stim.endsilence=-1;    stim.beginSilence=0;
rmeddis@0:     [audio, msg]=stimulusCreate(globalStimParams, stim);
rmeddis@0:     carrier=audio(:,1);
rmeddis@0:     x= (carrier.*envelope)';
rmeddis@0:     % base amplitude on peak of unmodulated carrier
rmeddis@0:     x=x/max(carrier);
rmeddis@0:     transposedStimulus=[transposedStimulus; x];
rmeddis@0: end
rmeddis@0: transposedStimulus=sum(transposedStimulus,1);
rmeddis@0: % clf,plot(transposedStimulus)
rmeddis@0: 
rmeddis@0: %--------------------------------------------------------------------makeClicks
rmeddis@0: function clickTrain=makeClicks(globalStimParams, stimComponents)
rmeddis@0: % makeClicks(F0, clickTimes, duration, FS);
rmeddis@0: % F0 specifies the repetition rate of the click sequence
rmeddis@0: % If F0=-1, a single click is generated at the start of the duration of the signal
rmeddis@0: %
rmeddis@0: % clickTimes a are fractions of the period
rmeddis@0: %  and specify when the click appears in the period
rmeddis@0: % A click time of 1 is reset to zero.
rmeddis@0: % if the clicktime plus the click width is greater than the period, no click is generated
rmeddis@0: % clicks are treated as 20 microPascal high before amplification
rmeddis@0: %  unless otherwise specified in stimComponents.clickHeight
rmeddis@0: % click width is dt unless specified in stimComponents.clickWidth
rmeddis@0: %
rmeddis@0: % for regular pulse train set clicktimes=1 or zero;
rmeddis@0: % FS is the sampling interval;
rmeddis@0: % CarlyonClickTrain(100, [.4 1], 40, 441000);
rmeddis@0: 
rmeddis@0: FS=globalStimParams.FS; % sample rate
rmeddis@0: dt=1/FS;
rmeddis@0: 
rmeddis@0: try,clickWidth=stimComponents.clickWidth;catch, clickWidth=dt; end
rmeddis@0: try,clickHeight=stimComponents.clickHeight; catch, clickHeight=28e-6 * 10^(stimComponents.amplitudesdB/20); end
rmeddis@0: try, clickTimes=stimComponents.clickTimes; catch, clickTimes=1; end
rmeddis@0: 
rmeddis@0: % clickTimes are the times in a cycle that the click
rmeddis@0: % occurs
rmeddis@0: checkClickTimes=clickTimes-1;
rmeddis@0: if max(checkClickTimes)>1
rmeddis@0:     msg= 'click times must be <= than 1 (period)';
rmeddis@0:     return
rmeddis@0: end
rmeddis@0: 
rmeddis@0: if clickWidth>=1/stimComponents.clickRepeatFrequency
rmeddis@0:     msg= 'click width is too great for frequency';
rmeddis@0:     return
rmeddis@0: end
rmeddis@0: 
rmeddis@0: duration=stimComponents.toneDuration;
rmeddis@0: totalLength=round(stimComponents.toneDuration/dt);
rmeddis@0: F0=stimComponents.clickRepeatFrequency;
rmeddis@0: F0=round(F0/dt)*dt;
rmeddis@0: if F0==-1 % single click required
rmeddis@0:     F0=1/duration;
rmeddis@0:     repetitions=1;
rmeddis@0:     clickStartTimes=1; %clicktimes are fractions of a period
rmeddis@0: else
rmeddis@0:     repetitions=round(F0*duration)-1;
rmeddis@0:     duration=repetitions/F0;
rmeddis@0:     clickStartTimes=clickTimes;
rmeddis@0: end
rmeddis@0: % if a clickTime=1 (end of duty period) set it to the beginning
rmeddis@0: clickStartTimes(clickStartTimes==1)=0;
rmeddis@0: 
rmeddis@0: period=1/F0;
rmeddis@0: time=dt:dt:period;
rmeddis@0: nPoints=length(time);
rmeddis@0: signal=zeros(1,nPoints);
rmeddis@0: dBSPL=stimComponents.amplitudesdB;
rmeddis@0: 
rmeddis@0: % compute click train for a single cycle
rmeddis@0: clickWidthNpoints=round(clickWidth*FS);
rmeddis@0: for i=1:length(clickStartTimes)
rmeddis@0:     %     if clickStartTimes(i)<clickWidth
rmeddis@0:     %         clickStartTimes(i)=dt;
rmeddis@0:     %     end
rmeddis@0:     clickTime=round(period*clickStartTimes(i)/dt -dt);
rmeddis@0:     % clicks are treated as 20 microPascal high
rmeddis@0:     if clickTime+clickWidthNpoints<length(signal)
rmeddis@0:         signal(clickTime+1:clickTime+clickWidthNpoints)=clickHeight;
rmeddis@0:     end
rmeddis@0: end
rmeddis@0: 
rmeddis@0: clickTrain=repmat(signal, 1, repetitions);
rmeddis@0: 
rmeddis@0: if length(clickTrain)>totalLength
rmeddis@0:     clickTrain=clickTrain(1:totalLength);
rmeddis@0: elseif length(clickTrain)<totalLength
rmeddis@0:     timeToAdd=zeros(1,round((totalLength-length(clickTrain))));
rmeddis@0:     clickTrain=[clickTrain timeToAdd];
rmeddis@0:     % figure (1), plot(clickTrain)
rmeddis@0: end
rmeddis@0: 
rmeddis@0: %----------------------------------------------------------------makePressnitzerClicks
rmeddis@0: function signal=makePressnitzerClicks(globalStimParams, stimComponents)
rmeddis@0: % PressnitzerClicks(k,duration,dt)
rmeddis@0: % Generates a sequence of clicks with intervals kxx
rmeddis@0: %  where x= rand*k/2
rmeddis@0: % This is not the same as Kaernbach and Demany clicks
rmeddis@0: 
rmeddis@0: FS=globalStimParams.FS; % sample rate
rmeddis@0: dt=1/FS;
rmeddis@0: 
rmeddis@0: % obligatory parameter
rmeddis@0: try
rmeddis@0:     k=stimComponents.k;
rmeddis@0: catch
rmeddis@0:     error('PressnitzerClicks: field ''k'' is missing from stimcomponent')
rmeddis@0: end
rmeddis@0: 
rmeddis@0: % optional parameters
rmeddis@0: if isfield(stimComponents,'clickWidth')
rmeddis@0:     clickWidth=stimComponents.clickWidth;
rmeddis@0: else
rmeddis@0:     clickWidth=dt;
rmeddis@0: end
rmeddis@0: clickWidthNpoints=round(clickWidth*FS);
rmeddis@0: 
rmeddis@0: if isfield(stimComponents,'clickHeight')
rmeddis@0:     clickHeight=stimComponents.clickHeight;
rmeddis@0: else
rmeddis@0:     clickHeight=28e-6 * 10^(stimComponents.amplitudesdB/20);
rmeddis@0: end
rmeddis@0: 
rmeddis@0: duration=stimComponents.toneDuration;
rmeddis@0: 
rmeddis@0: signalLength=round(duration/dt);
rmeddis@0: signal=zeros(1,signalLength);
rmeddis@0: kInterval=round(k/dt);
rmeddis@0: halfk=k/2;
rmeddis@0: signal(1)=clickHeight;
rmeddis@0: timeIdx=0;
rmeddis@0: while timeIdx<signalLength
rmeddis@0:     % first interval = k
rmeddis@0:     clickTime=timeIdx+kInterval;
rmeddis@0:     signal(clickTime:clickTime+clickWidthNpoints)=clickHeight;
rmeddis@0:     timeIdx=timeIdx+kInterval;
rmeddis@0: 
rmeddis@0:     % second interval = 0 : k/2
rmeddis@0:     intervalx1=round(rand*halfk/dt);
rmeddis@0:     clickTime=timeIdx+intervalx1;
rmeddis@0:     signal(clickTime:clickTime+clickWidthNpoints)=clickHeight;
rmeddis@0:     timeIdx=timeIdx+intervalx1;
rmeddis@0: 
rmeddis@0:     % third interval = 0 : k/2
rmeddis@0:     intervalx1=round(rand*halfk/dt);
rmeddis@0:     clickTime=timeIdx+intervalx1;
rmeddis@0:     signal(clickTime:clickTime+clickWidthNpoints)=clickHeight;
rmeddis@0:     timeIdx=timeIdx+intervalx1;
rmeddis@0: end
rmeddis@0: 
rmeddis@0: signal=signal(1:signalLength);
rmeddis@0: % figure(1),	plot(dt:dt:duration,signal)
rmeddis@0: 
rmeddis@0: %----------------------------------------------------------------makePressnitzerABXClicks
rmeddis@0: function signal=makePressnitzerABxClicks(globalStimParams, stimComponents)
rmeddis@0: % Generates a sequence of clicks with intervals ABx
rmeddis@0: % AB interval is 2*k
rmeddis@0: % where A= rand* k
rmeddis@0: %       B= k-A
rmeddis@0: %       x= k/2
rmeddis@0: % These are second order clicks
rmeddis@0: 
rmeddis@0: FS=globalStimParams.FS; % sample rate
rmeddis@0: dt=1/FS;
rmeddis@0: 
rmeddis@0: % obligatory parameter
rmeddis@0: try
rmeddis@0:     k=stimComponents.k;
rmeddis@0: catch
rmeddis@0:     error('PressnitzerClicks: field ''k'' is missing from stimcomponent')
rmeddis@0: end
rmeddis@0: 
rmeddis@0: % optional parameters
rmeddis@0: if isfield(stimComponents,'clickWidth')
rmeddis@0:     clickWidth=stimComponents.clickWidth;
rmeddis@0: else
rmeddis@0:     clickWidth=dt;
rmeddis@0: end
rmeddis@0: clickWidthNpoints=round(clickWidth*FS);
rmeddis@0: 
rmeddis@0: if isfield(stimComponents,'clickHeight')
rmeddis@0:     clickHeight=stimComponents.clickHeight;
rmeddis@0: else
rmeddis@0:     clickHeight=28e-6 * 10^(stimComponents.amplitudesdB/20);
rmeddis@0: end
rmeddis@0: 
rmeddis@0: duration=stimComponents.toneDuration;
rmeddis@0: 
rmeddis@0: signalLength=round(duration/dt);
rmeddis@0: signal=zeros(1,2*signalLength); % allow for overrun
rmeddis@0: ABinterval=k/dt;                % i.e. the number of dt steps
rmeddis@0: randomInterval=ABinterval/2;
rmeddis@0: signal(1)=clickHeight;
rmeddis@0: time=0;
rmeddis@0: while time<signalLength
rmeddis@0:     % first interval = A
rmeddis@0:     intervalA=rand*ABinterval;
rmeddis@0:     clickTime=round(time+intervalA)+1;   % +1 avoids zero index
rmeddis@0:     signal(clickTime:clickTime+clickWidthNpoints)=clickHeight;
rmeddis@0:     time=time+intervalA;
rmeddis@0: 
rmeddis@0:     % second interval = B
rmeddis@0:     intervalB=ABinterval-intervalA;
rmeddis@0:     clickTime=round(time+intervalB)+1;
rmeddis@0:     signal(clickTime:clickTime+clickWidthNpoints-1)=clickHeight;
rmeddis@0:     time=time+intervalB;
rmeddis@0: 
rmeddis@0:     % third interval = 0 : k/2
rmeddis@0:     intervalx1=rand*randomInterval;    % mean random interval=k
rmeddis@0:     clickTime=round(time+intervalx1)+1;
rmeddis@0:     signal(clickTime:clickTime+clickWidthNpoints-1)=clickHeight;
rmeddis@0:     time=time+intervalx1;
rmeddis@0: end
rmeddis@0: 
rmeddis@0: signal=signal(1:signalLength);
rmeddis@0: % figure(1),	plot(dt:dt:duration,signal)
rmeddis@0: 
rmeddis@0: %-----------------------------------------------------makeABxClicks
rmeddis@0: function signal=makeABxClicks(globalStimParams, stimComponents)
rmeddis@0: % Generates a sequence of clicks with intervals ABx
rmeddis@0: % AB interval is 2*k
rmeddis@0: % where A= rand* k
rmeddis@0: %       B= k-A
rmeddis@0: %       x= rand*2*k
rmeddis@0: % These are second order clicks
rmeddis@0: 
rmeddis@0: FS=globalStimParams.FS; % sample rate
rmeddis@0: dt=1/FS;
rmeddis@0: 
rmeddis@0: % obligatory parameter
rmeddis@0: try
rmeddis@0:     k=stimComponents.k;
rmeddis@0: catch
rmeddis@0:     error('PressnitzerClicks: field ''k'' is missing from stimcomponent')
rmeddis@0: end
rmeddis@0: 
rmeddis@0: % optional parameters
rmeddis@0: if isfield(stimComponents,'clickWidth')
rmeddis@0:     clickWidth=stimComponents.clickWidth;
rmeddis@0: else
rmeddis@0:     clickWidth=dt;
rmeddis@0: end
rmeddis@0: clickWidthNpoints=round(clickWidth*FS);
rmeddis@0: 
rmeddis@0: if isfield(stimComponents,'clickHeight')
rmeddis@0:     clickHeight=stimComponents.clickHeight;
rmeddis@0: else
rmeddis@0:     clickHeight=28e-6 * 10^(stimComponents.amplitudesdB/20);
rmeddis@0: end
rmeddis@0: 
rmeddis@0: duration=stimComponents.toneDuration;
rmeddis@0: 
rmeddis@0: signalLength=round(duration/dt);
rmeddis@0: signal=zeros(1,2*signalLength); % allow for overrun
rmeddis@0: ABinterval=2*k/dt;
rmeddis@0: randomInterval=ABinterval;
rmeddis@0: signal(1)=clickHeight;
rmeddis@0: timeIdx=0;
rmeddis@0: while timeIdx<signalLength
rmeddis@0:     % first interval = A
rmeddis@0:     intervalA=round(rand*ABinterval);
rmeddis@0:     clickTime=timeIdx+intervalA+1;
rmeddis@0:     signal(clickTime:clickTime+clickWidthNpoints-1)=clickHeight;
rmeddis@0:     timeIdx=timeIdx+intervalA;
rmeddis@0: 
rmeddis@0:     % second interval = B
rmeddis@0:     intervalB=round(ABinterval-intervalA);
rmeddis@0:     clickTime=timeIdx+intervalB;
rmeddis@0:     signal(clickTime:clickTime+clickWidthNpoints-1)=clickHeight;
rmeddis@0:     timeIdx=timeIdx+intervalB;
rmeddis@0: 
rmeddis@0:     % third interval = 0 : k
rmeddis@0:     intervalx1=round(rand*randomInterval);    % mean random interval=k
rmeddis@0:     clickTime=timeIdx+intervalx1;
rmeddis@0:     signal(clickTime:clickTime+clickWidthNpoints-1)=clickHeight;
rmeddis@0:     timeIdx=timeIdx+intervalx1;
rmeddis@0: end
rmeddis@0: 
rmeddis@0: signal=signal(1:signalLength);
rmeddis@0: % figure(1),	plot(dt:dt:duration,signal)
rmeddis@0: 
rmeddis@0: %----------------------------------------------------------------makeYostClicks
rmeddis@0: function signal=makeYostClicks(globalStimParams, stimComponents)
rmeddis@0: % Generates a shuffled sequence of clicks with intervals kxxxx
rmeddis@0: %  where max(x)= 2*k
rmeddis@0: %  and there are n occurrences of x
rmeddis@0: % this section requires:
rmeddis@0: %  stimComponents.k
rmeddis@0: %  stimComponents.nXs
rmeddis@0: %  stimComponents.toneDuration
rmeddis@0: % optional:
rmeddis@0: %  stimComponents.clickWidth	%useful because width depends on dt
rmeddis@0: %  stimComponents.clickHeight	%best left to amplitude rescaling later
rmeddis@0: 
rmeddis@0: FS=globalStimParams.FS; % sample rate
rmeddis@0: dt=1/FS;
rmeddis@0: 
rmeddis@0: % obligatory parameters
rmeddis@0: try
rmeddis@0:     k=stimComponents.k;
rmeddis@0: catch
rmeddis@0:     error('makeYostClicks: field ''k'' is missing from stimComponents')
rmeddis@0: end
rmeddis@0: 
rmeddis@0: try
rmeddis@0:     nXs=stimComponents.nXs;
rmeddis@0: catch
rmeddis@0:     error('makeYostClicks: field ''nXs'' is missing from stimComponents')
rmeddis@0: end
rmeddis@0: 
rmeddis@0: try
rmeddis@0:     shuffled=stimComponents.shuffled;
rmeddis@0: catch
rmeddis@0:     error('makeYostClicks: field ''shuffled'' is missing from stimComponents')
rmeddis@0: end
rmeddis@0: 
rmeddis@0: try
rmeddis@0:     duration=stimComponents.toneDuration;
rmeddis@0: catch
rmeddis@0:     error('makeYostClicks: field ''toneDuration'' is missing from stimComponents')
rmeddis@0: end
rmeddis@0: 
rmeddis@0: % optional parameters
rmeddis@0: if isfield(stimComponents,'clickWidth')
rmeddis@0:     clickWidth=stimComponents.clickWidth;
rmeddis@0: else
rmeddis@0:     clickWidth=dt;
rmeddis@0: end
rmeddis@0: clickWidthNpoints=round(clickWidth*FS);
rmeddis@0: 
rmeddis@0: if isfield(stimComponents,'clickHeight')
rmeddis@0:     clickHeight=stimComponents.clickHeight;
rmeddis@0: else
rmeddis@0:     clickHeight=28e-6 * 10^(stimComponents.amplitudesdB/20);
rmeddis@0: end
rmeddis@0: 
rmeddis@0: kInterval=round(k/dt);
rmeddis@0: twoK=k*2;				% max width of x interval
rmeddis@0: 
rmeddis@0: signalLength=round(duration/dt);
rmeddis@0: signal=zeros(1,signalLength);
rmeddis@0: 
rmeddis@0: timeIdx=0;
rmeddis@0: intervalCount=0;
rmeddis@0: while timeIdx<signalLength
rmeddis@0:     timeIdx=timeIdx+kInterval;
rmeddis@0:     if timeIdx>signalLength, break,end
rmeddis@0:     intervalCount=intervalCount+1;
rmeddis@0:     intervals(intervalCount)=kInterval;
rmeddis@0: 
rmeddis@0:     % repeat x intervals as required
rmeddis@0:     if nXs>0
rmeddis@0:         for nX=1:nXs
rmeddis@0:             xInterval=round(rand*twoK/dt);
rmeddis@0:             timeIdx=timeIdx+xInterval;
rmeddis@0:             if timeIdx>signalLength, break,end
rmeddis@0:             intervalCount=intervalCount+1;
rmeddis@0:             intervals(intervalCount)=xInterval;
rmeddis@0:         end
rmeddis@0:     end
rmeddis@0:     if timeIdx>signalLength, break,end
rmeddis@0: end
rmeddis@0: 
rmeddis@0: % shuffle intervals
rmeddis@0: if shuffled
rmeddis@0:     randomNumbers=rand(1,length(intervals));
rmeddis@0:     [randomNumbers idx]=sort(randomNumbers);
rmeddis@0:     intervals=intervals(idx);
rmeddis@0:     idx=intervals>0;
rmeddis@0:     intervals=intervals(idx);
rmeddis@0: end
rmeddis@0: 
rmeddis@0: intervalCount=length(intervals);
rmeddis@0: signal(1)=clickHeight;
rmeddis@0: clickTime=0;
rmeddis@0: for i=1:intervalCount
rmeddis@0:     clickTime=clickTime+intervals(i);
rmeddis@0:     signal(clickTime:clickTime+clickWidthNpoints)=clickHeight;
rmeddis@0: end
rmeddis@0: signal=signal(1:signalLength);
rmeddis@0: %  figure(1),	plot(dt:dt:duration,signal)
rmeddis@0: 
rmeddis@0: %--------------------------------------------------------------------makeKxxClicks
rmeddis@0: function signal=makeKxxClicks(globalStimParams, stimComponents)
rmeddis@0: % Click train consists of kkxxx.. sequences
rmeddis@0: % k is the duration of a fixed interval (seconds)
rmeddis@0: % random intervals are distributed 0 : 2* k (NB not like Pressnitzer clicks)
rmeddis@0: % nKs is the number of successive 'k' intervals
rmeddis@0: % nXs is the number of random intervals between k sequences
rmeddis@0: % sequences of 3 x intervals > k are replaced with new sequences
rmeddis@0: % shuffled causes all intervals to be reordered randomly
rmeddis@0: % NB 1/k is the mean click rate
rmeddis@0: 
rmeddis@0: FS=globalStimParams.FS; % sample rate
rmeddis@0: dt=1/FS;
rmeddis@0: 
rmeddis@0: try
rmeddis@0:     k=stimComponents.k;		% duration (s) of fixed intervals
rmeddis@0: catch
rmeddis@0:     error('makeYostClicks: field ''k'' is missing from stimComponents')
rmeddis@0: end
rmeddis@0: 
rmeddis@0: try
rmeddis@0:     duration=stimComponents.toneDuration;
rmeddis@0: catch
rmeddis@0:     error('makeYostClicks: field ''duration'' is missing from stimComponents')
rmeddis@0: end
rmeddis@0: 
rmeddis@0: if isfield(stimComponents,'clickWidth')
rmeddis@0:     clickWidth=stimComponents.clickWidth;
rmeddis@0: else
rmeddis@0:     clickWidth=dt;
rmeddis@0: end
rmeddis@0: 
rmeddis@0: if isfield(stimComponents,'clickHeight')
rmeddis@0:     clickHeight=stimComponents.clickHeight;
rmeddis@0: else
rmeddis@0:     clickHeight=28e-6 * 10^(stimComponents.amplitudesdB/20);
rmeddis@0: end
rmeddis@0: 
rmeddis@0: 
rmeddis@0: if isfield(stimComponents,'order')
rmeddis@0:     order=stimComponents.order;
rmeddis@0: else
rmeddis@0:     order=1;
rmeddis@0: end
rmeddis@0: 
rmeddis@0: if isfield(stimComponents,'nKs')
rmeddis@0:     nKs=stimComponents.nKs;
rmeddis@0: else
rmeddis@0:     nKs=1;
rmeddis@0: end
rmeddis@0: 
rmeddis@0: if isfield(stimComponents,'nXs')
rmeddis@0:     nXs=stimComponents.nXs;
rmeddis@0: else
rmeddis@0:     nXs=1;
rmeddis@0: end
rmeddis@0: 
rmeddis@0: if isfield(stimComponents,'shuffled')
rmeddis@0:     shuffled=stimComponents.shuffled;
rmeddis@0: else
rmeddis@0:     shuffled=1;
rmeddis@0: end
rmeddis@0: 
rmeddis@0: kLength=round(k/dt);    % fixed interval
rmeddis@0: xLength=2*kLength;      % maximum random interval
rmeddis@0: requiredSignalLength=round(duration/dt);
rmeddis@0: intervalsPerCycle=(nKs+nXs);
rmeddis@0: cycleLength=nKs*kLength+nXs*xLength;
rmeddis@0: % more cycles to allow for uncertainty
rmeddis@0: nCycles=5*round(requiredSignalLength/cycleLength);
rmeddis@0: nIntervals=nCycles*intervalsPerCycle;
rmeddis@0: 
rmeddis@0: % random intervals
rmeddis@0: if nXs>0
rmeddis@0:     xIntervals=floor(rand(1,nIntervals)*2*kLength);
rmeddis@0:     % weed out triple intervals > 2*k
rmeddis@0:     rogues=1;
rmeddis@0:     while sum(rogues)
rmeddis@0:         y=(xIntervals>kLength);
rmeddis@0:         rogues=(sum([y(1:end-2)' y(2:end-1)' y(3:end)'],2)>2);
rmeddis@0:         xIntervals(rogues)=floor(rand*2*kLength);
rmeddis@0:     end
rmeddis@0:     xIntervals=reshape(xIntervals,nCycles,intervalsPerCycle);
rmeddis@0: else
rmeddis@0:     xIntervals=[];
rmeddis@0: end
rmeddis@0: 
rmeddis@0: % insert constrained (k) intervals
rmeddis@0: if nKs>0
rmeddis@0:     switch order
rmeddis@0:         case 1
rmeddis@0:             kIntervals=floor(ones(nCycles,nKs)*kLength);
rmeddis@0:         case 2
rmeddis@0:             nKs=1; % force this
rmeddis@0:             kIntervals=floor(rand(nCycles,1)*kLength);
rmeddis@0:             kIntervals=[kIntervals kLength-kIntervals];
rmeddis@0:     end
rmeddis@0: else
rmeddis@0:     kIntervals=[];
rmeddis@0: end
rmeddis@0: 
rmeddis@0: % combine fixed and random
rmeddis@0: intervals=[kIntervals xIntervals(:,nKs+1:end)];
rmeddis@0: % make a single array;
rmeddis@0: [r c]=size(intervals);
rmeddis@0: intervals=reshape(intervals',1,r*c);
rmeddis@0: 
rmeddis@0: % shuffle intervals
rmeddis@0: if shuffled
rmeddis@0:     randomNumbers=rand(1,length(intervals));
rmeddis@0:     [randomNumbers idx]=sort(randomNumbers);
rmeddis@0:     intervals=intervals(idx);
rmeddis@0:     idx=intervals>0;
rmeddis@0:     intervals=intervals(idx);
rmeddis@0: end
rmeddis@0: 
rmeddis@0: % convert intervals to clicks
rmeddis@0: clickTimes=cumsum(intervals);
rmeddis@0: signal(clickTimes)=clickHeight;
rmeddis@0: signal=signal(1:requiredSignalLength);
rmeddis@0: % figure(1), clf, plot(signal)
rmeddis@0: 
rmeddis@0: 
rmeddis@0: 
rmeddis@0: %--------------------------------------------------------------------makeNoise
rmeddis@0: function noise=makeNoise(globalStimParams, stimComponents)
rmeddis@0: % FS in Hz, noiseDuration in s, delay in s;
rmeddis@0: % noise is returned with overall level dB(rms) = amplitudesdB
rmeddis@0: %
rmeddis@0: % % You need
rmeddis@0: %
rmeddis@0: % stim.type='noise'; % or 'IRN', or 'pinkNoise'
rmeddis@0: % stim.toneDuration=.05;
rmeddis@0: % stim.amplitudesdB=50;
rmeddis@0: % stim.beginSilence=.01;
rmeddis@0: % stim.endSilence=-1;
rmeddis@0: % stim.rampOnDur=.002;
rmeddis@0: % stim.rampOffDur=-1;
rmeddis@0: %
rmeddis@0: % % Mandatory structure fields
rmeddis@0: % globalStimParams.FS=100000;
rmeddis@0: % globalStimParams.overallDuration=.1;  % s
rmeddis@0: % globalStimParams.doPlot=1;
rmeddis@0: % globalStimParams.doPlay=1;
rmeddis@0: %
rmeddis@0: % [audio, msg]=stimulusCreate(globalStimParams, stim, );
rmeddis@0: %
rmeddis@0: % % local
rmeddis@0: % stim.type='noise'; % or 'IRN'
rmeddis@0: %
rmeddis@0: FS=globalStimParams.FS;
rmeddis@0: noiseDuration= stimComponents.toneDuration;
rmeddis@0: npts=round(noiseDuration*FS);
rmeddis@0: noise=randn(1,npts);    % NB randn (normally distributed)
rmeddis@0: 
rmeddis@0: switch stimComponents.type
rmeddis@0:     case  'pinkNoise'
rmeddis@0:         %         noise=UTIL_lowpassFilterFreq(noise, 100, 1/FS);
rmeddis@0:         noise=UTIL_bandPassFilter(noise, 1, 100, 200, 1/FS,[]);
rmeddis@0: end
rmeddis@0: 
rmeddis@0: rms=(mean(noise.^2)^.5);  %should be 20 microPascals for 0 dB SPL
rmeddis@0: adjust=20e-6/rms;
rmeddis@0: noise=noise*adjust;
rmeddis@0: rms=(mean(noise.^2)^.5);
rmeddis@0: amplitude=10.^(stimComponents.amplitudesdB/20);
rmeddis@0: noise=amplitude*noise;
rmeddis@0: % rms=(mean(noise.^2)^.5);
rmeddis@0: % dBnoise=20*log10(rms/20e-6)
rmeddis@0: 
rmeddis@0: 
rmeddis@0: %--------------------------------------------------------------------makeWhiteNoise
rmeddis@0: function noise=makeWhiteNoise(globalStimParams, stimComponents)
rmeddis@0: % FS in Hz, noiseDuration in s, delay in s;
rmeddis@0: % noise is bandpass filtered between 100 and 10000 Hz
rmeddis@0: % spectrum level (dB/Hz) is 40 dB below nominal level.
rmeddis@0: % noise is returned with dB(rms) = amplitudesdB
rmeddis@0: %
rmeddis@0: % % You need
rmeddis@0: %
rmeddis@0: % stim.type='noise'; % or 'IRN', or 'pinkNoise'
rmeddis@0: % stim.toneDuration=.05;
rmeddis@0: % stim.amplitudesdB=50;
rmeddis@0: % stim.beginSilence=.01;
rmeddis@0: % stim.endSilence=-1;
rmeddis@0: % stim.rampOnDur=.002;
rmeddis@0: % stim.rampOffDur=-1;
rmeddis@0: %
rmeddis@0: % % Mandatory structure fields
rmeddis@0: % globalStimParams.FS=100000;
rmeddis@0: % globalStimParams.overallDuration=.1;  % s
rmeddis@0: % globalStimParams.doPlot=1;
rmeddis@0: % globalStimParams.doPlay=1;
rmeddis@0: %
rmeddis@0: % [audio, msg]=stimulusCreate(globalStimParams, stim, );
rmeddis@0: %
rmeddis@0: % % local
rmeddis@0: % stim.type='noise'; % or 'IRN'
rmeddis@0: %
rmeddis@0: FS=globalStimParams.FS;
rmeddis@0: noiseDuration= stimComponents.toneDuration;
rmeddis@0: npts=round(noiseDuration*FS);
rmeddis@0: noise=randn(1,npts);
rmeddis@0: 
rmeddis@0: noise=UTIL_bandPassFilter (noise, 6, 100, 10000, 1/FS, []);
rmeddis@0: 
rmeddis@0: rms=(mean(noise.^2)^.5);  %should be 20 microPascals for 0 dB SPL
rmeddis@0: adjust=20e-6/rms;
rmeddis@0: noise=noise*adjust;
rmeddis@0: rms=(mean(noise.^2)^.5);
rmeddis@0: amplitude=10.^(stimComponents.amplitudesdB/20);
rmeddis@0: noise=amplitude*noise;
rmeddis@0: % rms=(mean(noise.^2)^.5);
rmeddis@0: % dBnoise=20*log10(rms/20e-6)
rmeddis@0: 
rmeddis@0: 
rmeddis@0: %-----------------------------------------------------------------makeIRN
rmeddis@0: function noise=makeIRN(globalStimParams, stimComponents)
rmeddis@0: % FS in Hz, noiseDuration in s, delay in s;
rmeddis@0: % noise is returned with dB(rms) = amplitudesdB
rmeddis@0: %
rmeddis@0: % % You need
rmeddis@0: %
rmeddis@0: % stim.type='noise'; % or 'IRN', or 'pinkNoise'
rmeddis@0: % stim.toneDuration=.05;
rmeddis@0: % stim.amplitudesdB=50;
rmeddis@0: % stim.beginSilence=.01;
rmeddis@0: % stim.endSilence=-1;
rmeddis@0: % stim.rampOnDur=.002;
rmeddis@0: % stim.rampOffDur=-1;
rmeddis@0: %
rmeddis@0: % % Mandatory structure fields
rmeddis@0: % globalStimParams.FS=100000;
rmeddis@0: % globalStimParams.overallDuration=.1;  % s
rmeddis@0: % globalStimParams.doPlot=1;
rmeddis@0: % globalStimParams.doPlay=1;
rmeddis@0: %
rmeddis@0: % [audio, msg]=stimulusCreate(globalStimParams, stim, );
rmeddis@0: %
rmeddis@0: % % local
rmeddis@0: % stim.type='noise'; % or 'IRN'
rmeddis@0: % % for IRN only
rmeddis@0: % stim.niterations = 8;   %0 for white noise
rmeddis@0: % stim.delay = 1/150;
rmeddis@0: % stim.irnGain = 1;
rmeddis@0: %
rmeddis@0: FS=globalStimParams.FS;
rmeddis@0: noiseDuration= stimComponents.toneDuration;
rmeddis@0: 
rmeddis@0: nIterations=stimComponents.niterations;
rmeddis@0: if nIterations==0
rmeddis@0:     % white noise is specified as nIterations=1
rmeddis@0:     nIterations=1;
rmeddis@0:     IRNgain=0;
rmeddis@0:     delay=0.01; % dummy
rmeddis@0: else
rmeddis@0:     % IRN
rmeddis@0:     delay=stimComponents.delay;
rmeddis@0:     IRNgain=stimComponents.irnGain;
rmeddis@0: end
rmeddis@0: 
rmeddis@0: npts=round(noiseDuration*FS);
rmeddis@0: dels=round(delay*FS);
rmeddis@0: noise=randn(1,npts);
rmeddis@0: 
rmeddis@0: %fringe=nIterations*dels;
rmeddis@0: %npts=npts+fringe;
rmeddis@0: 
rmeddis@0: for i=1:nIterations,
rmeddis@0:     dnoise=[noise(dels+1:npts) noise(1:dels)];
rmeddis@0:     dnoise=dnoise.*IRNgain;
rmeddis@0:     noise=noise+dnoise;
rmeddis@0: end;
rmeddis@0: 
rmeddis@0: switch stimComponents.type
rmeddis@0:     case  'pinkNoise'
rmeddis@0:         noise=UTIL_lowpassFilterFreq(noise, 10000, 1/FS);
rmeddis@0: end
rmeddis@0: 
rmeddis@0: rms=(mean(noise.^2)^.5);  %should be 20 microPascals for 0 dB SPL
rmeddis@0: adjust=20e-6/rms;
rmeddis@0: noise=noise*adjust;
rmeddis@0: rms=(mean(noise.^2)^.5);
rmeddis@0: amplitude=10.^(stimComponents.amplitudesdB/20);
rmeddis@0: noise=amplitude*noise;
rmeddis@0: % rms=(mean(noise.^2)^.5);
rmeddis@0: % dBnoise=20*log10(rms/20e-6)
rmeddis@0: 
rmeddis@0: %------------------------------------------------------------------ makeRPN
rmeddis@0: function RPN=makeRPN(globalStimParams, stim)
rmeddis@0: % 'period' is a collection of samples - AAABCD
rmeddis@0: % you need
rmeddis@0: %
rmeddis@0: % stim.type='RPN';
rmeddis@0: % stim.toneDuration=.2;
rmeddis@0: % stim.amplitudesdB=50;
rmeddis@0: % stim.beginSilence=.01;
rmeddis@0: % stim.rampOnDur=.002;
rmeddis@0: % stim.rampOffDur=-1;
rmeddis@0: %
rmeddis@0: % stim.sampleDuration=.005;  %200 Hz pitch
rmeddis@0: % stim.nSimilarSamples=5;   % pitch strength
rmeddis@0: % stim.nIndependentSamples=1% dilutes strength
rmeddis@0: %
rmeddis@0: % % Mandatory structure fields
rmeddis@0: %  globalStimParams.FS=44100;
rmeddis@0: %  globalStimParams.overallDuration=.21;  % s
rmeddis@0: %
rmeddis@0: % globalStimParams.doPlot=1;
rmeddis@0: % globalStimParams.doPlay=1;
rmeddis@0: % [audio, msg]=stimulusCreate(globalStimParams, stim);
rmeddis@0: 
rmeddis@0: FS=globalStimParams.FS;
rmeddis@0: ptsPerSample=floor(stim.sampleDuration*FS);
rmeddis@0: 
rmeddis@0: samplesPerPeriod=stim.nSimilarSamples+stim.nIndependentSamples;
rmeddis@0: periodDuration=samplesPerPeriod*stim.sampleDuration;
rmeddis@0: 
rmeddis@0: totalNumPeriods=2*floor(stim.toneDuration/periodDuration);  % longer than necessary
rmeddis@0: if totalNumPeriods<1
rmeddis@0:     error('stimulusCreate: RPN, stimulus duration needs to be longer')
rmeddis@0: end
rmeddis@0: 
rmeddis@0: RPN=[];
rmeddis@0: for j=1:totalNumPeriods
rmeddis@0:     noise=randn(1,ptsPerSample);
rmeddis@0:     for i=1:stim.nSimilarSamples
rmeddis@0:         RPN=[RPN noise];
rmeddis@0:     end
rmeddis@0: 
rmeddis@0:     for i=1:stim.nIndependentSamples
rmeddis@0:         noise=randn(1,ptsPerSample);
rmeddis@0:         RPN=[RPN noise];
rmeddis@0:     end
rmeddis@0: end
rmeddis@0: 
rmeddis@0: targetStimulusLength=round(stim.toneDuration/FS);
rmeddis@0: RPN=RPN(1:floor(stim.toneDuration*FS));     % take enough for stimulus
rmeddis@0: 
rmeddis@0: rms=(mean(RPN.^2)^.5);  %should be 20 microPascals for 0 dB SPL
rmeddis@0: adjust=20e-6/rms;
rmeddis@0: RPN=RPN*adjust;
rmeddis@0: rms=(mean(RPN.^2)^.5);
rmeddis@0: amplitude=10.^(stim.amplitudesdB/20);
rmeddis@0: RPN=amplitude*RPN;
rmeddis@0: % rms=(mean(noise.^2)^.5);
rmeddis@0: % dBnoise=20*log10(rms/20e-6)
rmeddis@0: 
rmeddis@0: %--------------------------------------------------------------------applyRampOn
rmeddis@0: function signal=applyRampOn(signal, rampDur, rampOnTime, sampleRate)
rmeddis@0: %applyRampOn applies raised cosine ramp
rmeddis@0: %rampOntime is the time at which the ramp begins
rmeddis@0: %At all other times the mask has a value of 1
rmeddis@0: %signal=applyRampOn(signal, rampDur, rampOnTime, sampleRate)
rmeddis@0: 
rmeddis@0: rampDurPoints=round(rampDur*sampleRate);
rmeddis@0: rampOn= (1+cos(pi:pi/(rampDurPoints-1): 2*pi))/2';
rmeddis@0: 
rmeddis@0: sigDurPoints=length(signal);
rmeddis@0: mask(1:sigDurPoints)=1;
rmeddis@0: rampOnStartIndex=round(rampOnTime*sampleRate+1);
rmeddis@0: mask(rampOnStartIndex: rampOnStartIndex+ rampDurPoints-1)=rampOn;
rmeddis@0: signal=signal.*mask;
rmeddis@0: %plot(mask)
rmeddis@0: 
rmeddis@0: %--------------------------------------------------------------------applyRampOff
rmeddis@0: function signal=applyRampOff(signal, rampDur, rampOffTime, sampleRate)
rmeddis@0: %applyRampOn applies raised cosine squared ramp
rmeddis@0: %rampOffTime is the time at which the ramp begins
rmeddis@0: %At all other times the mask has a value of 1
rmeddis@0: % signal=applyRampOff(signal, rampDur, rampOffTime, sampleRate)
rmeddis@0: 
rmeddis@0: rampDurPoints=round(rampDur*sampleRate);
rmeddis@0: rampOff= (1+cos(0:pi/(rampDurPoints-1): pi))/2';
rmeddis@0: 
rmeddis@0: sigDurPoints=length(signal);
rmeddis@0: mask=ones(1,sigDurPoints);
rmeddis@0: rampOffStartIndex=round(rampOffTime*sampleRate+1);
rmeddis@0: mask(rampOffStartIndex: round(rampOffStartIndex+ rampDurPoints-1))=rampOff;
rmeddis@0: if length(mask)>sigDurPoints, mask=mask(1:sigDurPoints); end
rmeddis@0: signal=signal.*mask;
rmeddis@0: %plot(mask)
rmeddis@0: 
rmeddis@0: function signal=applyGaussianRamps(signal, sigma, sampleRate)
rmeddis@0: dt=1/sampleRate;
rmeddis@0: time=dt:dt:dt*length(signal);
rmeddis@0: ramp=1-exp(-time.^2/(2*sigma^2));
rmeddis@0: % apply onset ramp
rmeddis@0: signal=signal.*ramp;
rmeddis@0: % apply offset ramp
rmeddis@0: ramp=fliplr(ramp);
rmeddis@0: signal=signal.*ramp;
rmeddis@0: 
rmeddis@0: 
rmeddis@0: 
rmeddis@0: %--------------------------------------------------------------------checkDescriptors
rmeddis@30: function [globalStimParams, stimComponents]=...
rmeddis@30:     checkDescriptors(globalStimParams, stimComponents)
rmeddis@0: 
rmeddis@0: try
rmeddis@0:     % if FS exists, it takes priority
rmeddis@0:     globalStimParams.dt=1/globalStimParams.FS;
rmeddis@0: catch
rmeddis@0:     % otherwise set FS using dt
rmeddis@0:     globalStimParams.FS=1/globalStimParams.dt;
rmeddis@0: end
rmeddis@0: 
rmeddis@30: sampleRate=globalStimParams.FS;
rmeddis@30: 
rmeddis@30: globalStimParams.nSignalPoints=...
rmeddis@30:     round(globalStimParams.overallDuration* sampleRate);
rmeddis@0: 
rmeddis@0: % optional field (ears)
rmeddis@0: try
rmeddis@0:     globalStimParams.ears;
rmeddis@0: catch
rmeddis@0:     % default: dichotic.
rmeddis@0:     globalStimParams.ears='dichotic';
rmeddis@0: end
rmeddis@0: 
rmeddis@0: % audioOutCorrection is optional
rmeddis@0: % audioOutCorrection is a scalar for reducing the sound
rmeddis@0: try
rmeddis@0:     globalStimParams.audioOutCorrection;
rmeddis@0: catch
rmeddis@0:     %      set to 1 if omitted
rmeddis@0:     globalStimParams.audioOutCorrection=1;
rmeddis@0: end
rmeddis@0: 
rmeddis@0: try
rmeddis@0:     globalStimParams.doPlay;
rmeddis@0: catch
rmeddis@0:     % default plays sound only if explicitly requested
rmeddis@0:     globalStimParams.doPlay=0;
rmeddis@0: end
rmeddis@0: 
rmeddis@0: try
rmeddis@0:     globalStimParams.doPlot;
rmeddis@0: catch
rmeddis@0:     % no plotting unless explicitly requested
rmeddis@0:     globalStimParams.doPlot=0;
rmeddis@0: end
rmeddis@0: 
rmeddis@0: [ears nComponentSounds]=size(stimComponents);
rmeddis@0: 
rmeddis@0: for ear=1:2 % 1=left/ 2=right
rmeddis@0: 
rmeddis@0:     % create a list of components whose type is specified
rmeddis@0:     % if no type is specified assume that it is an empty specification
rmeddis@0:     % this is allowed
rmeddis@0:     validStimComponents=[];
rmeddis@0:     for i=1:nComponentSounds
rmeddis@0:         try
rmeddis@0:             if ~isempty(stimComponents(ear,i).type)
rmeddis@0:                 validStimComponents=[validStimComponents i];
rmeddis@0:             end
rmeddis@0:         catch
rmeddis@0:         end
rmeddis@0:     end
rmeddis@0: 
rmeddis@0:     for componentNo=validStimComponents
rmeddis@0:         % If no AM filed is present, create it for completeness
rmeddis@0:         if ~isfield(stimComponents(ear,componentNo),'AMfrequency') |...
rmeddis@0:                 ~isfield(stimComponents(ear,componentNo),'AMdepth')
rmeddis@0:             stimComponents(ear,componentNo).AMfrequency=0;
rmeddis@0:             stimComponents(ear,componentNo).AMdepth=0;
rmeddis@0:         end
rmeddis@0: 
rmeddis@0:         % all signals must have durations, amplitudes and ramps
rmeddis@0:         if ...
rmeddis@0:                 isempty(stimComponents(ear,componentNo).type) |...
rmeddis@0:                 isempty(stimComponents(ear,componentNo).toneDuration) |...
rmeddis@0:                 isempty(stimComponents(ear,componentNo).amplitudesdB) |...
rmeddis@0:                 isempty(stimComponents(ear,componentNo).rampOnDur)
rmeddis@0:             descriptorError( 'missing stimComponent descriptor', stimComponents, ear, componentNo)
rmeddis@0:         end
rmeddis@0: 
rmeddis@0:         try, stimComponents(ear,componentNo).endSilence; catch, stimComponents(ear,componentNo).endSilence=-1; end
rmeddis@0: 
rmeddis@0:         % ramp checks do not apply to file input
rmeddis@0:         if ~strcmp(stimComponents(ear,componentNo).type, 'file')
rmeddis@0:             % match offset ramp to onset if not explicitly specified
rmeddis@0:             if stimComponents(ear,componentNo).rampOffDur==-1,
rmeddis@0:                 stimComponents(ear,componentNo).rampOffDur=stimComponents(ear,componentNo).rampOnDur;
rmeddis@0:             end
rmeddis@0:             % ramps must be shorter than the stimulus
rmeddis@0:             if stimComponents(ear,componentNo).rampOffDur> stimComponents(ear,componentNo).toneDuration | ...
rmeddis@0:                     stimComponents(ear,componentNo).rampOnDur> stimComponents(ear,componentNo).toneDuration
rmeddis@0:                 descriptorError( 'ramp longer than sound component', stimComponents, ear, componentNo)
rmeddis@0:             end
rmeddis@0:         end
rmeddis@0: 
rmeddis@0:         % end silence is measured to fit into the global duration
rmeddis@0:         if stimComponents(ear,componentNo).endSilence==-1,
rmeddis@30:             stimComponents(ear,componentNo).endSilence=...
rmeddis@30:                 globalStimParams.overallDuration-...
rmeddis@30:                 stimComponents(ear,componentNo).beginSilence -...
rmeddis@30:                 stimComponents(ear,componentNo).toneDuration;
rmeddis@30:             
rmeddis@30:             endSilenceNpoints=stimComponents(ear,componentNo).endSilence...
rmeddis@30:                 *sampleRate;
rmeddis@30:         end
rmeddis@30:         if stimComponents(ear,componentNo).endSilence<0
rmeddis@30:             globalStimParams
rmeddis@30:             descriptorError( 'component durations greater than overallDuration', stimComponents, ear, componentNo)
rmeddis@0:         end
rmeddis@0: 
rmeddis@0:         % check overall duration of this component against global duration
rmeddis@0:         totalDuration= ...
rmeddis@0:             stimComponents(ear,componentNo).beginSilence+...
rmeddis@0:             stimComponents(ear,componentNo).toneDuration+...
rmeddis@0:             stimComponents(ear,componentNo).endSilence;
rmeddis@0: 
rmeddis@0:         % avoid annoying error message for single stimulus component
rmeddis@0:         if ears==1 && nComponentSounds==1
rmeddis@0:             globalStimParams.overallDuration= totalDuration;
rmeddis@0:         end
rmeddis@0: 
rmeddis@0:         % check total duration
rmeddis@30:         totalSignalPoints= round(totalDuration* sampleRate);
rmeddis@0:         if totalSignalPoints  >globalStimParams.nSignalPoints
rmeddis@0:             descriptorError( 'Signal component duration does not match overall duration ', stimComponents, ear, componentNo)
rmeddis@0:         end
rmeddis@0: 
rmeddis@0:         if isfield(stimComponents(ear,componentNo), 'filter')
rmeddis@0:             if ~isequal(length(stimComponents(ear,componentNo).filter), 3)
rmeddis@0:                 descriptorError( 'Filter parameter must have three elements ', stimComponents, ear, componentNo)
rmeddis@0:             end
rmeddis@0:         end
rmeddis@0:     end % component
rmeddis@0:     % ??
rmeddis@0:     if strcmp(globalStimParams.ears,'monoticL') | strcmp(globalStimParams.ears, 'monoticR'), break, end
rmeddis@0: end		% ear
rmeddis@0: 
rmeddis@0: 
rmeddis@0: %-------------------------------------------------------------------- descriptorError
rmeddis@0: function descriptorError( msg, stimComponents, ear, componentNo)
rmeddis@0: stimComponents(ear, componentNo)
rmeddis@0: 
rmeddis@0: disp(' *********** **************** ************')
rmeddis@0: disp([ '...Error in stimComponents description: '])
rmeddis@0: disp([msg ])
rmeddis@0: disp(['Ear = ' num2str(ear) ' component No ' num2str(componentNo)])
rmeddis@0: disp(' *********** **************** ************')
rmeddis@0: error('myError ')
rmeddis@0: 
rmeddis@0: 
rmeddis@0: %-------------------------------------------------------------------- normalize
rmeddis@0: function [normalizedSignal, gain]= normalize(signal)
rmeddis@0: % normalize (signal)
rmeddis@0: maxSignal=max(max(signal));
rmeddis@0: minSignal=min(min(signal));
rmeddis@0: if -minSignal>maxSignal, normalizingFactor=-minSignal; else normalizingFactor=maxSignal; end
rmeddis@0: normalizingFactor=1.01*normalizingFactor;
rmeddis@0: gain= 20*log10(normalizingFactor);
rmeddis@0: normalizedSignal=signal/normalizingFactor;
rmeddis@0: 
rmeddis@0: 
rmeddis@0: %--------------------------------------------------------------------Butterworth
rmeddis@0: function y=Butterworth (x, dt, fl, fu, order)
rmeddis@0: % Butterworth (x, dt, fu, fl, order)
rmeddis@0: % Taken from Yuel and beauchamp page 261
rmeddis@0: % NB error in their table for K (see their text)
rmeddis@0: % x is original signal
rmeddis@0: % fu, fl upper and lower cutoff
rmeddis@0: % order is the number of times the filter is applied
rmeddis@0: 
rmeddis@0: q=(pi*dt*(fu-fl));
rmeddis@0: J=1/(1+ cot(q));
rmeddis@0: K= (2*cos(pi*dt*(fu+fl)))/(1+tan(q)*cos(q));
rmeddis@0: L= (tan(q)-1)/(tan(q)+1);
rmeddis@0: b=[J -J];
rmeddis@0: a=[1 -K  -L];
rmeddis@0: for i=1:order
rmeddis@0:     y=filter(b, a, x);
rmeddis@0:     x=y;
rmeddis@0: end
rmeddis@0: 
rmeddis@0: 
rmeddis@0: % -------------------------------------------------------- UTIL_amp2dB
rmeddis@0: function [y] = UTIL_amp2dB (x, ref)
rmeddis@0: % Calculates a dB (ref. ref) value 'y' from a peak amplitude number 'x'.
rmeddis@0: % if ref omitted treat as dB
rmeddis@0: % Check the number of arguments that are passed in.
rmeddis@0: if (nargin < 2)
rmeddis@0:     ref=1;
rmeddis@0: end
rmeddis@0: if (nargin > 2)
rmeddis@0:     error ('Too many arguments');
rmeddis@0: end
rmeddis@0: 
rmeddis@0: % Check arguments.
rmeddis@0: if x < 0.0
rmeddis@0:     error ('Can not calculate the log10 of a negative number');
rmeddis@0: elseif x == 0.0
rmeddis@0:     warning ('log10 of zero.  The result is set to -1000.0');
rmeddis@0:     y = -1000.0;
rmeddis@0: else
rmeddis@0:     % Do calculations.
rmeddis@0:     y = 20.0 * log10(x/(sqrt(2)*ref));
rmeddis@0: 
rmeddis@0: end
rmeddis@0: 
rmeddis@0: %-------------------------------------------------------------------- FFT
rmeddis@0: function showFFT (getFFTsignal, dt)
rmeddis@0: color='r';
rmeddis@0: figure(2), clf,
rmeddis@0: hold off
rmeddis@0: 
rmeddis@0: % trim initial silence
rmeddis@0: idx=find(getFFTsignal>0);
rmeddis@0: if ~isempty(idx)
rmeddis@0:     getFFTsignal=getFFTsignal(idx(1):end);
rmeddis@0: end
rmeddis@0: %trim final silence
rmeddis@0: getFFTsignal=getFFTsignal(end:-1:1);
rmeddis@0: idx=find(getFFTsignal>0);
rmeddis@0: if ~isempty(idx)
rmeddis@0:     getFFTsignal=getFFTsignal(idx(1):end);
rmeddis@0:     getFFTsignal=getFFTsignal(end:-1:1);
rmeddis@0: end
rmeddis@0: 
rmeddis@0: % Analyse make stimComponents length a power of 2
rmeddis@0: x=length(getFFTsignal);
rmeddis@0: squareLength=2;
rmeddis@0: while squareLength<=x
rmeddis@0:     squareLength=squareLength*2;
rmeddis@0: end
rmeddis@0: squareLength=round(squareLength/2);
rmeddis@0: getFFTsignal=getFFTsignal(1:squareLength);
rmeddis@0: n=length(getFFTsignal);
rmeddis@0: 
rmeddis@0: minf=100; maxf=20000;
rmeddis@0: 
rmeddis@0: fft_result = fft(getFFTsignal, n);				% Compute FFT of the input signal.
rmeddis@0: fft_power = fft_result .* conj(fft_result);% / n;	% Compute power spectrum.  Dividing by 'n' we get the power spectral density.
rmeddis@0: fft_phase = angle(fft_result);			% Compute the phase spectrum.
rmeddis@0: 
rmeddis@0: frequencies = (1/dt)*(1:n/2)/n;
rmeddis@0: fft_power=fft_power(1:length(fft_power)/2); % remove mirror frequencies
rmeddis@0: fft_phase=fft_phase(1:length(fft_phase)/2); % remove mirror frequencies
rmeddis@0: fft_powerdB = UTIL_amp2dB (fft_power, max(fft_power)); % convert to dB
rmeddis@0: %     jags=find(diff(fft_phase)>0); % unwrap phase
rmeddis@0: %     for i=jags, fft_phase(i+1:end)=fft_phase(i+1:end)-2*pi; end
rmeddis@0: 
rmeddis@0: xlim([minf maxf])
rmeddis@0: semilogx(frequencies, fft_powerdB-max(fft_powerdB), color)
rmeddis@0: ylim([ -20 5])
rmeddis@0: 
rmeddis@0: 
rmeddis@0: 
rmeddis@0: function y=UTIL_lowpassFilterFreq(x, cutOffFrequency, dt)
rmeddis@0: % UTIL_lowpassFilterFreq multi-channel filter
rmeddis@0: %
rmeddis@0: % Usage:
rmeddis@0: % output=UTIL_lowpassFilterFreq(input, cutOffFrequency, dt)
rmeddis@0: %
rmeddis@0: % cutoff should be around 100 Hz for audio work
rmeddis@0: % dt should be <1/50000 s for audio work
rmeddis@0: %
rmeddis@0: % Attenuation 	       is - 6 dB per octave above cutoff.
rmeddis@0: 
rmeddis@0: 
rmeddis@0: sampleRate=1/dt;
rmeddis@0: 
rmeddis@0: if 4*cutOffFrequency>sampleRate
rmeddis@0:     warning(['UTIL_lowpassFilterFreq: sample rate ' num2str(1/dt) ' is too low for this BF.  Sampling rate should be >' num2str(4*cutOffFrequency) 'or cutoff (' num2str(4*cutOffFrequency) ') should be lower' ])
rmeddis@0:     cutOffFrequency=sampleRate/4;
rmeddis@0: end
rmeddis@0: 
rmeddis@0: tau=1/(2*pi*cutOffFrequency);
rmeddis@0: 
rmeddis@0: y=zeros(size(x));
rmeddis@0: [numChannels numPoints]=size(x);
rmeddis@0: for i=1:numChannels
rmeddis@0:     y(i,:)=filter([dt/tau], [1 -(1-dt/tau)], x(i,:));
rmeddis@0: end
rmeddis@0: 
rmeddis@0: