Multitaper spectrogram - Prerau et al. (2017)

Hello there Brainstorm team, (sorry for the long post in advance)

For the past few days, I've been trying to generate a multitaper spectrogram (in the context of sleep EEG) using the method provided by Prerau et al. (2017). The authors of the article have thankfully shared the MATLAB code online that can generate their multitaper spectrogram.

I was able to successfully generate the below multitaper spectrogram (Fpz channel) using the codes provided by the authors. The spectrogram shows a participant during sleep, specifically stage N2 of NREM sleep, which is shown by the sigma band representing sleep spindles (EEG hallmark of stage N2 sleep).

The code works relatively well but I was wishing to implement this multitaper spectrogram within Brainstorm to make things easier for me. Then I realised that Brainstorm has the process "FieldTrip: ft_mtmconvol (Multitaper)" already. Upon inspection, this process is somewhat similar to the multitaper spectrogram provided by Prerau et al. (2017). In an ideal situation, both the multitaper spectrogram and the TF sigma peaks (TF sigma peaks are my research interests) would be implemented in Brainstorm and I really would like to ask the Brainstorm team to perhaps consider implementing these two codes but of course, I understand that this will be a long task and you probably have more important implementations to do.

Therefore I thought about comparing and contrasting "ft_mtmconvol (Multitaper)" versus the multitaper spectrogram as provided by Prerau et al. (2017) in the hopes that I can somehow closely replicate the above multitaper spectrogram but using "ft_mtmconvol (Multitaper)". To properly understand "ft_mtmconvol (Multitaper)" I have consulted these forum/tutorial posts first.

First of all there are some things that I don't really understand so I will begin with some questions.
The below is the correspondence with FieldTrip inputs that Francois included in the "ft_mtmconvol (Multitaper)" process

  • timeoi = time_window(1)+time_resolution/2 : time_step : time_window(2)-time_resolution/2-1/sfreq
  • tapsmofrq = frequencies / modulation_factor
  • timwin = repmat(time_resolution, 1, length(frequencies))

From the above statements:
a) what does time_window (1) and time_window (2) mean?
b) in this notation ": time_step :" what does the notation ":" indicate?
Is this saying "timeoi = time_window(1)+time_resolution/2 = time_step = time_window(2)-time_resolution/2-1/sfreq"?
Or is this saying "timeoi = time_window(1)+time_resolution/2 and time_step = time_window(2)-time_resolution/2-1/sfreq"?
The notation written seems a bit confusing for me.
c) what does "1/sfreq" mean?
d) in "tapsmofrq = frequencies / modulation_factor", what does "frequencies" mean exactly?
Is this the interval between the start and stop frequency? (i.e. stop minus start frequency)?

After playing with "ft_mtmconvol (Multitaper)" process I was able to generate the following TF decomposition from Fpz channel using the below inputs. If you compare Prerau et al. (2017)'s multitaper spectrogram versus "ft_mtmconvol (Multitaper)" TF decomposition they do look similar with respect to the overall TF activity pattern. However, it would be better if they looked even closer in similarity versus what I generated. But before I go any further, what are the units for the colour bar on the right side of the TF decomposition?

Therefore I did some further compare and contrasting...
a) Looking at MATLAB command window I discovered something interesting.
The command window states that only one taper was used. Unlike Prerau et al. (2017) FieldTrip does not let the user specify the number of tapers therefore I am uncertain as to how only one taper was used. Any thoughts would be appreciated.

b) The below is the raw code that I copied from Prerau et al. (2017)'s multitaper spectrogram MATLAB code which highlights all the inputs required to run the function.
If you compare the below code to the inputs required by "ft_mtmconvol (Multitaper)" one can make the following observations.

  • data and Fs are given by the data present - so this is of no concern
  • frequency_range (I think) is the FieldTrip equivalent of "frequencies (start:step:stop)" - except Prerau et al. (2017) does not specify the "step" frequency
  • taper_params are the inputs that FieldTrip does not require but Prerau et al. (2017) does.
  • window_params (I think) are the FieldTrip equivalent of "Time resolution" and "Time step".
  • FieldTrip does not have detrend_opt, min_NFFT and weighting of tapers.

