wellton/dreamento
1
0
Fork 0
Código Issues Pull Requests Ações Pacotes Projetos Versões Wiki Atividade
Arquivos
main
dreamento/periodogram-testing.ipynb
T
Mahdad Jafarzadehesfahani ac0139f2cd Dreamento class and models added
2022-07-03 14:40:53 +02:00

534 linhas
17 KiB
Plaintext
Original Link permanente Anotar Histórico

{
"cells": [
{
"cell_type": "markdown",
"id": "f668bf21",
"metadata": {},
"source": [
"# Base Code"
]
},
{
"cell_type": "code",
"execution_count": 38,
"id": "f4404bf7",
"metadata": {},
"outputs": [],
"source": [
"import matplotlib.pyplot as plt\n",
"from matplotlib.mlab import psd\n",
"import numpy as np\n",
"from scipy.integrate import simps\n",
"import scipy.signal as ssignal\n",
"\n",
"%matplotlib qt\n",
"\n",
"# %% Plot Power spectral density\n",
"def calculatePowerSpectralDensity(sig, fs, noverlap, NFFT=2 ** 11, scale_by_freq=True):\n",
"\n",
" # Compute power spectrums\n",
" try:\n",
" # psd_sig, f_psd_sig = psd(x=sig, Fs=fs, NFFT=NFFT, scale_by_freq=scale_by_freq, noverlap=noverlap)\n",
" f_psd_sig, psd_sig = ssignal.periodogram(x=sig, fs=fs)#, NFFT=NFFT, scale_by_freq=scale_by_freq, noverlap=noverlap)\n",
"\n",
" except ValueError:\n",
" f_psd_sig, psd_sig = ssignal.periodogram(x=sig, fs=fs)\n",
" # psd_sig, f_psd_sig = psd(x=np.ravel(sig), Fs=fs, NFFT=NFFT, scale_by_freq=scale_by_freq, noverlap=noverlap)\n",
"\n",
" # Compute log of power (optional)\n",
"\n",
" # ======================== Compute band power =========================== #\n",
"\n",
" # Defining EEG bands:\n",
" eeg_bands = {'Delta': (0.5, 4),\n",
" 'Theta': (4, 8),\n",
" 'Alpha': (8, 11),\n",
" 'Beta': (12, 24),\n",
" 'Sigma': (12, 15)}\n",
"\n",
" freq_ix = dict()\n",
"\n",
" fxx = f_psd_sig\n",
" pxx = psd_sig\n",
"\n",
" for band in eeg_bands:\n",
" freq_ix[band] = np.where((fxx >= eeg_bands[band][0]) &\n",
" (fxx <= eeg_bands[band][1]))[0]\n",
"\n",
" freq_resolu_per = fxx[1] - fxx[0]\n",
"\n",
" pow_total = simps(pxx[np.arange(freq_ix['Delta'][0], freq_ix['Beta'][-1])], dx=freq_resolu_per)\n",
" Pow_Delta = simps(pxx[freq_ix['Delta']], dx=freq_resolu_per)\n",
" Pow_Theta = simps(pxx[freq_ix['Theta']], dx=freq_resolu_per)\n",
" Pow_Alpha = simps(pxx[freq_ix['Alpha']], dx=freq_resolu_per)\n",
" Pow_Beta = simps(pxx[freq_ix['Beta']], dx=freq_resolu_per)\n",
" Pow_Sigma = simps(pxx[freq_ix['Sigma']], dx=freq_resolu_per)\n",
"\n",
" # Ratios\n",
" Pow_Delta_ratio = Pow_Delta / pow_total\n",
" Pow_Theta_ratio = Pow_Theta / pow_total\n",
" Pow_Alpha_ratio = Pow_Alpha / pow_total\n",
" Pow_Beta_ratio = Pow_Beta / pow_total\n",
" Pow_Sigma_ratio = Pow_Sigma / pow_total\n",
"\n",
" return psd_sig, f_psd_sig, Pow_Delta_ratio, Pow_Theta_ratio, Pow_Alpha_ratio, Pow_Beta_ratio, Pow_Sigma_ratio\n",
"\n",
"def plotPowerSpectralDensity(figure=None, axis=None, sig=None,\n",
" filtering_status=True,\n",
" lowcut=.3,\n",
" highcut=30,\n",
" ):\n",
" log_power = False\n",
" ylimit = 'auto' # [-50, 10]\n",
" label = 'psd'\n",
" freq_range = [0, 30]\n",
" f = 20\n",
" fs = 256\n",
" Ts = 1 / fs\n",
" t = np.arange(0, 30, Ts)\n",
" if sig is None:\n",
" sig = 1e-6* np.sin(2 * np.pi * f * t) + 2e-6* np.sin(2 * np.pi * 3 * t) + 10e-6*np.squeeze(np.random.rand(len(t), 1))\n",
" else:\n",
" if filtering_status:\n",
" lowcut = lowcut\n",
" highcut = highcut\n",
" nyquist_freq = fs / 2.\n",
" low = lowcut / nyquist_freq\n",
" high = highcut / nyquist_freq\n",
" # Req channel\n",
" print(\"filtering for periodogram\")\n",
" b, a = ssignal.butter(3, [low, high], btype='band')\n",
" sig = ssignal.filtfilt(b, a, sig)\n",
" psd_sig, f_psd_sig, Pow_Delta_ratio, Pow_Theta_ratio, Pow_Alpha_ratio, Pow_Beta_ratio, Pow_Sigma_ratio = \\\n",
" calculatePowerSpectralDensity(sig=sig, fs=256, noverlap= 0, NFFT=len(sig), scale_by_freq=True)\n",
"\n",
" if log_power:\n",
" psd_sig = 20 * np.log10(psd_sig)\n",
"\n",
" if figure == None or axis == None:\n",
" # Open a new fig\n",
" figure, axis = plt.subplots(1, 1, figsize=(20, 10))\n",
"\n",
" # Global setting for axes values size\n",
" plt.rc('xtick', labelsize=11)\n",
" plt.rc('ytick', labelsize=11)\n",
"\n",
" # Plot signals\n",
" axis.plot(f_psd_sig, psd_sig, label=label, color='blue')\n",
"\n",
" # Delta\n",
" axis.axvline(.5, linestyle='--', color='black')\n",
" axis.axvline(4, linestyle='--', color='black')\n",
"\n",
" # Theta\n",
" axis.axvline(8, linestyle='--', color='black')\n",
"\n",
" # Alpha\n",
" axis.axvline(12, linestyle='--', color='black')\n",
"\n",
" # Title and labels\n",
" axis.set(title='Power spectral density', \\\n",
" xlabel='Frequency (Hz)', \\\n",
" ylabel='Power spectral density (dB/ Hz)')\n",
"\n",
" # Legend\n",
" legend = axis.legend([f\"Delta: {round(Pow_Delta_ratio * 100, 2)}% \\n\\\n",
"Theta: {round(Pow_Theta_ratio * 100, 2)}% \\n\\\n",
"Alpha: {round(Pow_Alpha_ratio * 100, 2)}% \\n\\\n",
"Beta: {round(Pow_Beta_ratio * 100, 2)}%\"], prop={'size': 10}, frameon=False)\n",
" legend.set_draggable(state=True)\n",
"\n",
" # Deactivate grid\n",
" plt.grid(False)\n",
"\n",
" # Adding labels\n",
" loc_pow_vals = np.mean([np.min(psd_sig), np.max(psd_sig)])\n",
" loc_labels = loc_pow_vals / 4\n",
" axis.text(1.5, loc_labels, 'Delta', size=10)\n",
" axis.text(5, loc_labels, 'Theta', size=10)\n",
" axis.text(9, loc_labels, 'Alpha', size=10)\n",
" axis.text(13, loc_labels, 'Beta', size=10)\n",
"\n",
" # loc_percents = np.min(psd_sig) + 10\n",
" # # Write power percentages\n",
" # axis.text(1, loc_percents, f'{round(Pow_Delta_ratio * 100, 2)}%', size=8)\n",
" # axis.text(5, loc_percents, f'{round(Pow_Theta_ratio * 100, 2)}%', size=8)\n",
" # axis.text(9, loc_percents, f'{round(Pow_Alpha_ratio * 100, 2)}%', size=8)\n",
" # axis.text(13, loc_percents, f'{round(Pow_Beta_ratio * 100, 2)}%', size=8)\n",
"\n",
" # Limiting x-axis to 0-30 Hz\n",
" axis.set_xlim(freq_range)\n",
"\n",
"# plt.figtext(0.5, 0, f\"Delta: {round(Pow_Delta_ratio * 100, 2)}% | \\\n",
"# Theta: {round(Pow_Theta_ratio * 100, 2)}% | \\\n",
"# Alpha: {round(Pow_Alpha_ratio * 100, 2)}% | \\\n",
"# Beta: {round(Pow_Beta_ratio * 100, 2)}%\", ha=\"center\", fontsize=12, \\\n",
"# bbox={\"facecolor\": \"orange\", \"alpha\": 0.5, \"pad\": 5})\n",
"\n",
" if ylimit == 'auto':\n",
" axis.set_ylim([np.min(psd_sig), np.max(psd_sig)])\n",
"\n",
" else:\n",
" axis.set_ylim(ylimit)\n",
"\n",
" if figure == None or axis == None:\n",
" plt.show()\n",
" # DONT FORGET PLT.SHOW() IF YOU HAVE YOUR OWN FIGURE AND AXIS AS ARGUMENT\n",
"\n",
"if __name__ == \"__main__\":\n",
" fig, ax = plt.subplots(1, 1, figsize=(10, 5))\n",
" plotPowerSpectralDensity(figure=fig, axis=ax)\n",
" plt.tight_layout()\n",
" plt.show()"
]
},
{
"cell_type": "code",
"execution_count": 119,
"id": "6bd404a0",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0.39853406300861594"
]
},
"execution_count": 119,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"import matplotlib.pyplot as plt\n",
"from matplotlib.mlab import psd\n",
"import numpy as np\n",
"from scipy.integrate import simps\n",
"import scipy.signal as ssignal\n",
"\n",
"freq_range = [0, 30]\n",
"f = 20\n",
"fs = 256\n",
"Ts = 1 / fs\n",
"t = np.arange(0, 30, Ts)\n",
"sig = np.sin(2 * np.pi * f * t) + 2*np.sin(2 * np.pi * 3 * t) + 10*np.squeeze(np.random.rand(len(t), 1))\n",
"\n",
"# Compute power spectrums\n",
"f_psd_sig, psd_sig = ssignal.periodogram(x=sig, fs=fs)#, NFFT=NFFT, scale_by_freq=scale_by_freq, noverlap=noverlap)\n",
"\n",
"# ======================== Compute band power =========================== #\n",
"\n",
"# Defining EEG bands:\n",
"eeg_bands = {'Delta': (0.5, 4),\n",
" 'Theta': (4, 8),\n",
" 'Alpha': (8, 11),\n",
" 'Beta': (11, 25)}\n",
"\n",
"freq_ix = dict()\n",
"\n",
"fxx = f_psd_sig\n",
"pxx = psd_sig\n",
"\n",
"for band in eeg_bands:\n",
" freq_ix[band] = np.where((fxx >= eeg_bands[band][0]) &\n",
" (fxx <= eeg_bands[band][1]))[0]\n",
"\n",
"freq_resolu_per = fxx[1] - fxx[0]\n",
"\n",
"pow_total = simps(pxx[np.arange(freq_ix['Delta'][0], freq_ix['Beta'][-1]+1)], dx=freq_resolu_per)\n",
"pow_total = simps(pxx, dx=freq_resolu_per)\n",
"\n",
"Pow_Delta = simps(pxx[freq_ix['Delta']+1], dx=freq_resolu_per)\n",
"Pow_Theta = simps(pxx[freq_ix['Theta']+1], dx=freq_resolu_per)\n",
"Pow_Alpha = simps(pxx[freq_ix['Alpha']+1], dx=freq_resolu_per)\n",
"Pow_Beta = simps(pxx[freq_ix['Beta']+1], dx=freq_resolu_per)\n",
"\n",
"# Ratios\n",
"Pow_Delta_ratio = Pow_Delta / pow_total\n",
"Pow_Theta_ratio = Pow_Theta / pow_total\n",
"Pow_Alpha_ratio = Pow_Alpha / pow_total\n",
"Pow_Beta_ratio = Pow_Beta / pow_total\n",
"\n",
"sum([Pow_Delta_ratio, Pow_Theta_ratio, Pow_Alpha_ratio, Pow_Beta_ratio])"
]
},
{
"cell_type": "code",
"execution_count": 116,
"id": "c6c0da98",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"4.0"
]
},
"execution_count": 116,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"fxx[freq_ix['Delta'][-1]]"
]
},
{
"cell_type": "code",
"execution_count": 133,
"id": "1ef08fe1",
"metadata": {},
"outputs": [],
"source": [
"import yasa\n",
"ret = yasa.bandpower(sig, sf=fs)#, bands = [(0.5, 4, \"Delta\")])"
]
},
{
"cell_type": "code",
"execution_count": 141,
"id": "dd6a82e5",
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"<div>\n",
"<style scoped>\n",
" .dataframe tbody tr th:only-of-type {\n",
" vertical-align: middle;\n",
" }\n",
"\n",
" .dataframe tbody tr th {\n",
" vertical-align: top;\n",
" }\n",
"\n",
" .dataframe thead th {\n",
" text-align: right;\n",
" }\n",
"</style>\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
" <th></th>\n",
" <th>Delta</th>\n",
" <th>Theta</th>\n",
" <th>Alpha</th>\n",
" <th>Sigma</th>\n",
" <th>Beta</th>\n",
" <th>Gamma</th>\n",
" <th>TotalAbsPow</th>\n",
" <th>FreqRes</th>\n",
" <th>Relative</th>\n",
" </tr>\n",
" <tr>\n",
" <th>Chan</th>\n",
" <th></th>\n",
" <th></th>\n",
" <th></th>\n",
" <th></th>\n",
" <th></th>\n",
" <th></th>\n",
" <th></th>\n",
" <th></th>\n",
" <th></th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>CHAN000</th>\n",
" <td>0.437042</td>\n",
" <td>0.045732</td>\n",
" <td>0.042928</td>\n",
" <td>0.046666</td>\n",
" <td>0.312485</td>\n",
" <td>0.115148</td>\n",
" <td>5.090986</td>\n",
" <td>0.25</td>\n",
" <td>True</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
" Delta Theta Alpha Sigma Beta Gamma \\\n",
"Chan \n",
"CHAN000 0.437042 0.045732 0.042928 0.046666 0.312485 0.115148 \n",
"\n",
" TotalAbsPow FreqRes Relative \n",
"Chan \n",
"CHAN000 5.090986 0.25 True "
]
},
"execution_count": 141,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"ret"
]
},
{
"cell_type": "code",
"execution_count": 144,
"id": "1466144d",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"1.0"
]
},
"execution_count": 144,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"ret[\"Delta\"][0]+ret[\"Theta\"][0]+ret[\"Alpha\"][0]+ret[\"Beta\"][0]+ret[\"Gamma\"][0]+ret[\"Sigma\"][0]"
]
},
{
"cell_type": "code",
"execution_count": 156,
"id": "efed1cda",
"metadata": {},
"outputs": [],
"source": [
"import yasa\n",
"ret = yasa.bandpower(sig, sf=fs)#, win_sec = 8)#, bands = [(0.5, 4, \"Delta\")])"
]
},
{
"cell_type": "code",
"execution_count": 157,
"id": "dbbcd72d",
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"<div>\n",
"<style scoped>\n",
" .dataframe tbody tr th:only-of-type {\n",
" vertical-align: middle;\n",
" }\n",
"\n",
" .dataframe tbody tr th {\n",
" vertical-align: top;\n",
" }\n",
"\n",
" .dataframe thead th {\n",
" text-align: right;\n",
" }\n",
"</style>\n",
"<table border=\"1\" class=\"dataframe\">\n",
" <thead>\n",
" <tr style=\"text-align: right;\">\n",
" <th></th>\n",
" <th>Delta</th>\n",
" <th>Theta</th>\n",
" <th>Alpha</th>\n",
" <th>Sigma</th>\n",
" <th>Beta</th>\n",
" <th>Gamma</th>\n",
" <th>TotalAbsPow</th>\n",
" <th>FreqRes</th>\n",
" <th>Relative</th>\n",
" </tr>\n",
" <tr>\n",
" <th>Chan</th>\n",
" <th></th>\n",
" <th></th>\n",
" <th></th>\n",
" <th></th>\n",
" <th></th>\n",
" <th></th>\n",
" <th></th>\n",
" <th></th>\n",
" <th></th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>CHAN000</th>\n",
" <td>0.437042</td>\n",
" <td>0.045732</td>\n",
" <td>0.042928</td>\n",
" <td>0.046666</td>\n",
" <td>0.312485</td>\n",
" <td>0.115148</td>\n",
" <td>5.090986</td>\n",
" <td>0.25</td>\n",
" <td>True</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
" Delta Theta Alpha Sigma Beta Gamma \\\n",
"Chan \n",
"CHAN000 0.437042 0.045732 0.042928 0.046666 0.312485 0.115148 \n",
"\n",
" TotalAbsPow FreqRes Relative \n",
"Chan \n",
"CHAN000 5.090986 0.25 True "
]
},
"execution_count": 157,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"ret"
]
},
{
"cell_type": "code",
"execution_count": 154,
"id": "b65775d9",
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"C:\\Users\\asus\\anaconda3\\envs\\zmax\\lib\\site-packages\\scipy\\signal\\spectral.py:1963: UserWarning: nperseg = 1024 is greater than input length = 735, using nperseg = 735\n",
" .format(nperseg, input_length))\n"
]
}
],
"source": [
"req_sig = sig[np.arange(freq_ix['Delta'][0], freq_ix['Beta'][-1])]\n",
"ret = yasa.bandpower(req_sig, sf=fs)#, bands = [(0.5, 4, \"Delta\")])\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.6.13"
}
},
"nbformat": 4,
"nbformat_minor": 5
}
Referência em uma Nova Issue Ver Git Blame Copiar Link Permanente