rmeddis@38: function test_toneInNoise
rmeddis@38: % test_MAP1_14 is a general purpose test routine that can be adjusted to
rmeddis@38: % test a number of different applications of MAP1_14
rmeddis@38: %
rmeddis@38: % A range of options are supplied in the early part of the program
rmeddis@38: %
rmeddis@38: % One use of the function is to create demonstrations; filenames <demoxx>
rmeddis@38: %  to illustrate particular features
rmeddis@38: %
rmeddis@38: % #1
rmeddis@38: % Identify the file (in 'MAPparamsName') containing the model parameters
rmeddis@38: %
rmeddis@38: % #2
rmeddis@38: % Identify the kind of model required (in 'AN_spikesOrProbability').
rmeddis@38: %  A full brainstem model (spikes) can be computed or a shorter model
rmeddis@38: %  (probability) that computes only so far as the auditory nerve
rmeddis@38: %
rmeddis@38: % #3
rmeddis@38: % Choose between a tone signal or file input (in 'signalType')
rmeddis@38: %
rmeddis@38: % #4
rmeddis@38: % Set the signal rms level (in leveldBSPL)
rmeddis@38: %
rmeddis@38: % #5
rmeddis@38: % Identify the channels in terms of their best frequencies in the vector
rmeddis@38: %  BFlist.
rmeddis@38: %
rmeddis@38: % Last minute changes to the parameters fetched earlier can be made using
rmeddis@38: %  the cell array of strings 'paramChanges'.
rmeddis@38: %  Each string must have the same format as the corresponding line in the
rmeddis@38: %  file identified in 'MAPparamsName'
rmeddis@38: %
rmeddis@38: % When the demonstration is satisfactory, freeze it by renaming it <demoxx>
rmeddis@38: 
rmeddis@38: global dt dtSpikes  savedBFlist saveAN_spikesOrProbability saveMAPparamsName...
rmeddis@38:     savedInputSignal OMEextEarPressure TMoutput OMEoutput ARattenuation ...
rmeddis@38:     DRNLoutput IHC_cilia_output IHCrestingCiliaCond IHCrestingV...
rmeddis@38:     IHCoutput ANprobRateOutput ANoutput savePavailable ANtauCas  ...
rmeddis@38:     CNtauGk CNoutput  ICoutput ICmembraneOutput ICfiberTypeRates ...
rmeddis@38:     MOCattenuation
rmeddis@38: global OMEParams DRNLParams IHC_cilia_RPParams IHCpreSynapseParams
rmeddis@38: global AN_IHCsynapseParams MacGregorParams MacGregorMultiParams
rmeddis@38: global ICrate
rmeddis@38: 
rmeddis@38: 
rmeddis@38: dbstop if error
rmeddis@38: restorePath=path;
rmeddis@38: addpath (['..' filesep 'MAP'],    ['..' filesep 'wavFileStore'], ...
rmeddis@38:     ['..' filesep 'utilities'])
rmeddis@38: 
rmeddis@38: %%  #1 parameter file name
rmeddis@38: MAPparamsName='Normal';
rmeddis@38: 
rmeddis@38: 
rmeddis@38: %% #2 probability (fast) or spikes (slow) representation
rmeddis@38: AN_spikesOrProbability='spikes';
rmeddis@38: %   or
rmeddis@38: AN_spikesOrProbability='probability'; 
rmeddis@38: 
rmeddis@38: 
rmeddis@38: %% #3 pure tone, harmonic sequence or speech file input
rmeddis@38: signalType= 'tones';
rmeddis@38: sampleRate= 44100;          % must agree with noise
rmeddis@38: duration=0.010;                 % seconds
rmeddis@38: rampDuration=.001;              % raised cosine ramp (seconds)
rmeddis@38: beginSilence=0.010;               
rmeddis@38: endSilence=0.010;                  
rmeddis@38: toneFrequency= 2000;            % or a pure tone (Hz)
rmeddis@38: 
rmeddis@38: %   or
rmeddis@38: % harmonic sequence (Hz)
rmeddis@38: % F0=210;
rmeddis@38: % toneFrequency= F0:F0:8000;    
rmeddis@38: 
rmeddis@38: %   or
rmeddis@38: % signalType= 'file';
rmeddis@38: % fileName='twister_44kHz';
rmeddis@38: 
rmeddis@38: 
rmeddis@38: 
rmeddis@38: %% #4 rms level
rmeddis@38: % signal details
rmeddis@38: leveldBSPL= 70;                  % dB SPL (80 for Lieberman)
rmeddis@38: leveldBSPLNoise=30;
rmeddis@38: 
rmeddis@38: %% #5 number of channels in the model
rmeddis@38: %   21-channel model (log spacing)
rmeddis@38: numChannels=21;
rmeddis@38: lowestBF=250; 	highestBF= 8000;
rmeddis@38: BFlist=round(logspace(log10(lowestBF), log10(highestBF), numChannels));
rmeddis@38: 
rmeddis@38: % %   or specify your own channel BFs
rmeddis@38: % numChannels=1;
rmeddis@38: % BFlist=toneFrequency;
rmeddis@38: 
rmeddis@38: 
rmeddis@38: %% #6 change model parameters
rmeddis@38: 
rmeddis@38: paramChanges={};
rmeddis@38: 
rmeddis@38: % Parameter changes can be used to change one or more model parameters
rmeddis@38: %  *after* the MAPparams file has been read
rmeddis@38: % This example declares only one fiber type with a calcium clearance time
rmeddis@38: % constant of 80e-6 s (HSR fiber) when the probability option is selected.
rmeddis@38: % paramChanges={'AN_IHCsynapseParams.ANspeedUpFactor=5;', ...
rmeddis@38: %     'IHCpreSynapseParams.tauCa=86e-6; '};
rmeddis@38: paramChanges={'DRNLParams.MOCtauProb =.25;', ...
rmeddis@38:     'DRNLParams.rateToAttenuationFactorProb = 0.02; '};
rmeddis@38: % paramChanges={'DRNLParams.MOCtauProb =.15;', ...
rmeddis@38: %     'DRNLParams.rateToAttenuationFactorProb = 0.00; '};
rmeddis@38: 
rmeddis@38: 
rmeddis@38: %% delare 'showMap' options to control graphical output
rmeddis@38: showMapOptions.printModelParameters=1;   % prints all parameters
rmeddis@38: showMapOptions.showModelOutput=1;       % plot of all stages
rmeddis@38: showMapOptions.printFiringRates=1;      % prints stage activity levels
rmeddis@38: showMapOptions.showACF=0;               % shows SACF (probability only)
rmeddis@38: showMapOptions.showEfferent=1;          % tracks of AR and MOC
rmeddis@38: showMapOptions.surfProbability=1;       % 2D plot of HSR response 
rmeddis@38: showMapOptions.surfSpikes=1;            % 2D plot of spikes histogram
rmeddis@38: showMapOptions.ICrates=0;               % IC rates by CNtauGk
rmeddis@38: 
rmeddis@38: % disable certain silly options
rmeddis@38: if strcmp(AN_spikesOrProbability, 'spikes')
rmeddis@38:     % avoid nonsensical options
rmeddis@38:     showMapOptions.surfProbability=0;
rmeddis@38:     showMapOptions.showACF=0;
rmeddis@38: end
rmeddis@38: 
rmeddis@38: if strcmp(signalType, 'file')
rmeddis@38:     % needed for labeling plot
rmeddis@38:     showMapOptions.fileName=fileName;
rmeddis@38: else
rmeddis@38:     showMapOptions.fileName=[];
rmeddis@38: end
rmeddis@38: 
rmeddis@38: %% Generate stimuli
rmeddis@38: 
rmeddis@38: switch signalType
rmeddis@38:     case 'tones'
rmeddis@38:         % Create pure tone stimulus
rmeddis@38:         dt=1/sampleRate; % seconds
rmeddis@38:         time=dt: dt: duration;
rmeddis@38:         inputSignal=sum(sin(2*pi*toneFrequency'*time), 1);
rmeddis@38:         amp=10^(leveldBSPL/20)*28e-6;   % converts to Pascals (peak)
rmeddis@38:         inputSignal=amp*inputSignal;
rmeddis@38:         % apply ramps
rmeddis@38:         % catch rampTime error
rmeddis@38:         if rampDuration>0.5*duration, rampDuration=duration/2; end
rmeddis@38:         rampTime=dt:dt:rampDuration;
rmeddis@38:         ramp=[0.5*(1+cos(2*pi*rampTime/(2*rampDuration)+pi)) ...
rmeddis@38:             ones(1,length(time)-length(rampTime))];
rmeddis@38:         inputSignal=inputSignal.*ramp;
rmeddis@38:         ramp=fliplr(ramp);
rmeddis@38:         inputSignal=inputSignal.*ramp;
rmeddis@38:         % add silence
rmeddis@38:         intialSilence= zeros(1,round(beginSilence/dt));
rmeddis@38:         finalSilence= zeros(1,round(endSilence/dt));
rmeddis@38:         inputSignal= [intialSilence inputSignal finalSilence];
rmeddis@38:         
rmeddis@38: %         [inputNoise sampleRateN]=wavread('babble');
rmeddis@38:         [inputNoise sampleRateN]=wavread('white noise');
rmeddis@38:         inputNoise=inputNoise(1:length(inputSignal));
rmeddis@38:         inputNoise=inputNoise(:,1);
rmeddis@38:         targetRMS=20e-6*10^(leveldBSPLNoise/20);
rmeddis@38:         rms=(mean(inputNoise.^2))^0.5;
rmeddis@38:         amp=targetRMS/rms;
rmeddis@38:         inputNoise=inputNoise*amp;
rmeddis@38:         inputSignal=inputSignal+inputNoise';
rmeddis@38:         figure(2), subplot(2,1,1)
rmeddis@38:         time=dt:dt:dt*length(inputSignal);
rmeddis@38:         plot(time,inputSignal,'k')
rmeddis@38: 
rmeddis@38: 
rmeddis@38:     case 'file'
rmeddis@38:         %% file input simple or mixed
rmeddis@38:         [inputSignal sampleRate]=wavread(fileName);
rmeddis@38:         dt=1/sampleRate;
rmeddis@38:         inputSignal=inputSignal(:,1);
rmeddis@38:         targetRMS=20e-6*10^(leveldBSPL/20);
rmeddis@38:         rms=(mean(inputSignal.^2))^0.5;
rmeddis@38:         amp=targetRMS/rms;
rmeddis@38:         inputSignal=inputSignal*amp;
rmeddis@38:         intialSilence= zeros(1,round(0.1/dt));
rmeddis@38:         finalSilence= zeros(1,round(0.2/dt));
rmeddis@38:         inputSignal= [intialSilence inputSignal' finalSilence];
rmeddis@38:         
rmeddis@38: end
rmeddis@38: 
rmeddis@38: 
rmeddis@38: %% run the model
rmeddis@38: tic
rmeddis@38: fprintf('\n')
rmeddis@38: disp(['Signal duration= ' num2str(length(inputSignal)/sampleRate)])
rmeddis@38: disp([num2str(numChannels) ' channel model: ' AN_spikesOrProbability])
rmeddis@38: disp('Computing ...')
rmeddis@38: 
rmeddis@38: MAP1_14(inputSignal, sampleRate, BFlist, ...
rmeddis@38:     MAPparamsName, AN_spikesOrProbability, paramChanges);
rmeddis@38: 
rmeddis@38: 
rmeddis@38: %% the model run is now complete. Now display the results
rmeddis@38: UTIL_showMAP(showMapOptions, paramChanges)
rmeddis@38: 
rmeddis@38: if strcmp(signalType,'tones')
rmeddis@38:     disp(['duration=' num2str(duration)])
rmeddis@38:     disp(['level=' num2str(leveldBSPL)])
rmeddis@38:     disp(['toneFrequency=' num2str(toneFrequency)])
rmeddis@38:     global DRNLParams
rmeddis@38:     disp(['attenuation factor =' ...
rmeddis@38:         num2str(DRNLParams.rateToAttenuationFactor, '%5.3f') ])
rmeddis@38:     disp(['attenuation factor (probability)=' ...
rmeddis@38:         num2str(DRNLParams.rateToAttenuationFactorProb, '%5.3f') ])
rmeddis@38:     disp(AN_spikesOrProbability)
rmeddis@38: end
rmeddis@38: 
rmeddis@38: figure(2), subplot(2,1,2), plot(ANprobRateOutput(13+21,:))
rmeddis@38: disp([ 'peak channel 13: ' num2str(max(ANprobRateOutput(13+21,:)))])
rmeddis@38: 
rmeddis@38: for i=1:length(paramChanges)
rmeddis@38:     disp(paramChanges{i})
rmeddis@38: end
rmeddis@38: 
rmeddis@38: path(restorePath)
rmeddis@38: