The purpose of this page is just to serve as todo or scratch pad for the development project and to list and share some ideas.

The code development project mentioned on this page has been finished by now. Chances are that this page is considerably outdated and irrelevant. The notes here might not reflect the current state of the code, and you should not use this as serious documentation.

Implement support for CTF synthetic gradiometers

  1. Improve accuracy of source reconstructions for synthetic gradient data
  2. Be able to convert the data to and from synthetic gradients

The coefficients are stored in a site specific *.coef file (e.g. /opt/ctf/hardware/M016/M017_1706.coef). Those should not be used, instead they should be read from the res4 file (see CTF comments in the appendix).

The coefficients should be put into a Nchan X Nchan matrix, where Nchan = Nref+Nmeg. Assuming that the uncorrected data has been read in for all channels, the conversion of 0th order to 3rd order only requires a multiplication with this matrix.

TODO: This requires looking through the low-level code (ctf_read_res4.m).

There should be a matrix for

0th -> 1st  (T01)
0th -> 2st  (T02)
0th -> 3st  (T03)

From those, each matrix can be made by combining the matrices, i.e.

data1 = T01 * data0;
data2 = T02 * data0;
data3 = T03 * data0;

and hence

data0 =  inv(T01) * data1;
data0 =  inv(T02) * data2;
data0 =  inv(T03) * data3;

so for example

data1 = T01 * inv(T03) * data3;

TODO: This requires testing on real data using the CTF software as gold standard, i.e. use the newDs command to construct a 1st, 2nd and 3rd order dataset from an original 0th order dataset (only one trial required, e.g. average). Regarding Matlab code, fieldtrip/syntheticgradient.m performs the conversions of the data.

The gradiometer array is described by all coil (positions and orientations) and the weights (+1 or -1) to combine the coils into hardware channels. From the FAQ:

The gradiometer definition generally consists of multiple coils per channel, e.g. two coils for a 1st order gradiometer in which the orientation of the coils is opposite. Each coil is described separately and one large matrix (can be sparse) has to be given that defines how the forward computed field is combined over the coils to generate the output of each channel. The gradiometer definition constsis of the following fields

grad.pnt   % Mx3 matrix with the position of each coil
grad.ori   % Mx3 matrix with the orientation of each coil
grad.tra   % NxM matrix with the weight of each coil into each channel
grad.label % cell-array of length N with the label of each of the channels (magnetometers or gradiometers)

The low-level forward computation code computes the field on each coil for a dipole in the x, y, and z-direction. That field (Ncoils x 3) is multiplied with the grad.tra matrix. Hence, to compute the forward model for synthetic gradients, it is sufficient to multiply this with the T03 matrix, i.e.

grad.tra0 =       grad.tra; 
grad.tra1 = T01 * grad.tra; 
grad.tra2 = T02 * grad.tra; 
grad.tra3 = T03 * grad.tra; 

gives the linear combination of the simulated field at all coils to form the synthetic gradients at the channels.

TODO: This requires a smarter implementation for linking the local spheres (one per channel) to the single channels (each channel combines 2 coils in the original 0th order data, but much more coils in the synthetic 3rd order data). The meg_leadfield1() function requires one sphere per coil, i.e. the coils and spheres have to be linked. Currently that gives an (intentional) error when it detects a higher order synthetic gradiometer array (fieldtrip/private/prepare_vol_sens.m, line 248).

The res4 file should specify somewhere how the raw data has been written (i.e. with 0th, 1st, 2nd or 3rd order balancing). That information should be added to the dataset in Matlab memory, and the gradiometer structure (grad.tra) should be consistent.

When changing the balancing (e.g. using the fieldtrip/ft_denoise_synthetic.m function), the data structure should remain internally consistent and should describe the balancing.

TODO: think of a way of adding this information to the output of fieldtrip/preprocessing and related functions (e.g. have a field data.grad.order='G1BR').

The sensor and coefficient information should ALWAYS be read from
the dataset (i.e., '.res4' file). This will guarantee that the sensor and
coefficient information are the corrected ones at the time the data was
collected (or modified).

To iterate, never use the '.sens' and/or '.coef' files. There is no
guarantee that the current '.sens' and/or '.coef' files are appropriate for
existing dataset.

