diff MAP/filteredSACF.m @ 38:c2204b18f4a2 tip

End nov big change

--- a/MAP/filteredSACF.m	Thu Oct 06 15:43:20 2011 +0100
+++ b/MAP/filteredSACF.m	Mon Nov 28 13:34:28 2011 +0000
@@ -1,53 +1,52 @@
-function [P, BFlist, sacf, boundaryValue] = ...
-    filteredSACF(inputSignalMatrix, method, params)
-% UTIL_filteredSACF computes within-channel, running autocorrelations (acfs)
-%  and finds the sum across channels (sacf).
-%  The SACF is smoothed to give the 'p function' (P).
+function [LP_SACF, BFlist, SACF, ACFboundaryValue] = ...
+    filteredSACF(inputSignalMatrix, dt, BFlist, params)
+% UTIL_filteredSACF computes within-channel, running autocorrelations (ACFs)
+%  and finds the running sum across channels (SACF).
+%  The SACF is smoothed to give the 'p function' (LP_SACF).
+%  (Balaguer-Ballestera, E. Denham, S.L. and Meddis, R. (2008))
+%
 %
 % INPUT
 %  	inputSignalMatrix: a matrix (channel x time) of AN activity
-%  	method.dt:			the signal sampling interval in seconds
-%   method.segmentNo:
-%   method.nonlinCF
-%   
+%  	dt:		the signal sampling interval in seconds
+%   BFlist: channel BFs
+%
 % params: 	a list of parmeters to guide the computation:
-%   filteredSACFParams.lags: an array of lags to be computed (seconds)
-%   filteredSACFParams.acfTau: time constant (sec) of the running acf calculations
-%     if acfTau>1 it is assumed that Wiegrebe'sacf method
+%   params.lags: an array of lags to be computed (seconds)
+%   params.acfTau: time constant (sec) of the running ACF
+%     if acfTau>1 it is assumed that Wiegrebe'SACF method
 %   	for calculating tau is to be used (see below)
-%   filteredSACFParams.Lambda: time constant  for smoothing thsacf to make P
-%   filteredSACFParams.lagsProcedure identifies a strategy for omitting some lags.
-%    Options are: 'useAllLags', 'omitShortLags', or 'useBernsteinLagWeights'
-%   filteredSACFParams.usePressnitzer applies lower weights longer lags
-%   parafilteredSACFParamsms.plotACFs (=1) creates movie of acf matrix (optional)
-%
+%   params.Lambda: smoothing factor for the SACF
+%   params.lagsProcedure: strategies for omitting some lags.
+%    (Options: 'useAllLags' or 'omitShortLags')
+%   params.usePressnitzer applies lower weights longer lags
+%   params.plotACFs (=1) creates movie of ACF matrix (optional)
 %
 % OUTPUT
-%  	P:	P function 	(time x lags), a low pass filtered version of sacf.
-%  method: updated version of input method (to include lags used)
-%   sacf:  	(time x lags)
+%  	LP_SACF:	LP_SACF function 	(time x lags), a low pass filtered version of SACF.
+%  method: updated version of input 'method' (to include lags used)
+%   SACF:  	(time x lags)
+%
+% Notes:
+% ACFboundaryValue refers to segmented evaluation and is currently not
+%  supported. However the code may be useful later when this function
+%  is incorporated into MAP1_14.
 
 %%
-boundaryValue=[];
+global savedInputSignal
+
+ACFboundaryValue=[];
+
 [nChannels inputLength]= size(inputSignalMatrix);
