So, I know one of the major advantages of MATLAB is in using vector math instead of loops, but I'm having some issues figuring how to make a particular function run faster.

A little background (that can be skipped, if not necessary)

I collect neuronal data (spikes, or action potentials) from an auditory brain structure during the presentation of a sound. I'm attempting to calculate the latency of the neuronal response onset for each trial. Sometimes, it's not simply as easy as counting the first spike as there are spontaneous spikes that occur even in the absence of sound. The data get recorded as a time stamp. I have the timestamps of each stimulus presentation and I sort the spike timestamps relative to each stimulus presentation. From 100 ms before to 200 ms after the stimulus onset, I keep a running total of number of spikes occuring, which looks like this over about 100 trials (the function looks much more coarse for individual trials, but the concept is the same). The x-axis shows the time relative to stimulus onset (time 0). The y-axis shows the average number of spikes that have occured by time x in a given trial. In general, slope represents instantaneous firing rate. The shallow slope at the beginning represents a lower "spontaneous" firing rate. When the neuron becomes active due to the stimulus, the slope increases sharply.

There are plenty of ways that have been used to calculate response onset, but for various reasons, I've been left to kind of "invent" my own way. My idea is that the onset is represented here as the deflection from the shallow to steep slope. In order to automate this process, I decided to calculate the slope in a sliding window across this particular function. The point where the slope is maximal, should be the onset, like this. I have to calculate this slope many times during each trial. Each trial has data from 301 ms. I use a window of 100 ms and slide it every 1 ms, so I calculate it 201 times per trial. My blocks of data can have 500-600 trials. I exclude "catch trials" where there are no stimuli to save time. Still, just to analyze the latency for each block of trials takes about a minute. This is just too long considering all the other analysis I do takes time. I know this is because on some block, I'm calling polyfit() ~12,000 times. I'll be happy to answer any conceptual questions you might have.

So here's my code

I know this might be confusing to follow, but I've included the entire function for completeness.

function latency = findLatency(spkcum,cdfbin,stimtrial,window)

% spkcum = cumulative spike counts (trial number x peristimulus time)
% cdfbin = peristimulus time corresponding to dimension 2 of spkcum
% stimtrial = logical where 1 = stimulus presented
% window = length of stimulus and correlation/slope calculation window.

for ii = 1:size(spkcum,1) % Loops through each trial
    if stimtrial(ii) % Skips catch trials
        for jj = 1:size(spkcum,2)-window % Loops through each time point in a given trial ii
            tempS = polyfit(cdfbin(jj:jj window),spkcum(ii,jj:jj window),1); % Fits a line to the function within the window (jj:jj window)
            S(jj) = tempS(1); % extracts the slope, discards intercept
        end
        try % Skips errors on trials where there were no spikes
            ts = cdfbin(find(diff(spkcum(ii,:)))); % Calculates the spike times per trial
            [~,peakS] = max(S.*(cdfbin(1:length(S))>0)); % Finds peak of the slope AFTER stimulus onset
            peakT = cdfbin(round(peakS)); % Uses the index to calculate the time of the peak
            latency(ii) = min(ts(ts>peakT)); % finds the first spike after the peak of the slope
        catch
            latency(ii) = NaN; % No response onset calculated when no spikes
        end      
    else
        latency(ii) = NaN; % No response onset calculated when catch trial
    end
end

Here is the relevant portion of the function, which I think takes the most time:

for ii = 1:size(spkcum,1)
    for jj = 1:size(spkcum,2)-window
        tempS = polyfit(cdfbin(jj:jj window),spkcum(ii,jj:jj window),1);
        S(jj) = tempS(1); % extracts the slope, discards intercept
    end
end

Help!

First off, yes I use "cum" as part of a variable. One of my favorite functions is cumtrapz(). It's the little things that keep you sane in science.

Seeing the first graph I posted, one might think I could just use diff() or something similar. That graph is a very pretty curve because it averages 100 trials. Single trial cumulative spike counts are often not very pretty and the diff(), well, sucks. Using a sliding window to calculate slope has been my most effective way that I have come up with to calculate onset. I've considered calculating slope only every 2 or 5 ms to speed up the process, but I'm wondering if it is possible to get rid one or both of the loops entirely. Additionally, I'd be open to a different method entirely to generate the slope function.

If you need the data that I'm using, it can be found in a mat file here.

Thanks in advance! I'm happy to answer any question or clarify anything.