Beamforming across trials not time

[Robert Oostenveld]

Dear Marie,

On 26 Feb 2007, at 18:49, Marie Smith wrote:
> In order to achieve a better temporal resolution with a beamformer
> analysis it has been suggested to use single trial data in place of
> single time point data. Eg. In an experiment with 1000trials,
> 400time points - rather than compute the covariance across the 400
> time points (averaged over trials), compute it over 1000 points
> (for a very short fixed time interval e.g. 10ms). I have tried to
> implement this and would like some advice on whether or not the
> steps i have taken make sense in terms of fieldtrip functions and
> in general. Any advice would be welcome.

I suspect that what you have implemented is largely supported by the
existing code (see below). Regardless of that, it would be
interestnig if you would share your experiences on this with us.

> Effectively what I am doing is:
>
> Step 1) Choose a time point of interest
> Step 2) Smooth the single trial data with a 30point moving average
> to increase SNR
> Step 3) Generate a fake "single trial" which is the concatenation
> of each of these points across all trials e.g. the matrix has dims
> (nchans, (ntrials*ntimes));
> Step 4) Perform time-lock analysis on this "single trial" to
> compute the covariance  - avg1
> Step 5) Perform time-lock analysis on the actual data over the same
> short time range with keeptrials = 'yes'  -avg2
> Step 6) Perform lcmv sourceanalysis on avg1 (output of step 4)
> Step 7) Apply this filter to avg2 (output of step 5 and
> sourceanalysis) to estimate the single trial time courses at each
> grid location, across the short time interval.

The timelockanalysis function computes covariances (and also
averages, but the averaging is more or less independent from the
covariance computation. Timelockanalysis does the following:

a- it filters the data in each trial (optional, see private/preproc),
one of the options is also to apply a boxcar smoothing kernel
b- it selects the timewindow for averaging the ERF and computing the
single sample variance over trials
c- it selects the timewindow for covariance computation
d- for each trial, it selects the avg-timewindow and sums the data
e- for each trial, it selects the cov-timewindow and computes the
covariance within that timewindow
f- in case of keeptrials=yes, the single trial covariance (and the
data trials) are returned
g- in case of keeptrials=no, the covariance matrices are averaged
over trials

When beaming the data in sourceanalysis, the single trial covariance
matrices are averaged over trials. This is a bit messy in the
documentation but also in the code (e.g. the sourceanalysis options
keepfilters, rawtrial and singletrial interact with each other)

Regarding your 2 -> that can be done with step "a". The smoothing can
be done with a boxcar, and the selection of samples of interest is
done in our step "e".

Regarding your 3 and 4 -> averaging covariance matrixes (step "g" or
alternatively in sourceanalysis) is the same as your steps:
consider
   for i=1:40
     a{i} = randn(151,1000);
   end
   for i=1:10
     ca{i} = cov(a{i});
   end
   % concatenate trials
   b  = cat(2, a{:});
   cb = cov(b);
in this case the average over ca is the same as cb.

However, I now realise what the main difference appears to be between
your approach and the obvious use of the existing code: it relates to
baseline correction. Covariance in general is computed like this
   a = randn(151, 1000);
   for i=1:151
     a(i,:) = a(i,:) - mean(a(i,:));  % ensure zero mean
   end
   c = a * a';     % compute covariance
   c = c ./ 1000;  % divide by number of samples
That means that in the current implementation, the data within each
trial is baseline subtracted (over the whole timewindow) prior to
computing the covariance within that trial. It is possible to specify
in timelockanalysis removemean=no, but that does not achieve exactly
the same as your implementation, since your implementation removes
the mean OVER trials and then computes the covariance by multiplying
the matrix with its transpose.

So part of your method should already be possible with the existing
code (although admittedly not well documented). But your different
approach in baseline correcting is not possible with existing code.
Please have a look at the timelockanalysis code together with the
suggestions above and let me know whether you concur. Of course we
could extend timelockanalysis with your approach.

best regards,
Robert