Tutorial
The following tutorial provides an overview of the major steps of the eBOSC code. There are two files in the eBOSC repository that you can use to get an impression of the workings of the code.
- eBOSC_example_empirical: 1 unsystematic eyes closed resting state dataset from a young adult
- eBOSC_example_simulation: pink noise with systematic additions of 10 Hz rhythms of varying duration and amplitude
The tutorial has been written with the empirical example in mind, but the major steps of the eBOSC code remain the same (and are contained in the eBOSC_wrapper function) and are the focus here. Specifics of the example code can be found in the comments of the example scripts.
The eBOSC code roughly entails 5 steps:
- Time-to-frequency domain transform (wavelet-based)
- Define static background, derive power threshold and (optional) duration threshold
- Apply thresholds on continuous signals --> binary 'detected' matrix
- (Optional) continuous episode creation via post-processing of detected matrix
- (Optional) refinement of continuous episodes
Step 4 is required to produce a table of individual rhythmic episodes.
| Parameter | Description |
|---|---|
| cfg.eBOSC.F | to-be-sampled frequencies |
| cfg.eBOSC.wavenumber | wavelet family parameter |
| cfg.eBOSC.fsample | sampling frequency of input data |
The following entries concern the power and duration thresholds that are applied to detect the presence of a rhythm in time.
| Parameter | Description |
|---|---|
| cfg.eBOSC.threshold.excludePeak | lower and upper bound of frequencies to be excluded during background fit (Hz); multiple ranges can be encoded as separate rows |
| cfg.eBOSC.threshold.duration | vector of duration thresholds at each frequency |
| cfg.eBOSC.threshold.percentile | percentile of background fit for power threshold |
Note that you can set the duration parameter to zero if you wish to characterize the duration of rhythmic episodes afterwards.
The following entries concern signal padding to exclude edge artifacts and allow for the duration threshold. All padding is applied prior and after the signal of interest. Note that this is important for the calculation. Input data should be exactly padded according to the entries below such that after padding removal, only signal of interest remains.
| Parameter | Description |
|---|---|
| cfg.eBOSC.pad.tfr_s | padding following wavelet transform to avoid edge artifacts in seconds |
| cfg.eBOSC.pad.detection_s | padding following rhythm detection in seconds |
| cfg.eBOSC.pad.total_s | complete padding (WL + shoulder) |
| cfg.eBOSC.pad.background_s | padding of segments for background; default: same as used for TFR |
The following entries describe options for post-processing. If cfg.eBOSC.postproc.use is set to no, all other entries are irrelevant. These options mainly attempt to reduce temporal smearing due to the wavelet, although they should be used with care and with sanity checks (e.g., visual inspection of detection results).
| Parameter | Description |
|---|---|
| cfg.eBOSC.postproc.use | Post-processing of rhythmic eBOSC.episodes, i.e., wavelet 'deconvolution' (default = 'no') |
| cfg.eBOSC.postproc.method | Deconvolution method (default = 'MaxBias', FWHM: 'FWHM') |
| cfg.eBOSC.postproc.edgeOnly | Deconvolution only at on- and offsets of eBOSC.episodes? (default = 'yes') |
| cfg.eBOSC.postproc.effSignal | Amplitude deconvolution on whole signal or signal above power threshold? (default = 'PT') |
The following section provides a breakdown of the individual steps involved in rhythm detection that can be replicated via the 'eBOSC_example' script.
After defining the parameters, we first perfom a wavelet transform of all individual trials, which will be (a) used for the background fit and (b) for detecting rhythms from a single trial.
TFR = [];
for indTrial = 1:eBOSC.Ntrial
TFR.trial{indTrial} = BOSC_tf(data.trial{indTrial}(e,:),cfg.eBOSC.F,cfg.eBOSC.fsample,cfg.eBOSC.wavenumber);
clear tmp_dat
end; clear indTrial
Next, we perform a robust linear fit of the background. To reduce the bias of any rhythmic peaks on the background fit, users can enter frequency ranges, within which a rhythmic peak will be excluded prior to the background fit.
[eBOSC, pt, dt] = eBOSC_getThresholds(cfg, TFR, eBOSC);
Running this produces some summary outputs, as well as the static power threshold and duration threshold (even if that is set to zero) at each frequency. At this point, it makes sense to visually assess the fit. Here is an example fit in the empirical example.
Now that the thresholds have been defined, we can go ahead and apply them on the continuous signal. At each frequency, we apply the power and duration threshold to check which time points exceed those criteria. Those time points get a label of '1' in the binary matrix 'detected'.
But before we do that, we exclude the padding from the wavelet transform to reduce potential edge artifacts from the transform.
TFR_ = TFR.trial{1,indTrial}(:,cfg.eBOSC.pad.tfr_sample+1:end-cfg.eBOSC.pad.tfr_sample);
Then we apply the thresholds:
detected = zeros(size(TFR_));
for f = 1:length(cfg.eBOSC.F)
detected(f,:) = BOSC_detect(TFR_(f,:),pt(f),dt(f),cfg.eBOSC.fsample);
end; clear f
If we are interested in this detected matrix, we afterwards have to additionally strip off padding that was used to allow for a given duration threshold (imagine that we have a segment that would be 3c long, but occurs right at the edge, rendering it shorter than the target cycle duration).
% remove padding for detection (matrix with padding required for refinement)
eBOSC.detected(indTrial,:,:) = detected(:,cfg.eBOSC.pad.detection_sample+1:end-cfg.eBOSC.pad.detection_sample);
curDetected = squeeze(eBOSC.detected(indTrial,:,:));
The original Pepisode index, described by Caplan et al. is now simply the average of detected time points (i.e., the relative duration that a rhythm has occured in the analyzed segement).
eBOSC.pepisode(indTrial,:) = mean(curDetected,2);
At this point, it also makes sense to sanity-check the detected segments, to see whether something is off. Here is how an example trace looks like in the current example.
The previous sections reproduce the functionality of the original BOSC algorithm, with the slight adjustment ot the background fit (i.e., allowing for the exclusion of rhythmic spectral peaks). A novelty of the eBOSC algorithm is the possibility to create continuous episodes as an optional step, which will result in a table of individual episodes (see) below.
% get rhythmic episodes
eBOSC.episodes = [];
[detected_ep,eBOSC.episodes] = eBOSC_episode_create(TFR_,eBOSC, cfg);
The main problem this routine is trying to solve is the identification of episode on- and offsets, which is not yet achieved by the binary 'detected' matrix.
Let's take a look at the individual steps that are executed in this function:
- The original 'detected' matrix is sparsified across frequencies. The aim is to determine the local frequency at which power is maximal for any detected time point. This sparsification will next allow us to form continuous rhythmic episodes.
detected = eBOSC_episode_sparsefreq(cfg, detected, TFR);
-
The following code section 'traces' continuous sections of the sparsified detected matrix to create individual rhythmic episodes. These are subsequently characterized within a table (see below for info on the table contents).
-
(Optional) eBOSC provides the choice of either of two post-processing steps on these individual episodes. The main idea here is to try and alleviate the impact of the wavelet smooting on the time-frequency estimates(i.e., intrinsic time-frequency tradeoff) on especially the on- and offset estimates of the episodes. Two choices are available via the parameter cfg.eBOSC.postproc.method: (a) FWHM; this choice calcualtes the full-width-at-half-maximum of the wavelet at each power value of the episode, and checks whether the empirical power values can simply be explained by the extension of the wavelet; and (b) the 'MaxBias' method, that simulated the empirical smoothing produced by the wavelet on adjacent time points. Both options can be run on the entire episode, or only on the edges via the parameter cfg.eBOSC.postproc.edgeOnly. Setting this parameter to yes is advised if longer episodes (esp. of low rhythmic SNR) are otherwise broken up into a lot of short episodes, that may even fall under the duration threshold. A final parameter, cfg.eBOSC.postproc.effSignal, defined whether the 'bias' computation is based on the 'raw' power estimates, or only the part that falls above the power threshold.
For Steps 1-3 see the following Figure, reproduced with permission from Kosciessa et al. (2020):
Panels C-D indicate the spectral sparsification and creation of continuous episodes (Steps 2-3). Panels E-F illustrate the post-processing in Step 3.
- Finally, the next function accounts for the bilateral padding that is used for detecting rhythms (see padding removal for 'detected' matrix above), now in the domain of rhythmic episodes.
[episodesTable] = eBOSC_episode_rm_shoulder(cfg,detected_new,episodesTable);
The final eBOSC.episodes table contains information regarding all detected rhythmic episodes.
| Parameter | Description |
|---|---|
| Trial | Trial number (# as it appears in the loop; may need to be converted to external trial ID if some trials were dropped during preproc) |
| Channel | Channel ID as provided in loop |
| FrequencyMean | mean frequency of the episode in Hz |
| DurationS | episode duration in seconds |
| DurationC | episode duration in cycles |
| AmplitudeMean | mean power of episode |
| Onset | onset of episode w.r.t. timing vector in input data sturcture |
| Offset | offset of episode w.r.t. timing vector in input data sturcture |
| Amplitude | cell containing wavelet-based power values for every time point in episode |
| Frequency | cell containing wavelet-based frequency values for every time point in episode |
| RowID | index of 2D row (frequency) index in the detected matrix (after accounting for all padding) |
| ColID | index of 2D column (time) index in the detected matrix (after accounting for all padding) |
| SNR | SNR index: power during episode/ power of background estimate [for each episode time point] |
| SNRMean | mean SNR of episode |
Currently, the code is designed to save as much information as possible directly within episodes. If this is too memory intensive, users may wish to delete select output fields before saving, as especially the cells in the table can take up a lot of memory if many episodes are detected.
With this, we conclude a quick tutorial on how rhythmic signals can be detected using the functions provided in the eBOSC toolbox. We hope that you could successfully follow threough. If there are any bugs in the code, feature requests etc. feel free to open up an issue on the GitHub page.
eBOSC allows for multiple implementation streams, depending on the required fidelity of output and efficiency. Here is a quick overview on common use cases.
extended BOSC, no episodes:
- Closely related to original BOSC implementation (Caplan et al.), but with reduced bias of the initial background fit if spectral peak is excluded
- Pepisode index
- Fast to compute, requires a priori duration threshold
extended BOSC, episodes
- Additional post-processing of original detection matrix
- ‘abundance’ index
- Creates individual episodes with table of characteristics
- Necessary for on- and offset identification
- Powerful w.r.t. indices, but time-consuming
- Different parameter variants and choices for optimization:
- FWHM
- Faster, but potentially less empirical benefit
- MaxBias
- Slow, simulation-based, better performance in simulations
- In general, these post-processing choices should be sanity-checked and should be considered as in-development.
- FWHM