As my goal is to replicate (as close as possible) the Prerau et al. (2017) multitaper spectrogram using "ft_mtmconvol (Multitaper)", my question here is what FieldTrip input do you think I should modify to make the FieldTrip TF decomposition as smooth as possible like Prerau et al. (2017)'s multitaper spectrogram?

For reference I used taper_params = [5 9]; and window_params = [6 0.25]. That is...

  • time-halfbandwidth product = 5
  • number of tapers = 9
  • window size (secs) = 6 s
  • step size (secs) = 0.25 s
%   Input:
%       data: <number of samples> x 1  vector - time series data-- required
%       Fs: double - sampling frequency in Hz  -- required
%       frequency_range: 1x2 vector - [<min frequency>, <max frequency>] (default: [0 nyquist])
%       taper_params: 1x2 vector - [<time-halfbandwidth product>, <number of tapers>] (default: [5 9])
%       window_params: 1x2 vector - [window size (seconds), step size (seconds)] (default: [5 1])
%       detrend_opt: string - detrend data window ('linear' (default), 'constant', 'off');
%       min_NFFT: double - minimum allowable NFFT size, adds zero padding for interpolation (closest 2^x) (default: 0)
%       weighting: string - weighting of tapers ('unity' (default), 'eigen', 'adapt');
%       plot_on: boolean to plot results (default: true)
%       verbose: boolean to display spectrogram properties (default: true)
1 Like

Hi @mcp0228, thank you for the detailed information.

In brief: timeoi is the time axis for the spectrogram.
timeoi(1) is the center of the first window, timeoi(2) of the second and so on...

You can refer to the current implementation in Brainstorm here:

In detatil:
a) These indicate the time window where the multitaper will be computed. Index (1) is the start and index (2) is the end. These are the values of the Time window field in the process GUI

b) This is a standard way to generate an array in Matlab, [start: step: stop], thus [1:2:7] = [1, 3, 5, 7]

In the spectrogram code you shared the window parameters are [window size (seconds), step size (seconds)], there are equivalent to the fields Time resolution and Time step fields in the process GUI. Which in the code are saved as mt.timeres and timestep respectively. Thus:

The time axis (timeoi) as an array with:

  • Start at StartTime plus half of a multitaper window,
    OPTIONS.TimeVector(1) + mt.timeres/2
  • With a step of window step,
  • Up to EndTime minus half of a multitaper window minus one sample,
    OPTIONS.TimeVector(end) - mt.timeres/2 - 1/fsample

c) 1/fsample, is the sampling period (time of one sample) as fsample is the sampling frequency.

d) mt.frequencies is an array with the frequencies in the GUI, indicated as start:step:stop

The implementation of multitaper is different on both approaches. This is why the input parameters do not match one to one. As such, there is not a set of parameters that will give you the exactly same answer in both implementations, but both can provide useful spectrograms.

The approach on Prerau (2017) follows the short-time Fourier transform approach, where the number of tappers is fixed, thus the frequency smooth changes for each frequency. While multitaper in FieldTrip, time-halfbandwidth is dependent on the frequency under analysis and the smoothing factor, thus the number of tapers is not the same for all frequencies, with the goal of control the frequency smooth across frequencies.

@Raymundo.Cassani thank you very much for these comments - appreciate your time!
Just following on from the above:

  1. What are the units for the colour bar on the right side of the FieldTrip TF decomposition representing power?
  2. If I wanted to make a custom process to implement Prerau (2017) - any ideas of how I might begin?
    I am guessing it will look similar to "process_ft_mtmconvol.m" what are your thoughts?

I've found the answer to this one.

  • Power (\muV^2/Hz)
  • Magnitude (\muV/sqrt(Hz))
  • Log power (dB)

What are your thoughts for the below question?

In your process, the output of the Prerau (2017) code, will need to be reorganized to match the recreate a the TimeFreq structure that is shared by all the time-frequency representations in Brainstorm.

You take process_ft_mtmconvol as inspiration for the input parameters. However internally calls bst_timefreq which is the code that handles the call to the FieldTrip code and organizes its output in the Brainstorm TimeFreq structure.