MAP

# ignore windows autosaves

*.asv

# ignore mac autosaves

*~

# ignore MSword

#*.doc

#*.docx

~*

# Mac OS dir file

*.DS_Store

# We probably should not redistribute 3rd party USB drivers

*/CedrusUSB/*

% example:

% <navigate to 'MAP1_14\MAP'>

%  [inputSignal FS] = wavread('../wavFileStore/twister_44kHz');

%  MAP1_14(inputSignal, FS, -1, 'Normal', 'probability', [])

%

%  BFlist is a vector of BFs but can be '-1' to allow MAPparams to choose

%  MAPparamsName='Normal';          % source of model parameters

%  e.g.  paramChanges{1}= 'DRNLParams.a=0;'; % disable OHCs

%  e.g.  paramchanges={};                    % no changes

global  savedParamChanges savedBFlist saveAN_spikesOrProbability ...

    saveMAPparamsName savedInputSignal dt dtSpikes ...

    OMEextEarPressure TMoutput OMEoutput DRNLoutput...

    IHC_cilia_output IHCrestingCiliaCond IHCrestingV...

    IHCoutput ANprobRateOutput ANoutput savePavailable saveNavailable ...

    ANtauCas CNtauGk CNoutput ICoutput ICmembraneOutput ICfiberTypeRates...

    MOCattenuation ARattenuation 

%   global dtSpikes ICoutput ANprobRateOutput

%  'spike' matrices using a lower sample rate (see dtSpikes).

%  They are represented with dt=1/sampleRate (not dtSpikes)

if nargin<1

    error(' MAP1_14 is not a script but a function that must be called')

end

if nargin<6

    paramChanges=[];

end

% Read parameters from MAPparams<***> file in 'parameterStore' folder

% Beware, 'BFlist=-1' is a legitimate argument for MAPparams<>

%  It means that the calling program allows MAPparams to specify the list

cmd=['method=MAPparams' MAPparamsName ...

    '(BFlist, sampleRate, 0, paramChanges);'];

eval(cmd);

BFlist=DRNLParams.nonlinCFs;

savedParamChanges=paramChanges;

ANspeedUpFactor=ceil(sampleRate/AN_IHCsynapseParams.spikesTargetSampleRate);

% ANspeedUpFactor=AN_IHCsynapseParams.ANspeedUpFactor;  % e.g.5 times

% spiking activity is computed at longer sampling intervals (dtSpikes)

%  so it has a smaller number of epochs per segment(see 'ANspeeUpFactor' above)

dtSpikes=dt*ANspeedUpFactor;

OMEhighPassHighCutOff=OMEParams.OMEstapesHPcutoff;

MOCrateThresholdProb=DRNLParams.MOCrateThresholdProb;

minMOCattenuation=10^(DRNLParams.minMOCattenuationdB/20);

a1=dt/DRNLParams.MOCtauProb-1; a0=1;

b0=1+ a1;

MOCfiltProb_b=b0; MOCfiltProb_a=[a0 a1];

% if DRNLParams.a>0

%     DRNLcompressionThreshold=10^((1/(1-DRNLParams.c))* ...

%     log10(DRNLParams.b/DRNLParams.a));

% else

%     DRNLcompressionThreshold=inf;

% end

% DRNLcompressionThreshold=DRNLParams.cTh;

% DRNLa2=DRNLParams.a2;

% DRNLb=DRNLParams.b;

ctBM=10e-9*10^(DRNLParams.ctBMdB/20);

CtS=ctBM/DRNLa;

% i.e. b= [1  dt/tc-1]  and a= dt/IHC_cilia_RPParams.tc

IHCGa= IHC_cilia_RPParams.Ga; % (leakage)

IHCGu0 = IHCGa+IHCGmax./(1+exp(IHCu0/IHCs0).*(1+exp(IHCu1/IHCs1)));

IHCrestingCiliaCond=IHCGu0;

IHCrestingV= (IHC_Gk*IHC_Ekp+IHCGu0*IHC_Et)/(IHCGu0+IHC_Gk);

ANtauCas= IHCpreSynapseParams.tauCa;

nANchannels= nANfiberTypes*nBFs;

synapticCa= zeros(nANchannels,segmentLength);

mICa=zeros(nANchannels,segmentLength);

%  (one tauCa per channel). 

%  Do not confuse with ANtauCas vector (one per fiber type)

tauCa=repmat(ANtauCas, nBFs,1);

tauCa=reshape(tauCa, nANchannels, 1);

% presynapse startup values (vectors, length:nANchannels)

nANfibers= nANchannels*nFibersPerChannel;

AN_refractory_period= AN_IHCsynapseParams.refractory_period;

PAN_ydt = repmat(AN_IHCsynapseParams.y*dt, nANchannels,1);

PAN_ldt = repmat(AN_IHCsynapseParams.l*dt, nANchannels,1);

PAN_xdt = repmat(AN_IHCsynapseParams.x*dt, nANchannels,1);

PAN_rdt = repmat(AN_IHCsynapseParams.r*dt, nANchannels,1);

ANprobability=zeros(nANchannels,segmentLength);

ANprobRateOutput=zeros(nANchannels,signalLength);

lengthAbsRefractoryP= round(AN_refractory_period/dt);

cumANnotFireProb=ones(nANchannels,signalLength);

savePavailableSeg=zeros(nANchannels,segmentLength);

savePavailable=zeros(nANchannels,signalLength);

% only one stream of available transmitter will be saved

saveNavailableSeg=zeros(1,reducedSegmentLength);

saveNavailable=zeros(1,reducedSignalLength);

lengthAbsRefractory= round(AN_refractory_period/dtSpikes);

AN_ydt= repmat(AN_IHCsynapseParams.y*dtSpikes, nANfibers,1);

AN_ldt= repmat(AN_IHCsynapseParams.l*dtSpikes, nANfibers,1);

AN_xdt= repmat(AN_IHCsynapseParams.x*dtSpikes, nANfibers,1);

AN_rdt= repmat(AN_IHCsynapseParams.r*dtSpikes, nANfibers,1);

% If there is more than one value, everything is replicated accordingly

ANavailableFibersPerChan=AN_IHCsynapseParams.numFibers;

ANfibersFanInToCN=MacGregorMultiParams.fibersPerNeuron;

CNtauGk=MacGregorMultiParams.tauGk; % row vector of CN types (by tauGk)

nCNtauGk=length(CNtauGk);

% the total number of 'channels' is now greater

% 'channel' is defined as collections of units with the same parameters

%  i.e. same BF, same ANtau, same CNtauGk

nCNchannels=nANchannels*nCNtauGk;

nCNneuronsPerChannel=MacGregorMultiParams.nNeuronsPerBF;

tauGk=repmat(CNtauGk, nCNneuronsPerChannel,1);

tauGk=reshape(tauGk,nCNneuronsPerChannel*nCNtauGk,1);

% Now the number of neurons has been increased

nCNneurons=nCNneuronsPerChannel*nCNchannels;

CNinputfiberLists=zeros(nANchannels*nCNneuronsPerChannel, ANfibersFanInToCN);

for ch=1:nANchannels

        for idx2=1:nCNtauGk

            fibersUsed=(ch-1)*ANavailableFibersPerChan + ...

                ceil(rand(1,ANfibersFanInToCN)* ANavailableFibersPerChan);

            CNinputfiberLists(unitNo,:)=fibersUsed;

            unitNo=unitNo+1;

        end

t=dtSpikes:dtSpikes:5*CNspikeToCurrentTau;

CNalphaFunction= (1 / ...

    CNspikeToCurrentTau)*t.*exp(-t /CNspikeToCurrentTau);

CNalphaFunction=CNalphaFunction*CNcurrentPerSpike;

% figure(98), plot(t,CNalphaFunction), xlim([0 .020]), xlabel('time (s)'), ylabel('I')

tauGk=repmat(tauGk,nANchannels,1);

nICcells=nANchannels*nCNtauGk;  % one cell per channel

CN_PSTH=zeros(nICcells ,reducedSegmentLength);

% ICspikeWidth=0.00015;   % this may need revisiting

% epochsPerSpike=round(ICspikeWidth/dtSpikes);

% if epochsPerSpike<1

%     error(['MacGregorMulti: sample rate too low to support ' ...

%         num2str(ICspikeWidth*1e6) '  microsec spikes']);

% end

t=dtSpikes:dtSpikes:3*ICspikeToCurrentTau;

ICoutput=false(nICcells,reducedSignalLength);

    inputPressureSegment=inputSignal...

    % plot(segmentTime,inputPressureSegment), title('signalSegment')

    y=inputPressureSegment;

            inputPressureSegment, OMEExtFilterBndry{n});

    inputPressureSegment=y;

    OMEextEarPressure(segmentStartPTR:segmentEndPTR)= inputPressureSegment;

        filter(TMdisp_b,TMdisp_a,inputPressureSegment, ...

                filter(GTlin_b(BFno,:), GTlin_a(BFno,:), linOutput, ...

                GTlinBdry{BFno,order});

        % apply MOC to nonlinear input function      

%         figure(88), plot(MOC)

        nonlinOutput=stapesDisplacement.* MOC;

        %       first gammatone filter (nonlin path)

                nonlinOutput, GTnonlinBdry1{BFno,order});

        % Nick's compression algorithm

        abs_x= abs(nonlinOutput);

        signs= sign(nonlinOutput);

        belowThreshold= abs_x<CtS;

        nonlinOutput(belowThreshold)= DRNLa *nonlinOutput(belowThreshold);

        aboveThreshold=~belowThreshold;

        nonlinOutput(aboveThreshold)= signs(aboveThreshold) *ctBM .* ...

            exp(DRNLc *log( DRNLa*abs_x(aboveThreshold)/ctBM ));

%         %    original   broken stick instantaneous compression

%         holdY=nonlinOutput;

%         abs_x = abs(nonlinOutput);

%         nonlinOutput=sign(nonlinOutput).*min(DRNLa*abs_x, DRNLb*abs_x.^DRNLc);

% 

%         %   new broken stick instantaneous compression

%         y= nonlinOutput.* DRNLa;  % linear section attenuation/gain.

%         % compress parts of the signal above the compression threshold

%         %         holdY=y;

%         abs_y = abs(y);

%         idx=find(abs_y>DRNLcompressionThreshold);

%         if ~isempty(idx)>0

%             %             y(idx)=sign(y(idx)).* (DRNLcompressionThreshold + ...

%             %                 (abs_y(idx)-DRNLcompressionThreshold).^DRNLc);

%             y(idx)=sign(y(idx)).* (DRNLcompressionThreshold + ...

%                 (abs_y(idx)-DRNLcompressionThreshold)*DRNLa2);

%         end

%         nonlinOutput=y;

% if segmentStartPTR==10*segmentLength+1

%     figure(90)

%     plot(holdY,'b'), hold on

%     plot(nonlinOutput, 'r'), hold off

%     ylim([-1e-5 1e-5])

%     pause(1)

% end

%       second filter removes distortion products

                filter(GTnonlin_b(BFno,:), GTnonlin_a(BFno,:), ...

                nonlinOutput, GTnonlinBdry2{BFno,order});

%         disp(num2str(max(linOutput)))

    %         title('DRNLresponse'); 

    G=IHCGmax./(1+exp(-(IHCciliaDisplacement-IHCu0)/IHCs0).*...

    Gu=G + IHCGa;

    % replicate IHC_RP for each fiber type to obtain the driving voltage

    % fraction of channels open - apply time membrane constant

    % calcium current

    % apply calcium channel time constant

    synapticCa=-synapticCa; % treat synapticCa as positive substance

    %% AN -------------------------------

            %% Apply refractory effect

               % the probability of a spike's occurring in the preceding

               %  refractory window: t= (tnow-refractory period) :dt: tnow

               %    pFired= 1 - II(1-p(t)),

               % we need a running account of cumProb=II(1-p(t)) in order

               % not to have to recompute this for each value of t

               %   cumProb(t)= cumProb(t-1)*(1-p(t))/(1-p(t-refracPeriod))

               %   cumProb(0)=0

               %   pFired(t)= 1-cumProb(t)

               % This gives the fraction of firing events that must be 

               %  discounted because of a firing event in the refractory

               %  period

               %   p(t)= ANprobOutput(t) * pFired(t)

               % where ANprobOutput is the uncorrected firing probability

               %  based on vesicle release rate

               % NB this covers only the absoute refractory period

               % not the relative refractory period. To approximate this it

               % is necessary to extend the refractory period by 50%

                for t = segmentStartPTR:segmentEndPTR;

                    if t>1

                    ANprobRateOutput(:,t)= ANprobRateOutput(:,t)...

                        .* cumANnotFireProb(:,t-1);

                    end

                    % add recent and remove distant probabilities

                    refrac=round(lengthAbsRefractoryP * 1.5);

                    if t>refrac

                        cumANnotFireProb(:,t)= cumANnotFireProb(:,t-1)...

                            .*(1-ANprobRateOutput(:,t)*dt)...

                            ./(1-ANprobRateOutput(:,t-refrac)*dt);

                    end

                end

%                 figure(88), plot(cumANnotFireProb'), title('cumNotFire')

%                 figure(89), plot(ANprobRateOutput'), title('ANprobRateOutput')

            %% Estimate acoustic reflex efferent effect:  0 < ARattenuation > 1

            [r c]=size(ANrate);

            if r>nBFs % Only if LSR fibers are computed

                ARAttSeg=mean(ANrate(1:nBFs,:),1); %LSR channels are 1:nBF

                % smooth

                [ARAttSeg, ARboundaryProb] = ...

                    filter(ARfilt_b, ARfilt_a, ARAttSeg, ARboundaryProb);

                ARAttSeg=ARAttSeg-ARrateThreshold;

                ARAttSeg(ARAttSeg<0)=0;   % prevent negative strengths

                ARattenuation(segmentStartPTR:segmentEndPTR)=...

                    (1-ARrateToAttenuationFactorProb.* ARAttSeg);

            end

            %             plot(ARattenuation)

            % MOC attenuation based on within-channel HSR fiber activity

            if rateToAttenuationFactorProb<0

                % negative factor implies a fixed attenuation

                MOCattenuation(:,segmentStartPTR:segmentEndPTR)= ...

                    ones(size(rates))* -rateToAttenuationFactorProb;

            else

                for idx=1:nBFs

                    [smoothedRates, MOCprobBoundary{idx}] = ...

                        filter(MOCfiltProb_b, MOCfiltProb_a, rates(idx,:), ...

                        MOCprobBoundary{idx});

                    smoothedRates=smoothedRates-MOCrateThresholdProb;

                    smoothedRates=max(smoothedRates, 0);

                    x=(1- smoothedRates* rateToAttenuationFactorProb);

                    MOCattenuation(idx,segmentStartPTR:segmentEndPTR)= x;

                end

            end

            MOCattenuation(MOCattenuation<minMOCattenuation)= minMOCattenuation;

            % plot(MOCattenuation)            

                releaseProb=vesicleReleaseRate(:,t)*dtSpikes;

                releaseProb=reshape(releaseProb, nFibersPerChannel*nANchannels,1);

%                 probabilities= 1-(1-releaseProb).^AN_available; % slow

%                 probabilities= 1-(1-AN_reprocess.*AN_xdt).^M_q; % slow

%                 probabilities= 1-(1-AN_ydt).^M_q; %slow

                saveNavailableSeg(:,ANtimeCount)=AN_available(end,:); % only last channel

            % plotInstructions.displaydt=dtSpikes;

            % monitor synapse contents (only sometimes used)

            saveNavailable(reducedSegmentPTR:reducedSegmentPTR+reducedSegmentLength-1)=...

                saveNavailableSeg;

            % input is from AN so dtSpikes is used throughout

            for ch=1:nCNchannels

                    % ANpsth has a bin width of dtSpikes

                    conv2(AN_PSTH(unitNo,:),CNalphaFunction);

%             disp(['sum(AN_PSTH)= ' num2str(sum(AN_PSTH(1,:)))])

%                 disp(['mean(CNcurrentInput)= ' num2str(mean(CNcurrentInput(1,:)))])

                    CNtimeSinceLastSpike=CNtimeSinceLastSpike-dtSpikes;

                        CNcurrentInput(:,t)/CN_cap+(...

                        CN_Gk/CN_cap).*(CN_Ek-CN_E))*dtSpikes;

                    dGk=-CN_Gk*dtSpikes./tauGk + CN_b*s;

                    dTh=-(CN_Th-CN_Th0)*dtSpikes/CN_tauTh + CN_c*s;

                E=zeros(1,reducedSegmentLength);

                Gk=zeros(1,reducedSegmentLength);

                ss=zeros(1,reducedSegmentLength);

                    % time of previous spike moves back in time

                    CNtimeSinceLastSpike=CNtimeSinceLastSpike-dtSpikes;

                    % action potential if E>threshold

                    %  allow time for s to reset between events

                    s=CN_E>CN_Th0 & CNtimeSinceLastSpike<0 ;  

                    ss(t)=s(1);

                    CNtimeSinceLastSpike(s)=0.0005; % 0.5 ms for sodium spike

                        CNcurrentInput(:,t)/CN_cap +...

                        (CN_Gk/CN_cap).*(CN_Ek-CN_E))*dtSpikes;

                    dGk=-CN_Gk*dtSpikes./tauGk +CN_b*s;

                    E(t)=CN_E(1);

                    Gk(t)=CN_Gk(1);

                    CNmembranePotential(:,t)=CN_E +s.*(CN_Eb-CN_E)+CN_Er;

%             disp(['CN_E= ' num2str(sum(CN_E(1,:)))])

%             disp(['CN_Gk= ' num2str(sum(CN_Gk(1,:)))])

%             disp(['CNmembranePotential= ' num2str(sum(CNmembranePotential(1,:)))])

%             plot(CNmembranePotential(1,:))

            CN_spikes=CNmembranePotential> -0.02;

%             disp(['CNspikesbefore= ' num2str(sum(sum(CN_spikes)))])

            % for each spike put a zero in the next epoch

%             disp(['CNspikes= ' num2str(sum(sum(CN_spikes)))])

            % plotInstructions.displaydt=dtSpikes;

                % combine CN neurons in same channel (BF x AN tau x CNtau) 

                %  i.e. same BF, same tauCa, same CNtau

                for idx=1:nCNneuronsPerChannel: ...

                        nCNneurons-nCNneuronsPerChannel+1;

                        conv2(CN_PSTH(ICneuronNo,:),  IC_CNalphaFunction);

                    % faster computation when threshold is stable (c==0)

                            (IC_Gk/IC_cap).*(IC_Ek-IC_E))*dtSpikes;

                        dGk=-IC_Gk*dtSpikes/IC_tauGk +IC_b*s;

                            .*(IC_Ek-IC_E))*dtSpikes;

                        dGk=-IC_Gk*dtSpikes/IC_tauGk +IC_b*s;

                        dTh=-(IC_Th-Th0)*dtSpikes/IC_tauTh +s*IC_c;

                %figure(2),plot(ICmembranePotential(2,:))

                    % currently only the last segment is saved

                        /(nCellsPerTau*dtSpikes);

                ICoutput(:,reducedSegmentPTR:shorterSegmentEndPTR)=ICspikes;

                % figure(3),plot(ICoutput(2,:))

                % store membrane output on original dt scale

                % do this for single channel models only 

                % and only for the HSR-driven IC cells

                if round(nICcells/length(ANtauCas))==1  % single channel

                    % select HSR driven cells

                    x= ICmembranePotential(length(ANtauCas),:);

                    % restore original dt

                    x= repmat(x, ANspeedUpFactor,1);

                        y=repmat(ICmembranePotential(end,:),...

                            ANspeedUpFactor,1);

                % figure(4),plot(ICmembraneOutput(2,:))

                % AR is based on LSR units. LSR channels are 1:nBF

                if nANfiberTypes>1  % use only if model is multi-fiber

                    ARAttSeg=mean(ICspikes(1:nBFs,:),1)/dtSpikes;

%                    ARAttSeg(ARAttSeg<0)=0;   % prevent negative strengths

                    % scale up to dt from dtSpikes

                    x= repmat(ARAttSeg, ANspeedUpFactor,1);

                    x= reshape(x,1,segmentLength);

                    % max 60 dB attenuation

                    ARattenuation(ARattenuation<0)=0.01;

                    % single fiber type; disable AR because no LSR fibers

                %  separate MOC effect for each BF

                %  there is only one unit per channel

                rates=ICspikes(HSRbegins:end,:)/dtSpikes;

                % figure(4),plot(rates(1,:))

                    % spont 'rates' is zero for IC

                    % expand timescale back to model dt from dtSpikes

                    % figure(4),plot(x)

                % max attenuation is 30 dB

                MOCattenuation(MOCattenuation<minMOCattenuation)=...

                    minMOCattenuation;

                % figure(4),plot(MOCattenuation)

                % plotInstructions.displaydt=dtSpikes;

path(restorePath)

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))

%

%  	dt:		the signal sampling interval in seconds

%   BFlist: channel BFs

% params: 	a list of parmeters to guide the computation:

%   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)

%   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)

%  	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.

global savedInputSignal

ACFboundaryValue=[];

[nChannels inputLength]= size(inputSignalMatrix);

% create ACF movie

    signalTime=dt:dt:dt*length(savedInputSignal);

params.plotACFsInterval=round(params.plotACFsInterval/dt);

        % no action required lagWeights set above

        error ('ACF: params.lagProcedure not recognised')

% Establish matrix of lag time constants

else                  % use Wiegrebe method: tau= 2*lag (for example)

% LP_SACF function lowpass filter decay (only one value needed)

LP_SACF=zeros(nLags,inputLength+1);   % LP_SACF must match segment length +1

SACF=zeros(nLags,inputLength);

method=[];  % legacy programming

    ACF=zeros(nChannels,nLags);

    % 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};

    % ACF is a continuously changing channels x lags matrix

    % 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

    % E.g. if max(lags) is .04 then this is when the ACf will begin (?).

    ACF= (repmat(inputSignalMatrix(:, timePointer), 1, nLags) .* ...

        inputSignalMatrix(:, timePointer-lagPointers)).*...

        lagWeights *dt + ACF.* acfDecayFactors;

    x=(mean(ACF,1)./acfTaus)';

    SACF(:,timeCounter)=x;

    % 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

        plot([time time], [a(1) a(1)+(a(2)-a(1))/4], 'r', 'linewidth', 5)

LP_SACF=LP_SACF(:,1:end-1);  % correction for speed up above

    [a nTimePoints]=size(LP_SACF);

    LP_SACF=LP_SACF.*PressnitzerWeights';

    SACF=SACF.*PressnitzerWeights';

% BoundaryValue is legacy programming used in segmented model (keep)

ACFboundaryValue{1}=ACF(:,end-nLags+1:end);

ACFboundaryValue{2}=LP_SACF(:,end);

ACFboundaryValue{3} = inputSignalMatrix(:, end-max(lagPointers)+1 : end);