preprocessing.Rmd

---
title: "Preprocessing"
date: "Last updated: `r format(Sys.time(), '%B %d, %Y')`"
bibliography: bibliography.bib
csl: american-medical-association.csl
link-citations: yes
output:
  html_document:
    toc: TRUE
    toc_depth: 2
    toc_float:
      collapsed: FALSE
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

```{r, child = "js/back-to-top.js"}
```

```{css, echo = FALSE}
table {
  margin: auto;
}

table thead th {
  border-bottom: 1px solid #ddd;
}

th, td {
  padding: 5px;
}
```

<br>

Listed below are the preprocessing routines that are run in `GannetLoad.m`. Whether a particular routine or subroutine is run will depend on the format of the inputted data and the options set in `GannetPreInitialise.m`.

## RF coil combination

Certain raw MRS data formats store data without coil combination; specifically, GE P-file (.7), NIfTI-MRS (if the source data were raw), Philips .raw, and Siemens TWIX (.dat) data. Gannet uses generalized least squares [@An2013] to optimally combine the signal from the multiple RF channels. If water files are provided, these data will be used as references for signal weighting and phasing of the coil data.

## Eddy-current correction

::: info
<i class="fa fa-info-circle" style="color: white"></i>&nbsp; Eddy-current correction can only be applied if water reference data are provided.
:::

In `GannetPreInitialise.m`, users have the option to apply eddy-current correction (ECC) to metabolite and water data. If applied, Gannet uses the method described by Klose (1990) [@Klose1990]. The code for the ECC routine can be found in `EddyCurrentCorrection.m`.

## Phase correction

It is common for unprocessed spectra to be out of phase. Gannet applies a global zero-order phase correction to all transients by fitting the real-valued 3 ppm Cr and 3.2 ppm Cho signals in the frequency domain and correcting the phase to be 0&deg; phase (i.e., to make all the Cr and Cho peaks have positive phase). This is performed in `PhaseCorrection.m`.

## Line-broadending (apodization)

FID data are multiplied by a time-varying exponential weighting function where the weighting constant is set in `GannetPreInitialise.m` (3 Hz is the default).

## Zero-filling

Gannet zero-pads all raw FIDs (i.e., adds zeros to the end of each FID) to obtain a nominal spectral resolution (the resolution between each frequency-domain data point) of 0.061 Hz/point. Differences in spectral width and number of complex data points of the raw FIDs are accounted for to obtain this nominal spectral resolution.

## Frequency and phase alignment

During acquisition, spectral data are affected by frequency and phase offsets as a result of biophysical, electronic, and participant factors. Gannet has several algorithms to correct for these errors during preprocessing. Users can choose which method to use in `GannetPreInitialise.m`:

| <u>Method</u> | <u>Option</u> | <u>Description</u> | <u>Function</u> |
| :---- | :---- | :--------- | :---- |
| Spectral registration [@Near2015] | `SpecReg` | A time-domain-based alignment method that uses nonlinear least-squares optimization to align each individual transient to a reference transient. | `SpectralRegistration.m` |
| Multi-step frequency and phase correction [@Mikkelsen2018] | `SpecRegHERMES` | A method originally developed to align multiplexed edited HERMES data. This approach is based on spectral registration. | `SpectralRegistrationHERMES.m` |
| Robust spectral registration [@Mikkelsen2020] (the default and recommended) | `RobustSpecReg` | A method based on spectral registration that is robust against spectral distortions caused by unstable residual water peaks and lipid contamination. | `RobustSpectralRegistration.m` | 
| Peak alignment | `Cr`, `Cho`, `NAA`, `H2O` | Frequency-domain-based alignment using one of the following peaks in each transient as the target signal: Cr, Cho, NAA, or residual H2O. | `AlignUsingPeak.m`, `AlignUsingH2O.m` |
| No alignment | `none` | Do not perform any frequency or phase alignment between transients. | n/a |

## Signal averaging

Gannet provides two methods for averaging individual transients (selected in `GannetPreInitialise.m`): arithmetic averaging (with outlier rejection) and weighted averaging (the default). The code for the signal averaging routines can be found in `SignalAveraging.m`.

### Arithmetic averaging

Arithmetic averaging is straightforward. All sequentially acquired $n$ pairs of subspectra $x_i$ (e.g., all edit-ON and edit-OFF subspectra) are averaged using the arithmetic mean: $\bar{x} = \frac{1}{n}\sum_{i=1}^nx_i$.

Note that before the arithmetic averaging of subspectra, individual transients are excluded based on the outlier rejection algorithm used during frequency and phase alignment.

### Weighted averaging

Weighted averaging down-weights individual difference subspectra that are corrupted by signal artifacts — this is an important distinction from traditional signal averaging. First, the difference between sequentially acquired pairs (e.g., all edit-ON and edit-OFF subspectra) is calculated. A similarity matrix $\mathbf{D}\in\mathbb{R}^{P{\times}P}$ is obtained by calculating the mean squared error between each real-valued difference subspectrum $p$ and every other real-valued difference subspectrum (in the range 1.8 to 3.4 ppm). A similarity metric $d_{p}$ is calculated as the column-wise median of $\mathbf{D}$. Normalized weights $w_{p}$ are then derived, $w_{p} = d^{-2}_p/\sum{d^{-2}_p}$, and applied to the difference pairs before summation.

Other algorithms to calculate weights can be found in `SignalAveraging.m`. (Note: These have not been tested extensively and should be considered experimental.)

## Residual water removal

If `water_removal` is set to 1 in `GannetPreInitialise.m` (the default), the residual water peak is removed from all difference spectra using an HSVD filter [@Barkhuijsen1987].

<br>

### References