Less For, more MATLABesque? I have written a piece of MATLAB code that (hopefully) calculates CECT. I have two matrices of psychophysiological measures - one for brain data, one for heart rate; both are of the format epochs * times (i.e., multiple epochs of 10-second long intervals). The supposed outcome is to calculate the correlation between EEG in a short window using narrow bins, and heart rate in a longer window using wider bins. The basic idea is that an early brain effect predicts a heart rate speedup a few seconds later. In the words of the author (from the paywalled paper in the link):


  To allow CECT computation across the 350–420 trials available for analysis per individual after artifact screening, EEG segments were divided into 50 bins of 15.625 ms each, ranging from 0 ms (stimulus onset) to 781.25 ms. HP segments were divided into 10 bins of 500 ms, ranging from 0–5 s. Next, CECTs were computed by correlating (Pearson’s r) each EEG bin with each HP bin across trials within subjects, resulting in a 10 × 50 cross-correlation CECT matrix.


My code differs from that a bit in that I calculate for overlapping sliding windows instead of non-overlapping consecutive bins. Also, my windows are slightly different. However, what I'm mostly concerned with is that I have the hunch my code is extremely suboptimal.

Currently, the procedure is that I For-loop over brain data bins, collect a vector of heart rate data using a For loop, do a correlation between the heart bin vector and the single brain bin, store the resulting vector of correlations, and proceed to the next vector.

Importantly, the search windows and the bin windows for heart and brain data are of different length (heart rate variations are slower than brain effects). 
Also, this will in the end be calculated for multiple subjects; I have omitted this dimension of the data.

% data is a matrix of format (channels * times * epochs)
% times is a matrix of time points corresponding to sampling points

% CECT is supposed to be a matrix of correlation coefficients,
% dimensions (braintime * hearttime)

scalefac = 1000/samplingrate;

brainchannel = 12;  % channel # for brain data
heartchannel = 62;  % channel # for heart data

brainbin = 10;      % bin length for brain data (symmetric, in msec)
heartbin = 250;     % bin length for heart data (= 500 msec window)

heart = 5500;       % max point (from zero, in msec) for heart data
brain = 2500;       % max point (from zero, in msec) for brain data
first = -500;       % first point (from zero, in msec)

yy = 1;             % counter for looping over brain data
for x = find(times&gt;first,1):find(times&gt;brain,1)
    t0 = (x-brainbin/scalefac);
    t1 = (x+brainbin/scalefac);
    eeg = mean(squeeze(data(brainchannel,t0:t1,:)),1);
    p=1;
    clear hps
    for y = find(times&gt;first,1):find(times&gt;heart,1)
        t0 = (y-heartbin/scalefac);
        t1 = (y+heartbin/scalefac);

        hps(p,:) = mean(squeeze(data(heartchannel,t0:t1,:)),1);
        p = p+1;    % counter for looping over heart data
    end
    CECT(yy,:) = corr(squeeze(eeg).',hps.');
    yy = yy+1;
end


How can I make this code less awful (faster, more structured)?
With this random data and a few tweeks I got the time reduced about 30 times. The inner loop can be taken outside which is the biggest difference.

clear all

n=800;
data = randn(2,n,30);
times=linspace(-1900,8600,n);

tic;

samplingrate = 200;
scalefac = 1000/samplingrate;

brainchannel = 12;  % channel # for brain data
heartchannel = 62;  % channel # for heart data

brainchannel = 1;
heartchannel = 2; %adjusted for test data

brainbin = 10;      % bin length for brain data (symmetric, in msec)
heartbin = 250;     % bin length for heart data (= 500 msec window)

heart = 5500;       % max point (from zero, in msec) for heart data
brain = 2500;       % max point (from zero, in msec) for brain data
first = -500;       % first point (from zero, in msec)


xvec=find(times&gt;first,1):find(times&gt;brain,1);
yvec=find(times&gt;first,1):find(times&gt;heart,1);

CECT=zeros(length(xvec),length(yvec));  %initialize

braindata=squeeze(data(brainchannel,(xvec(1)-heartbin/scalefac):(xvec(end)+brainbin/scalefac),:));
bn=length((xvec(1)-brainbin/scalefac):(xvec(1)+brainbin/scalefac));

heartdata=squeeze(data(heartchannel,(yvec(1)-heartbin/scalefac):(yvec(end)+heartbin/scalefac),:));
hn=length((yvec(1)-heartbin/scalefac):(yvec(1)+heartbin/scalefac));
hps=zeros(length(yvec),size(data,3));

p=1;
for y = yvec
    hps(p,:) = mean(heartdata(p:p-1+hn,:),1);
    p = p+1;    % counter for looping over heart data
end

yy = 1;             % counter for looping over brain data
for x = xvec
    eeg = mean(braindata(yy:yy-1+bn,:),1);
    CECT(yy,:) = corr(squeeze(eeg).',hps.');
    %%CECT(yy,:) = eeg * hps.'; %for testing
    yy = yy+1;
end

toc