rmeddis@9: function demoTwisterProbability
rmeddis@9: 
rmeddis@9: % MAPdemo runs the MATLAB auditory periphery model (MAP1_14) as far as
rmeddis@9: %  the AN (probabilities) or IC (spikes) with graphical output
rmeddis@9: 
rmeddis@25: dbstop if error
rmeddis@24: restorePath=path;
rmeddis@24: addpath (['..' filesep 'MAP'],    ['..' filesep 'wavFileStore'], ...
rmeddis@24:     ['..' filesep 'utilities'])
rmeddis@9: 
rmeddis@9: %%  #1 parameter file name
rmeddis@9: MAPparamsName='Normal';
rmeddis@9: 
rmeddis@9: 
rmeddis@24: %% #2 probability (fast) 
rmeddis@9: AN_spikesOrProbability='probability';
rmeddis@9: 
rmeddis@9: 
rmeddis@24: %% #3  speech file input
rmeddis@9: signalType= 'file';
rmeddis@9: fileName='twister_44kHz';
rmeddis@9: 
rmeddis@9: 
rmeddis@9: %% #4 rms level
rmeddis@24: leveldBSPL=60;        % dB SPL
rmeddis@9: 
rmeddis@9: 
rmeddis@9: %% #5 number of channels in the model
rmeddis@9: %   21-channel model (log spacing)
rmeddis@9: numChannels=21;
rmeddis@24: lowestBF=250; 	highestBF= 8000;
rmeddis@9: BFlist=round(logspace(log10(lowestBF), log10(highestBF), numChannels));
rmeddis@9: 
rmeddis@24: %% #6 no change to model parameters
rmeddis@9: paramChanges=[];
rmeddis@9: 
rmeddis@9: 
rmeddis@9: %% delare showMap options
rmeddis@24: showMapOptions.printModelParameters=1;
rmeddis@24: showMapOptions.showModelOutput=1;
rmeddis@24: showMapOptions.printFiringRates=1;
rmeddis@24: showMapOptions.showACF=0;
rmeddis@24: showMapOptions.showEfferent=0;
rmeddis@24: showMapOptions.surfProbability=1;       % 2D plot of HSR response 
rmeddis@9: 
rmeddis@9: %% Generate stimuli
rmeddis@9: 
rmeddis@9: switch signalType
rmeddis@9:     case 'tones'
rmeddis@9:         inputSignal=createMultiTone(sampleRate, toneFrequency, ...
rmeddis@9:             leveldBSPL, duration, rampDuration);
rmeddis@24: 
rmeddis@9:     case 'file'
rmeddis@24:         [inputSignal sampleRate]=wavread(fileName);
rmeddis@9:         inputSignal(:,1);
rmeddis@9:         targetRMS=20e-6*10^(leveldBSPL/20);
rmeddis@9:         rms=(mean(inputSignal.^2))^0.5;
rmeddis@9:         amp=targetRMS/rms;
rmeddis@9:         inputSignal=inputSignal*amp;
rmeddis@9: end
rmeddis@9: 
rmeddis@9: 
rmeddis@9: %% run the model
rmeddis@9: tic
rmeddis@9: 
rmeddis@26: fprintf('\n')
rmeddis@26: disp(['Signal duration= ' num2str(length(inputSignal)/sampleRate)])
rmeddis@26: disp([num2str(numChannels) ' channel model'])
rmeddis@26: disp('Computing ...')
rmeddis@26: 
rmeddis@9: MAP1_14(inputSignal, sampleRate, BFlist, ...
rmeddis@9:     MAPparamsName, AN_spikesOrProbability, paramChanges);
rmeddis@24: 
rmeddis@9: 
rmeddis@9: % the model run is now complete. Now display the results
rmeddis@26: UTIL_showMAP(showMapOptions, paramChanges)
rmeddis@9: 
rmeddis@9: toc
rmeddis@9: path(restorePath)
rmeddis@9: 
rmeddis@9: 
rmeddis@26: 
rmeddis@9: function inputSignal=createMultiTone(sampleRate, toneFrequency, ...
rmeddis@9:     leveldBSPL, duration, rampDuration)
rmeddis@9: % Create pure tone stimulus
rmeddis@9: dt=1/sampleRate; % seconds
rmeddis@9: time=dt: dt: duration;
rmeddis@9: inputSignal=sum(sin(2*pi*toneFrequency'*time), 1);
rmeddis@9: amp=10^(leveldBSPL/20)*28e-6;   % converts to Pascals (peak)
rmeddis@9: inputSignal=amp*inputSignal;
rmeddis@9: 
rmeddis@9: % apply ramps
rmeddis@9: % catch rampTime error
rmeddis@9: if rampDuration>0.5*duration, rampDuration=duration/2; end
rmeddis@9: rampTime=dt:dt:rampDuration;
rmeddis@9: ramp=[0.5*(1+cos(2*pi*rampTime/(2*rampDuration)+pi)) ...
rmeddis@9:     ones(1,length(time)-length(rampTime))];
rmeddis@9: inputSignal=inputSignal.*ramp;
rmeddis@9: ramp=fliplr(ramp);
rmeddis@9: inputSignal=inputSignal.*ramp;
rmeddis@9: 
rmeddis@9: % add 10 ms silence
rmeddis@9: silence= zeros(1,round(0.03/dt));
rmeddis@9: % inputSignal= [silence inputSignal silence];
rmeddis@9: