langouSNR.ipynb

{
 "metadata": {
  "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.9.5"
  },
  "orig_nbformat": 4,
  "kernelspec": {
   "name": "python3",
   "display_name": "Python 3.9.5 64-bit ('mne': conda)"
  },
  "interpreter": {
   "hash": "e3970cf2a056ed71c6703e76baeca137440641ad1b48b7e7a1b0e3e35696df1b"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2,
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "# License: BSD (3-clause)\n",
    "\n",
    "import numpy as np\n",
    "import os.path as op\n",
    "import matplotlib.pyplot as plt\n",
    "from langouEEG import *\n",
    "\n",
    "import mne\n",
    "import pickle\n",
    "from mne.datasets import sample\n",
    "from mne.minimum_norm import apply_inverse_epochs, read_inverse_operator\n",
    "from mne.connectivity import spectral_connectivity\n",
    "from mne.viz import circular_layout, plot_connectivity_circle\n",
    "import mne\n",
    "from mne.datasets import eegbci\n",
    "from mne.datasets import fetch_fsaverage\n",
    "from mne.datasets import sample\n",
    "from mne.minimum_norm import make_inverse_operator, apply_inverse\n",
    "from mne.minimum_norm import write_inverse_operator\n",
    "import matplotlib.pyplot as plt\n",
    "import mne\n",
    "import numpy as np\n",
    "from scipy.stats import ttest_rel\n",
    "dataRoot = \"/data/home/viscent/Light\""
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "epochs_4F = []\n",
    "epochs_RF = []\n",
    "for subject_name in range(1,21):\n",
    "    if subject_name<10:\n",
    "        subject_name='S0'+str(subject_name)\n",
    "    else:\n",
    "        subject_name='S'+str(subject_name)\n",
    "    with open(dataRoot+'/clean_data/'+subject_name+'_clean.lgeeg','rb') as f:\n",
    "        raw=pickle.load(f)\n",
    "    raw.set_channel_types({'Trigger':'stim','VEO':'eog'})\n",
    "    raw.set_eeg_reference(projection=True)\n",
    "    events, event_dict=extractEvents(raw)\n",
    "    picks = mne.pick_types(raw.info, meg=False, eeg=True, stim=False, eog=True,\n",
    "                        exclude='bads')\n",
    "    epoch_RR,epoch_RF,epoch_4R,epoch_4F = extractEpochs(raw,events,picks)\n",
    "    epochs_4F.append(epoch_4F)\n",
    "    epochs_RF.append(epoch_RF)\n",
    "epochs_4F=mne.concatenate_epochs(epochs_4F)\n",
    "epochs_RF=mne.concatenate_epochs(epochs_RF)"
   ],
   "outputs": [],
   "metadata": {
    "tags": [
     "outputPrepend"
    ]
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "tmin = 10.\n",
    "tmax = 20.\n",
    "fmin = 1.\n",
    "fmax = 120.\n",
    "sfreq = epoch_4F.info['sfreq']\n",
    "\n",
    "# psds, freqs = mne.time_frequency.psd_welch(\n",
    "#     mne.concatenate_epochs([epoch_RF,epoch_RF]),\n",
    "#     n_fft=int(sfreq * (tmax - tmin)),\n",
    "#     n_overlap=0, n_per_seg=None,\n",
    "#     tmin=tmin, tmax=tmax,\n",
    "#     fmin=fmin, fmax=fmax,\n",
    "#     window='boxcar',\n",
    "#     verbose=False)\n",
    "psds, freqs = mne.time_frequency.psd_welch(\n",
    "    epochs_4F,\n",
    "    n_fft=int(sfreq * (tmax - tmin)),\n",
    "    n_overlap=0, n_per_seg=None,\n",
    "    tmin=tmin, tmax=tmax,\n",
    "    fmin=fmin, fmax=fmax,\n",
    "    window='boxcar',\n",
    "    verbose=False)"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "snrs = snr_spectrum(psds, noise_n_neighbor_freqs=3,\n",
    "                    noise_skip_neighbor_freqs=1)"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "fig, axes = plt.subplots(2, 1, sharex='all', sharey='none', figsize=(8, 5))\n",
    "freq_range = range(np.where(np.floor(freqs) == 1.)[0][0],\n",
    "                   np.where(np.ceil(freqs) == fmax - 1)[0][0])\n",
    "\n",
    "psds_plot = 10 * np.log10(psds)\n",
    "psds_mean = psds_plot.mean(axis=(0, 1))[freq_range]\n",
    "psds_std = psds_plot.std(axis=(0, 1))[freq_range]\n",
    "axes[0].plot(freqs[freq_range], psds_mean, color='b')\n",
    "axes[0].fill_between(\n",
    "    freqs[freq_range], psds_mean - psds_std, psds_mean + psds_std,\n",
    "    color='b', alpha=.2)\n",
    "axes[0].set(title=\"PSD spectrum\", ylabel='Power Spectral Density [dB]')\n",
    "\n",
    "# SNR spectrum\n",
    "snr_mean = snrs.mean(axis=(0, 1))[freq_range]\n",
    "snr_std = snrs.std(axis=(0, 1))[freq_range]\n",
    "\n",
    "axes[1].plot(freqs[freq_range], snr_mean, color='r')\n",
    "axes[1].fill_between(\n",
    "    freqs[freq_range], snr_mean - snr_std, snr_mean + snr_std,\n",
    "    color='r', alpha=.2)\n",
    "axes[1].set(\n",
    "    title=\"SNR spectrum\", xlabel='Frequency [Hz]',\n",
    "    ylabel='SNR', ylim=[-2, 30], xlim=[fmin, fmax])\n",
    "fig.show()\n",
    "fig.savefig(op.join(dataRoot,'result','spectrum.svg'))\n",
    "fig.savefig(op.join(dataRoot,'result','spectrum.png'))"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "# define stimulation frequency\n",
    "stim_freq = 40.\n",
    "# find index of frequency bin closest to stimulation frequency\n",
    "i_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.\n",
    "i_bin_40hz = np.argmin(abs(freqs - 40))\n",
    "i_bin_80hz = np.argmin(abs(freqs - 80))"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "i_trial_40 = np.where(epochs.events[:, 2] == 4)[0]\n",
    "i_trial_rand = np.where(epochs.events[:, 2] == 1)[0]"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "roi_vis = ['POZ', 'OZ', 'O1', 'O2', 'PO3', 'PO4', 'PO7',\n",
    "           'PO8']\n",
    "picks_roi_vis = mne.pick_types(epoch_4F.info, eeg=True, stim=False,\n",
    "                               exclude='bads', selection=roi_vis)"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "snrs_target = snrs[i_trial_40, :, i_bin_40hz][:, picks_roi_vis]\n",
    "print(\"40 Hz trials, SNR at 40 Hz\")\n",
    "print(f'average SNR (occipital ROI): {snrs_target.mean()}')"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "snrs_roi = snrs[:, picks_roi_vis, :].mean(axis=1)\n",
    "\n",
    "freq_plot = [40, 80]\n",
    "color_plot = [\n",
    "    'darkblue',  'blue'\n",
    "]\n",
    "xpos_plot = [-5. / 12, -3. / 12, -1. / 12, 1. / 12, 3. / 12, 5. / 12]\n",
    "fig, ax = plt.subplots()\n",
    "labels = ['40 Hz trials', 'random trials']\n",
    "x = np.arange(len(labels))  # the label locations\n",
    "width = 0.6  # the width of the bars\n",
    "res = dict()\n",
    "\n",
    "# loop to plot SNRs at stimulation frequencies and harmonics\n",
    "for i, f in enumerate(freq_plot):\n",
    "    # extract snrs\n",
    "    stim_12hz_tmp = \\\n",
    "        snrs_roi[i_trial_40, np.argmin(abs(freqs - f))]\n",
    "    stim_15hz_tmp = \\\n",
    "        snrs_roi[i_trial_rand, 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], SNR_tmp, width / len(freq_plot),\n",
    "        yerr=np.std(SNR_tmp),\n",
    "        label='%i Hz SNR' % f, color=color_plot[i])\n",
    "    # store results for statistical comparison\n",
    "    res['stim_40hz_snrs_%ihz' % f] = stim_12hz_tmp\n",
    "    res['stim_random_snrs_%ihz' % f] = stim_15hz_tmp\n",
    "\n",
    "# Add some text for labels, title and custom x-axis tick labels, etc.\n",
    "ax.set_ylabel('SNR')\n",
    "ax.set_title('Average SNR at target frequencies')\n",
    "ax.set_xticks(x)\n",
    "ax.set_xticklabels(labels)\n",
    "ax.legend(['%i Hz' % f for f in freq_plot], title='SNR at:')\n",
    "ax.set_ylim([0, 3])\n",
    "ax.axhline(1, ls='--', c='r')\n",
    "fig.show()"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "res_pd=pd.DataFrame(dict([(k, pd.Series(v)) for k, v in res.items()]))\n"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "res_pd"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "res_pd.to_excel(op.join(dataRoot,'result','SNR_data.xlsx'), index=False)"
   ],
   "outputs": [],
   "metadata": {}
  }
 ]
}