-% list of BFs must be repeated is many fiber types used
-BFlist=method.nonlinCF;
-nfibertypes=nChannels/length(BFlist);
-BFlist=repmat(BFlist',2,1)';
 
-dt=method.dt;
-% adjust sample rate, if required
-if isfield(params,'dt')
-    [inputSignalMatrix dt]=UTIL_adjustDT(params.dt, method.dt, inputSignalMatrix);
-    method.dt=dt;
-end
-
-% create acf movie
+% create ACF movie
 if isfield(params, 'plotACFs') && params.plotACFs==1
     plotACF=1;
+    signalTime=dt:dt:dt*length(savedInputSignal);
 else
     plotACF=0;  % default
 end
+params.plotACFsInterval=round(params.plotACFsInterval/dt);
 
 if isfield(params,'lags')
     lags=params.lags;
@@ -64,11 +63,11 @@
 else
     lagsProcedure='useAllLags';  % default
 end
-% disp(['lag procedure= ''' lagsProcedure ''''])
+
 lagWeights=ones(nChannels,nLags);
 switch lagsProcedure
     case 'useAllLags'
-       % no action required lagWeights set above
+        % no action required lagWeights set above
     case 'omitShortLags'
         % remove lags that are short relative to CF
         allLags=repmat(lags,nChannels,1);
@@ -76,14 +75,12 @@
         criterionForOmittingLags=1./(params.criterionForOmittingLags*allCFs);
         idx= allLags < criterionForOmittingLags;	% ignore these lags
         lagWeights(idx)=0;
-    case 'useBernsteinLagWeights'
-        lagWeights=BernsteinLags(BFlist, lags)';
     otherwise
-        error ('acf: params.lagProcedure not recognised')
+        error ('ACF: params.lagProcedure not recognised')
 end
 
 
-% Establish matrix of lag time constants 
+% Establish matrix of lag time constants
 %   these are all the same if tau<1
 % if acfTau>1, it is assumed that we are using the Wiegrebe method
 %   and a different decay factor is applied to each lag
@@ -92,7 +89,7 @@
 if acfTau<1          % all taus are the same
     acfTaus=repmat(acfTau, 1, nLags);
     acfDecayFactors=ones(size(lags)).*exp(-dt/acfTau);
-else                      % use Wiegrebe method: tau= 2*lag (for example)
+else                  % use Wiegrebe method: tau= 2*lag (for example)
     WiegrebeFactor=acfTau;
     acfTaus=WiegrebeFactor*lags;
     idx= acfTaus<0.0025; acfTaus(idx)=0.0025;
@@ -101,7 +98,7 @@
 % make acfDecayFactors into a (channels x lags) matrix for speedy computation
 acfDecayFactors=repmat(acfDecayFactors,nChannels, 1);
 
-% P function lowpass filter decay (only one value needed)
+% LP_SACF function lowpass filter decay (only one value needed)
 pDecayFactor=exp(-dt/params.lambda);
 
 % ACF
@@ -112,124 +109,93 @@
 end
 
 
-P=zeros(nLags,inputLength+1);   % P must match segment length +1
-sacf=zeros(nLags,inputLength);
+LP_SACF=zeros(nLags,inputLength+1);   % LP_SACF must match segment length +1
+SACF=zeros(nLags,inputLength);
+method=[];  % legacy programming
 if   ~isfield(method,'segmentNumber') || method.segmentNumber==1
-    acf=zeros(nChannels,nLags);
+    ACF=zeros(nChannels,nLags);
     % create a runup buffer of signal
     buffer= zeros(nChannels, max(lagPointers));
 else
-    % boundaryValue picks up from a previous calculation
-    acf=params.boundaryValue{1};
-    P(: , 1)=params.boundaryValue{2}; % NB first value is last value of previous segment
-    buffer=params.boundaryValue{3};
+    % ACFboundaryValue picks up from a previous calculation
+    ACF=params.ACFboundaryValue{1};
+    % NB first value is last value of previous segment
+    LP_SACF(: , 1)=params.ACFboundaryValue{2};
+    buffer=params.ACFboundaryValue{3};
 end
 inputSignalMatrix=[buffer inputSignalMatrix];
 [nChannels inputLength]= size(inputSignalMatrix);
 
 timeCounter=0; biggestSACF=0;
 for timePointer= max(lagPointers)+1:inputLength
-    % acf is a continuously changing channels x lags matrix
+    % ACF is a continuously changing channels x lags matrix
     %   Only the current value is stored
-    % sacf is the vertical summary of acf ( a vector) and all values are kept and returned
-    % P is the smoothed version of sacf and all values are kept and returned
+    % SACF is the vertical sum of ACFs; all values are kept and returned
+    % LP_SACF is the smoothed version of SACF:all values are kept and returned
     % lagWeights emphasise some BF/lag combinations and ignore others
     % NB time now begins at the longest lag.
-    % E.g. if max(lags) is .04 then this is when the ACf will begin.
-    %            AN                                       AN delayed                             weights               filtering
-    
+    % E.g. if max(lags) is .04 then this is when the ACf will begin (?).
+
     % This is the ACF calculation
     timeCounter=timeCounter+1;
-    acf= (repmat(inputSignalMatrix(:, timePointer), 1, nLags) .* inputSignalMatrix(:, timePointer-lagPointers)).*lagWeights *dt + acf.* acfDecayFactors;
-    x=(mean(acf,1)./acfTaus)';
-%     disp(num2str(x'))
-    sacf(:,timeCounter)=x;
-    P(:,timeCounter+1)=sacf(:,timeCounter)*(1-pDecayFactor)+P(:,timeCounter)*pDecayFactor;
+    ACF= (repmat(inputSignalMatrix(:, timePointer), 1, nLags) .* ...
+        inputSignalMatrix(:, timePointer-lagPointers)).*...
+        lagWeights *dt + ACF.* acfDecayFactors;
+    x=(mean(ACF,1)./acfTaus)';
+    SACF(:,timeCounter)=x;
     
-    % plot at intervals of 200 points
-    if plotACF && ~mod(timePointer,params.plotACFsInterval)
-        %       mark cursor on chart to signal progress
-        % this assumes that the user has already plotted
-        % the signal in subplot(2,1,1) of figure (13)
-        figure(13)
-        hold on
-        subplot(4,1,1)
+    % smoothed version of SACF
+    LP_SACF(:,timeCounter+1)=SACF(:,timeCounter)* (1-pDecayFactor) + ...
+        LP_SACF(:,timeCounter)* pDecayFactor;
+
+    % plot ACF at intervals if requested to do so
+    if plotACF && ~mod(timePointer,params.plotACFsInterval) && ...
+            timePointer*dt>3*max(lags)
+        figure(89), clf
+        %           plot ACFs one per channel
+        subplot(2,1,1)
+        UTIL_cascadePlot(ACF, lags)
+        title(['running ACF  at ' num2str(dt*timePointer,'%5.3f') ' s'])
+        ylabel('channel BF'), xlabel('period (lag, ms)')
+        set(gca,'ytickLabel',[])
+
+        %           plot SACF
+        subplot(4,1,3), cla
+        plotSACF=SACF(:,timeCounter)-min(SACF(:,timeCounter));
+        plot(lags*1000, plotSACF, 'k')
+        biggestSACF=max(biggestSACF, max(plotSACF));
+        if biggestSACF>0, ylim([0 biggestSACF]), else ylim([0 1]), end
+        ylabel('SACF'), set(gca,'ytickLabel',[])
+%         xlim([min(lags*1000) max(lags*1000)])
+
+        %           plot signal
+        subplot(4,1,4)
+        plot(signalTime, savedInputSignal, 'k'), hold on
+        xlim([0 max(signalTime)])
+        a= ylim;
+        %       mark cursor on chart to indicate progress
         time=timePointer*dt;
-        a =ylim;
-        plot([time time], [a(1) a(1)+(a(2)-a(1))/4]) % current signal point marker
-        
-        %         plot ACFs one per channel
-        subplot(2,1,2), cla
-        cascadePlot(acf, lags, BFlist)
-        xlim([min(lags) max(lags)])
-        %         set(gca,'xscale','log')
-        title(num2str(method.dt*timePointer))
-        ylabel('BF'), xlabel('period (lag)')
-        
-        %         plot SACF
-        subplot(4,1,2), hold off
-        plot(lags,sacf(:,timeCounter)-min(sacf(:,timeCounter)))
-        biggestSACF=max(biggestSACF, max(sacf(:,timeCounter)));
-        if biggestSACF>0, ylim([0 biggestSACF]), end
-        %         set(gca,'xscale','log')
-        title('SACF')
+        plot([time time], [a(1) a(1)+(a(2)-a(1))/4], 'r', 'linewidth', 5)
         pause(params.plotMoviePauses)
     end
 end
-P=P(:,1:end-1);  % correction for speed up above
+LP_SACF=LP_SACF(:,1:end-1);  % correction for speed up above
 
 % Pressnitzer weights
 if ~isfield(params, 'usePressnitzer'),     params.usePressnitzer=0; end
 if params.usePressnitzer
-    [a nTimePoints]=size(P);
+    [a nTimePoints]=size(LP_SACF);
     % higher pitches get higher weights
     %     PressnitzerWeights=repmat(min(lags)./lags,nTimePoints,1);
     % weaker weighting
     PressnitzerWeights=repmat((min(lags)./lags).^0.5, nTimePoints,1);
-    P=P.*PressnitzerWeights';
-    sacf=sacf.*PressnitzerWeights';
+    LP_SACF=LP_SACF.*PressnitzerWeights';
+    SACF=SACF.*PressnitzerWeights';
 end
 
 % wrap up
-method.acfLags=lags;
-method.filteredSACFdt=dt;
-
-boundaryValue{1}=acf(:,end-nLags+1:end);
-boundaryValue{2}=P(:,end);
+% BoundaryValue is legacy programming used in segmented model (keep)
+ACFboundaryValue{1}=ACF(:,end-nLags+1:end);
+ACFboundaryValue{2}=LP_SACF(:,end);
 % save signal buffer for next segment
-boundaryValue{3} = inputSignalMatrix(:, end-max(lagPointers)+1 : end);
-
-method.displaydt=method.filteredSACFdt;
-
-% if ~isfield(params, 'plotUnflteredSACF'), params.plotUnflteredSACF=0; end
-% if method.plotGraphs
-%     method.plotUnflteredSACF=params.plotUnflteredSACF;
-%     if ~method.plotUnflteredSACF
-%         method=filteredSACFPlot(P,method);
-%     else
-%         method=filteredSACFPlot(SACF,method);
-%     end
-% end
-
-
-% ------------------------------------------------  plotting ACFs
-function cascadePlot(toPlot, lags, BFs)
-% % useful code
-[nChannels nLags]=size(toPlot);
-
-% cunning code to represent channels as parallel lines
-[nRows nCols]=size(toPlot);
-if nChannels>1
-    % max(toPlot) defines the spacing between lines
-    a=max(max(toPlot))*(0:nRows-1)';
-    % a is the height to be added to each channel
-    peakGain=10;
-    % peakGain emphasises the peak height
-    x=peakGain*toPlot+repmat(a,1,nCols);
-    x=nRows*x/max(max(x));
-else
-    x=toPlot;                            % used when only the stimulus is returned
-end
-plot(lags, x','k')
-ylim([0 nRows])
-
+ACFboundaryValue{3} = inputSignalMatrix(:, end-max(lagPointers)+1 : end);