The sensor positions in the '.sens' file is ALWAYS relative to the dewar
coordinates. The sensor positions in the dataset's '.res4' file are
specified in both dewar coordinates and head coordinates. You should use the
one in head coordinates. (As noted above, DON'T USE THE '.sens' FILE).
Our coils are 1st order gradiometers. Synthetic 2nd and 3rd order
gradiometers are formed using information >from the references. See below for
details on synthetic gradiometers.
To simulate (forward compute) the 1st gradiometer, the order of summing over
each coil doesn't matter (if you follow the steps prescribed below).

Simulating data / forward solution computation
Each sensor is described by:
- one or more coils
- zero or more baselines
- signed gain

Each coil is described by:
- position of centre of coil (use the position relative to head coordinate system)
- number of turns (N)
- coil area (A)
- unit normal vector (p) (use the position relative to head coordinate system)
- baseline from previous coil (for coils other than the first)

To compute the "field" picked up by a sensor:
1. Compute and sum the flux picked up by each coil.
2. To convert from flux to "field", divide the total flux by the effective
area (N * A) of the sensing coil (i.e., the first coil).
3. Here is the confusing part.... Our internal polarity definition is
opposite of the general convention, thus to get the polarity of the
simulated data consistent with the measured data, you have to apply the
opposite polarity (of the signed gain) to the simulated data. Therefore if
the sensor's gain is positive, you have to invert the simulated data; if the
sensor's gain is negative, leave the simulated data alone.

Note: If you are working with data with 'Tesla', you don't need to use the
gain value, just it's sign (as per step 3 above).
To compute the flux picked up by a coil, take the projection of the field
vector onto the coil's unit normal vector. That is,
Flux(at coil i) = N[i] * A[i] * Dot(B,p[i]),
where B is the field vector of your simulated signal source at the center of
coil i.
N[i] = number of turns of coil i
A[i] = area of coil i
p[i] = unit vector normal to area of coil i.
Total flux(of sensor) = Sum{ flux[coil i], i = 1 to number of coils}.
If you are integrating over the coil area, the procedure is similar.

Reading and using coefficients
If directly reading the coefficients from the dataset's '.res4'
file (or from the '.coef' file too), then the coefficients are Phi0 data.
Thus if the data is in Tesla, then you will have to convert the coefficients 
to be relevant to Tesla as well. To do this, use the formula
CoefOfRefInTesla = CoefOfRefInPhi0 * gainOfRef / gainOfSensor
For example, if the 3rd gradient coefficients for sensor MLC11 are (for Phi0
BG1: cBG1
BG2: cBG2
BG3: cBG3
G11: cG11
and the gains of the references and sensors are:
BG1: gBG1
BG2: gBG2
MLC11: gMLC11
Then the coefficients for MLC11 corresponding to Tesla data are:
BG1: cBG1 * gBG1 / gMLC11
BG2: cBG2 * gBG2 / gMLC11
BG3: cBG3 * gBG3 / gMLC11
G11: cG11 * gG11 / gMLC11
G3OI coefficient for data in phi0 for MLC12:
G11 (G11-1105): 0.143393
G12 (G12-1105): -0.00166001
G13 (G13-1105): -0.129106
G22 (G22-1105): -0.136122
G23 (G23-1105): 0.0653427
P11 (P11-1105): 0.00202835
Q11 (Q11-1105): -0.00211531
Q13 (Q13-1105): 0.265981
R11 (R11-1105): -0.127228
R12 (R12-1105): -0.00188768
R13 (R13-1105): 0.13007
R22 (R22-1105): 0.0777415
R23 (R23-1105): 0.330706
G3OI coefficient for data in Tesla for MLC12:
G11 (G11-1105): -0.285225
G12 (G12-1105): 0.00313698
G13 (G13-1105): -0.241933
G22 (G22-1105): -0.270226
G23 (G23-1105): -0.125866
P11 (P11-1105): 0.00393464
Q11 (Q11-1105): -0.0039391
Q13 (Q13-1105): 0.477763
R11 (R11-1105): -0.238921
R12 (R12-1105): -0.00330739
R13 (R13-1105): -0.234992
R22 (R22-1105): -0.145796
R23 (R23-1105): 0.602411

The coefficients are applied by SUBTRACTING the linear combination of
reference signals from the sensor signal. E.g.,
MLC11(3rd) = MLC11(raw) - Sum{ ref[i] * coef[i], i=1 to num of references}

Ideal vs Real coefficients (G3OI vs G3BR).
The ideal coefficients are based on geometry of the sensor configuration.
The real coefficients accounts for geometrical and common mode errors; these
are determined experimentally.

FOR DipoleFit AND dfit VERSION 4.12 and up

We use the real coefficients for both processing real data and forward solution
computation. Recent studies at CTF showed that using the real coefficients for the
forward solution computation gave a more accurate match with the real data compared
to the using the ideal coefficients.