rmeddis@32: function pitchModel_RM
rmeddis@32: % Modification of testMAP_14 to replicate the pitch model published
rmeddis@32: %  in JASA 2006.
rmeddis@32: %
rmeddis@32: % test_MAP1_14 is a general purpose test routine that can be adjusted to
rmeddis@32: % test a number of different applications of MAP1_14
rmeddis@32: %
rmeddis@32: % A range of options are supplied in the early part of the program
rmeddis@32: %
rmeddis@32: % One use of the function is to create demonstrations; filenames <demoxx>
rmeddis@32: %  to illustrate particular features
rmeddis@32: %
rmeddis@32: % #1
rmeddis@32: % Identify the file (in 'MAPparamsName') containing the model parameters
rmeddis@32: %
rmeddis@32: % #2
rmeddis@32: % Identify the kind of model required (in 'AN_spikesOrProbability').
rmeddis@32: %  A full brainstem model (spikes) can be computed or a shorter model
rmeddis@32: %  (probability) that computes only so far as the auditory nerve
rmeddis@32: %
rmeddis@32: % #3
rmeddis@32: % Choose between a tone signal or file input (in 'signalType')
rmeddis@32: %
rmeddis@32: % #4
rmeddis@32: % Set the signal rms level (in leveldBSPL)
rmeddis@32: %
rmeddis@32: % #5
rmeddis@32: % Identify the channels in terms of their best frequencies in the vector
rmeddis@32: %  BFlist.
rmeddis@32: %
rmeddis@32: % Last minute changes to the parameters fetched earlier can be made using
rmeddis@32: %  the cell array of strings 'paramChanges'.
rmeddis@32: %  Each string must have the same format as the corresponding line in the
rmeddis@32: %  file identified in 'MAPparamsName'
rmeddis@32: %
rmeddis@32: % When the demonstration is satisfactory, freeze it by renaming it <demoxx>
rmeddis@32: 
rmeddis@32: dbstop if error
rmeddis@32: restorePath=path;
rmeddis@32: addpath (['..' filesep 'MAP'],    ['..' filesep 'wavFileStore'], ...
rmeddis@32:     ['..' filesep 'utilities'])
rmeddis@32: 
rmeddis@32: % Pitch model modification here
rmeddis@32: global ICrate           % used to collect rate profile from showMAP temporary
rmeddis@32: rates=[]; F0count=0;
rmeddis@32: 
rmeddis@32: % F0s=[150 200 250];          % fundamental frequency
rmeddis@32: % harmonics= 3:5;
rmeddis@32: 
rmeddis@32: % F0s=[3000];          % fundamental frequency
rmeddis@32: F0s=50:5:1000;
rmeddis@32: harmonics= 1;
rmeddis@32: % F0s=150;
rmeddis@32: for F0=F0s
rmeddis@32:     F0count=F0count+1;
rmeddis@32:     
rmeddis@32: 
rmeddis@32: %%  #1 parameter file name
rmeddis@32: MAPparamsName='Normal';
rmeddis@32: 
rmeddis@32: 
rmeddis@32: %% #2 probability (fast) or spikes (slow) representation
rmeddis@32: AN_spikesOrProbability='spikes';
rmeddis@32: 
rmeddis@32: % or
rmeddis@32: % NB probabilities are not corrected for refractory effects
rmeddis@32: % AN_spikesOrProbability='probability';
rmeddis@32: 
rmeddis@32: 
rmeddis@32: %% #3 pure tone, harmonic sequence or speech file input
rmeddis@32: signalType= 'tones';
rmeddis@32: sampleRate= 50000;
rmeddis@32: duration=0.50;                 % seconds
rmeddis@32: % toneFrequency= 1000;            % or a pure tone (Hz8
rmeddis@32: 
rmeddis@32: % F0=210;
rmeddis@32: toneFrequency= F0*harmonics;  % harmonic sequence (Hz)
rmeddis@32: 
rmeddis@32: rampDuration=.005;              % raised cosine ramp (seconds)
rmeddis@32: 
rmeddis@32: % or
rmeddis@32: 
rmeddis@32: % signalType= 'file';
rmeddis@32: % fileName='twister_44kHz';
rmeddis@32: 
rmeddis@32: 
rmeddis@32: %% #4 rms level
rmeddis@32: % signal details
rmeddis@32: leveldBSPL= 50;                  % dB SPL
rmeddis@32: 
rmeddis@32: 
rmeddis@32: %% #5 number of channels in the model
rmeddis@32: %   21-channel model (log spacing)
rmeddis@32: numChannels=21;
rmeddis@32: lowestBF=250; 	highestBF= 8000;
rmeddis@32: BFlist=round(logspace(log10(lowestBF), log10(highestBF), numChannels));
rmeddis@32: 
rmeddis@32: %   or specify your own channel BFs
rmeddis@32: % numChannels=1;
rmeddis@32: BFlist=toneFrequency;
rmeddis@32: % BFlist=500;
rmeddis@32: 
rmeddis@32: 
rmeddis@32: %% #6 change model parameters
rmeddis@32: paramChanges=[];
rmeddis@32: 
rmeddis@32: % or
rmeddis@32: % Parameter changes can be used to change one or more model parameters
rmeddis@32: %  *after* the MAPparams file has been read
rmeddis@32: % This example declares only one fiber type with a calcium clearance time
rmeddis@32: % constant of 80e-6 s (HSR fiber) when the probability option is selected.
rmeddis@32: 
rmeddis@32: % paramChanges={'AN_IHCsynapseParams.ANspeedUpFactor=5;', ...
rmeddis@32: %     'IHCpreSynapseParams.tauCa=86e-6;'};
rmeddis@32: 
rmeddis@32: % paramChanges={'DRNLParams.rateToAttenuationFactorProb = 0;'};
rmeddis@32: 
rmeddis@32: % paramChanges={'IHCpreSynapseParams.tauCa=86e-6;',
rmeddis@32: %     'AN_IHCsynapseParams.numFibers=	1000;'};
rmeddis@32: 
rmeddis@32: % fixed MOC attenuation(using negative factor)
rmeddis@32: % paramChanges={'DRNLParams.rateToAttenuationFactorProb=-0.005;'};
rmeddis@32: 
rmeddis@32: % slow the CN chopping rate
rmeddis@32: % paramChanges={'IHCpreSynapseParams.tauCa= 70e-6;'...'
rmeddis@32: %     'MacGregorMultiParams.tauGk=	[0.75e-3:.0001 : 3e-3];'...
rmeddis@32: %     ' MacGregorParams.dendriteLPfreq=4000;'...
rmeddis@32: %     'MacGregorParams.tauGk=	1e-4;'...
rmeddis@32: %     'MacGregorParams.currentPerSpike=220e-8;'...
rmeddis@32: % };
rmeddis@32: paramChanges={...
rmeddis@32:     'MacGregorMultiParams.currentPerSpike=25e-9;'...
rmeddis@32:     'MacGregorMultiParams.tauGk=	[0.1e-3:.00005 : 1e-3];'...
rmeddis@32: 'MacGregorParams.currentPerSpike=40e-9;'...
rmeddis@32: };
rmeddis@32: 
rmeddis@32: %% delare 'showMap' options to control graphical output
rmeddis@32: 
rmeddis@32: showMapOptions.printModelParameters=0;   % prints all parameters
rmeddis@32: showMapOptions.showModelOutput=1;       % plot of all stages
rmeddis@32: showMapOptions.printFiringRates=1;      % prints stage activity levels
rmeddis@32: showMapOptions.showACF=0;               % shows SACF (probability only)
rmeddis@32: showMapOptions.showEfferent=0;          % tracks of AR and MOC
rmeddis@32: showMapOptions.surfProbability=0;       % 2D plot of HSR response 
rmeddis@32: showMapOptions.surfSpikes=0;            % 2D plot of spikes histogram
rmeddis@32: showMapOptions.ICrates=1;               % IC rates by CNtauGk
rmeddis@32: 
rmeddis@32: % disable certain silly options
rmeddis@32: if strcmp(AN_spikesOrProbability, 'spikes')
rmeddis@32:     % avoid nonsensical options
rmeddis@32:     showMapOptions.surfProbability=0;
rmeddis@32:     showMapOptions.showACF=0;
rmeddis@32: else
rmeddis@32:     showMapOptions.surfSpikes=0;
rmeddis@32: end
rmeddis@32: if strcmp(signalType, 'file')
rmeddis@32:     % needed for labeling plot
rmeddis@32:     showMapOptions.fileName=fileName;
rmeddis@32: else
rmeddis@32:     showMapOptions.fileName=[];
rmeddis@32: end
rmeddis@32: 
rmeddis@32: %% Generate stimuli
rmeddis@32: 
rmeddis@32: switch signalType
rmeddis@32:     case 'tones'
rmeddis@32:         inputSignal=createMultiTone(sampleRate, toneFrequency, ...
rmeddis@32:             leveldBSPL, duration, rampDuration);
rmeddis@32: 
rmeddis@32:     case 'file'
rmeddis@32:         %% file input simple or mixed
rmeddis@32:         [inputSignal sampleRate]=wavread(fileName);
rmeddis@32:         dt=1/sampleRate;
rmeddis@32:         inputSignal=inputSignal(:,1);
rmeddis@32:         targetRMS=20e-6*10^(leveldBSPL/20);
rmeddis@32:         rms=(mean(inputSignal.^2))^0.5;
rmeddis@32:         amp=targetRMS/rms;
rmeddis@32:         inputSignal=inputSignal*amp;
rmeddis@32:         silence= zeros(1,round(0.1/dt));
rmeddis@32:         inputSignal= [silence inputSignal' silence];
rmeddis@32: end
rmeddis@32: 
rmeddis@32: 
rmeddis@32: %% run the model
rmeddis@32: tic
rmeddis@32: 
rmeddis@32: fprintf('\n')
rmeddis@32: disp(['Signal duration= ' num2str(length(inputSignal)/sampleRate)])
rmeddis@32: disp([num2str(numChannels) ' channel model'])
rmeddis@32: disp([num2str(F0) ' F0'])
rmeddis@32: disp('Computing ...')
rmeddis@32: 
rmeddis@32: MAP1_14(inputSignal, sampleRate, BFlist, ...
rmeddis@32:     MAPparamsName, AN_spikesOrProbability, paramChanges);
rmeddis@32: 
rmeddis@32: 
rmeddis@32: % the model run is now complete. Now display the results
rmeddis@32: UTIL_showMAP(showMapOptions, paramChanges)
rmeddis@32: 
rmeddis@32: %% pitch model Collect and analyse data
rmeddis@32: % ICrate is global and computed in showMAP
rmeddis@32: %  a vector of 'stage4' rates; one value for each tauCNGk
rmeddis@32: rates=[rates; ICrate];
rmeddis@32: figure(92), imagesc(rates)
rmeddis@32: ylabel ('F0 no'), xlabel('tauGk')
rmeddis@32: % figure(92), plot(rates), ylim([0 inf])
rmeddis@32: 
rmeddis@32: h=figure(99); CNmovie(F0count)=getframe(h);
rmeddis@32: figure(91), plot(rates'),ylim([0 inf])
rmeddis@32: pause (0.1)
rmeddis@32: path(restorePath)
rmeddis@32: 
rmeddis@32: end
rmeddis@32: %% show results
rmeddis@32: toc
rmeddis@32: figure(91), plot(F0s,rates'), xlabel('F0'), ylabel('rate'),ylim([0 inf])
rmeddis@32: % figure(99),clf,movie(CNmovie,1,4)
rmeddis@32: 
rmeddis@32: 
rmeddis@32: function inputSignal=createMultiTone(sampleRate, toneFrequency, ...
rmeddis@32:     leveldBSPL, duration, rampDuration)
rmeddis@32: % Create pure tone stimulus
rmeddis@32: dt=1/sampleRate; % seconds
rmeddis@32: time=dt: dt: duration;
rmeddis@32: inputSignal=sum(sin(2*pi*toneFrequency'*time), 1);
rmeddis@32: amp=10^(leveldBSPL/20)*28e-6;   % converts to Pascals (peak)
rmeddis@32: inputSignal=amp*inputSignal;
rmeddis@32: 
rmeddis@32: % apply ramps
rmeddis@32: % catch rampTime error
rmeddis@32: if rampDuration>0.5*duration, rampDuration=duration/2; end
rmeddis@32: rampTime=dt:dt:rampDuration;
rmeddis@32: ramp=[0.5*(1+cos(2*pi*rampTime/(2*rampDuration)+pi)) ...
rmeddis@32:     ones(1,length(time)-length(rampTime))];
rmeddis@32: inputSignal=inputSignal.*ramp;
rmeddis@32: ramp=fliplr(ramp);
rmeddis@32: inputSignal=inputSignal.*ramp;
rmeddis@32: 
rmeddis@32: % add 10 ms silence
rmeddis@32: silence= zeros(1,round(0.005/dt));
rmeddis@32: inputSignal= [silence inputSignal silence];
rmeddis@32: