rmeddis@38: % function [P dt lags SACF]= testACF
rmeddis@38: % testACF is a *script* to demonstrate the smoothed ACF of 
rmeddis@38: % Balaguer-Ballestera, E. Denham, S.L. and Meddis, R. (2008). 
rmeddis@38: %
rmeddis@38: % Convert this to a *function* by uncommenting the first line
rmeddis@38: %  The function returns the P matrix plotted in Figure 96.
rmeddis@38: %  If a function is used, the following outputs are returned:
rmeddis@38: %   P:  smoothed SACF (lags x time matrix)
rmeddis@38: %   dt: time interval between successive columns of P
rmeddis@38: %   lags: lags used in computing P
rmeddis@38: %   SACF: unsmoothed SACFs
rmeddis@38: %  
rmeddis@38: % A range of options are supplied in the early part of the program
rmeddis@38: %
rmeddis@38: % #1
rmeddis@38: % Identify the model parameter file (in 'MAPparamsName') 
rmeddis@38: %
rmeddis@38: % #2
rmeddis@38: % Identify the kind of model required (in 'AN_spikesOrProbability')
rmeddis@38: %  'probability' is recommended for ACF work
rmeddis@38: %
rmeddis@38: % #3
rmeddis@38: % Choose between a harmonic complex or file input
rmeddis@38: %  by commenting out unwanted code
rmeddis@38: %
rmeddis@38: % #4
rmeddis@38: % Set the signal rms level (in leveldBSPL)
rmeddis@38: %
rmeddis@38: % #5
rmeddis@38: % Identify the model channel BFs in the vector 'BFlist'.
rmeddis@38: %
rmeddis@38: % #6
rmeddis@38: % Last minute changes to the model parameters can be made using
rmeddis@38: %  the cell array of strings 'paramChanges'.
rmeddis@38: %  This is used here to control the details of the ACF computations
rmeddis@38: %  Read the notes in this section for more information
rmeddis@38: %
rmeddis@38: % displays:
rmeddis@38: % Figure 97 shows the AN response to the stimulus. this is a channel x time
rmeddis@38: %  display. The z-axis (and colour) is the AN fiber firing rate
rmeddis@38: %
rmeddis@38: % Figure 96 shows the P-matrix, the smoothed SACF.
rmeddis@38: %
rmeddis@38: % Figure 89 shows a summary of the evolution of the unsmoothed SACF 
rmeddis@38: %  over time. If you wish to take a snapshot of the P-matrix at a
rmeddis@38: %  particular time, this figure can help identify when to take it. 
rmeddis@38: %  The index on the y-axis, identifies the required row numbers 
rmeddis@38: %    of the P or SACF matrix, e.g. P(:,2000)
rmeddis@38: %
rmeddis@38: %  On request, (filteredSACFParams.plotACFs=1) Figure 89 shows the channel 
rmeddis@38: %   by channel ACFs at intervals during the computation as a movie.
rmeddis@38: %  The number of ACF displays is controlled by 'plotACFsInterval'
rmeddis@38: %   and the movie can be slowed or speeded up using 'plotMoviePauses' 
rmeddis@38: %   (see paramChanges section below).
rmeddis@38: 
rmeddis@38: % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
rmeddis@38: % This global will find results from MAP1_14
rmeddis@38: global savedInputSignal ANprobRateOutput ANoutput dt dtSpikes savedBFlist 
rmeddis@38:  % This global,from model parameter file
rmeddis@38: global filteredSACFParams                      
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: % User sets up requirements
rmeddis@38: %%  #1 parameter file name
rmeddis@38: MAPparamsName='Normal';                 % recommended
rmeddis@38: 
rmeddis@38: 
rmeddis@38: %% #2 probability (fast) or spikes (slow) representation: select one
rmeddis@38: % AN_spikesOrProbability='spikes';
rmeddis@38: %   or
rmeddis@38: AN_spikesOrProbability='probability';   % recommended
rmeddis@38: 
rmeddis@38: %% #3 A. harmonic sequence or B. speech file input
rmeddis@38: % Comment out unwanted code
rmeddis@38: % A. harmonic tone (Hz) - useful to demonstrate a broadband sound
rmeddis@38: sampleRate= 44100;              % recommended 44100
rmeddis@38: signalType= 'tones';
rmeddis@38: duration=0.100;                 % seconds
rmeddis@38: beginSilence=0.020;
rmeddis@38: endSilence=0.020;
rmeddis@38: rampDuration=.005;              % raised cosine ramp (seconds)
rmeddis@38: 
rmeddis@38: % toneFrequency is a vector of component frequencies
rmeddis@38: F0=120;
rmeddis@38: toneFrequency= [3*F0 4*F0 5*F0];
rmeddis@38: 
rmeddis@38: %   or
rmeddis@38: % B. file input
rmeddis@38: signalType= 'file';
rmeddis@38: fileName='Oh No';
rmeddis@38: fileName='1z67931a_44kHz';
rmeddis@38: 
rmeddis@38: 
rmeddis@38: %% #4 rms level
rmeddis@38: leveldBSPL= 60;                  % dB SPL (80 for Lieberman)
rmeddis@38: 
rmeddis@38: %% #5 number of channels in the model
rmeddis@38: %   21-channel model (log spacing of BFs)
rmeddis@38: numChannels=21;
rmeddis@38: lowestBF=150; 	highestBF= 6000;
rmeddis@38: BFlist=round(logspace(log10(lowestBF), log10(highestBF), numChannels));
rmeddis@38: 
rmeddis@38: %% #6 change model parameters
rmeddis@38: % Parameter changes can be used to change one or more model parameters
rmeddis@38: %  *after* the MAPparams file has been read (see manual)
rmeddis@38: 
rmeddis@38: % Take control of ACF parameters
rmeddis@38: %  The filteredACF parameters are set in the MAPparamsNormal file
rmeddis@38: %  However, it is convenient to change them here leving the file intacta
rmeddis@38:     minPitch=	80; maxPitch=	500; numPitches=50;  
rmeddis@38:     minPitch=	200; maxPitch=	4000; numPitches=40;  
rmeddis@38:     maxLag=1/minPitch; minLag=1/maxPitch;
rmeddis@38:     lags= linspace(minLag, maxLag, numPitches);
rmeddis@38:     
rmeddis@38: paramChanges={...
rmeddis@38:     'filteredSACFParams.lags=lags;     % autocorrelation lags vector;',...
rmeddis@38:     'filteredSACFParams.acfTau=	2;     % (Wiegrebe) time constant ACF;',...
rmeddis@38:     'filteredSACFParams.lambda=	0.12;  % slower filter to smooth ACF;',...
rmeddis@38:     'filteredSACFParams.plotACFs=1;    % plot ACFs while computing;',...
rmeddis@38:     'filteredSACFParams.plotACFsInterval=0.01;',...
rmeddis@38:     'filteredSACFParams.plotMoviePauses=.1;  ',...
rmeddis@38:     'filteredSACFParams.usePressnitzer=0; % attenuates ACF at  long lags;',...
rmeddis@38:     'filteredSACFParams.lagsProcedure=  ''useAllLags'';',...
rmeddis@38:     };
rmeddis@38: 
rmeddis@38: % Notes:
rmeddis@38: % acfTau:   time constant of unsmoothed ACF
rmeddis@38: % lambda:   time constant of smoothed ACFS
rmeddis@38: % plotACFs: plot ACFs during computation (0 to switch off, for speed)
rmeddis@38: % plotACFsInterval: sampling interval for plots
rmeddis@38: % plotMoviePauses:  pause duration between frames to allow viewing
rmeddis@38: % usePressnitzer:   gives low weights to long lags
rmeddis@38: % lagsProcedure:    used to fiddle with output (ignore)
rmeddis@38: 
rmeddis@38: %% delare 'showMap' options to control graphical output
rmeddis@38: % see UTIL_showMAP for more options
rmeddis@38: showMapOptions=[];
rmeddis@38: % showMapOptions.showModelOutput=0;     % plot of all stages
rmeddis@38: showMapOptions.surfAN=1;                % surface plot of HSR response
rmeddis@38: showMapOptions.PSTHbinwidth=0.001;      % smoothing for PSTH
rmeddis@38: 
rmeddis@38: if exist('fileName','var')
rmeddis@38:     % needed for labeling plot
rmeddis@38:     showMapOptions.fileName=fileName;
rmeddis@38: end
rmeddis@38: 
rmeddis@38: %% Generate stimuli
rmeddis@38: switch signalType
rmeddis@38:     case 'tones'
rmeddis@38:         % Create 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:     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: end
rmeddis@38: 
rmeddis@38: wavplay(inputSignal, sampleRate)
rmeddis@38: 
rmeddis@38: %% run the model
rmeddis@38: 
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: % display the AN response
rmeddis@38: UTIL_showMAP(showMapOptions)
rmeddis@38: 
rmeddis@38: % compute ACF
rmeddis@38: switch AN_spikesOrProbability
rmeddis@38:     case 'probability'
rmeddis@38:         % use only HSR fibers
rmeddis@38:         inputToACF=ANprobRateOutput(end-length(savedBFlist)+1:end,:);
rmeddis@38:     otherwise
rmeddis@38:         inputToACF=ANoutput;
rmeddis@38:         dt=dtSpikes;
rmeddis@38: end
rmeddis@38: 
rmeddis@38: disp ('computing ACF...')
rmeddis@38: 
rmeddis@38: % read paramChanges to get new filteredSACFParams
rmeddis@38: for i=1:length(paramChanges)
rmeddis@38: eval(paramChanges{i});
rmeddis@38: end
rmeddis@38: 
rmeddis@38: [P BFlist SACF]= filteredSACF(inputToACF, dt, savedBFlist, filteredSACFParams);
rmeddis@38: disp(' ACF done.')
rmeddis@38: 
rmeddis@38: %% plot original waveform on summary/smoothed ACF plot
rmeddis@38: figure(96), clf
rmeddis@38: subplot(3,1,3)
rmeddis@38: t=dt*(1:length(savedInputSignal));
rmeddis@38: plot(t,savedInputSignal, 'k')
rmeddis@38: xlim([0 t(end)])
rmeddis@38: title(['stimulus: ' num2str(leveldBSPL, '%4.0f') ' dB SPL']);
rmeddis@38: 
rmeddis@38: % plot SACF
rmeddis@38: figure(96)
rmeddis@38: subplot(2,1,1)
rmeddis@38: imagesc(P.^2)
rmeddis@38: colormap bone
rmeddis@38: ylabel('periodicities (Hz)'), xlabel('time (s)')
rmeddis@38: title('smoothed SACF. (periodicity x time)')
rmeddis@38: % y-axis specifies pitches (1/lags)
rmeddis@38: % Force MATLAB to show the lowest pitch
rmeddis@38: postedYvalues=[1 get(gca,'ytick')]; set(gca,'ytick',postedYvalues)
rmeddis@38: pitches=1./filteredSACFParams.lags;
rmeddis@38: set(gca,'ytickLabel', round(pitches(postedYvalues)))
rmeddis@38: % x-axis is time at which P is samples
rmeddis@38: [nCH nTimes]=size(P);
rmeddis@38: t=dt:dt:dt*nTimes;
rmeddis@38: tt=get(gca,'xtick');
rmeddis@38: set(gca,'xtickLabel', round(100*t(tt))/100)
rmeddis@38: 
rmeddis@38: %% On a new figure show a cascade of SACFs
rmeddis@38: figure(89), clf
rmeddis@38: % select 100 samples;
rmeddis@38: [r c]=size(SACF);
rmeddis@38: step=round(c/100);
rmeddis@38: idx=step:step:c;
rmeddis@38: 
rmeddis@38: UTIL_cascadePlot(SACF(:,idx)', 1./pitches)
rmeddis@38: 
rmeddis@38: xlabel('lag (s)'), ylabel('time pointer -->')
rmeddis@38: title(' SACF summary over time')
rmeddis@38: yValues=get(gca,'yTick');
rmeddis@38: set(gca,'yTickLabel', num2str(yValues'*100))
rmeddis@38: 
rmeddis@38: path(restorePath)
rmeddis@38: