{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n\n# Frequency-tagging: Basic analysis of an SSVEP/vSSR dataset\n\nIn this tutorial we compute the frequency spectrum and quantify signal-to-noise\nratio (SNR) at a target frequency in EEG data recorded during fast periodic\nvisual stimulation (FPVS) at 12 Hz and 15 Hz in different trials.\nExtracting SNR at stimulation frequency is a simple way to quantify frequency\ntagged responses in MEEG (a.k.a. steady state visually evoked potentials,\nSSVEP, or visual steady-state responses, vSSR in the visual domain,\nor auditory steady-state responses, ASSR in the auditory domain).\n\nFor a general introduction to the method see\n[Norcia et al. (2015)](https://doi.org/10.1167/15.6.4) for the visual domain,\nand [Picton et al. (2003)](https://doi.org/10.3109/14992020309101316) for\nthe auditory domain.\n\n**Data and outline:**\n\nWe use a simple example dataset with frequency tagged visual stimulation:\nN=2 participants observed checkerboard patterns inverting with a constant\nfrequency of either 12.0 Hz of 15.0 Hz.\n32 channels wet EEG was recorded.\n(see `ssvep-dataset` for more information).\n\nWe will visualize both the power-spectral density (PSD) and the SNR\nspectrum of the epoched data,\nextract SNR at stimulation frequency,\nplot the topography of the response,\nand statistically separate 12 Hz and 15 Hz responses in the different trials.\nSince the evoked response is mainly generated in early visual areas of the\nbrain the statistical analysis will be carried out on an occipital\nROI.\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# Authors: Dominik Welke \n# Evgenii Kalenkovich \n#\n# License: BSD-3-Clause\n# Copyright the MNE-Python contributors."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"import matplotlib.pyplot as plt\nimport numpy as np\nfrom scipy.stats import ttest_rel\n\nimport mne"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Data preprocessing\nDue to a generally high SNR in SSVEP/vSSR, typical preprocessing steps\nare considered optional. This doesn't mean, that a proper cleaning would not\nincrease your signal quality!\n\n* Raw data have FCz reference, so we will apply common-average rereferencing.\n\n* We will apply a 0.1 highpass filter.\n\n* Lastly, we will cut the data in 20 s epochs corresponding to the trials.\n\n\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# Load raw data\ndata_path = mne.datasets.ssvep.data_path()\nbids_fname = (\n data_path / \"sub-02\" / \"ses-01\" / \"eeg\" / \"sub-02_ses-01_task-ssvep_eeg.vhdr\"\n)\n\nraw = mne.io.read_raw_brainvision(bids_fname, preload=True, verbose=False)\nraw.info[\"line_freq\"] = 50.0\n\n# Set montage\nmontage = mne.channels.make_standard_montage(\"easycap-M1\")\nraw.set_montage(montage, verbose=False)\n\n# Set common average reference\nraw.set_eeg_reference(\"average\", projection=False, verbose=False)\n\n# Apply bandpass filter\nraw.filter(l_freq=0.1, h_freq=None, fir_design=\"firwin\", verbose=False)\n\n# Construct epochs\nraw.annotations.rename({\"Stimulus/S255\": \"12hz\", \"Stimulus/S155\": \"15hz\"})\ntmin, tmax = -1.0, 20.0 # in s\nbaseline = None\nepochs = mne.Epochs(\n raw,\n event_id=[\"12hz\", \"15hz\"],\n tmin=tmin,\n tmax=tmax,\n baseline=baseline,\n verbose=False,\n)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Frequency analysis\nNow we compute the frequency spectrum of the EEG data.\nYou will already see the peaks at the stimulation frequencies and some of\ntheir harmonics, without any further processing.\n\nThe 'classical' PSD plot will be compared to a plot of the SNR spectrum.\nSNR will be computed as a ratio of the power in a given frequency bin\nto the average power in its neighboring bins.\nThis procedure has two advantages over using the raw PSD:\n\n* it normalizes the spectrum and accounts for 1/f power decay.\n\n* power modulations which are not very narrow band will disappear.\n\n### Calculate power spectral density (PSD)\nThe frequency spectrum will be computed using Fast Fourier transform (FFT).\nThis seems to be common practice in the steady-state literature and is\nbased on the exact knowledge of the stimulus and the assumed response -\nespecially in terms of it's stability over time.\nFor a discussion see e.g.\n[Bach & Meigen (1999)](https://doi.org/10.1023/A:1002648202420)\n\nWe will exclude the first second of each trial from the analysis:\n\n* steady-state response often take a while to stabilize, and the\n transient phase in the beginning can distort the signal estimate.\n\n* this section of data is expected to be dominated by responses related to\n the stimulus onset, and we are not interested in this.\n\nIn MNE we call plain FFT as a special case of Welch's method, with only a\nsingle Welch window spanning the entire trial and no specific windowing\nfunction (i.e. applying a boxcar window).\n\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"tmin = 1.0\ntmax = 20.0\nfmin = 1.0\nfmax = 90.0\nsfreq = epochs.info[\"sfreq\"]\n\nspectrum = epochs.compute_psd(\n \"welch\",\n n_fft=int(sfreq * (tmax - tmin)),\n n_overlap=0,\n n_per_seg=None,\n tmin=tmin,\n tmax=tmax,\n fmin=fmin,\n fmax=fmax,\n window=\"boxcar\",\n verbose=False,\n)\npsds, freqs = spectrum.get_data(return_freqs=True)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Calculate signal to noise ratio (SNR)\n\nSNR - as we define it here - is a measure of relative power:\nit's the ratio of power in a given frequency bin - the 'signal' -\nto a 'noise' baseline - the average power in the surrounding frequency bins.\nThis approach was initially proposed by\n[Meigen & Bach (1999)](https://doi.org/10.1023/A:1002097208337)\n\nHence, we need to set some parameters for this baseline - how many\nneighboring bins should be taken for this computation, and do we want to skip\nthe direct neighbors (this can make sense if the stimulation frequency is not\nsuper constant, or frequency bands are very narrow).\n\nThe function below does what we want.\n\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"def snr_spectrum(psd, noise_n_neighbor_freqs=1, noise_skip_neighbor_freqs=1):\n \"\"\"Compute SNR spectrum from PSD spectrum using convolution.\n\n Parameters\n ----------\n psd : ndarray, shape ([n_trials, n_channels,] n_frequency_bins)\n Data object containing PSD values. Works with arrays as produced by\n MNE's PSD functions or channel/trial subsets.\n noise_n_neighbor_freqs : int\n Number of neighboring frequencies used to compute noise level.\n increment by one to add one frequency bin ON BOTH SIDES\n noise_skip_neighbor_freqs : int\n set this >=1 if you want to exclude the immediately neighboring\n frequency bins in noise level calculation\n\n Returns\n -------\n snr : ndarray, shape ([n_trials, n_channels,] n_frequency_bins)\n Array containing SNR for all epochs, channels, frequency bins.\n NaN for frequencies on the edges, that do not have enough neighbors on\n one side to calculate SNR.\n \"\"\"\n # Construct a kernel that calculates the mean of the neighboring\n # frequencies\n averaging_kernel = np.concatenate(\n (\n np.ones(noise_n_neighbor_freqs),\n np.zeros(2 * noise_skip_neighbor_freqs + 1),\n np.ones(noise_n_neighbor_freqs),\n )\n )\n averaging_kernel /= averaging_kernel.sum()\n\n # Calculate the mean of the neighboring frequencies by convolving with the\n # averaging kernel.\n mean_noise = np.apply_along_axis(\n lambda psd_: np.convolve(psd_, averaging_kernel, mode=\"valid\"), axis=-1, arr=psd\n )\n\n # The mean is not defined on the edges so we will pad it with nas. The\n # padding needs to be done for the last dimension only so we set it to\n # (0, 0) for the other ones.\n edge_width = noise_n_neighbor_freqs + noise_skip_neighbor_freqs\n pad_width = [(0, 0)] * (mean_noise.ndim - 1) + [(edge_width, edge_width)]\n mean_noise = np.pad(mean_noise, pad_width=pad_width, constant_values=np.nan)\n\n return psd / mean_noise"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now we call the function to compute our SNR spectrum.\n\nAs described above, we have to define two parameters.\n\n* how many noise bins do we want?\n\n* do we want to skip the n bins directly next to the target bin?\n\n\nTweaking these parameters *can* drastically impact the resulting spectrum,\nbut mainly if you choose extremes.\nE.g. if you'd skip very many neighboring bins, broad band power modulations\n(such as the alpha peak) should reappear in the SNR spectrum.\nOn the other hand, if you skip none you might miss or smear peaks if the\ninduced power is distributed over two or more frequency bins (e.g. if the\nstimulation frequency isn't perfectly constant, or you have very narrow\nbins).\n\nHere, we want to compare power at each bin with average power of the\n**three neighboring bins** (on each side) and **skip one bin** directly next\nto it.\n\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"snrs = snr_spectrum(psds, noise_n_neighbor_freqs=3, noise_skip_neighbor_freqs=1)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Plot PSD and SNR spectra\nNow we will plot grand average PSD (in blue) and SNR (in red) \u00b1 sd\nfor every frequency bin.\nPSD is plotted on a log scale.\n\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"fig, axes = plt.subplots(2, 1, sharex=\"all\", sharey=\"none\", figsize=(8, 5))\nfreq_range = range(\n np.where(np.floor(freqs) == 1.0)[0][0], np.where(np.ceil(freqs) == fmax - 1)[0][0]\n)\n\npsds_plot = 10 * np.log10(psds)\npsds_mean = psds_plot.mean(axis=(0, 1))[freq_range]\npsds_std = psds_plot.std(axis=(0, 1))[freq_range]\naxes[0].plot(freqs[freq_range], psds_mean, color=\"b\")\naxes[0].fill_between(\n freqs[freq_range], psds_mean - psds_std, psds_mean + psds_std, color=\"b\", alpha=0.2\n)\naxes[0].set(title=\"PSD spectrum\", ylabel=\"Power Spectral Density [dB]\")\n\n# SNR spectrum\nsnr_mean = snrs.mean(axis=(0, 1))[freq_range]\nsnr_std = snrs.std(axis=(0, 1))[freq_range]\n\naxes[1].plot(freqs[freq_range], snr_mean, color=\"r\")\naxes[1].fill_between(\n freqs[freq_range], snr_mean - snr_std, snr_mean + snr_std, color=\"r\", alpha=0.2\n)\naxes[1].set(\n title=\"SNR spectrum\",\n xlabel=\"Frequency [Hz]\",\n ylabel=\"SNR\",\n ylim=[-2, 30],\n xlim=[fmin, fmax],\n)\nfig.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"You can see that the peaks at the stimulation frequencies (12 Hz, 15 Hz)\nand their harmonics are visible in both plots (just as the line noise at\n50 Hz).\nYet, the SNR spectrum shows them more prominently as peaks from a\nnoisy but more or less constant baseline of SNR = 1.\nYou can further see that the SNR processing removes any broad-band power\ndifferences (such as the increased power in alpha band around 10 Hz),\nand also removes the 1/f decay in the PSD.\n\nNote, that while the SNR plot implies the possibility of values below 0\n(mean minus sd) such values do not make sense.\nEach SNR value is a ratio of positive PSD values, and the lowest possible PSD\nvalue is 0 (negative Y-axis values in the upper panel only result from\nplotting PSD on a log scale).\nHence SNR values must be positive and can minimally go towards 0.\n\n## Extract SNR values at the stimulation frequency\n\nOur processing yielded a large array of many SNR values for each trial \u00d7\nchannel \u00d7 frequency-bin of the PSD array.\n\nFor statistical analysis we obviously need to define specific subsets of this\narray. First of all, we are only interested in SNR at the stimulation\nfrequency, but we also want to restrict the analysis to a spatial ROI.\nLastly, answering your interesting research questions will probably rely on\ncomparing SNR in different trials.\n\nTherefore we will have to find the indices of trials, channels, etc.\nAlternatively, one could subselect the trials already at the epoching step,\nusing MNE's event information, and process different epoch structures\nseparately.\n\nLet's only have a look at the trials with 12 Hz stimulation, for now.\n\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# define stimulation frequency\nstim_freq = 12.0"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Get index for the stimulation frequency (12 Hz)\nIdeally, there would be a bin with the stimulation frequency exactly in its\ncenter. However, depending on your Spectral decomposition this is not\nalways the case. We will find the bin closest to it - this one should contain\nour frequency tagged response.\n\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# find index of frequency bin closest to stimulation frequency\ni_bin_12hz = np.argmin(abs(freqs - stim_freq))\n# could be updated to support multiple frequencies\n\n# for later, we will already find the 15 Hz bin and the 1st and 2nd harmonic\n# for both.\ni_bin_24hz = np.argmin(abs(freqs - 24))\ni_bin_36hz = np.argmin(abs(freqs - 36))\ni_bin_15hz = np.argmin(abs(freqs - 15))\ni_bin_30hz = np.argmin(abs(freqs - 30))\ni_bin_45hz = np.argmin(abs(freqs - 45))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Get indices for the different trial types\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"i_trial_12hz = np.where(epochs.annotations.description == \"12hz\")[0]\ni_trial_15hz = np.where(epochs.annotations.description == \"15hz\")[0]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Get indices of EEG channels forming the ROI\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# Define different ROIs\nroi_vis = [\n \"POz\",\n \"Oz\",\n \"O1\",\n \"O2\",\n \"PO3\",\n \"PO4\",\n \"PO7\",\n \"PO8\",\n \"PO9\",\n \"PO10\",\n \"O9\",\n \"O10\",\n] # visual roi\n\n# Find corresponding indices using mne.pick_types()\npicks_roi_vis = mne.pick_types(\n epochs.info, eeg=True, stim=False, exclude=\"bads\", selection=roi_vis\n)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Apply the subset, and check the result\nNow we simply need to apply our selection and yield a result. Therefore,\nwe typically report grand average SNR over the subselection.\n\nIn this tutorial we don't verify the presence of a neural response.\nThis is commonly done in the ASSR literature where SNR is\noften lower. An F-test or Hotelling T\u00b2 would be\nappropriate for this purpose.\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"snrs_target = snrs[i_trial_12hz, :, i_bin_12hz][:, picks_roi_vis]\nprint(\"sub 2, 12 Hz trials, SNR at 12 Hz\")\nprint(f\"average SNR (occipital ROI): {snrs_target.mean()}\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Topography of the vSSR\nBut wait...\nAs described in the intro, we have decided *a priori* to work with average\nSNR over a subset of occipital channels - a visual region of interest (ROI)\n- because we expect SNR to be higher on these channels than in other\nchannels.\n\nLet's check out, whether this was a good decision!\n\nHere we will plot average SNR for each channel location as a topoplot.\nThen we will do a simple paired T-test to check, whether average SNRs over\nthe two sets of channels are significantly different.\n\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# get average SNR at 12 Hz for ALL channels\nsnrs_12hz = snrs[i_trial_12hz, :, i_bin_12hz]\nsnrs_12hz_chaverage = snrs_12hz.mean(axis=0)\n\n# plot SNR topography\nfig, ax = plt.subplots(1)\nmne.viz.plot_topomap(snrs_12hz_chaverage, epochs.info, vlim=(1, None), axes=ax)\n\nprint(\"sub 2, 12 Hz trials, SNR at 12 Hz\")\nprint(f\"average SNR (all channels): {snrs_12hz_chaverage.mean()}\")\nprint(f\"average SNR (occipital ROI): {snrs_target.mean()}\")\n\ntstat_roi_vs_scalp = ttest_rel(snrs_target.mean(axis=1), snrs_12hz.mean(axis=1))\nprint(\n \"12 Hz SNR in occipital ROI is significantly larger than 12 Hz SNR over all \"\n f\"channels: t = {tstat_roi_vs_scalp[0]:.3f}, p = {tstat_roi_vs_scalp[1]}\"\n)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can see, that 1) this participant indeed exhibits a cluster of channels\nwith high SNR in the occipital region and 2) that the average SNR over all\nchannels is smaller than the average of the visual ROI computed above.\nThe difference is statistically significant. Great!\n\nSuch a topography plot can be a nice tool to explore and play with your data\n- e.g. you could try how changing the reference will affect the spatial\ndistribution of SNR values.\n\nHowever, we also wanted to show this plot to point at a potential\nproblem with frequency-tagged (or any other brain imaging) data:\nthere are many channels and somewhere you will likely find some\nstatistically significant effect.\nIt is very easy - even unintended - to end up double-dipping or p-hacking.\nSo if you want to work with an ROI or individual channels, ideally select\nthem *a priori* - before collecting or looking at the data - and preregister\nthis decision so people will believe you.\nIf you end up selecting an ROI or individual channel for reporting *because\nthis channel or ROI shows an effect*, e.g. in an explorative analysis, this\nis also fine but make it transparently and correct for multiple comparison.\n\n## Statistical separation of 12 Hz and 15 Hz vSSR\nAfter this little detour into open science, let's move on and\ndo the analyses we actually wanted to do:\n\nWe will show that we can easily detect and discriminate the brains responses\nin the trials with different stimulation frequencies.\n\nIn the frequency and SNR spectrum plot above, we had all trials mixed up.\nNow we will extract 12 and 15 Hz SNR in both types of trials individually,\nand compare the values with a simple t-test.\nWe will also extract SNR of the 1st and 2nd harmonic for both stimulation\nfrequencies. These are often reported as well and can show interesting\ninteractions.\n\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"snrs_roi = snrs[:, picks_roi_vis, :].mean(axis=1)\n\nfreq_plot = [12, 15, 24, 30, 36, 45]\ncolor_plot = [\"darkblue\", \"darkgreen\", \"mediumblue\", \"green\", \"blue\", \"seagreen\"]\nxpos_plot = [-5.0 / 12, -3.0 / 12, -1.0 / 12, 1.0 / 12, 3.0 / 12, 5.0 / 12]\nfig, ax = plt.subplots()\nlabels = [\"12 Hz trials\", \"15 Hz trials\"]\nx = np.arange(len(labels)) # the label locations\nwidth = 0.6 # the width of the bars\nres = dict()\n\n# loop to plot SNRs at stimulation frequencies and harmonics\nfor i, f in enumerate(freq_plot):\n # extract snrs\n stim_12hz_tmp = snrs_roi[i_trial_12hz, np.argmin(abs(freqs - f))]\n stim_15hz_tmp = snrs_roi[i_trial_15hz, np.argmin(abs(freqs - f))]\n SNR_tmp = [stim_12hz_tmp.mean(), stim_15hz_tmp.mean()]\n # plot (with std)\n ax.bar(\n x + width * xpos_plot[i],\n SNR_tmp,\n width / len(freq_plot),\n yerr=np.std(SNR_tmp),\n label=f\"{f} Hz SNR\",\n color=color_plot[i],\n )\n # store results for statistical comparison\n res[f\"stim_12hz_snrs_{f}hz\"] = stim_12hz_tmp\n res[f\"stim_15hz_snrs_{f}hz\"] = stim_15hz_tmp\n\n# Add some text for labels, title and custom x-axis tick labels, etc.\nax.set_ylabel(\"SNR\")\nax.set_title(\"Average SNR at target frequencies\")\nax.set_xticks(x)\nax.set_xticklabels(labels)\nax.legend([f\"{f} Hz\" for f in freq_plot], title=\"SNR at:\")\nax.set_ylim([0, 70])\nax.axhline(1, ls=\"--\", c=\"r\")\nfig.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"As you can easily see there are striking differences between the trials.\nLet's verify this using a series of two-tailed paired T-Tests.\n\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# Compare 12 Hz and 15 Hz SNR in trials after averaging over channels\n\ntstat_12hz_trial_stim = ttest_rel(\n res[\"stim_12hz_snrs_12hz\"], res[\"stim_12hz_snrs_15hz\"]\n)\nprint(\n \"12 Hz Trials: 12 Hz SNR is significantly higher than 15 Hz SNR: t = \"\n f\"{tstat_12hz_trial_stim[0]:.3f}, p = {tstat_12hz_trial_stim[1]}\"\n)\n\ntstat_12hz_trial_1st_harmonic = ttest_rel(\n res[\"stim_12hz_snrs_24hz\"], res[\"stim_12hz_snrs_30hz\"]\n)\nprint(\n \"12 Hz Trials: 24 Hz SNR is significantly higher than 30 Hz SNR: t = \"\n f\"{tstat_12hz_trial_1st_harmonic[0]:.3f}, p = {tstat_12hz_trial_1st_harmonic[1]}\"\n)\n\ntstat_12hz_trial_2nd_harmonic = ttest_rel(\n res[\"stim_12hz_snrs_36hz\"], res[\"stim_12hz_snrs_45hz\"]\n)\nprint(\n \"12 Hz Trials: 36 Hz SNR is significantly higher than 45 Hz SNR: t = \"\n f\"{tstat_12hz_trial_2nd_harmonic[0]:.3f}, p = {tstat_12hz_trial_2nd_harmonic[1]}\"\n)\n\nprint()\ntstat_15hz_trial_stim = ttest_rel(\n res[\"stim_15hz_snrs_12hz\"], res[\"stim_15hz_snrs_15hz\"]\n)\nprint(\n \"15 Hz trials: 12 Hz SNR is significantly lower than 15 Hz SNR: t = \"\n f\"{tstat_15hz_trial_stim[0]:.3f}, p = {tstat_15hz_trial_stim[1]}\"\n)\n\ntstat_15hz_trial_1st_harmonic = ttest_rel(\n res[\"stim_15hz_snrs_24hz\"], res[\"stim_15hz_snrs_30hz\"]\n)\nprint(\n \"15 Hz trials: 24 Hz SNR is significantly lower than 30 Hz SNR: t = \"\n f\"{tstat_15hz_trial_1st_harmonic[0]:.3f}, p = {tstat_15hz_trial_1st_harmonic[1]}\"\n)\n\ntstat_15hz_trial_2nd_harmonic = ttest_rel(\n res[\"stim_15hz_snrs_36hz\"], res[\"stim_15hz_snrs_45hz\"]\n)\nprint(\n \"15 Hz trials: 36 Hz SNR is significantly lower than 45 Hz SNR: t = \"\n f\"{tstat_15hz_trial_2nd_harmonic[0]:.3f}, p = {tstat_15hz_trial_2nd_harmonic[1]}\"\n)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Debriefing\nSo that's it, we hope you enjoyed our little tour through this example\ndataset.\n\nAs you could see, frequency-tagging is a very powerful tool that can yield\nvery high signal to noise ratios and effect sizes that enable you to detect\nbrain responses even within a single participant and single trials of only\na few seconds duration.\n\n## Bonus exercises\nFor the overly motivated amongst you, let's see what else we can show with\nthese data.\n\nUsing the PSD function as implemented in MNE makes it very easy to change\nthe amount of data that is actually used in the spectrum\nestimation.\n\nHere we employ this to show you some features of frequency\ntagging data that you might or might not have already intuitively expected:\n\n### Effect of trial duration on SNR\nFirst we will simulate shorter trials by taking only the first x s of our 20s\ntrials (2, 4, 6, 8, ..., 20 s), and compute the SNR using a FFT window\nthat covers the entire epoch:\n\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"stim_bandwidth = 0.5\n\n# shorten data and welch window\nwindow_lengths = [i for i in range(2, 21, 2)]\nwindow_snrs = [[]] * len(window_lengths)\nfor i_win, win in enumerate(window_lengths):\n # compute spectrogram\n this_spectrum = epochs[\"12hz\"].compute_psd(\n \"welch\",\n n_fft=int(sfreq * win),\n n_overlap=0,\n n_per_seg=None,\n tmin=0,\n tmax=win,\n window=\"boxcar\",\n fmin=fmin,\n fmax=fmax,\n verbose=False,\n )\n windowed_psd, windowed_freqs = this_spectrum.get_data(return_freqs=True)\n # define a bandwidth of 1 Hz around stimfreq for SNR computation\n bin_width = windowed_freqs[1] - windowed_freqs[0]\n skip_neighbor_freqs = (\n round((stim_bandwidth / 2) / bin_width - bin_width / 2.0 - 0.5)\n if (bin_width < stim_bandwidth)\n else 0\n )\n n_neighbor_freqs = int(\n (\n sum((windowed_freqs <= 13) & (windowed_freqs >= 11))\n - 1\n - 2 * skip_neighbor_freqs\n )\n / 2\n )\n # compute snr\n windowed_snrs = snr_spectrum(\n windowed_psd,\n noise_n_neighbor_freqs=n_neighbor_freqs if (n_neighbor_freqs > 0) else 1,\n noise_skip_neighbor_freqs=skip_neighbor_freqs,\n )\n window_snrs[i_win] = windowed_snrs[\n :, picks_roi_vis, np.argmin(abs(windowed_freqs - 12.0))\n ].mean(axis=1)\n\nfig, ax = plt.subplots(1)\nax.boxplot(window_snrs, tick_labels=window_lengths, vert=True)\nax.set(\n title=\"Effect of trial duration on 12 Hz SNR\",\n ylabel=\"Average SNR\",\n xlabel=\"Trial duration [s]\",\n)\nax.axhline(1, ls=\"--\", c=\"r\")\nfig.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"You can see that the signal estimate / our SNR measure increases with the\ntrial duration.\n\nThis should be easy to understand: in longer recordings there is simply\nmore signal (one second of additional stimulation adds, in our case, 12\ncycles of signal) while the noise is (hopefully) stochastic and not locked\nto the stimulation frequency.\nIn other words: with more data the signal term grows faster than the noise\nterm.\n\nWe can further see that the very short trials with FFT windows < 2-3s are not\ngreat - here we've either hit the noise floor and/or the transient response\nat the trial onset covers too much of the trial.\n\nAgain, this tutorial doesn't statistically test for the presence of a neural\nresponse, but an F-test or Hotelling T\u00b2 would be appropriate for this\npurpose.\n\n### Time resolved SNR\n..and finally we can trick MNE's PSD implementation to make it a\nsliding window analysis and come up with a time resolved SNR measure.\nThis will reveal whether a participant blinked or scratched their head..\n\nEach of the ten trials is coded with a different color in the plot below.\n\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# 3s sliding window\nwindow_length = 4\nwindow_starts = [i for i in range(20 - window_length)]\nwindow_snrs = [[]] * len(window_starts)\n\nfor i_win, win in enumerate(window_starts):\n # compute spectrogram\n this_spectrum = epochs[\"12hz\"].compute_psd(\n \"welch\",\n n_fft=int(sfreq * window_length) - 1,\n n_overlap=0,\n n_per_seg=None,\n window=\"boxcar\",\n tmin=win,\n tmax=win + window_length,\n fmin=fmin,\n fmax=fmax,\n verbose=False,\n )\n windowed_psd, windowed_freqs = this_spectrum.get_data(return_freqs=True)\n # define a bandwidth of 1 Hz around stimfreq for SNR computation\n bin_width = windowed_freqs[1] - windowed_freqs[0]\n skip_neighbor_freqs = (\n round((stim_bandwidth / 2) / bin_width - bin_width / 2.0 - 0.5)\n if (bin_width < stim_bandwidth)\n else 0\n )\n n_neighbor_freqs = int(\n (\n sum((windowed_freqs <= 13) & (windowed_freqs >= 11))\n - 1\n - 2 * skip_neighbor_freqs\n )\n / 2\n )\n # compute snr\n windowed_snrs = snr_spectrum(\n windowed_psd,\n noise_n_neighbor_freqs=n_neighbor_freqs if (n_neighbor_freqs > 0) else 1,\n noise_skip_neighbor_freqs=skip_neighbor_freqs,\n )\n window_snrs[i_win] = windowed_snrs[\n :, picks_roi_vis, np.argmin(abs(windowed_freqs - 12.0))\n ].mean(axis=1)\n\nfig, ax = plt.subplots(1)\ncolors = plt.colormaps[\"Greys\"](np.linspace(0, 1, 10))\nfor i in range(10):\n ax.plot(window_starts, np.array(window_snrs)[:, i], color=colors[i])\nax.set(\n title=f\"Time resolved 12 Hz SNR - {window_length}s sliding window\",\n ylabel=\"Average SNR\",\n xlabel=\"t0 of analysis window [s]\",\n)\nax.axhline(1, ls=\"--\", c=\"r\")\nax.legend([\"individual trials in greyscale\"])\nfig.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Well.. turns out this was a bit too optimistic ;)\n\nBut seriously: this was a nice idea, but we've reached the limit of\nwhat's possible with this single-subject example dataset.\nHowever, there might be data, applications, or research questions\nwhere such an analysis makes sense.\n\n\n"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.12.2"
}
},
"nbformat": 4,
"nbformat_minor": 0
}