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 350420 trials available for analysis per individual after artifact screening, EEG segments were divided into 50 bins of 15.625ms each, ranging from 0ms (stimulus onset) to 781.25ms. HP segments were divided into 10 bins of 500ms, ranging from 05s. Next, CECTs were computed by correlating (Pearsons 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 dataheartchannel = 62;  % channel # for heart databrainbin = 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 databrain = 2500;       % max point (from zero, in msec) for brain datafirst = -500;       % first point (from zero, in msec)yy = 1;             % counter for looping over brain datafor x = find(times>first,1):find(times>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>first,1):find(times>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;endHow can I make this code less awful (faster, more structured)?
Less For, more MATLABesque?
performance;matrix;matlab
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 alln=800;data = randn(2,n,30);times=linspace(-1900,8600,n);tic;samplingrate = 200;scalefac = 1000/samplingrate;brainchannel = 12;  % channel # for brain dataheartchannel = 62;  % channel # for heart databrainchannel = 1;heartchannel = 2; %adjusted for test databrainbin = 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 databrain = 2500;       % max point (from zero, in msec) for brain datafirst = -500;       % first point (from zero, in msec)xvec=find(times>first,1):find(times>brain,1);yvec=find(times>first,1):find(times>heart,1);CECT=zeros(length(xvec),length(yvec));  %initializebraindata=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 dataendyy = 1;             % counter for looping over brain datafor 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;endtoc