MAP

function [audio, msg, time]=stimulusCreate(globalStimParams,...

    stimComponents, doPlot)

% updated June 2007

% the only authorised version of stimulus create is the version to be found

% in MAP1_6.  Other versions are copies!!

%

% for a simple tone you need

%

% % Mandatory structure fields

%  globalStimParams.FS=100000;

%  globalStimParams.overallDuration=.1;  % s

% doPlot=1;

%

% stim.type='tone';

% stim.phases='sin';

% stim.toneDuration=.05;;

% stim.frequencies=500;

% stim.amplitudesdB=50;

% stim.beginSilence=.01;

% stim.endSilence=-1;

% stim.rampOnDur=.002;

% stim.rampOffDur=-1;

% [audio, msg]=stimulusCreate(globalStimParams, stim, doPlot);

%

% % or for multi-component stimuli

% % Mandatory structure fields

%  globalStimParams.FS=100000;

%  globalStimParams.overallDuration=.1;  % s

% ear=1;

% componentNo=1;

%

% stimComponents(ear, componentNo).type='tone';

% stimComponents(ear, componentNo).phases='sin';

% stimComponents(ear, componentNo).toneDuration=.05;;

% stimComponents(ear, componentNo).frequencies=500;

% stimComponents(ear, componentNo).amplitudesdB=50;

% stimComponents(ear, componentNo).beginSilence=.01;

% stimComponents(ear, componentNo).endSilence=-1;

% stimComponents(ear, componentNo).rampOnDur=.002;

% stimComponents(ear, componentNo).rampOffDur=-1;

%

% % All components are forced to have the same overall duration and sample rate

%

%

% [audio, msg]=stimulusCreate(globalStimParams, stimComponents);

%

%  Optional fields

%  .ears overides ear setting by component

%  globalStimParams.ears='monoticL'; % 'monoticL', 'monoticR', 'diotic', 'dichotic'

%

%  globalStimParams.filter = [leftfl leftfu rightfl right fu]

%    filter is applied separately to left and right combined sounds

%

%  correction factor is applied to all signals to compensate for differences in output devices.

% audioOutCorrection is a scalar

%  globalStimParams.audioOutCorrection=2;

%

%  globalStimParams.FFT= 1; % {0, 1}

%  globalStimParams.doPlot=1; % {0, 1}

%  globalStimParams.doPlay=1; % {0, 1}

%

%  stimComponents(ear, componentNo).filter=[100 10000 2] % [lower, upper, order] applies to an

% individual component

%

% Output arguments:

%  audio is a stereo signal, a 2-column vector

%

%

% stim.type='noise';    % {'IRN', 'irn', 'noise', 'pinkNoise'}

% stim.toneDuration=.05;;

% stim.amplitudesdB=50;

% stim.beginSilence=.01;

% stim.endSilence=-1;

% stim.rampOnDur=.002;

% stim.rampOffDur=-1;

% [audio, msg]=stimulusCreate(globalStimParams, stim);

%

% % for IRN only

% stim.type='IRN';

% stim.toneDuration=.05;;

% stim.amplitudesdB=50;

% stim.beginSilence=.01;

% stim.endSilence=-1;

% stim.rampOnDur=.002;

% stim.rampOffDur=-1;

% stim.niterations = 8;   %0 for white noise

% stim.delay = 1/150;

% stim.irnGain = 1;

% [audio, msg]=stimulusCreate(globalStimParams, stim);stimComponents.clickWidth;

%

% stim.type='clicks';    % uses frequencies for duty cycle

% stim.toneDuration=.05;;

% stim.frequencies=500;

% stim.amplitudesdB=50;

% stim.beginSilence=.01;

% stim.endSilence=-1;

% stim.rampOnDur=.002;

% stim.rampOffDur=-1;

% [audio, msg]=stimulusCreate(globalStimParams, stim);

% stimComponents.clickWidth=.0001;    % default= dt

% stimComponents.clickHeight= 20e-6;  % default= 28e-6 * 10^(stimComponents.amplitudesdB/20);

% stim.clickTimes=[.4 .6];            % times in cylce when clicks occur, default = 1

%

%

msg=''; %error messages; no message is a good message

% plotting can be set as a separate argument or as a globalstimParams

% variable. this is for backwards compatibility only

if nargin>2,

    globalStimParams.doPlot=doPlot;

end

% stimComponents(1,1).endSilence=-1;  % end silence is always computed

% perform checks and set defaults

[globalStimParams, stimComponents]=checkDescriptors(globalStimParams, stimComponents);

% create empty stereo audio of appropriate length

audio=zeros(globalStimParams.nSignalPoints, 2);

% signal=zeros(globalStimParams.nSignalPoints, 1);

dt=globalStimParams.dt;

[Nears nComponentSounds]=size(stimComponents);

