Mercurial > hg > map
view utilities/stimulusCreate.m @ 38:c2204b18f4a2 tip
End nov big change
author | Ray Meddis <rmeddis@essex.ac.uk> |
---|---|
date | Mon, 28 Nov 2011 13:34:28 +0000 |
parents | 25d53244d5c8 |
children |
line wrap: on
line source
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 % begin 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 % The length of the list of frequencies determines the numberof tones % Tones are either presented at 10-ms intervals or simultaneously % all tones are individually ramped % Each tone has its own amplitude and its own ramp. frequencies=stimComponents.frequencies; amplitudesdB=stimComponents.amplitudesdB; nFrequencies=length(frequencies); if amplitudesdB==-100 % catch trial amplitudesdB=repmat(-100,1,nFrequencies); end dt=globalStimParams.dt; toneDuration=.010; time=dt:dt:toneDuration; % fixed ramp, silenceDuration, toneDuration rampDuration=0.005; 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); silenceDuration=0.010; silenceDurationLength=round(silenceDuration/dt); initialSilence=zeros(1,silenceDurationLength); silenceToneDuration=toneDuration + silenceDuration; silenceToneDurationLength=round(silenceToneDuration/dt); % OHIO spect produces simultaneous tones if strcmp(stimComponents.OHIOtype,'OHIOspect') totalDuration=silenceToneDuration; else totalDuration=silenceToneDuration*nFrequencies; end totalDurationLength=round(totalDuration/dt); stimulus=zeros(1,totalDurationLength); toneBeginPTR=1; for i=1:nFrequencies frequency=frequencies(i); dBSPL=amplitudesdB(i); amplitude=28e-6* 10.^(dBSPL/20); tone=amplitude*sin(2*pi*frequency*time); tone=tone.*ramp; % stimulus is normally zeros except for OHIOspect stimulus(toneBeginPTR:toneBeginPTR+silenceToneDurationLength-1)=... [initialSilence tone]+... stimulus(toneBeginPTR:toneBeginPTR+silenceToneDurationLength-1); if ~strcmp(stimComponents.OHIOtype,'OHIOspect') toneBeginPTR=toneBeginPTR+silenceToneDurationLength; end 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; else nKs=1; end if isfield(stimComponents,'nXs') nXs=stimComponents.nXs; else nXs=1; end if isfield(stimComponents,'shuffled') shuffled=stimComponents.shuffled; else shuffled=1; end kLength=round(k/dt); % fixed interval xLength=2*kLength; % maximum random interval requiredSignalLength=round(duration/dt); intervalsPerCycle=(nKs+nXs); cycleLength=nKs*kLength+nXs*xLength; % more cycles to allow for uncertainty nCycles=5*round(requiredSignalLength/cycleLength); nIntervals=nCycles*intervalsPerCycle; % random intervals if nXs>0 xIntervals=floor(rand(1,nIntervals)*2*kLength); % weed out triple intervals > 2*k rogues=1; while sum(rogues) y=(xIntervals>kLength); rogues=(sum([y(1:end-2)' y(2:end-1)' y(3:end)'],2)>2); xIntervals(rogues)=floor(rand*2*kLength); end xIntervals=reshape(xIntervals,nCycles,intervalsPerCycle); else xIntervals=[]; end % insert constrained (k) intervals if nKs>0 switch order case 1 kIntervals=floor(ones(nCycles,nKs)*kLength); case 2 nKs=1; % force this kIntervals=floor(rand(nCycles,1)*kLength); kIntervals=[kIntervals kLength-kIntervals]; end else kIntervals=[]; end % combine fixed and random intervals=[kIntervals xIntervals(:,nKs+1:end)]; % make a single array; [r c]=size(intervals); intervals=reshape(intervals',1,r*c); % shuffle intervals if shuffled randomNumbers=rand(1,length(intervals)); [randomNumbers idx]=sort(randomNumbers); intervals=intervals(idx); idx=intervals>0; intervals=intervals(idx); end % convert intervals to clicks clickTimes=cumsum(intervals); signal(clickTimes)=clickHeight; signal=signal(1:requiredSignalLength); % figure(1), clf, plot(signal) %--------------------------------------------------------------------makeNoise function noise=makeNoise(globalStimParams, stimComponents) % FS in Hz, noiseDuration in s, delay in s; % noise is returned with overall level dB(rms) = amplitudesdB % % % You need % % stim.type='noise'; % or 'IRN', or 'pinkNoise' % stim.toneDuration=.05; % stim.amplitudesdB=50; % stim.beginSilence=.01; % stim.endSilence=-1; % stim.rampOnDur=.002; % stim.rampOffDur=-1; % % % Mandatory structure fields % globalStimParams.FS=100000; % globalStimParams.overallDuration=.1; % s % globalStimParams.doPlot=1; % globalStimParams.doPlay=1; % % [audio, msg]=stimulusCreate(globalStimParams, stim, ); % % % local % stim.type='noise'; % or 'IRN' % FS=globalStimParams.FS; noiseDuration= stimComponents.toneDuration; npts=round(noiseDuration*FS); noise=randn(1,npts); % NB randn (normally distributed) switch stimComponents.type case 'pinkNoise' % noise=UTIL_lowpassFilterFreq(noise, 100, 1/FS); noise=UTIL_bandPassFilter(noise, 1, 100, 200, 1/FS,[]); end rms=(mean(noise.^2)^.5); %should be 20 microPascals for 0 dB SPL adjust=20e-6/rms; noise=noise*adjust; rms=(mean(noise.^2)^.5); amplitude=10.^(stimComponents.amplitudesdB/20); noise=amplitude*noise; % rms=(mean(noise.^2)^.5); % dBnoise=20*log10(rms/20e-6) %--------------------------------------------------------------------makeWhiteNoise function noise=makeWhiteNoise(globalStimParams, stimComponents) % FS in Hz, noiseDuration in s, delay in s; % noise is bandpass filtered between 100 and 10000 Hz % spectrum level (dB/Hz) is 40 dB below nominal level. % noise is returned with dB(rms) = amplitudesdB % % % You need % % stim.type='noise'; % or 'IRN', or 'pinkNoise' % stim.toneDuration=.05; % stim.amplitudesdB=50; % stim.beginSilence=.01; % stim.endSilence=-1; % stim.rampOnDur=.002; % stim.rampOffDur=-1; % % % Mandatory structure fields % globalStimParams.FS=100000; % globalStimParams.overallDuration=.1; % s % globalStimParams.doPlot=1; % globalStimParams.doPlay=1; % % [audio, msg]=stimulusCreate(globalStimParams, stim, ); % % % local % stim.type='noise'; % or 'IRN' % FS=globalStimParams.FS; noiseDuration= stimComponents.toneDuration; npts=round(noiseDuration*FS); noise=randn(1,npts); noise=UTIL_bandPassFilter (noise, 6, 100, 10000, 1/FS, []); rms=(mean(noise.^2)^.5); %should be 20 microPascals for 0 dB SPL adjust=20e-6/rms; noise=noise*adjust; rms=(mean(noise.^2)^.5); amplitude=10.^(stimComponents.amplitudesdB/20); noise=amplitude*noise; % rms=(mean(noise.^2)^.5); % dBnoise=20*log10(rms/20e-6) %-----------------------------------------------------------------makeIRN function noise=makeIRN(globalStimParams, stimComponents) % FS in Hz, noiseDuration in s, delay in s; % noise is returned with dB(rms) = amplitudesdB % % % You need % % stim.type='noise'; % or 'IRN', or 'pinkNoise' % stim.toneDuration=.05; % stim.amplitudesdB=50; % stim.beginSilence=.01; % stim.endSilence=-1; % stim.rampOnDur=.002; % stim.rampOffDur=-1; % % % Mandatory structure fields % globalStimParams.FS=100000; % globalStimParams.overallDuration=.1; % s % globalStimParams.doPlot=1; % globalStimParams.doPlay=1; % % [audio, msg]=stimulusCreate(globalStimParams, stim, ); % % % local % stim.type='noise'; % or 'IRN' % % for IRN only % stim.niterations = 8; %0 for white noise % stim.delay = 1/150; % stim.irnGain = 1; % FS=globalStimParams.FS; noiseDuration= stimComponents.toneDuration; nIterations=stimComponents.niterations; if nIterations==0 % white noise is specified as nIterations=1 nIterations=1; IRNgain=0; delay=0.01; % dummy else % IRN delay=stimComponents.delay; IRNgain=stimComponents.irnGain; end npts=round(noiseDuration*FS); dels=round(delay*FS); noise=randn(1,npts); %fringe=nIterations*dels; %npts=npts+fringe; for i=1:nIterations, dnoise=[noise(dels+1:npts) noise(1:dels)]; dnoise=dnoise.*IRNgain; noise=noise+dnoise; end; switch stimComponents.type case 'pinkNoise' noise=UTIL_lowpassFilterFreq(noise, 10000, 1/FS); end rms=(mean(noise.^2)^.5); %should be 20 microPascals for 0 dB SPL adjust=20e-6/rms; noise=noise*adjust; rms=(mean(noise.^2)^.5); amplitude=10.^(stimComponents.amplitudesdB/20); noise=amplitude*noise; % rms=(mean(noise.^2)^.5); % dBnoise=20*log10(rms/20e-6) %------------------------------------------------------------------ makeRPN function RPN=makeRPN(globalStimParams, stim) % 'period' is a collection of samples - AAABCD % you need % % stim.type='RPN'; % stim.toneDuration=.2; % stim.amplitudesdB=50; % stim.beginSilence=.01; % stim.rampOnDur=.002; % stim.rampOffDur=-1; % % stim.sampleDuration=.005; %200 Hz pitch % stim.nSimilarSamples=5; % pitch strength % stim.nIndependentSamples=1% dilutes strength % % % Mandatory structure fields % globalStimParams.FS=44100; % globalStimParams.overallDuration=.21; % s % % globalStimParams.doPlot=1; % globalStimParams.doPlay=1; % [audio, msg]=stimulusCreate(globalStimParams, stim); FS=globalStimParams.FS; ptsPerSample=floor(stim.sampleDuration*FS); samplesPerPeriod=stim.nSimilarSamples+stim.nIndependentSamples; periodDuration=samplesPerPeriod*stim.sampleDuration; totalNumPeriods=2*floor(stim.toneDuration/periodDuration); % longer than necessary if totalNumPeriods<1 error('stimulusCreate: RPN, stimulus duration needs to be longer') end RPN=[]; for j=1:totalNumPeriods noise=randn(1,ptsPerSample); for i=1:stim.nSimilarSamples RPN=[RPN noise]; end for i=1:stim.nIndependentSamples noise=randn(1,ptsPerSample); RPN=[RPN noise]; end end targetStimulusLength=round(stim.toneDuration/FS); RPN=RPN(1:floor(stim.toneDuration*FS)); % take enough for stimulus rms=(mean(RPN.^2)^.5); %should be 20 microPascals for 0 dB SPL adjust=20e-6/rms; RPN=RPN*adjust; rms=(mean(RPN.^2)^.5); amplitude=10.^(stim.amplitudesdB/20); RPN=amplitude*RPN; % rms=(mean(noise.^2)^.5); % dBnoise=20*log10(rms/20e-6) %--------------------------------------------------------------------applyRampOn function signal=applyRampOn(signal, rampDur, rampOnTime, sampleRate) %applyRampOn applies raised cosine ramp %rampOntime is the time at which the ramp begins %At all other times the mask has a value of 1 %signal=applyRampOn(signal, rampDur, rampOnTime, sampleRate) rampDurPoints=round(rampDur*sampleRate); rampOn= (1+cos(pi:pi/(rampDurPoints-1): 2*pi))/2'; sigDurPoints=length(signal); mask(1:sigDurPoints)=1; rampOnStartIndex=round(rampOnTime*sampleRate+1); mask(rampOnStartIndex: rampOnStartIndex+ rampDurPoints-1)=rampOn; signal=signal.*mask; %plot(mask) %--------------------------------------------------------------------applyRampOff function signal=applyRampOff(signal, rampDur, rampOffTime, sampleRate) %applyRampOn applies raised cosine squared ramp %rampOffTime is the time at which the ramp begins %At all other times the mask has a value of 1 % signal=applyRampOff(signal, rampDur, rampOffTime, sampleRate) rampDurPoints=round(rampDur*sampleRate); rampOff= (1+cos(0:pi/(rampDurPoints-1): pi))/2'; sigDurPoints=length(signal); mask=ones(1,sigDurPoints); rampOffStartIndex=round(rampOffTime*sampleRate+1); mask(rampOffStartIndex: round(rampOffStartIndex+ rampDurPoints-1))=rampOff; if length(mask)>sigDurPoints, mask=mask(1:sigDurPoints); end signal=signal.*mask; %plot(mask) function signal=applyGaussianRamps(signal, sigma, sampleRate) dt=1/sampleRate; time=dt:dt:dt*length(signal); ramp=1-exp(-time.^2/(2*sigma^2)); % apply onset ramp signal=signal.*ramp; % apply offset ramp ramp=fliplr(ramp); signal=signal.*ramp; %--------------------------------------------------------------------checkDescriptors function [globalStimParams, stimComponents]=... checkDescriptors(globalStimParams, stimComponents) try % if FS exists, it takes priority globalStimParams.dt=1/globalStimParams.FS; catch % otherwise set FS using dt globalStimParams.FS=1/globalStimParams.dt; end sampleRate=globalStimParams.FS; globalStimParams.nSignalPoints=... round(globalStimParams.overallDuration* sampleRate); % optional field (ears) try globalStimParams.ears; catch % default: dichotic. globalStimParams.ears='dichotic'; end % audioOutCorrection is optional % audioOutCorrection is a scalar for reducing the sound try globalStimParams.audioOutCorrection; catch % set to 1 if omitted globalStimParams.audioOutCorrection=1; end try globalStimParams.doPlay; catch % default plays sound only if explicitly requested globalStimParams.doPlay=0; end try globalStimParams.doPlot; catch % no plotting unless explicitly requested globalStimParams.doPlot=0; end [ears nComponentSounds]=size(stimComponents); for ear=1:2 % 1=left/ 2=right % create a list of components whose type is specified % if no type is specified assume that it is an empty specification % this is allowed validStimComponents=[]; for i=1:nComponentSounds try if ~isempty(stimComponents(ear,i).type) validStimComponents=[validStimComponents i]; end catch end end for componentNo=validStimComponents % If no AM filed is present, create it for completeness if ~isfield(stimComponents(ear,componentNo),'AMfrequency') |... ~isfield(stimComponents(ear,componentNo),'AMdepth') stimComponents(ear,componentNo).AMfrequency=0; stimComponents(ear,componentNo).AMdepth=0; end % all signals must have durations, amplitudes and ramps if ... isempty(stimComponents(ear,componentNo).type) |... isempty(stimComponents(ear,componentNo).toneDuration) |... isempty(stimComponents(ear,componentNo).amplitudesdB) |... isempty(stimComponents(ear,componentNo).rampOnDur) descriptorError( 'missing stimComponent descriptor', stimComponents, ear, componentNo) end try, stimComponents(ear,componentNo).endSilence; catch, stimComponents(ear,componentNo).endSilence=-1; end % ramp checks do not apply to file input if ~strcmp(stimComponents(ear,componentNo).type, 'file') % match offset ramp to onset if not explicitly specified if stimComponents(ear,componentNo).rampOffDur==-1, stimComponents(ear,componentNo).rampOffDur=stimComponents(ear,componentNo).rampOnDur; end % ramps must be shorter than the stimulus if stimComponents(ear,componentNo).rampOffDur> stimComponents(ear,componentNo).toneDuration | ... stimComponents(ear,componentNo).rampOnDur> stimComponents(ear,componentNo).toneDuration descriptorError( 'ramp longer than sound component', stimComponents, ear, componentNo) end end % end silence is measured to fit into the global duration if stimComponents(ear,componentNo).endSilence==-1, stimComponents(ear,componentNo).endSilence=... globalStimParams.overallDuration-... stimComponents(ear,componentNo).beginSilence -... stimComponents(ear,componentNo).toneDuration; endSilenceNpoints=stimComponents(ear,componentNo).endSilence... *sampleRate; end if stimComponents(ear,componentNo).endSilence<0 globalStimParams descriptorError( 'component durations greater than overallDuration', stimComponents, ear, componentNo) end % check overall duration of this component against global duration totalDuration= ... stimComponents(ear,componentNo).beginSilence+... stimComponents(ear,componentNo).toneDuration+... stimComponents(ear,componentNo).endSilence; % avoid annoying error message for single stimulus component if ears==1 && nComponentSounds==1 globalStimParams.overallDuration= totalDuration; end % check total duration totalSignalPoints= round(totalDuration* sampleRate); if totalSignalPoints >globalStimParams.nSignalPoints descriptorError( 'Signal component duration does not match overall duration ', stimComponents, ear, componentNo) end if isfield(stimComponents(ear,componentNo), 'filter') if ~isequal(length(stimComponents(ear,componentNo).filter), 3) descriptorError( 'Filter parameter must have three elements ', stimComponents, ear, componentNo) end end end % component % ?? if strcmp(globalStimParams.ears,'monoticL') | strcmp(globalStimParams.ears, 'monoticR'), break, end end % ear %-------------------------------------------------------------------- descriptorError function descriptorError( msg, stimComponents, ear, componentNo) stimComponents(ear, componentNo) disp(' *********** **************** ************') disp([ '...Error in stimComponents description: ']) disp([msg ]) disp(['Ear = ' num2str(ear) ' component No ' num2str(componentNo)]) disp(' *********** **************** ************') error('myError ') %-------------------------------------------------------------------- normalize function [normalizedSignal, gain]= normalize(signal) % normalize (signal) maxSignal=max(max(signal)); minSignal=min(min(signal)); if -minSignal>maxSignal, normalizingFactor=-minSignal; else normalizingFactor=maxSignal; end normalizingFactor=1.01*normalizingFactor; gain= 20*log10(normalizingFactor); normalizedSignal=signal/normalizingFactor; %--------------------------------------------------------------------Butterworth function y=Butterworth (x, dt, fl, fu, order) % Butterworth (x, dt, fu, fl, order) % Taken from Yuel and beauchamp page 261 % NB error in their table for K (see their text) % x is original signal % fu, fl upper and lower cutoff % order is the number of times the filter is applied q=(pi*dt*(fu-fl)); J=1/(1+ cot(q)); K= (2*cos(pi*dt*(fu+fl)))/(1+tan(q)*cos(q)); L= (tan(q)-1)/(tan(q)+1); b=[J -J]; a=[1 -K -L]; for i=1:order y=filter(b, a, x); x=y; end % -------------------------------------------------------- UTIL_amp2dB function [y] = UTIL_amp2dB (x, ref) % Calculates a dB (ref. ref) value 'y' from a peak amplitude number 'x'. % if ref omitted treat as dB % Check the number of arguments that are passed in. if (nargin < 2) ref=1; end if (nargin > 2) error ('Too many arguments'); end % Check arguments. if x < 0.0 error ('Can not calculate the log10 of a negative number'); elseif x == 0.0 warning ('log10 of zero. The result is set to -1000.0'); y = -1000.0; else % Do calculations. y = 20.0 * log10(x/(sqrt(2)*ref)); end %-------------------------------------------------------------------- FFT function showFFT (getFFTsignal, dt) color='r'; figure(2), clf, hold off % trim initial silence idx=find(getFFTsignal>0); if ~isempty(idx) getFFTsignal=getFFTsignal(idx(1):end); end %trim final silence getFFTsignal=getFFTsignal(end:-1:1); idx=find(getFFTsignal>0); if ~isempty(idx) getFFTsignal=getFFTsignal(idx(1):end); getFFTsignal=getFFTsignal(end:-1:1); end % Analyse make stimComponents length a power of 2 x=length(getFFTsignal); squareLength=2; while squareLength<=x squareLength=squareLength*2; end squareLength=round(squareLength/2); getFFTsignal=getFFTsignal(1:squareLength); n=length(getFFTsignal); minf=100; maxf=20000; fft_result = fft(getFFTsignal, n); % Compute FFT of the input signal. fft_power = fft_result .* conj(fft_result);% / n; % Compute power spectrum. Dividing by 'n' we get the power spectral density. fft_phase = angle(fft_result); % Compute the phase spectrum. frequencies = (1/dt)*(1:n/2)/n; fft_power=fft_power(1:length(fft_power)/2); % remove mirror frequencies fft_phase=fft_phase(1:length(fft_phase)/2); % remove mirror frequencies fft_powerdB = UTIL_amp2dB (fft_power, max(fft_power)); % convert to dB % jags=find(diff(fft_phase)>0); % unwrap phase % for i=jags, fft_phase(i+1:end)=fft_phase(i+1:end)-2*pi; end xlim([minf maxf]) semilogx(frequencies, fft_powerdB-max(fft_powerdB), color) ylim([ -20 5]) function y=UTIL_lowpassFilterFreq(x, cutOffFrequency, dt) % UTIL_lowpassFilterFreq multi-channel filter % % Usage: % output=UTIL_lowpassFilterFreq(input, cutOffFrequency, dt) % % cutoff should be around 100 Hz for audio work % dt should be <1/50000 s for audio work % % Attenuation is - 6 dB per octave above cutoff. sampleRate=1/dt; if 4*cutOffFrequency>sampleRate 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' ]) cutOffFrequency=sampleRate/4; end tau=1/(2*pi*cutOffFrequency); y=zeros(size(x)); [numChannels numPoints]=size(x); for i=1:numChannels y(i,:)=filter([dt/tau], [1 -(1-dt/tau)], x(i,:)); end