for ear=1:Nears % left/ right

    % combinedSignal is the sum of all sounds in one ear

    % it is a column vector

    combinedSignal=zeros(globalStimParams.nSignalPoints,1);

    % find valid components

    % if a component is empty, it is not a validComponent and is ignored

    validStimComponents=[];

    for i=1:nComponentSounds

        if ~isempty(stimComponents(ear,i).type)

            validStimComponents=[validStimComponents i];

        end

    end

    for componentNo=validStimComponents

        % compute individual components before adding

        stim=stimComponents(ear,componentNo);

        switch stim.type

            case 'tone'

                stimulus=maketone(globalStimParams, stim);

            case 'fmTone'

                stimulus=makeFMtone(globalStimParams, stim);

            case 'OHIO'

                stim.beginSilence=0;

                stimulus=makeOHIOtones(globalStimParams, stim);

            case 'transposedStimulus'

                stim.beginSilence=0;    % necessary because of recursion

                stimulus=makeTransposedStimulus(globalStimParams, stim);

            case { 'noise', 'pinkNoise'}

                stimulus=makeNoise(globalStimParams, stim);

            case { 'whiteNoise'}

                stimulus=makeWhiteNoise(globalStimParams, stim);

            case {'IRN', 'irn'}

                stimulus=makeIRN(globalStimParams, stim);

            case {'RPN'}

                stimulus=makeRPN(globalStimParams, stim);

            case 'clicks'

                stimulus=makeClicks(globalStimParams, stim);

            case 'PressnitzerClicks'

                % kxx clicks

                % k is 1/clickRepeatFrequency

                stimulus=makePressnitzerClicks(globalStimParams, stimComponents);

            case 'PressnitzerABxClicks'

                % kxx clicks

                % k is 1/clickRepeatFrequency

                stimulus=makePressnitzerABxClicks(globalStimParams, stimComponents);

            case 'ABxClicks'

                % A=rand*k,  B=k-A, x=rand*k.

                stimulus=makeABxClicks(globalStimParams, stimComponents);

            case 'YostClicks'

                % kxx clicks

                % k is 1/clickRepeatFrequency

                stimulus=makeYostClicks(globalStimParams, stimComponents);

            case 'kxxClicks'

                % kxx clicks

                % k is 1/clickRepeatFrequency

                stimulus=makeKxxClicks(globalStimParams, stimComponents);

%             case 'babble'

%                 % NB random start in a long file

%                 [stimulus sampleRate]= wavread('babble');

%                 nPoints=round(sampleRate*...

%                     stimComponents(ear,componentNo).toneDuration);

%                 start=round(rand(1,1)*(length(stimulus)-nPoints));

%                 stimulus=stimulus(start:start+nPoints-1);

%                 rms= 20*log10((mean(stimulus.^2).^0.5)/20e-6);

%                 dBSPLrms=stimComponents(ear,componentNo).amplitudesdB;

%                 gain=10.^((dBSPLrms-rms)/20);

%                 stimulus=stimulus'*gain;

            case 'speech'

                [stimulus sampleRate]= wavread('speech');

                stimulus=stimulus';

                nPoints=sampleRate*stimComponents(ear,componentNo).toneDuration;

                if nPoints > length(stimulus)

                    initialSilence=zeros(1,nPoints-length(stimulus));

                else

                    initialSilence=[];

                    start=round(rand(1,1)*(length(stimulus)-nPoints));

                    stimulus=stimulus(start:start+nPoints-1);

                end

                rms= 20*log10((mean(stimulus.^2).^0.5)/20e-6);

                dBSPLrms=stimComponents(ear,componentNo).amplitudesdB;

                gain=10.^((dBSPLrms-rms)/20);

                stimulus=stimulus*gain;

                stimulus=[stimulus initialSilence ];

            case 'file'

                % component already read from file and stored in stimulus. Insert it here

                % additional code for establishing signal rms level

                % NB signal is always mono at this stage

                stimulus=stim.stimulus;

                dBSPL=stim.amplitudesdB;

                nPoints=round(stim.toneDuration/dt);

                [r c]=size(stimulus);

                if r>c, stimulus=stimulus'; end    % secure horizontal vector

                stimulus=stimulus(1,1:nPoints); % only mono taken from file

                try

                    % dB rms

                    rms= 20*log10((mean(stimulus.^2).^0.5)/20e-6);

                    % special request to set fixed rms for stimulus

                    dBSPLrms=stimComponents(ear,componentNo).amplitudesdBrms;

                    if ~(dBSPLrms==-1)

                        gain=10.^((dBSPLrms-rms)/20);

                        stimulus=stimulus*gain;

                    end

                catch

                    % If no request for rms level is made

                    % set dB as peak amp

                    [stimulus gain]= normalize(stimulus);

                    dBSPL=stimComponents(ear,componentNo).amplitudesdB;

                    if ~(dBSPL==-1)

                        amplitude=28e-6*10.^(dBSPL/20);

                        stimulus=stimulus*amplitude;

                    end

                end

            case 'none'

                % no stimulus

                stimulus=zeros(1,round(stim.toneDuration/dt));

            case 'digitStrings'

                % select a digit string at random anduse as target

                files=dir(['..' filesep '..' filesep 'multithresholdResources\digitStrings']);

                files=files(3:end);

                nFiles=length(files);

                fileNo=ceil(nFiles*rand);

                fileData=files(fileNo);

                fileName=['..\..\multithresholdResources\digitStrings\' fileData.name];

                [stimulus sampleRate]=wavread(fileName);

                stimulus=stimulus(:,1)';  % make into a row vector

                % estimate the extend of endsilence padding

                nPoints=sampleRate*...

                    stimComponents(ear,componentNo).toneDuration;

                if nPoints > length(stimulus)

                    endSilence=zeros(1,nPoints-length(stimulus));

                else

                    % or truncate the file

                    endSilence=[];

                    stimulus=stimulus(1:nPoints);

                end

                % compute rms before adding silence

                rms= 20*log10((mean(stimulus.^2).^0.5)/20e-6);

                dBSPLrms=stimComponents(ear,componentNo).amplitudesdB;

                gain=10.^((dBSPLrms-rms)/20);

                stimulus=stimulus*gain;

                stimulus=[stimulus endSilence];

                global stimulusParameters

                stimulusParameters.digitString=fileName(end-7:end-5);

            otherwise

                switch stim.type(end-5:end)

                    % any file name with 'Babble' at the end is a

                    % multiThreshold file

                    case 'Babble'

                        % one of the many babbles is selected.

                        % NB random start in a long file

                        %  stim.type should contain the name of the babble file

                        fileName=['..' filesep '..' filesep ...

                            'multithresholdResources' filesep ...

                            'backgrounds and maskers'...

                             filesep stim.type];

                        [stimulus sampleRate]= wavread(fileName);

                        if ~isequal(sampleRate, globalStimParams.FS)

                            % NB the file will read but will disagree with

                            % tone stimuli or other files read

                            msg= ['error: file sample rate disagrees ' ...

                                'with sample rate requested in paradigm'...

                                ' file (' ...

                            num2str(globalStimParams.FS) ').'];

                            error(msg);

                        end

                        nPoints=round(sampleRate*...

                            stimComponents(ear,componentNo).toneDuration);

                        start=round(rand(1,1)*(length(stimulus)-nPoints));

                        stimulus=stimulus(start:start+nPoints-1);

                        rms= 20*log10((mean(stimulus.^2).^0.5)/20e-6);

                        dBSPLrms=stimComponents(ear,componentNo).amplitudesdB;

                        gain=10.^((dBSPLrms-rms)/20);

                        stimulus=stimulus'*gain;

                    otherwise

                        % stim.type may be the name of a file to be read

                        % play from beginning for stimulus duration

                        try

                            [stimulus sampleRate]= wavread([stim.type '.wav']);

                        catch

                            error(['stimulusCreate: unrecognised stimulus type -> '...

                                stim.type])

                        end

                        if ~isequal(sampleRate, globalStimParams.FS)

                            % NB the file will read but will disagree with

                            % tone stimuli or other files read

                            msg= ['error: file sample rate disagrees ' ...

                                'with sample rate requested in paradigm'...

                                ' file (' ...

                            num2str(globalStimParams.FS) ').'];

                            error(msg);

                        end

                        stimulus=stimulus';  % make into a row vector

                        % estimate the extend of endsilence padding

                        nPoints=sampleRate*...

                            stimComponents(ear,componentNo).toneDuration;

                        if nPoints > length(stimulus)

                            endSilence=zeros(1,nPoints-length(stimulus));

                        else

                            % or truncate the file

                            endSilence=[];

                            stimulus=stimulus(1:nPoints);

                        end

                        % compute rms before adding silence

                        rms= 20*log10((mean(stimulus.^2).^0.5)/20e-6);

                        dBSPLrms=stimComponents(ear,componentNo).amplitudesdB;

                        gain=10.^((dBSPLrms-rms)/20);

                        stimulus=stimulus*gain;

                        stimulus=[stimulus endSilence];

                end

        end

        % NB stimulus is a row vector now!

        % audio and combinedSignal were column vectors

        % signal will be a row vector

        % filter stimulus

        try

            % if filter field is present, [lower upper order]

            if stim.filter(1)>0	% 0 means don't filter

                stimulus=Butterworth (stimulus, dt, stim.filter(1), ...

                    stim.filter(2), stim.filter(3));

            end

        catch

        end

        % apply amplitude modulation

        if isfield(stim,'AMfrequency') & isfield(stim,'AMdepth')

            if stim.AMfrequency>0 & stim.AMdepth>0

                time=dt:dt:dt*length(stimulus);

                modulator=sin(2*pi*stim.AMfrequency*time);

                modulator=modulator*stim.AMdepth/100 + 1; % depth is percent

                stimulus=stimulus.*modulator/2;

            end

        end

        % Basic stimulus is now created.

        % Add trappings, ramps, silences to main stimulus

        rampOnTime=0;		 %ramp starts at the beginning of the stimulus

        rampOffTime=stim.toneDuration-stim.rampOffDur;

        if stim.rampOnDur>0.0001

            stimulus=applyRampOn(stimulus, stim.rampOnDur, rampOnTime, 1/dt);

            stimulus=applyRampOff(stimulus, stim.rampOffDur, rampOffTime, 1/dt);

        end

        if stim.rampOnDur<-0.0001   % apply Gaussian ramp

            stimulus=applyGaussianRamps(stimulus, -stim.rampOnDur, 1/dt);

        end

        % Initial silence

        % start with a signal of the right length consisting of zeros

        signal=zeros(1, globalStimParams.nSignalPoints);

        % identify start of stimulus

        insertStimulusAt=round(stim.beginSilence/dt)+1;

        % add stimulus

        endOfStimulus=insertStimulusAt+length(stimulus)-1;

        if endOfStimulus<=globalStimParams.nSignalPoints

            signal(1, insertStimulusAt: endOfStimulus)=stimulus;

        else

            error('stimulusCreate: tone too long to fit into the overall duration')

        end

        % time signature

        time=dt:dt:dt*length(signal);

        %     figure(22), plot(signal), title([num2str(ear) ' - ' num2str(componentNo)]),pause (1)

        try

            % create a column vector and trap if no signal has been created

            signal=signal';

            % also traps if signals are not the same length

            combinedSignal=combinedSignal+signal;

            %             figure(21), plot(combinedSignal), title([num2str(ear) ' - ' num2str(componentNo)]),pause (1)

        catch

            % dump everything because there is a problem

            globalStimParams

            [ear  componentNo]

            stim

            size(combinedSignal)

            size(signal)

            [   size(initialSilence)            size(signal)            size(endSilence)]

            [ ear  componentNo...

                round(globalStimParams.overallDuration*globalStimParams.FS)...

                round(stim.beginSilence*globalStimParams.FS)...

                round(stim.toneDuration*globalStimParams.FS) ...

                round(stim.endSilence*globalStimParams.FS)...

                (                       round(globalStimParams.overallDuration*globalStimParams.FS)...

                -round(stim.beginSilence*globalStimParams.FS)...

                -round(stim.toneDuration*globalStimParams.FS) ...

                -round(stim.endSilence*globalStimParams.FS))...

                ]

            error(' trouble in stimulusCreate: signals are the wrong length ')

        end

    end % component no

    audio(:,ear)= combinedSignal;

    % FFT

    try

        if globalStimParams.FFT

            showFFT (audio, dt)

        end

    catch

    end

end % ear

switch globalStimParams.ears

    % normally the signals are created in appropriate ears but .ears can

    % overide this to produce a mono signal.

    case 'monoticL';

        % combine left and right ears to make a mono signal in the left ear

        audio(:,1)=audio(:,1)+audio(:,2);

        audio(:,2)=zeros(globalStimParams.nSignalPoints, 1);

    case 'monoticR';

        % combine left and right ears to make a mono signal in the right ear

        audio(:,2)=audio(:,1)+audio(:,2);

        audio(:,1)=zeros(globalStimParams.nSignalPoints, 1);

    case 'diotic';

        % combine left and right ears to make a mono signal in both ears

        bothEarsCombined=audio(:,1)+audio(:,2);

        audio(:,2)=bothEarsCombined;

        audio(:,1)=bothEarsCombined;

    otherwise

        % Any other denomination produces no effect here

end

% Plotting as required

if globalStimParams.doPlot

    figure(9), clf

    plot(time,audio(:,1)'), hold on,

    if Nears>1

        offSet=(max(audio(:,1))+max(audio(:,2)))/10;

        offSet=2*offSet+max(audio(:,1))+max(audio(:,2));

        plot(time,audio(:,2)'+offSet,'r'), hold off

    end

    ylabel('left                             right')

end

% Calibration

% all signals are divided by this correction factor

% peakAmp=globalStimParams.audioOutCorrection; % microPascals = 100 dB SPL

audio=audio/globalStimParams.audioOutCorrection;

if globalStimParams.doPlay

        if ispc

            wavplay(audio,globalStimParams.FS)

        else

            sound(audio,globalStimParams.FS)

        end

    wavplay(audio,globalStimParams.FS)

end

% all Done

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%---------------------------------------------------- maketone

function tone=maketone(globalStimParams, stimComponents)

% maketone generates a stimComponents tone

% tone=maketone(dt, frequencies, toneDuration, dBSPL, phases)

% tone is returned in Pascals

% frequencies is a list of frequencies

% phase is list of phases or 'sin', 'cos', 'alt'

%

% defaults:

%  phase = sin

%  dBSPL=56 dB SPL

dt=globalStimParams.dt;

frequencies=stimComponents.frequencies;

toneDuration=stimComponents.toneDuration;

dBSPL=stimComponents.amplitudesdB;

phases=stimComponents.phases;

if ischar(phases)

    switch phases

        case 'sin'

            phases= zeros(1,length(frequencies));

        case 'cos'

            phases= pi/2*ones(1,length(frequencies));

        case 'alt'

            phases= repmat([0 pi/2], 1, floor(length(frequencies)/2));

            if length(phases)<length(frequencies)

                phases=[phases 0];

            end

        case {'ran', 'rand'}

            phases= 2*pi*rand(1,length(frequencies));

    end

end

if length(phases)==1, phases=repmat(phases(1), 1, length(frequencies)); end

if length(phases)<length(frequencies)

    error('makeTone:phase specification must have the same length as frequencies')

end

if length(dBSPL)==1, dBSPL=repmat(dBSPL(1), 1, length(frequencies)); end

if length(dBSPL)<length(dBSPL)

    error('makeTone:dBSPL specification must have the same length as frequencies')

end

time=dt:dt:toneDuration;

amplitudes=28e-6* 10.^(dBSPL/20);

tone=zeros(size(time));

for i=1:length(frequencies)

    frequency=frequencies(i);

    phase=phases(i);

    tone=tone+amplitudes(i)*sin(2*pi*frequency*time+phase);

end

% figure(1), plot(time, signal)

%---------------------------------------------------- makeOHIOtones

function stimulus=makeOHIOtones(globalStimParams, stimComponents)

% Generates a stimulus consisting of one or more 10-ms tones

% Tones are presented at 10-ms intervals

% Each tone has its own amplitude and its own ramp.

frequencies=stimComponents.frequencies;

amplitudesdB=stimComponents.amplitudesdB;

nFrequencies=length(frequencies);

dt=globalStimParams.dt;

toneDuration=.010;

time=dt:dt:toneDuration;

rampDuration=stimComponents.rampOnDur;

rampTime=dt:dt:rampDuration;

ramp=[0.5*(1+cos(2*pi*rampTime/(2*rampDuration)+pi)) ones(1,length(time)-length(rampTime))];

ramp=ramp.*fliplr(ramp);

stimulus= zeros(1,round((toneDuration+globalStimParams.beginSilences(end))/dt)+1);

for i=1:nFrequencies

    toneBeginPTR=round(globalStimParams.beginSilences(i)/dt)+1;

    frequency=frequencies(i);

    dBSPL=amplitudesdB(i);

    amplitude=28e-6* 10.^(dBSPL/20);

    tone=amplitude*sin(2*pi*frequency*time);

    tone=tone.*ramp;

    stimulus(toneBeginPTR:toneBeginPTR+length(tone)-1)=stimulus(toneBeginPTR:toneBeginPTR+length(tone)-1)+tone;

end

% figure(2), plot( stimulus')

%---------------------------------------------------- makeFMtone

function tone=makeFMtone(globalStimParams, stimComponents)

% maketone generates a stimComponents tone

% tone=maketone(dt, frequencies, toneDuration, dBSPL, phases)

% tone is returned in Pascals

% frequencies is a list of frequencies

% phase is list of phases or 'sin', 'cos', 'alt'

%

% defaults:

%  phase = sin

%  dBSPL=56 dB SPL

dt=globalStimParams.dt;

frequencies=stimComponents.frequencies;

toneDuration=stimComponents.toneDuration;

dBSPL=stimComponents.amplitudesdB;

phases=stimComponents.phases;

fmDepth=stimComponents.fmDepth;

fmFrequency=stimComponents.fmFrequency;

if ischar(phases)

    switch phases

        case 'sin'

            phases= zeros(1,length(frequencies));

        case 'cos'

            phases= pi/2*ones(1,length(frequencies));

        case 'alt'

            phases= repmat([0 pi/2], 1, floor(length(frequencies)/2));

            if length(phases)<length(frequencies)

                phases=[phases 0];

            end

    end

end

if length(phases)==1, phases=repmat(phases(1), 1, length(frequencies)); end

if length(phases)<length(frequencies)

    error('makeTone:phase specification must have the same length as frequencies')

end

if length(dBSPL)==1, dBSPL=repmat(dBSPL(1), 1, length(frequencies)); end

if length(dBSPL)<length(dBSPL)

    error('makeTone:dBSPL specification must have the same length as frequencies')

end

time=dt:dt:toneDuration;

amplitudes=28e-6* 10.^(dBSPL/20);

tone=zeros(size(time));

for i=1:length(frequencies)

    frequency=frequencies(i);

    phase=phases(i);

    tone=tone+amplitudes(i)*sin(2*pi*frequency*time+phase + fmDepth*sin(2*pi*fmFrequency*time));

end

% figure(1), plot(time, signal)

%----------------------------------------------------makeTransposedStimulus

function [transposedStimulus]=makeTransposedStimulus(globalStimParams, stim)

% globalStimParams.FS=100000;

% globalStimParams.overallDuration=.1;  % s

% globalStimParams.doPlot=1;

%

% stim.type='transposedStimulus';

% stim.phases='sin';

% stim.toneDuration=.05;;

% stim.frequencies=500;

% stim.amplitudesdB=50;

% stim.beginSilence=.01;

% stim.endSilence=-1;

% stim.rampOnDur=.002;

% stim.rampOffDur=-1;

%

% stim.transCarrierFreq=4000;

% stim.transModFreq=100;

transposedStimulus=[];

% make envelope of transposed tone

for i=1:length(stim.transModFreq)

    stim.type='tone';

    stim.frequencies=stim.transModFreq(i);

    stim.endsilence=-1;     stim.beginSilence=0;

    [envelope, msg]=stimulusCreate(globalStimParams, stim);    % NB recursive

    envelope=envelope(:,1);   % mono

    % HW rectify

    envelope(find(envelope<0))=0;

    % LP filter

    maxEnvelope=max(envelope);

    envelope=UTIL_Butterworth (envelope, globalStimParams.dt, 10, .2*stim.transModFreq(i), 2);

    envelope=envelope*maxEnvelope/max(envelope);

    % make the carrier

    stim.frequencies=stim.transCarrierFreq(i);

    stim.endsilence=-1;    stim.beginSilence=0;

    [audio, msg]=stimulusCreate(globalStimParams, stim);

    carrier=audio(:,1);

    x= (carrier.*envelope)';

    % base amplitude on peak of unmodulated carrier

    x=x/max(carrier);

    transposedStimulus=[transposedStimulus; x];

end

transposedStimulus=sum(transposedStimulus,1);

% clf,plot(transposedStimulus)

%--------------------------------------------------------------------makeClicks

function clickTrain=makeClicks(globalStimParams, stimComponents)

% makeClicks(F0, clickTimes, duration, FS);

% F0 specifies the repetition rate of the click sequence

% If F0=-1, a single click is generated at the start of the duration of the signal

%

% clickTimes a are fractions of the period

%  and specify when the click appears in the period

% A click time of 1 is reset to zero.

% if the clicktime plus the click width is greater than the period, no click is generated

% clicks are treated as 20 microPascal high before amplification

%  unless otherwise specified in stimComponents.clickHeight

% click width is dt unless specified in stimComponents.clickWidth

%

% for regular pulse train set clicktimes=1 or zero;

% FS is the sampling interval;

% CarlyonClickTrain(100, [.4 1], 40, 441000);

FS=globalStimParams.FS; % sample rate

dt=1/FS;

try,clickWidth=stimComponents.clickWidth;catch, clickWidth=dt; end

try,clickHeight=stimComponents.clickHeight; catch, clickHeight=28e-6 * 10^(stimComponents.amplitudesdB/20); end

try, clickTimes=stimComponents.clickTimes; catch, clickTimes=1; end

% clickTimes are the times in a cycle that the click

% occurs

checkClickTimes=clickTimes-1;

if max(checkClickTimes)>1

    msg= 'click times must be <= than 1 (period)';

    return

end

if clickWidth>=1/stimComponents.clickRepeatFrequency

    msg= 'click width is too great for frequency';

    return

end

duration=stimComponents.toneDuration;

totalLength=round(stimComponents.toneDuration/dt);

F0=stimComponents.clickRepeatFrequency;

F0=round(F0/dt)*dt;

if F0==-1 % single click required

    F0=1/duration;

    repetitions=1;

    clickStartTimes=1; %clicktimes are fractions of a period

else

    repetitions=round(F0*duration)-1;

    duration=repetitions/F0;

    clickStartTimes=clickTimes;

end

% if a clickTime=1 (end of duty period) set it to the beginning

clickStartTimes(clickStartTimes==1)=0;

period=1/F0;

time=dt:dt:period;

nPoints=length(time);

signal=zeros(1,nPoints);

dBSPL=stimComponents.amplitudesdB;

% compute click train for a single cycle

clickWidthNpoints=round(clickWidth*FS);

for i=1:length(clickStartTimes)

    %     if clickStartTimes(i)<clickWidth

    %         clickStartTimes(i)=dt;

    %     end

    clickTime=round(period*clickStartTimes(i)/dt -dt);

    % clicks are treated as 20 microPascal high

    if clickTime+clickWidthNpoints<length(signal)

        signal(clickTime+1:clickTime+clickWidthNpoints)=clickHeight;

    end

end

clickTrain=repmat(signal, 1, repetitions);

if length(clickTrain)>totalLength

    clickTrain=clickTrain(1:totalLength);

elseif length(clickTrain)<totalLength

    timeToAdd=zeros(1,round((totalLength-length(clickTrain))));

    clickTrain=[clickTrain timeToAdd];

    % figure (1), plot(clickTrain)

end

%----------------------------------------------------------------makePressnitzerClicks

function signal=makePressnitzerClicks(globalStimParams, stimComponents)

% PressnitzerClicks(k,duration,dt)

% Generates a sequence of clicks with intervals kxx

%  where x= rand*k/2

% This is not the same as Kaernbach and Demany clicks

FS=globalStimParams.FS; % sample rate

dt=1/FS;

% obligatory parameter

try

    k=stimComponents.k;

catch

    error('PressnitzerClicks: field ''k'' is missing from stimcomponent')

end

% optional parameters

if isfield(stimComponents,'clickWidth')

    clickWidth=stimComponents.clickWidth;

else

    clickWidth=dt;

end

clickWidthNpoints=round(clickWidth*FS);

if isfield(stimComponents,'clickHeight')

    clickHeight=stimComponents.clickHeight;

else

    clickHeight=28e-6 * 10^(stimComponents.amplitudesdB/20);

end

duration=stimComponents.toneDuration;

signalLength=round(duration/dt);

signal=zeros(1,signalLength);

kInterval=round(k/dt);

halfk=k/2;

signal(1)=clickHeight;

timeIdx=0;

while timeIdx<signalLength

    % first interval = k

    clickTime=timeIdx+kInterval;

    signal(clickTime:clickTime+clickWidthNpoints)=clickHeight;

    timeIdx=timeIdx+kInterval;

    % second interval = 0 : k/2

    intervalx1=round(rand*halfk/dt);

    clickTime=timeIdx+intervalx1;

    signal(clickTime:clickTime+clickWidthNpoints)=clickHeight;

    timeIdx=timeIdx+intervalx1;

    % third interval = 0 : k/2

    intervalx1=round(rand*halfk/dt);

    clickTime=timeIdx+intervalx1;

    signal(clickTime:clickTime+clickWidthNpoints)=clickHeight;

    timeIdx=timeIdx+intervalx1;

end

signal=signal(1:signalLength);

% figure(1),	plot(dt:dt:duration,signal)

%----------------------------------------------------------------makePressnitzerABXClicks

function signal=makePressnitzerABxClicks(globalStimParams, stimComponents)

% Generates a sequence of clicks with intervals ABx

% AB interval is 2*k

% where A= rand* k

%       B= k-A

%       x= k/2

% These are second order clicks

FS=globalStimParams.FS; % sample rate

dt=1/FS;

% obligatory parameter

try

    k=stimComponents.k;

catch

    error('PressnitzerClicks: field ''k'' is missing from stimcomponent')

end

% optional parameters

if isfield(stimComponents,'clickWidth')

    clickWidth=stimComponents.clickWidth;

else

    clickWidth=dt;

end

clickWidthNpoints=round(clickWidth*FS);

if isfield(stimComponents,'clickHeight')

    clickHeight=stimComponents.clickHeight;

else

    clickHeight=28e-6 * 10^(stimComponents.amplitudesdB/20);

end

duration=stimComponents.toneDuration;

signalLength=round(duration/dt);

signal=zeros(1,2*signalLength); % allow for overrun

ABinterval=k/dt;                % i.e. the number of dt steps

randomInterval=ABinterval/2;

signal(1)=clickHeight;

time=0;

while time<signalLength

    % first interval = A

    intervalA=rand*ABinterval;

    clickTime=round(time+intervalA)+1;   % +1 avoids zero index

    signal(clickTime:clickTime+clickWidthNpoints)=clickHeight;

    time=time+intervalA;

    % second interval = B

    intervalB=ABinterval-intervalA;

    clickTime=round(time+intervalB)+1;

    signal(clickTime:clickTime+clickWidthNpoints-1)=clickHeight;

    time=time+intervalB;

    % third interval = 0 : k/2

    intervalx1=rand*randomInterval;    % mean random interval=k

    clickTime=round(time+intervalx1)+1;

    signal(clickTime:clickTime+clickWidthNpoints-1)=clickHeight;

    time=time+intervalx1;

end

signal=signal(1:signalLength);

% figure(1),	plot(dt:dt:duration,signal)

%-----------------------------------------------------makeABxClicks

function signal=makeABxClicks(globalStimParams, stimComponents)

% Generates a sequence of clicks with intervals ABx

% AB interval is 2*k

% where A= rand* k

%       B= k-A

%       x= rand*2*k

% These are second order clicks

FS=globalStimParams.FS; % sample rate

dt=1/FS;

% obligatory parameter

try

    k=stimComponents.k;

catch

    error('PressnitzerClicks: field ''k'' is missing from stimcomponent')

end

% optional parameters

if isfield(stimComponents,'clickWidth')

    clickWidth=stimComponents.clickWidth;

else

    clickWidth=dt;

end

clickWidthNpoints=round(clickWidth*FS);

if isfield(stimComponents,'clickHeight')

    clickHeight=stimComponents.clickHeight;

else

    clickHeight=28e-6 * 10^(stimComponents.amplitudesdB/20);

end

duration=stimComponents.toneDuration;

signalLength=round(duration/dt);

signal=zeros(1,2*signalLength); % allow for overrun

ABinterval=2*k/dt;

randomInterval=ABinterval;

signal(1)=clickHeight;

timeIdx=0;

while timeIdx<signalLength

    % first interval = A

    intervalA=round(rand*ABinterval);

    clickTime=timeIdx+intervalA+1;

    signal(clickTime:clickTime+clickWidthNpoints-1)=clickHeight;

    timeIdx=timeIdx+intervalA;

    % second interval = B

    intervalB=round(ABinterval-intervalA);

    clickTime=timeIdx+intervalB;

    signal(clickTime:clickTime+clickWidthNpoints-1)=clickHeight;

    timeIdx=timeIdx+intervalB;

    % third interval = 0 : k

    intervalx1=round(rand*randomInterval);    % mean random interval=k

    clickTime=timeIdx+intervalx1;

    signal(clickTime:clickTime+clickWidthNpoints-1)=clickHeight;

    timeIdx=timeIdx+intervalx1;

end

signal=signal(1:signalLength);

% figure(1),	plot(dt:dt:duration,signal)

%----------------------------------------------------------------makeYostClicks

function signal=makeYostClicks(globalStimParams, stimComponents)

% Generates a shuffled sequence of clicks with intervals kxxxx

%  where max(x)= 2*k

%  and there are n occurrences of x

% this section requires:

%  stimComponents.k

%  stimComponents.nXs

%  stimComponents.toneDuration

% optional:

%  stimComponents.clickWidth	%useful because width depends on dt

%  stimComponents.clickHeight	%best left to amplitude rescaling later

FS=globalStimParams.FS; % sample rate

dt=1/FS;

% obligatory parameters

try

    k=stimComponents.k;

catch

    error('makeYostClicks: field ''k'' is missing from stimComponents')

end

try

    nXs=stimComponents.nXs;

catch

    error('makeYostClicks: field ''nXs'' is missing from stimComponents')

end

try

    shuffled=stimComponents.shuffled;

catch

    error('makeYostClicks: field ''shuffled'' is missing from stimComponents')

end

try

    duration=stimComponents.toneDuration;

catch

    error('makeYostClicks: field ''toneDuration'' is missing from stimComponents')

end

% optional parameters

if isfield(stimComponents,'clickWidth')

    clickWidth=stimComponents.clickWidth;

else

    clickWidth=dt;

end

clickWidthNpoints=round(clickWidth*FS);

if isfield(stimComponents,'clickHeight')

    clickHeight=stimComponents.clickHeight;

else

    clickHeight=28e-6 * 10^(stimComponents.amplitudesdB/20);

end

kInterval=round(k/dt);

twoK=k*2;				% max width of x interval

signalLength=round(duration/dt);

signal=zeros(1,signalLength);

timeIdx=0;

intervalCount=0;

while timeIdx<signalLength

    timeIdx=timeIdx+kInterval;

    if timeIdx>signalLength, break,end

    intervalCount=intervalCount+1;

    intervals(intervalCount)=kInterval;

    % repeat x intervals as required

    if nXs>0

        for nX=1:nXs

            xInterval=round(rand*twoK/dt);

            timeIdx=timeIdx+xInterval;

            if timeIdx>signalLength, break,end

            intervalCount=intervalCount+1;

            intervals(intervalCount)=xInterval;

        end

    end

    if timeIdx>signalLength, break,end

end

% shuffle intervals

if shuffled

    randomNumbers=rand(1,length(intervals));

    [randomNumbers idx]=sort(randomNumbers);

    intervals=intervals(idx);

    idx=intervals>0;

    intervals=intervals(idx);

end

intervalCount=length(intervals);

signal(1)=clickHeight;

clickTime=0;

for i=1:intervalCount

    clickTime=clickTime+intervals(i);

    signal(clickTime:clickTime+clickWidthNpoints)=clickHeight;

end

signal=signal(1:signalLength);

%  figure(1),	plot(dt:dt:duration,signal)

%--------------------------------------------------------------------makeKxxClicks

function signal=makeKxxClicks(globalStimParams, stimComponents)

% Click train consists of kkxxx.. sequences

% k is the duration of a fixed interval (seconds)

% random intervals are distributed 0 : 2* k (NB not like Pressnitzer clicks)

% nKs is the number of successive 'k' intervals

% nXs is the number of random intervals between k sequences

% sequences of 3 x intervals > k are replaced with new sequences

% shuffled causes all intervals to be reordered randomly

% NB 1/k is the mean click rate

FS=globalStimParams.FS; % sample rate

dt=1/FS;

try

    k=stimComponents.k;		% duration (s) of fixed intervals

catch

    error('makeYostClicks: field ''k'' is missing from stimComponents')

end

try

    duration=stimComponents.toneDuration;

catch

    error('makeYostClicks: field ''duration'' is missing from stimComponents')

end

if isfield(stimComponents,'clickWidth')

    clickWidth=stimComponents.clickWidth;

else

    clickWidth=dt;

end

if isfield(stimComponents,'clickHeight')

    clickHeight=stimComponents.clickHeight;

else

    clickHeight=28e-6 * 10^(stimComponents.amplitudesdB/20);

end

if isfield(stimComponents,'order')

    order=stimComponents.order;

else

    order=1;

end

if isfield(stimComponents,'nKs')

    nKs=stimComponents.nKs;