%matplotlib inline
%reload_ext autoreload
%autoreload 2
from IPython.display import HTML
HTML('''<script>
code_show=true;
function code_toggle() {
if (code_show){
$('div.input').hide();
$('div.output_stderr').hide();
} else {
$('div.input').show();
$('div.output_stderr').show();
}
code_show = !code_show
}
$( document ).ready(code_toggle);
</script>
<form action='javascript:code_toggle()'><input STYLE='color: #4286f4' type='submit' value='Click here to toggle on/off the raw code.'></form>''')
import numpy as np
import matplotlib.pyplot as plt
from spectral_connectivity import Multitaper
from spectral_connectivity import Connectivity
from spectral_connectivity import multitaper_connectivity
frequency_of_interest = 200
sampling_frequency = 1500
time_extent = (0, 50)
n_time_samples = ((time_extent[1] - time_extent[0]) * sampling_frequency) + 1
time = np.linspace(time_extent[0], time_extent[1], num=n_time_samples, endpoint=True)
signal = np.sin(2 * np.pi * time * frequency_of_interest)
noise = np.random.normal(0, 4, len(signal))
plt.figure(figsize=(15, 6))
plt.subplot(2, 2, 1)
plt.plot(time, signal)
plt.xlabel('Time')
plt.ylabel('Amplitude')
plt.title('Signal', fontweight='bold')
plt.xlim((0.95, 1.05))
plt.ylim((-10, 10))
plt.subplot(2, 2, 2)
plt.plot(time, signal + noise)
plt.xlabel('Time')
plt.ylabel('Amplitude')
plt.title('Signal + Noise', fontweight='bold')
plt.xlim((0.95, 1.05))
plt.ylim((-10, 10))
plt.subplot(2, 2, 3)
m = Multitaper(signal,
sampling_frequency=sampling_frequency,
time_halfbandwidth_product=3,
start_time=time[0])
c = Connectivity(fourier_coefficients=m.fft(),
frequencies=m.frequencies)
plt.plot(c.frequencies, c.power().squeeze());
plt.subplot(2, 2, 4)
m = Multitaper(signal + noise,
sampling_frequency=sampling_frequency,
time_halfbandwidth_product=3,
start_time=time[0])
c = Connectivity(fourier_coefficients=m.fft(),
frequencies=m.frequencies)
plt.plot(c.frequencies, c.power().squeeze());
frequency_of_interest = 30
sampling_frequency = 1500
time_extent = (0, 50)
n_time_samples = ((time_extent[1] - time_extent[0]) * sampling_frequency) + 1
time = np.linspace(time_extent[0], time_extent[1], num=n_time_samples, endpoint=True)
signal = np.sin(2 * np.pi * time * frequency_of_interest)
noise = np.random.normal(0, 4, len(signal))
plt.figure(figsize=(15, 6))
plt.subplot(2, 2, 1)
plt.plot(time, signal)
plt.xlabel('Time')
plt.ylabel('Amplitude')
plt.title('Signal', fontweight='bold')
plt.xlim((0.95, 1.05))
plt.ylim((-10, 10))
plt.subplot(2, 2, 2)
plt.plot(time, signal + noise)
plt.xlabel('Time')
plt.ylabel('Amplitude')
plt.title('Signal + Noise', fontweight='bold')
plt.xlim((0.95, 1.05))
plt.ylim((-10, 10))
plt.subplot(2, 2, 3)
m = Multitaper(signal,
sampling_frequency=sampling_frequency,
time_halfbandwidth_product=3,
start_time=time[0])
c = Connectivity(fourier_coefficients=m.fft(),
frequencies=m.frequencies)
plt.plot(c.frequencies, c.power().squeeze());
plt.subplot(2, 2, 4)
m = Multitaper(signal + noise,
sampling_frequency=sampling_frequency,
time_halfbandwidth_product=3,
start_time=time[0])
c = Connectivity(fourier_coefficients=m.fft(),
frequencies=m.frequencies)
plt.plot(c.frequencies, c.power().squeeze());
frequency_of_interest = [200, 50]
sampling_frequency = 1500
time_extent = (0, 50)
n_time_samples = ((time_extent[1] - time_extent[0]) * sampling_frequency) + 1
time = np.linspace(time_extent[0], time_extent[1], num=n_time_samples, endpoint=True)
signal = np.sin(2 * np.pi * time[:, np.newaxis] * frequency_of_interest)
signal[:n_time_samples // 2, 1] = 0
signal = signal.sum(axis=1)
noise = np.random.normal(0, 4, signal.shape)
fig, axes = plt.subplots(nrows=3, ncols=2, figsize=(15, 9))
axes[0, 0].plot(time, signal)
axes[0, 0].set_xlabel('Time')
axes[0, 0].set_ylabel('Amplitude')
axes[0, 0].set_title('Signal', fontweight='bold')
axes[0, 0].set_xlim((24.90, 25.10))
axes[0, 0].set_ylim((-10, 10))
axes[0, 1].plot(time, signal + noise)
axes[0, 1].set_xlabel('Time')
axes[0, 1].set_ylabel('Amplitude')
axes[0, 1].set_title('Signal + Noise', fontweight='bold')
axes[0, 1].set_xlim((24.90, 25.10))
axes[0, 1].set_ylim((-10, 10))
m = Multitaper(signal,
sampling_frequency=sampling_frequency,
time_halfbandwidth_product=3,
start_time=time[0])
c = Connectivity(fourier_coefficients=m.fft(),
frequencies=m.frequencies)
axes[1, 0].plot(c.frequencies, c.power().squeeze())
m = Multitaper(signal + noise,
sampling_frequency=sampling_frequency,
time_halfbandwidth_product=3,
start_time=time[0])
c = Connectivity(fourier_coefficients=m.fft(),
frequencies=m.frequencies)
axes[1, 1].plot(c.frequencies, c.power().squeeze())
m = Multitaper(signal,
sampling_frequency=sampling_frequency,
time_halfbandwidth_product=3,
time_window_duration=0.600,
time_window_step=0.300,
start_time=time[0])
c = Connectivity(fourier_coefficients=m.fft(),
frequencies=m.frequencies,
time=m.time)
mesh = axes[2, 0].pcolormesh(c.time, c.frequencies, c.power().squeeze().T,
vmin=0.0, vmax=0.03, cmap='viridis')
axes[2, 0].set_ylim((0, 300))
axes[2, 0].axvline(time[int(np.fix(n_time_samples / 2))], color='black')
m = Multitaper(signal + noise,
sampling_frequency=sampling_frequency,
time_halfbandwidth_product=3,
time_window_duration=0.600,
time_window_step=0.300,
start_time=time[0])
c = Connectivity(fourier_coefficients=m.fft(),
frequencies=m.frequencies,
time=m.time)
mesh = axes[2, 1].pcolormesh(c.time, c.frequencies, c.power().squeeze().T,
vmin=0.0, vmax=0.03, cmap='viridis')
axes[2, 1].set_ylim((0, 300))
axes[2, 1].axvline(time[int(np.fix(n_time_samples / 2))], color='black')
plt.tight_layout()
cb = fig.colorbar(mesh, ax=axes.ravel().tolist(), orientation='horizontal',
shrink=.5, aspect=15, pad=0.1, label='Power')
cb.outline.set_linewidth(0)
/Users/edeno/miniconda3/envs/spectral_connectivity/lib/python3.7/site-packages/ipykernel_launcher.py:53: MatplotlibDeprecationWarning: shading='flat' when X and Y have the same dimensions as C is deprecated since 3.3. Either specify the corners of the quadrilaterals with X and Y, or pass shading='auto', 'nearest' or 'gouraud', or set rcParams['pcolor.shading']. This will become an error two minor releases later. /Users/edeno/miniconda3/envs/spectral_connectivity/lib/python3.7/site-packages/ipykernel_launcher.py:67: MatplotlibDeprecationWarning: shading='flat' when X and Y have the same dimensions as C is deprecated since 3.3. Either specify the corners of the quadrilaterals with X and Y, or pass shading='auto', 'nearest' or 'gouraud', or set rcParams['pcolor.shading']. This will become an error two minor releases later.
time_halfbandwidth_product = 1
frequency_of_interest = [200, 50]
time_extent = (0, 0.600)
n_trials = 100
sampling_frequency = 1500
n_time_samples = int(((time_extent[1] - time_extent[0]) * sampling_frequency) + 1)
time = np.linspace(time_extent[0], time_extent[1], num=n_time_samples, endpoint=True)
signal = np.sin(2 * np.pi * time[:, np.newaxis] * frequency_of_interest)
signal[:n_time_samples // 2, 1] = 0
signal = signal.sum(axis=1)[:, np.newaxis, np.newaxis]
noise = np.random.normal(0, 2, size=(n_time_samples, n_trials, 1))
fig, axes = plt.subplots(nrows=3, ncols=2, figsize=(15, 9))
axes[0, 0].plot(time, signal.squeeze())
axes[0, 0].set_xlim(time_extent)
axes[0, 0].set_xlabel('Time')
axes[0, 0].set_ylabel('Amplitude')
axes[0, 0].set_title('Signal', fontweight='bold')
axes[0, 0].set_ylim((-10, 10))
axes[0, 1].plot(time, signal[:, 0, 0] + noise[:, 0, 0])
axes[0, 1].set_xlabel('Time')
axes[0, 1].set_ylabel('Amplitude')
axes[0, 1].set_title('Signal + Noise', fontweight='bold')
axes[0, 1].set_xlim(time_extent)
axes[0, 1].set_ylim((-10, 10))
m = Multitaper(signal,
sampling_frequency=sampling_frequency,
time_halfbandwidth_product=3,
start_time=time[0])
c = Connectivity(fourier_coefficients=m.fft(),
frequencies=m.frequencies)
axes[1, 0].plot(c.frequencies, c.power().squeeze())
m = Multitaper(signal + noise,
sampling_frequency=sampling_frequency,
time_halfbandwidth_product=3,
start_time=time[0])
c = Connectivity(fourier_coefficients=m.fft(),
frequencies=m.frequencies)
axes[1, 1].plot(c.frequencies, c.power().squeeze())
m = Multitaper(signal,
sampling_frequency=sampling_frequency,
time_halfbandwidth_product=time_halfbandwidth_product,
time_window_duration=0.060,
time_window_step=0.060,
start_time=time[0])
c = Connectivity(fourier_coefficients=m.fft(),
frequencies=m.frequencies,
time=m.time)
mesh = axes[2, 0].pcolormesh(c.time, c.frequencies, c.power().squeeze().T,
vmin=0.0, vmax=0.03, cmap='viridis')
axes[2, 0].set_ylim((0, 300))
axes[2, 0].axvline(time[int(np.fix(n_time_samples / 2))], color='black')
m = Multitaper(signal + noise,
sampling_frequency=sampling_frequency,
time_halfbandwidth_product=time_halfbandwidth_product,
time_window_duration=0.060,
time_window_step=0.060,
start_time=time[0])
c = Connectivity(fourier_coefficients=m.fft(),
frequencies=m.frequencies,
time=m.time)
mesh = axes[2, 1].pcolormesh(c.time, c.frequencies, c.power().squeeze().T,
vmin=0.0, vmax=0.03, cmap='viridis')
axes[2, 1].set_ylim((0, 300))
axes[2, 1].axvline(time[int(np.fix(n_time_samples / 2))], color='black')
plt.tight_layout()
cb = fig.colorbar(mesh, ax=axes.ravel().tolist(), orientation='horizontal',
shrink=.5, aspect=15, pad=0.1, label='Power')
cb.outline.set_linewidth(0)
print('frequency resolution: {}'.format(m.frequency_resolution))
/Users/edeno/miniconda3/envs/spectral_connectivity/lib/python3.7/site-packages/ipykernel_launcher.py:56: MatplotlibDeprecationWarning: shading='flat' when X and Y have the same dimensions as C is deprecated since 3.3. Either specify the corners of the quadrilaterals with X and Y, or pass shading='auto', 'nearest' or 'gouraud', or set rcParams['pcolor.shading']. This will become an error two minor releases later. /Users/edeno/miniconda3/envs/spectral_connectivity/lib/python3.7/site-packages/ipykernel_launcher.py:70: MatplotlibDeprecationWarning: shading='flat' when X and Y have the same dimensions as C is deprecated since 3.3. Either specify the corners of the quadrilaterals with X and Y, or pass shading='auto', 'nearest' or 'gouraud', or set rcParams['pcolor.shading']. This will become an error two minor releases later.
frequency resolution: 16.666666666666668
time_halfbandwidth_product = 3
frequency_of_interest = [200, 50]
time_extent = (0, 0.600)
n_trials = 100
sampling_frequency = 1500
n_time_samples = int(((time_extent[1] - time_extent[0]) * sampling_frequency) + 1)
time = np.linspace(time_extent[0], time_extent[1], num=n_time_samples, endpoint=True)
signal = np.sin(2 * np.pi * time[:, np.newaxis] * frequency_of_interest)
signal[:n_time_samples // 2, 1] = 0
signal = signal.sum(axis=1)[:, np.newaxis, np.newaxis]
noise = np.random.normal(0, 2, size=(n_time_samples, n_trials, 1))
fig, axes = plt.subplots(nrows=3, ncols=2, figsize=(15, 9))
axes[0, 0].plot(time, signal.squeeze())
axes[0, 0].set_xlim(time_extent)
axes[0, 0].set_xlabel('Time')
axes[0, 0].set_ylabel('Amplitude')
axes[0, 0].set_title('Signal', fontweight='bold')
axes[0, 0].set_ylim((-10, 10))
axes[0, 1].plot(time, signal[:, 0, 0] + noise[:, 0, 0])
axes[0, 1].set_xlabel('Time')
axes[0, 1].set_ylabel('Amplitude')
axes[0, 1].set_title('Signal + Noise', fontweight='bold')
axes[0, 1].set_xlim(time_extent)
axes[0, 1].set_ylim((-10, 10))
m = Multitaper(signal,
sampling_frequency=sampling_frequency,
time_halfbandwidth_product=time_halfbandwidth_product,
start_time=time[0])
c = Connectivity(fourier_coefficients=m.fft(),
frequencies=m.frequencies)
axes[1, 0].plot(c.frequencies, c.power().squeeze())
m = Multitaper(signal + noise,
sampling_frequency=sampling_frequency,
time_halfbandwidth_product=time_halfbandwidth_product,
start_time=time[0])
c = Connectivity(fourier_coefficients=m.fft(),
frequencies=m.frequencies)
axes[1, 1].plot(c.frequencies, c.power().squeeze())
m = Multitaper(signal,
sampling_frequency=sampling_frequency,
time_halfbandwidth_product=time_halfbandwidth_product,
time_window_duration=0.060,
time_window_step=0.060,
start_time=time[0])
c = Connectivity(fourier_coefficients=m.fft(),
frequencies=m.frequencies,
time=m.time)
mesh = axes[2, 0].pcolormesh(c.time, c.frequencies, c.power().squeeze().T,
vmin=0.0, vmax=0.03, cmap='viridis')
axes[2, 0].set_ylim((0, 300))
axes[2, 0].axvline(time[int(np.fix(n_time_samples / 2))], color='black')
m = Multitaper(signal + noise,
sampling_frequency=sampling_frequency,
time_halfbandwidth_product=time_halfbandwidth_product,
time_window_duration=0.060,
time_window_step=0.060,
start_time=time[0])
c = Connectivity(fourier_coefficients=m.fft(),
frequencies=m.frequencies,
time=m.time)
mesh = axes[2, 1].pcolormesh(c.time, c.frequencies, c.power().squeeze().T,
vmin=0.0, vmax=0.03, cmap='viridis')
axes[2, 1].set_ylim((0, 300))
axes[2, 1].axvline(time[int(np.fix(n_time_samples / 2))], color='black')
plt.tight_layout()
cb = fig.colorbar(mesh, ax=axes.ravel().tolist(), orientation='horizontal',
shrink=.5, aspect=15, pad=0.1, label='Power')
cb.outline.set_linewidth(0)
print('frequency resolution: {}'.format(m.frequency_resolution))
/Users/edeno/miniconda3/envs/spectral_connectivity/lib/python3.7/site-packages/ipykernel_launcher.py:56: MatplotlibDeprecationWarning: shading='flat' when X and Y have the same dimensions as C is deprecated since 3.3. Either specify the corners of the quadrilaterals with X and Y, or pass shading='auto', 'nearest' or 'gouraud', or set rcParams['pcolor.shading']. This will become an error two minor releases later. /Users/edeno/miniconda3/envs/spectral_connectivity/lib/python3.7/site-packages/ipykernel_launcher.py:70: MatplotlibDeprecationWarning: shading='flat' when X and Y have the same dimensions as C is deprecated since 3.3. Either specify the corners of the quadrilaterals with X and Y, or pass shading='auto', 'nearest' or 'gouraud', or set rcParams['pcolor.shading']. This will become an error two minor releases later.
frequency resolution: 50.0
time_halfbandwidth_product = 5
frequency_of_interest = 200
sampling_frequency = 1500
time_extent = (0, 50)
n_signals = 2
n_time_samples = ((time_extent[1] - time_extent[0]) * sampling_frequency) + 1
time = np.linspace(time_extent[0], time_extent[1], num=n_time_samples, endpoint=True)
signal = np.zeros((n_time_samples, n_signals))
signal[:, 0] = np.sin(2 * np.pi * time * frequency_of_interest)
phase_offset = np.pi / 2
signal[:, 1] = np.sin((2 * np.pi * time * frequency_of_interest) + phase_offset)
noise = np.random.normal(0, 4, signal.shape)
plt.figure(figsize=(15, 6))
plt.subplot(2, 2, 1)
plt.title('Signal', fontweight='bold')
plt.plot(time, signal[:, 0], label='Signal1')
plt.plot(time, signal[:, 1], label='Signal2')
plt.xlabel('Time')
plt.ylabel('Amplitude')
plt.xlim((0.95, 1.05))
plt.ylim((-10, 10))
plt.legend()
plt.subplot(2, 2, 2)
plt.title('Signal + Noise', fontweight='bold')
plt.plot(time, signal + noise)
plt.xlabel('Time')
plt.ylabel('Amplitude')
plt.xlim((0.95, 1.05))
plt.ylim((-10, 10))
plt.legend()
m = Multitaper(signal,
sampling_frequency=sampling_frequency,
time_halfbandwidth_product=time_halfbandwidth_product,
start_time=time[0])
c = Connectivity(fourier_coefficients=m.fft(),
frequencies=m.frequencies)
plt.subplot(2, 2, 3)
plt.plot(c.frequencies, c.coherence_magnitude()[0, :, 0, 1])
m = Multitaper(signal + noise,
sampling_frequency=sampling_frequency,
time_halfbandwidth_product=time_halfbandwidth_product,
start_time=time[0])
c = Connectivity(fourier_coefficients=m.fft(),
frequencies=m.frequencies)
plt.subplot(2, 2, 4)
plt.plot(c.frequencies, c.coherence_magnitude()[0, :, 0, 1]);
No handles with labels found to put in legend.
time_halfbandwidth_product = 5
frequency_of_interest = 200
sampling_frequency = 1500
time_extent = (0, 0.600)
n_trials = 100
n_signals = 2
n_time_samples = int(((time_extent[1] - time_extent[0]) * sampling_frequency) + 1)
time = np.linspace(time_extent[0], time_extent[1], num=n_time_samples, endpoint=True)[:, np.newaxis]
signal = np.zeros((n_time_samples, n_trials, n_signals))
signal[:, :, 0] = np.sin(2 * np.pi * time * frequency_of_interest)
phase_offset = np.pi / 2
signal[:, :, 1] = np.sin((2 * np.pi * time * frequency_of_interest) + phase_offset)
noise = np.random.normal(0, 4, signal.shape)
plt.figure(figsize=(15, 6))
plt.subplot(2, 2, 1)
plt.title('Signal', fontweight='bold')
plt.plot(time, signal[:, 0, 0], label='Signal1')
plt.plot(time, signal[:, 0, 1], label='Signal2')
plt.xlabel('Time')
plt.ylabel('Amplitude')
plt.xlim(time_extent)
plt.ylim((-2, 2))
plt.legend()
plt.subplot(2, 2, 2)
plt.title('Signal + Noise', fontweight='bold')
plt.plot(time, signal[:, 0, :] + noise[:, 0, :])
plt.xlabel('Time')
plt.ylabel('Amplitude')
plt.xlim(time_extent)
plt.ylim((-10, 10))
plt.legend()
m = Multitaper(signal,
sampling_frequency=sampling_frequency,
time_halfbandwidth_product=time_halfbandwidth_product,
start_time=time[0])
c = Connectivity(fourier_coefficients=m.fft(),
frequencies=m.frequencies)
plt.subplot(2, 2, 3)
plt.plot(c.frequencies, c.coherence_magnitude()[0, :, 0, 1])
m = Multitaper(signal + noise,
sampling_frequency=sampling_frequency,
time_halfbandwidth_product=time_halfbandwidth_product,
start_time=time[0])
c = Connectivity(fourier_coefficients=m.fft(),
frequencies=m.frequencies)
plt.subplot(2, 2, 4)
plt.plot(c.frequencies, c.coherence_magnitude()[0, :, 0, 1]);
No handles with labels found to put in legend.
time_halfbandwidth_product = 2
frequency_of_interest = 200
sampling_frequency = 1500
time_extent = (0, 2.400)
n_trials = 100
n_signals = 2
n_time_samples = int(((time_extent[1] - time_extent[0]) * sampling_frequency) + 1)
time = np.linspace(time_extent[0], time_extent[1], num=n_time_samples, endpoint=True)[:, np.newaxis]
signal = np.zeros((n_time_samples, n_trials, n_signals))
signal[:, :, 0] = np.sin(2 * np.pi * time * frequency_of_interest)
phase_offset = np.random.uniform(-np.pi, np.pi, size=(n_time_samples, 1))
phase_offset[np.where(time > 1.5), :] = np.pi / 2
signal[:, :, 1] = np.sin((2 * np.pi * time * frequency_of_interest) + phase_offset)
noise = np.random.normal(0, 1, signal.shape)
fig, axes = plt.subplots(nrows=3, ncols=2, figsize=(15, 9))
axes[0, 0].set_title('Signal', fontweight='bold')
axes[0, 0].plot(time, signal[:, 0, 0], label='Signal1')
axes[0, 0].plot(time, signal[:, 0, 1], label='Signal2')
axes[0, 0].set_xlabel('Time')
axes[0, 0].set_ylabel('Amplitude')
axes[0, 0].set_xlim(time_extent)
axes[0, 0].set_ylim((-2, 2))
plt.legend()
axes[0, 1].set_title('Signal + Noise', fontweight='bold')
axes[0, 1].plot(time, signal[:, 0, :] + noise[:, 0, :])
axes[0, 1].set_xlabel('Time')
axes[0, 1].set_ylabel('Amplitude')
axes[0, 1].set_xlim(time_extent)
axes[0, 1].set_ylim((-10, 10))
plt.legend()
m = Multitaper(signal,
sampling_frequency=sampling_frequency,
time_halfbandwidth_product=time_halfbandwidth_product,
start_time=time[0])
c = Connectivity(fourier_coefficients=m.fft(),
frequencies=m.frequencies)
axes[1, 0].plot(c.frequencies, c.coherence_magnitude()[..., 0, 1].squeeze())
axes[1, 0].set_xlim((0, m.nyquist_frequency))
m = Multitaper(signal + noise,
sampling_frequency=sampling_frequency,
time_halfbandwidth_product=time_halfbandwidth_product,
start_time=time[0])
c = Connectivity(fourier_coefficients=m.fft(),
frequencies=m.frequencies)
axes[1, 1].plot(c.frequencies, c.coherence_magnitude()[..., 0, 1].squeeze())
axes[1, 1].set_xlim((0, m.nyquist_frequency))
m = Multitaper(signal,
sampling_frequency=sampling_frequency,
time_halfbandwidth_product=time_halfbandwidth_product,
time_window_duration=0.080,
time_window_step=0.080,
start_time=time[0])
c = Connectivity(fourier_coefficients=m.fft(),
frequencies=m.frequencies,
time=m.time)
time_grid, freq_grid = np.meshgrid(
np.append(c.time, time_extent[-1]),
np.append(c.frequencies, m.nyquist_frequency))
mesh = axes[2, 0].pcolormesh(time_grid, freq_grid, c.coherence_magnitude()[..., 0, 1].squeeze().T,
vmin=0.0, vmax=1.0, cmap='viridis')
axes[2, 0].set_ylim((0, 300))
axes[2, 0].set_xlim(time_extent)
axes[2, 0].axvline(1.5, color='black')
m = Multitaper(signal + noise,
sampling_frequency=sampling_frequency,
time_halfbandwidth_product=time_halfbandwidth_product,
time_window_duration=0.080,
time_window_step=0.080,
start_time=time[0])
c = Connectivity(fourier_coefficients=m.fft(),
frequencies=m.frequencies,
time=m.time)
time_grid, freq_grid = np.meshgrid(
np.append(c.time, time_extent[-1]),
np.append(c.frequencies, m.nyquist_frequency))
mesh = axes[2, 1].pcolormesh(time_grid, freq_grid, c.coherence_magnitude()[..., 0, 1].squeeze().T,
vmin=0.0, vmax=1.0, cmap='viridis')
axes[2, 1].set_ylim((0, 300))
axes[2, 1].set_xlim(time_extent)
axes[2, 1].axvline(1.5, color='black')
plt.tight_layout()
cb = fig.colorbar(mesh, ax=axes.ravel().tolist(), orientation='horizontal',
shrink=.5, aspect=15, pad=0.1, label='Coherence')
cb.outline.set_linewidth(0)
print('frequency resolution: {}'.format(m.frequency_resolution))
No handles with labels found to put in legend. No handles with labels found to put in legend.
frequency resolution: 25.0
time_halfbandwidth_product = 2
frequency_of_interest = 200
sampling_frequency = 1500
time_extent = (0, 2.400)
n_trials = 100
n_signals = 2
n_time_samples = int(((time_extent[1] - time_extent[0]) * sampling_frequency) + 1)
time = np.linspace(time_extent[0], time_extent[1], num=n_time_samples, endpoint=True)[:, np.newaxis]
signal = np.zeros((n_time_samples, n_trials, n_signals))
signal[:, :, 0] = np.sin(2 * np.pi * time * frequency_of_interest)
phase_offset = np.random.uniform(-np.pi, np.pi, size=(n_time_samples, 1))
phase_offset[np.where(time > 1.5), :] = np.pi / 2
signal[:, :, 1] = np.sin((2 * np.pi * time * frequency_of_interest) + phase_offset)
noise = np.random.normal(0, 1, signal.shape)
fig, axes = plt.subplots(nrows=3, ncols=2, figsize=(15, 9))
axes[0, 0].set_title('Signal', fontweight='bold')
axes[0, 0].plot(time, signal[:, 0, 0], label='Signal1')
axes[0, 0].plot(time, signal[:, 0, 1], label='Signal2')
axes[0, 0].set_xlabel('Time')
axes[0, 0].set_ylabel('Amplitude')
axes[0, 0].set_xlim(time_extent)
axes[0, 0].set_ylim((-2, 2))
plt.legend()
axes[0, 1].set_title('Signal + Noise', fontweight='bold')
axes[0, 1].plot(time, signal[:, 0, :] + noise[:, 0, :])
axes[0, 1].set_xlabel('Time')
axes[0, 1].set_ylabel('Amplitude')
axes[0, 1].set_xlim(time_extent)
axes[0, 1].set_ylim((-10, 10))
plt.legend()
m = Multitaper(signal,
sampling_frequency=sampling_frequency,
time_halfbandwidth_product=time_halfbandwidth_product,
start_time=time[0])
c = Connectivity(fourier_coefficients=m.fft(),
frequencies=m.frequencies)
axes[1, 0].plot(c.frequencies, c.imaginary_coherence()[..., 0, 1].squeeze())
axes[1, 0].set_xlim((0, m.nyquist_frequency))
m = Multitaper(signal + noise,
sampling_frequency=sampling_frequency,
time_halfbandwidth_product=time_halfbandwidth_product,
start_time=time[0])
c = Connectivity(fourier_coefficients=m.fft(),
frequencies=m.frequencies)
axes[1, 1].plot(c.frequencies, c.imaginary_coherence()[..., 0, 1].squeeze())
axes[1, 1].set_xlim((0, m.nyquist_frequency))
m = Multitaper(signal,
sampling_frequency=sampling_frequency,
time_halfbandwidth_product=time_halfbandwidth_product,
time_window_duration=0.080,
time_window_step=0.080,
start_time=time[0])
c = Connectivity(fourier_coefficients=m.fft(),
frequencies=m.frequencies,
time=m.time)
time_grid, freq_grid = np.meshgrid(
np.append(c.time, time_extent[-1]),
np.append(c.frequencies, m.nyquist_frequency))
mesh = axes[2, 0].pcolormesh(time_grid, freq_grid, c.imaginary_coherence()[..., 0, 1].squeeze().T,
vmin=0.0, vmax=1.0, cmap='viridis')
axes[2, 0].set_ylim((0, 300))
axes[2, 0].set_xlim(time_extent)
axes[2, 0].axvline(1.5, color='black')
m = Multitaper(signal + noise,
sampling_frequency=sampling_frequency,
time_halfbandwidth_product=time_halfbandwidth_product,
time_window_duration=0.080,
time_window_step=0.080,
start_time=time[0])
c = Connectivity(fourier_coefficients=m.fft(),
frequencies=m.frequencies,
time=m.time)
time_grid, freq_grid = np.meshgrid(
np.append(c.time, time_extent[-1]),
np.append(c.frequencies, m.nyquist_frequency))
mesh = axes[2, 1].pcolormesh(time_grid, freq_grid, c.imaginary_coherence()[..., 0, 1].squeeze().T,
vmin=0.0, vmax=1.0, cmap='viridis')
axes[2, 1].set_ylim((0, 300))
axes[2, 1].set_xlim(time_extent)
axes[2, 1].axvline(1.5, color='black')
plt.tight_layout()
cb = fig.colorbar(mesh, ax=axes.ravel().tolist(), orientation='horizontal',
shrink=.5, aspect=15, pad=0.1, label='Imaginary Coherence')
cb.outline.set_linewidth(0)
print('frequency resolution: {}'.format(m.frequency_resolution))
No handles with labels found to put in legend. No handles with labels found to put in legend.
frequency resolution: 25.0
time_halfbandwidth_product = 2
frequency_of_interest = 200
sampling_frequency = 1500
time_extent = (0, 2.400)
n_trials = 100
n_signals = 2
n_time_samples = int(((time_extent[1] - time_extent[0]) * sampling_frequency) + 1)
time = np.linspace(time_extent[0], time_extent[1], num=n_time_samples, endpoint=True)[:, np.newaxis]
signal = np.zeros((n_time_samples, n_trials, n_signals))
signal[:, :, 0] = np.sin(2 * np.pi * time * frequency_of_interest)
phase_offset = np.random.uniform(-np.pi, np.pi, size=(n_time_samples, 1))
phase_offset[np.where(time > 1.5), :] = np.pi / 2
signal[:, :, 1] = np.sin((2 * np.pi * time * frequency_of_interest) + phase_offset)
noise = np.random.normal(0, 1, signal.shape)
fig, axes = plt.subplots(nrows=3, ncols=2, figsize=(15, 9))
axes[0, 0].set_title('Signal', fontweight='bold')
axes[0, 0].plot(time, signal[:, 0, 0], label='Signal1')
axes[0, 0].plot(time, signal[:, 0, 1], label='Signal2')
axes[0, 0].set_xlabel('Time')
axes[0, 0].set_ylabel('Amplitude')
axes[0, 0].set_xlim(time_extent)
axes[0, 0].set_ylim((-2, 2))
plt.legend()
axes[0, 1].set_title('Signal + Noise', fontweight='bold')
axes[0, 1].plot(time, signal[:, 0, :] + noise[:, 0, :])
axes[0, 1].set_xlabel('Time')
axes[0, 1].set_ylabel('Amplitude')
axes[0, 1].set_xlim(time_extent)
axes[0, 1].set_ylim((-10, 10))
plt.legend()
m = Multitaper(signal,
sampling_frequency=sampling_frequency,
time_halfbandwidth_product=time_halfbandwidth_product,
start_time=time[0])
c = Connectivity(fourier_coefficients=m.fft(),
frequencies=m.frequencies)
axes[1, 0].plot(c.frequencies, abs(c.phase_locking_value())[..., 0, 1].squeeze())
axes[1, 0].set_xlim((0, m.nyquist_frequency))
m = Multitaper(signal + noise,
sampling_frequency=sampling_frequency,
time_halfbandwidth_product=time_halfbandwidth_product,
start_time=time[0])
c = Connectivity(fourier_coefficients=m.fft(),
frequencies=m.frequencies)
axes[1, 1].plot(c.frequencies, abs(c.phase_locking_value())[..., 0, 1].squeeze())
axes[1, 1].set_xlim((0, m.nyquist_frequency))
m = Multitaper(signal,
sampling_frequency=sampling_frequency,
time_halfbandwidth_product=time_halfbandwidth_product,
time_window_duration=0.080,
time_window_step=0.080,
start_time=time[0])
c = Connectivity(fourier_coefficients=m.fft(),
frequencies=m.frequencies,
time=m.time)
time_grid, freq_grid = np.meshgrid(
np.append(c.time, time_extent[-1]),
np.append(c.frequencies, m.nyquist_frequency))
mesh = axes[2, 0].pcolormesh(time_grid, freq_grid, abs(c.phase_locking_value())[..., 0, 1].squeeze().T,
vmin=0.0, vmax=1.0, cmap='viridis')
axes[2, 0].set_ylim((0, 300))
axes[2, 0].set_xlim(time_extent)
axes[2, 0].axvline(1.5, color='black')
m = Multitaper(signal + noise,
sampling_frequency=sampling_frequency,
time_halfbandwidth_product=time_halfbandwidth_product,
time_window_duration=0.080,
time_window_step=0.080,
start_time=time[0])
c = Connectivity(fourier_coefficients=m.fft(),
frequencies=m.frequencies,
time=m.time)
time_grid, freq_grid = np.meshgrid(
np.append(c.time, time_extent[-1]),
np.append(c.frequencies, m.nyquist_frequency))
mesh = axes[2, 1].pcolormesh(time_grid, freq_grid, abs(c.phase_locking_value())[..., 0, 1].squeeze().T,
vmin=0.0, vmax=1.0, cmap='viridis')
axes[2, 1].set_ylim((0, 300))
axes[2, 1].set_xlim(time_extent)
axes[2, 1].axvline(1.5, color='black')
plt.tight_layout()
cb = fig.colorbar(mesh, ax=axes.ravel().tolist(), orientation='horizontal',
shrink=.5, aspect=15, pad=0.1, label='Phase Locking Value')
cb.outline.set_linewidth(0)
print('frequency resolution: {}'.format(m.frequency_resolution))
No handles with labels found to put in legend. No handles with labels found to put in legend.
frequency resolution: 25.0
time_halfbandwidth_product = 2
frequency_of_interest = 200
sampling_frequency = 1500
time_extent = (0, 2.400)
n_trials = 100
n_signals = 2
n_time_samples = int(((time_extent[1] - time_extent[0]) * sampling_frequency) + 1)
time = np.linspace(time_extent[0], time_extent[1], num=n_time_samples, endpoint=True)[:, np.newaxis]
signal = np.zeros((n_time_samples, n_trials, n_signals))
signal[:, :, 0] = np.sin(2 * np.pi * time * frequency_of_interest)
phase_offset = np.random.uniform(-np.pi, np.pi, size=(n_time_samples, 1))
phase_offset[np.where(time > 1.5), :] = np.pi / 2
signal[:, :, 1] = np.sin((2 * np.pi * time * frequency_of_interest) + phase_offset)
noise = np.random.normal(0, 1, signal.shape)
fig, axes = plt.subplots(nrows=3, ncols=2, figsize=(15, 9))
axes[0, 0].set_title('Signal', fontweight='bold')
axes[0, 0].plot(time, signal[:, 0, 0], label='Signal1')
axes[0, 0].plot(time, signal[:, 0, 1], label='Signal2')
axes[0, 0].set_xlabel('Time')
axes[0, 0].set_ylabel('Amplitude')
axes[0, 0].set_xlim(time_extent)
axes[0, 0].set_ylim((-2, 2))
plt.legend()
axes[0, 1].set_title('Signal + Noise', fontweight='bold')
axes[0, 1].plot(time, signal[:, 0, :] + noise[:, 0, :])
axes[0, 1].set_xlabel('Time')
axes[0, 1].set_ylabel('Amplitude')
axes[0, 1].set_xlim(time_extent)
axes[0, 1].set_ylim((-10, 10))
plt.legend()
m = Multitaper(signal,
sampling_frequency=sampling_frequency,
time_halfbandwidth_product=time_halfbandwidth_product,
start_time=time[0])
c = Connectivity(fourier_coefficients=m.fft(),
frequencies=m.frequencies)
axes[1, 0].plot(c.frequencies, c.phase_lag_index()[..., 0, 1].squeeze())
axes[1, 0].set_xlim((0, m.nyquist_frequency))
m = Multitaper(signal + noise,
sampling_frequency=sampling_frequency,
time_halfbandwidth_product=time_halfbandwidth_product,
start_time=time[0])
c = Connectivity(fourier_coefficients=m.fft(),
frequencies=m.frequencies)
axes[1, 1].plot(c.frequencies, c.phase_lag_index()[..., 0, 1].squeeze())
axes[1, 1].set_xlim((0, m.nyquist_frequency))
m = Multitaper(signal,
sampling_frequency=sampling_frequency,
time_halfbandwidth_product=time_halfbandwidth_product,
time_window_duration=0.080,
time_window_step=0.080,
start_time=time[0])
c = Connectivity(fourier_coefficients=m.fft(),
frequencies=m.frequencies,
time=m.time)
time_grid, freq_grid = np.meshgrid(
np.append(c.time, time_extent[-1]),
np.append(c.frequencies, m.nyquist_frequency))
mesh = axes[2, 0].pcolormesh(time_grid, freq_grid, c.phase_lag_index()[..., 0, 1].squeeze().T,
vmin=-1.0, vmax=1.0, cmap='RdBu_r')
axes[2, 0].set_ylim((0, 300))
axes[2, 0].set_xlim(time_extent)
axes[2, 0].axvline(1.5, color='black')
m = Multitaper(signal + noise,
sampling_frequency=sampling_frequency,
time_halfbandwidth_product=time_halfbandwidth_product,
time_window_duration=0.080,
time_window_step=0.080,
start_time=time[0])
c = Connectivity(fourier_coefficients=m.fft(),
frequencies=m.frequencies,
time=m.time)
time_grid, freq_grid = np.meshgrid(
np.append(c.time, time_extent[-1]),
np.append(c.frequencies, m.nyquist_frequency))
mesh = axes[2, 1].pcolormesh(time_grid, freq_grid, c.phase_lag_index()[..., 0, 1].squeeze().T,
vmin=-1.0, vmax=1.0, cmap='RdBu_r')
axes[2, 1].set_ylim((0, 300))
axes[2, 1].set_xlim(time_extent)
axes[2, 1].axvline(1.5, color='black')
plt.tight_layout()
cb = fig.colorbar(mesh, ax=axes.ravel().tolist(), orientation='horizontal',
shrink=.5, aspect=15, pad=0.1, label='Phase Lag Index')
cb.outline.set_linewidth(0)
print('frequency resolution: {}'.format(m.frequency_resolution))
No handles with labels found to put in legend. No handles with labels found to put in legend.
frequency resolution: 25.0
time_halfbandwidth_product = 2
frequency_of_interest = 200
sampling_frequency = 1500
time_extent = (0, 2.400)
n_trials = 100
n_signals = 2
n_time_samples = int(((time_extent[1] - time_extent[0]) * sampling_frequency) + 1)
time = np.linspace(time_extent[0], time_extent[1], num=n_time_samples, endpoint=True)[:, np.newaxis]
signal = np.zeros((n_time_samples, n_trials, n_signals))
signal[:, :, 0] = np.sin(2 * np.pi * time * frequency_of_interest)
phase_offset = np.random.uniform(-np.pi, np.pi, size=(n_time_samples, 1))
phase_offset[np.where(time > 1.5), :] = np.pi / 2
signal[:, :, 1] = np.sin((2 * np.pi * time * frequency_of_interest) + phase_offset)
noise = np.random.normal(0, 1, signal.shape)
fig, axes = plt.subplots(nrows=3, ncols=2, figsize=(15, 9))
axes[0, 0].set_title('Signal', fontweight='bold')
axes[0, 0].plot(time, signal[:, 0, 0], label='Signal1')
axes[0, 0].plot(time, signal[:, 0, 1], label='Signal2')
axes[0, 0].set_xlabel('Time')
axes[0, 0].set_ylabel('Amplitude')
axes[0, 0].set_xlim(time_extent)
axes[0, 0].set_ylim((-2, 2))
plt.legend()
axes[0, 1].set_title('Signal + Noise', fontweight='bold')
axes[0, 1].plot(time, signal[:, 0, :] + noise[:, 0, :])
axes[0, 1].set_xlabel('Time')
axes[0, 1].set_ylabel('Amplitude')
axes[0, 1].set_xlim(time_extent)
axes[0, 1].set_ylim((-10, 10))
plt.legend()
m = Multitaper(signal,
sampling_frequency=sampling_frequency,
time_halfbandwidth_product=time_halfbandwidth_product,
start_time=time[0])
c = Connectivity(fourier_coefficients=m.fft(),
frequencies=m.frequencies)
axes[1, 0].plot(c.frequencies, c.weighted_phase_lag_index()[..., 0, 1].squeeze())
axes[1, 0].set_xlim((0, m.nyquist_frequency))
m = Multitaper(signal + noise,
sampling_frequency=sampling_frequency,
time_halfbandwidth_product=time_halfbandwidth_product,
start_time=time[0])
c = Connectivity(fourier_coefficients=m.fft(),
frequencies=m.frequencies)
axes[1, 1].plot(c.frequencies, c.weighted_phase_lag_index()[..., 0, 1].squeeze())
axes[1, 1].set_xlim((0, m.nyquist_frequency))
m = Multitaper(signal,
sampling_frequency=sampling_frequency,
time_halfbandwidth_product=time_halfbandwidth_product,
time_window_duration=0.080,
time_window_step=0.080,
start_time=time[0])
c = Connectivity(fourier_coefficients=m.fft(),
frequencies=m.frequencies,
time=m.time)
time_grid, freq_grid = np.meshgrid(
np.append(c.time, time_extent[-1]),
np.append(c.frequencies, m.nyquist_frequency))
mesh = axes[2, 0].pcolormesh(time_grid, freq_grid, c.weighted_phase_lag_index()[..., 0, 1].squeeze().T,
vmin=-1.0, vmax=1.0, cmap='RdBu_r')
axes[2, 0].set_ylim((0, 300))
axes[2, 0].set_xlim(time_extent)
axes[2, 0].axvline(1.5, color='black')
m = Multitaper(signal + noise,
sampling_frequency=sampling_frequency,
time_halfbandwidth_product=time_halfbandwidth_product,
time_window_duration=0.080,
time_window_step=0.080,
start_time=time[0])
c = Connectivity(fourier_coefficients=m.fft(),
frequencies=m.frequencies,
time=m.time)
time_grid, freq_grid = np.meshgrid(
np.append(c.time, time_extent[-1]),
np.append(c.frequencies, m.nyquist_frequency))
mesh = axes[2, 1].pcolormesh(time_grid, freq_grid, c.weighted_phase_lag_index()[..., 0, 1].squeeze().T,
vmin=-1.0, vmax=1.0, cmap='RdBu_r')
axes[2, 1].set_ylim((0, 300))
axes[2, 1].set_xlim(time_extent)
axes[2, 1].axvline(1.5, color='black')
plt.tight_layout()
cb = fig.colorbar(mesh, ax=axes.ravel().tolist(), orientation='horizontal',
shrink=.5, aspect=15, pad=0.1, label='Weighted Phase Lag Index')
cb.outline.set_linewidth(0)
print('frequency resolution: {}'.format(m.frequency_resolution))
No handles with labels found to put in legend. No handles with labels found to put in legend.
frequency resolution: 25.0
time_halfbandwidth_product = 2
frequency_of_interest = 200
sampling_frequency = 1500
time_extent = (0, 2.400)
n_trials = 100
n_signals = 2
n_time_samples = int(((time_extent[1] - time_extent[0]) * sampling_frequency) + 1)
time = np.linspace(time_extent[0], time_extent[1], num=n_time_samples, endpoint=True)[:, np.newaxis]
signal = np.zeros((n_time_samples, n_trials, n_signals))
signal[:, :, 0] = np.sin(2 * np.pi * time * frequency_of_interest)
phase_offset = np.random.uniform(-np.pi, np.pi, size=(n_time_samples, 1))
phase_offset[np.where(time > 1.5), :] = np.pi / 2
signal[:, :, 1] = np.sin((2 * np.pi * time * frequency_of_interest) + phase_offset)
noise = np.random.normal(0, 1, signal.shape)
fig, axes = plt.subplots(nrows=3, ncols=2, figsize=(15, 9))
axes[0, 0].set_title('Signal', fontweight='bold')
axes[0, 0].plot(time, signal[:, 0, 0], label='Signal1')
axes[0, 0].plot(time, signal[:, 0, 1], label='Signal2')
axes[0, 0].set_xlabel('Time')
axes[0, 0].set_ylabel('Amplitude')
axes[0, 0].set_xlim(time_extent)
axes[0, 0].set_ylim((-2, 2))
plt.legend()
axes[0, 1].set_title('Signal + Noise', fontweight='bold')
axes[0, 1].plot(time, signal[:, 0, :] + noise[:, 0, :])
axes[0, 1].set_xlabel('Time')
axes[0, 1].set_ylabel('Amplitude')
axes[0, 1].set_xlim(time_extent)
axes[0, 1].set_ylim((-10, 10))
plt.legend()
m = Multitaper(signal,
sampling_frequency=sampling_frequency,
time_halfbandwidth_product=time_halfbandwidth_product,
start_time=time[0])
c = Connectivity(fourier_coefficients=m.fft(),
frequencies=m.frequencies)
axes[1, 0].plot(c.frequencies, c.debiased_squared_phase_lag_index()[..., 0, 1].squeeze())
axes[1, 0].set_xlim((0, m.nyquist_frequency))
m = Multitaper(signal + noise,
sampling_frequency=sampling_frequency,
time_halfbandwidth_product=time_halfbandwidth_product,
start_time=time[0])
c = Connectivity(fourier_coefficients=m.fft(),
frequencies=m.frequencies)
axes[1, 1].plot(c.frequencies, c.debiased_squared_phase_lag_index()[..., 0, 1].squeeze())
axes[1, 1].set_xlim((0, m.nyquist_frequency))
m = Multitaper(signal,
sampling_frequency=sampling_frequency,
time_halfbandwidth_product=time_halfbandwidth_product,
time_window_duration=0.080,
time_window_step=0.080,
start_time=time[0])
c = Connectivity(fourier_coefficients=m.fft(),
frequencies=m.frequencies,
time=m.time)
time_grid, freq_grid = np.meshgrid(
np.append(c.time, time_extent[-1]),
np.append(c.frequencies, m.nyquist_frequency))
mesh = axes[2, 0].pcolormesh(time_grid, freq_grid, c.debiased_squared_phase_lag_index()[..., 0, 1].squeeze().T,
vmin=0.0, vmax=1.0, cmap='viridis')
axes[2, 0].set_ylim((0, 300))
axes[2, 0].set_xlim(time_extent)
axes[2, 0].axvline(1.5, color='black')
m = Multitaper(signal + noise,
sampling_frequency=sampling_frequency,
time_halfbandwidth_product=time_halfbandwidth_product,
time_window_duration=0.080,
time_window_step=0.080,
start_time=time[0])
c = Connectivity(fourier_coefficients=m.fft(),
frequencies=m.frequencies,
time=m.time)
time_grid, freq_grid = np.meshgrid(
np.append(c.time, time_extent[-1]),
np.append(c.frequencies, m.nyquist_frequency))
mesh = axes[2, 1].pcolormesh(time_grid, freq_grid, c.debiased_squared_phase_lag_index()[..., 0, 1].squeeze().T,
vmin=0.0, vmax=1.0, cmap='viridis')
axes[2, 1].set_ylim((0, 300))
axes[2, 1].set_xlim(time_extent)
axes[2, 1].axvline(1.5, color='black')
plt.tight_layout()
cb = fig.colorbar(mesh, ax=axes.ravel().tolist(), orientation='horizontal',
shrink=.5, aspect=15, pad=0.1, label='Debiased Squared Phase Lag Index')
cb.outline.set_linewidth(0)
print('frequency resolution: {}'.format(m.frequency_resolution))
No handles with labels found to put in legend. No handles with labels found to put in legend.
frequency resolution: 25.0
time_halfbandwidth_product = 2
frequency_of_interest = 200
sampling_frequency = 1500
time_extent = (0, 2.400)
n_trials = 100
n_signals = 2
n_time_samples = int(((time_extent[1] - time_extent[0]) * sampling_frequency) + 1)
time = np.linspace(time_extent[0], time_extent[1], num=n_time_samples, endpoint=True)[:, np.newaxis]
signal = np.zeros((n_time_samples, n_trials, n_signals))
signal[:, :, 0] = np.sin(2 * np.pi * time * frequency_of_interest)
phase_offset = np.random.uniform(-np.pi, np.pi, size=(n_time_samples, 1))
phase_offset[np.where(time > 1.5), :] = np.pi / 2
signal[:, :, 1] = np.sin((2 * np.pi * time * frequency_of_interest) + phase_offset)
noise = np.random.normal(0, 1, signal.shape)
fig, axes = plt.subplots(nrows=3, ncols=2, figsize=(15, 9))
axes[0, 0].set_title('Signal', fontweight='bold')
axes[0, 0].plot(time, signal[:, 0, 0], label='Signal1')
axes[0, 0].plot(time, signal[:, 0, 1], label='Signal2')
axes[0, 0].set_xlabel('Time')
axes[0, 0].set_ylabel('Amplitude')
axes[0, 0].set_xlim(time_extent)
axes[0, 0].set_ylim((-2, 2))
plt.legend()
axes[0, 1].set_title('Signal + Noise', fontweight='bold')
axes[0, 1].plot(time, signal[:, 0, :] + noise[:, 0, :])
axes[0, 1].set_xlabel('Time')
axes[0, 1].set_ylabel('Amplitude')
axes[0, 1].set_xlim(time_extent)
axes[0, 1].set_ylim((-10, 10))
plt.legend()
m = Multitaper(signal,
sampling_frequency=sampling_frequency,
time_halfbandwidth_product=time_halfbandwidth_product,
start_time=time[0])
c = Connectivity(fourier_coefficients=m.fft(),
frequencies=m.frequencies)
axes[1, 0].plot(c.frequencies, c.debiased_squared_weighted_phase_lag_index()[..., 0, 1].squeeze())
axes[1, 0].set_xlim((0, m.nyquist_frequency))
m = Multitaper(signal + noise,
sampling_frequency=sampling_frequency,
time_halfbandwidth_product=time_halfbandwidth_product,
start_time=time[0])
c = Connectivity(fourier_coefficients=m.fft(),
frequencies=m.frequencies)
axes[1, 1].plot(c.frequencies, c.debiased_squared_weighted_phase_lag_index()[..., 0, 1].squeeze())
axes[1, 1].set_xlim((0, m.nyquist_frequency))
m = Multitaper(signal,
sampling_frequency=sampling_frequency,
time_halfbandwidth_product=time_halfbandwidth_product,
time_window_duration=0.080,
time_window_step=0.080,
start_time=time[0])
c = Connectivity(fourier_coefficients=m.fft(),
frequencies=m.frequencies,
time=m.time)
time_grid, freq_grid = np.meshgrid(
np.append(c.time, time_extent[-1]),
np.append(c.frequencies, m.nyquist_frequency))
mesh = axes[2, 0].pcolormesh(time_grid, freq_grid, c.debiased_squared_weighted_phase_lag_index()[..., 0, 1].squeeze().T,
vmin=0.0, vmax=1.0, cmap='viridis')
axes[2, 0].set_ylim((0, 300))
axes[2, 0].set_xlim(time_extent)
axes[2, 0].axvline(1.5, color='black')
m = Multitaper(signal + noise,
sampling_frequency=sampling_frequency,
time_halfbandwidth_product=time_halfbandwidth_product,
time_window_duration=0.080,
time_window_step=0.080,
start_time=time[0])
c = Connectivity(fourier_coefficients=m.fft(),
frequencies=m.frequencies,
time=m.time)
time_grid, freq_grid = np.meshgrid(
np.append(c.time, time_extent[-1]),
np.append(c.frequencies, m.nyquist_frequency))
mesh = axes[2, 1].pcolormesh(time_grid, freq_grid, c.debiased_squared_weighted_phase_lag_index()[..., 0, 1].squeeze().T,
vmin=0.0, vmax=1.0, cmap='viridis')
axes[2, 1].set_ylim((0, 300))
axes[2, 1].set_xlim(time_extent)
axes[2, 1].axvline(1.5, color='black')
plt.tight_layout()
cb = fig.colorbar(mesh, ax=axes.ravel().tolist(), orientation='horizontal',
shrink=.5, aspect=15, pad=0.1, label='Debiased Weighted Squared Phase Lag Index')
cb.outline.set_linewidth(0)
print('frequency resolution: {}'.format(m.frequency_resolution))
No handles with labels found to put in legend. No handles with labels found to put in legend.
frequency resolution: 25.0
time_halfbandwidth_product = 2
frequency_of_interest = 200
sampling_frequency = 1500
time_extent = (0, 2.400)
n_trials = 100
n_signals = 2
n_time_samples = int(((time_extent[1] - time_extent[0]) * sampling_frequency) + 1)
time = np.linspace(time_extent[0], time_extent[1], num=n_time_samples, endpoint=True)[:, np.newaxis]
signal = np.zeros((n_time_samples, n_trials, n_signals))
signal[:, :, 0] = np.sin(2 * np.pi * time * frequency_of_interest)
phase_offset = np.random.uniform(-np.pi, np.pi, size=(n_time_samples, 1))
phase_offset[np.where(time > 1.5), :] = np.pi / 2
signal[:, :, 1] = np.sin((2 * np.pi * time * frequency_of_interest) + phase_offset)
noise = np.random.normal(0, 1, signal.shape)
fig, axes = plt.subplots(nrows=3, ncols=2, figsize=(15, 9))
axes[0, 0].set_title('Signal', fontweight='bold')
axes[0, 0].plot(time, signal[:, 0, 0], label='Signal1')
axes[0, 0].plot(time, signal[:, 0, 1], label='Signal2')
axes[0, 0].set_xlabel('Time')
axes[0, 0].set_ylabel('Amplitude')
axes[0, 0].set_xlim(time_extent)
axes[0, 0].set_ylim((-2, 2))
plt.legend()
axes[0, 1].set_title('Signal + Noise', fontweight='bold')
axes[0, 1].plot(time, signal[:, 0, :] + noise[:, 0, :])
axes[0, 1].set_xlabel('Time')
axes[0, 1].set_ylabel('Amplitude')
axes[0, 1].set_xlim(time_extent)
axes[0, 1].set_ylim((-10, 10))
plt.legend()
m = Multitaper(signal,
sampling_frequency=sampling_frequency,
time_halfbandwidth_product=time_halfbandwidth_product,
start_time=time[0])
c = Connectivity(fourier_coefficients=m.fft(),
frequencies=m.frequencies)
axes[1, 0].plot(c.frequencies, c.pairwise_phase_consistency()[..., 0, 1].squeeze())
axes[1, 0].set_xlim((0, m.nyquist_frequency))
m = Multitaper(signal + noise,
sampling_frequency=sampling_frequency,
time_halfbandwidth_product=time_halfbandwidth_product,
start_time=time[0])
c = Connectivity(fourier_coefficients=m.fft(),
frequencies=m.frequencies)
axes[1, 1].plot(c.frequencies, c.pairwise_phase_consistency()[..., 0, 1].squeeze())
axes[1, 1].set_xlim((0, m.nyquist_frequency))
m = Multitaper(signal,
sampling_frequency=sampling_frequency,
time_halfbandwidth_product=time_halfbandwidth_product,
time_window_duration=0.080,
time_window_step=0.080,
start_time=time[0])
c = Connectivity(fourier_coefficients=m.fft(),
frequencies=m.frequencies,
time=m.time)
time_grid, freq_grid = np.meshgrid(
np.append(c.time, time_extent[-1]),
np.append(c.frequencies, m.nyquist_frequency))
mesh = axes[2, 0].pcolormesh(time_grid, freq_grid, c.pairwise_phase_consistency()[..., 0, 1].squeeze().T,
vmin=0.0, vmax=1.0, cmap='viridis')
axes[2, 0].set_ylim((0, 300))
axes[2, 0].set_xlim(time_extent)
axes[2, 0].axvline(1.5, color='black')
m = Multitaper(signal + noise,
sampling_frequency=sampling_frequency,
time_halfbandwidth_product=time_halfbandwidth_product,
time_window_duration=0.080,
time_window_step=0.080,
start_time=time[0])
c = Connectivity(fourier_coefficients=m.fft(),
frequencies=m.frequencies,
time=m.time)
time_grid, freq_grid = np.meshgrid(
np.append(c.time, time_extent[-1]),
np.append(c.frequencies, m.nyquist_frequency))
mesh = axes[2, 1].pcolormesh(time_grid, freq_grid, c.pairwise_phase_consistency()[..., 0, 1].squeeze().T,
vmin=0.0, vmax=1.0, cmap='viridis')
axes[2, 1].set_ylim((0, 300))
axes[2, 1].set_xlim(time_extent)
axes[2, 1].axvline(1.5, color='black')
plt.tight_layout()
cb = fig.colorbar(mesh, ax=axes.ravel().tolist(), orientation='horizontal',
shrink=.5, aspect=15, pad=0.1, label='Pairwise Phase Consistency')
cb.outline.set_linewidth(0)
print('frequency resolution: {}'.format(m.frequency_resolution))
No handles with labels found to put in legend. No handles with labels found to put in legend.
frequency resolution: 25.0
import scipy
sampling_frequency = 1000
time_extent = (0, 1)
n_trials = 500
time_halfbandwidth_product = 1
n_time_samples = ((time_extent[1] - time_extent[0]) * sampling_frequency) + 1
time = np.linspace(time_extent[0], time_extent[1], num=n_time_samples)
signal1 = (scipy.stats.norm.pdf(time, .40, .025) - scipy.stats.norm.pdf(time, .45, .025)) / 10
signal1 = signal1[:, np.newaxis] * np.ones((len(time), n_trials))
signal2 = (scipy.stats.norm.pdf(time, .43, .025) - scipy.stats.norm.pdf(time, .48, .025)) / 10
signal2 = signal2[:, np.newaxis] * np.ones((len(time), n_trials))
noise1 = np.random.normal(0, .2, size=(len(time), n_trials))
noise2 = np.random.normal(0, .1, size=(len(time), n_trials))
data1 = signal1 + noise1
data2 = signal2 + noise2
signals = np.stack((signal1, signal2), axis=-1)
data = np.stack((data1, data2), axis=-1)
fig, axis_handles = plt.subplots(5, 2, figsize=(12, 9))
axis_handles[0, 0].plot(time, signal1, color='blue');
axis_handles[0, 0].plot(time, signal2, color='green');
axis_handles[0, 0].set_xlabel('Time')
axis_handles[0, 1].plot(time, data1, color='blue');
axis_handles[0, 1].plot(time, data2, color='green');
axis_handles[0, 1].set_xlabel('Time')
m = Multitaper(signals,
sampling_frequency=sampling_frequency,
time_halfbandwidth_product=time_halfbandwidth_product,
start_time=time[0])
c = Connectivity(fourier_coefficients=m.fft(),
frequencies=m.frequencies,
time=m.time)
axis_handles[1, 0].plot(c.frequencies, c.power()[..., 0].squeeze())
axis_handles[1, 0].plot(c.frequencies, c.power()[..., 1].squeeze())
axis_handles[2, 0].plot(c.frequencies, c.coherence_magnitude()[..., 0, 1].squeeze())
axis_handles[3, 0].plot(c.frequencies, c.coherence_phase()[..., 0, 1].squeeze(), linestyle='None', marker='8')
delay, slope, r_value = c.group_delay(frequencies_of_interest=[2, 10], frequency_resolution=m.frequency_resolution)
axis_handles[4, 0].bar([1, 2], [delay[..., 0, 1].squeeze(), delay[..., 1, 0].squeeze()], color=['b', 'g']);
axis_handles[4, 0].set_xlim((0.5, 2.5))
axis_handles[4, 0].axhline(0, color='black');
axis_handles[4, 0].set_xticks([1])
axis_handles[4, 0].set_xticklabels(['x1 → x2'])
m = Multitaper(data,
sampling_frequency=sampling_frequency,
time_halfbandwidth_product=time_halfbandwidth_product,
start_time=time[0])
c = Connectivity(fourier_coefficients=m.fft(),
frequencies=m.frequencies,
time=m.time)
axis_handles[1, 1].plot(c.frequencies, c.power()[..., 0].squeeze())
axis_handles[1, 1].plot(c.frequencies, c.power()[..., 1].squeeze())
axis_handles[2, 1].plot(c.frequencies, c.coherence_magnitude()[..., 0, 1].squeeze())
axis_handles[3, 1].plot(c.frequencies, c.coherence_phase()[..., 0, 1].squeeze(), linestyle='None', marker='8')
delay, slope, r_value = c.group_delay(frequencies_of_interest=[2, 10], frequency_resolution=m.frequency_resolution)
axis_handles[4, 1].bar([1, 2], [delay[..., 0, 1].squeeze(), delay[..., 1, 0].squeeze()], color=['b', 'g']);
axis_handles[4, 1].set_xlim((0.5, 2.5))
axis_handles[4, 1].axhline(0, color='black');
axis_handles[4, 1].set_xticks([1])
axis_handles[4, 1].set_xticklabels(['x1 → x2']);
sampling_frequency = 1000
time_extent = (0, 1)
n_trials = 500
time_halfbandwidth_product = 1
n_time_samples = ((time_extent[1] - time_extent[0]) * sampling_frequency) + 1
time = np.linspace(time_extent[0], time_extent[1], num=n_time_samples)
signal1 = (scipy.stats.norm.pdf(time, .43, .025) - scipy.stats.norm.pdf(time, .48, .025)) / 10
signal1 = signal1[:, np.newaxis] * np.ones((len(time), n_trials))
signal2 = (scipy.stats.norm.pdf(time, .40, .025) - scipy.stats.norm.pdf(time, .45, .025)) / 10
signal2 = signal2[:, np.newaxis] * np.ones((len(time), n_trials))
noise1 = np.random.normal(0, .2, size=(len(time), n_trials))
noise2 = np.random.normal(0, .1, size=(len(time), n_trials))
data1 = signal1 + noise1
data2 = signal2 + noise2
signals = np.stack((signal1, signal2), axis=-1)
data = np.stack((data1, data2), axis=-1)
fig, axis_handles = plt.subplots(5, 2, figsize=(12, 9))
axis_handles[0, 0].plot(time, signal1, color='blue');
axis_handles[0, 0].plot(time, signal2, color='green');
axis_handles[0, 0].set_xlabel('Time')
axis_handles[0, 1].plot(time, data1, color='blue');
axis_handles[0, 1].plot(time, data2, color='green');
axis_handles[0, 1].set_xlabel('Time')
m = Multitaper(signals,
sampling_frequency=sampling_frequency,
time_halfbandwidth_product=time_halfbandwidth_product,
start_time=time[0])
c = Connectivity(fourier_coefficients=m.fft(),
frequencies=m.frequencies,
time=m.time)
axis_handles[1, 0].plot(c.frequencies, c.power()[..., 0].squeeze())
axis_handles[1, 0].plot(c.frequencies, c.power()[..., 1].squeeze())
axis_handles[2, 0].plot(c.frequencies, c.coherence_magnitude()[..., 0, 1].squeeze())
axis_handles[3, 0].plot(c.frequencies, c.coherence_phase()[..., 0, 1].squeeze(), linestyle='None', marker='8')
delay, slope, r_value = c.group_delay(frequencies_of_interest=[2, 10], frequency_resolution=m.frequency_resolution)
axis_handles[4, 0].bar([1, 2], [delay[..., 0, 1].squeeze(), delay[..., 1, 0].squeeze()], color=['b', 'g']);
axis_handles[4, 0].set_xlim((0.5, 2.5))
axis_handles[4, 0].axhline(0, color='black');
axis_handles[4, 0].set_xticks([1])
axis_handles[4, 0].set_xticklabels(['x1 → x2']);
m = Multitaper(data,
sampling_frequency=sampling_frequency,
time_halfbandwidth_product=time_halfbandwidth_product,
start_time=time[0])
c = Connectivity(fourier_coefficients=m.fft(),
frequencies=m.frequencies,
time=m.time)
axis_handles[1, 1].plot(c.frequencies, c.power()[..., 0].squeeze())
axis_handles[1, 1].plot(c.frequencies, c.power()[..., 1].squeeze())
axis_handles[2, 1].plot(c.frequencies, c.coherence_magnitude()[..., 0, 1].squeeze())
axis_handles[3, 1].plot(c.frequencies, c.coherence_phase()[..., 0, 1].squeeze(), linestyle='None', marker='8')
delay, slope, r_value = c.group_delay(frequencies_of_interest=[2, 10], frequency_resolution=m.frequency_resolution)
axis_handles[4, 1].bar([1, 2], [delay[..., 0, 1].squeeze(), delay[..., 1, 0].squeeze()], color=['b', 'g']);
axis_handles[4, 1].set_xlim((0.5, 2.5))
axis_handles[4, 1].axhline(0, color='black');
axis_handles[4, 1].set_xticks([1])
axis_handles[4, 1].set_xticklabels(['x1 → x2']);
sampling_frequency = 1000
time_extent = (0, 2)
n_trials = 500
time_halfbandwidth_product = 1
n_time_samples = ((time_extent[1] - time_extent[0]) * sampling_frequency) + 1
time = np.linspace(time_extent[0], time_extent[1], num=n_time_samples)
signal1 = (scipy.stats.norm.pdf(time, .43, .025) - scipy.stats.norm.pdf(time, .48, .025)) / 10
signal1 = signal1[:, np.newaxis] * np.ones((len(time), n_trials))
signal2 = (scipy.stats.norm.pdf(time, .40, .025) - scipy.stats.norm.pdf(time, .45, .025)) / 10
signal2 = signal2[:, np.newaxis] * np.ones((len(time), n_trials))
noise1 = np.random.normal(0, .2, size=(len(time), n_trials))
noise2 = np.random.normal(0, .1, size=(len(time), n_trials))
data1 = signal1 + noise1
data2 = signal2 + noise2
signals = np.stack((signal1, signal2), axis=-1)
data = np.stack((data1, data2), axis=-1)
fig, axis_handles = plt.subplots(2, 2, figsize=(12, 9))
axis_handles[0, 0].plot(time, signal1, color='blue');
axis_handles[0, 0].plot(time, signal2, color='green');
axis_handles[0, 0].set_xlabel('Time')
axis_handles[0, 0].set_xlim(time_extent);
axis_handles[0, 1].plot(time, data1, color='blue');
axis_handles[0, 1].plot(time, data2, color='green');
axis_handles[0, 1].set_xlabel('Time')
axis_handles[0, 1].set_xlim(time_extent);
m = Multitaper(signals,
sampling_frequency=sampling_frequency,
time_halfbandwidth_product=time_halfbandwidth_product,
time_window_duration=0.500,
time_window_step=0.100,
start_time=time[0])
c = Connectivity(fourier_coefficients=m.fft(),
frequencies=m.frequencies,
time=m.time)
delay, slope, r_value = c.group_delay(frequencies_of_interest=[0, 15], frequency_resolution=m.frequency_resolution)
axis_handles[1, 0].plot(c.time + m.time_window_duration/2, delay[..., 0, 1]);
axis_handles[1, 0].set_xlim(time_extent)
m = Multitaper(data,
sampling_frequency=sampling_frequency,
time_halfbandwidth_product=time_halfbandwidth_product,
time_window_duration=0.500,
time_window_step=0.100,
start_time=time[0])
c = Connectivity(fourier_coefficients=m.fft(),
frequencies=m.frequencies,
time=m.time)
delay, slope, r_value = c.group_delay(frequencies_of_interest=[0, 15], frequency_resolution=m.frequency_resolution)
axis_handles[1, 1].plot(c.time + m.time_window_duration/2, delay[..., 0, 1]);
axis_handles[1, 1].set_xlim(time_extent);
sampling_frequency = 1000
time_extent = (0, 1)
n_trials = 500
time_halfbandwidth_product = 1
n_time_samples = ((time_extent[1] - time_extent[0]) * sampling_frequency) + 1
time = np.linspace(time_extent[0], time_extent[1], num=n_time_samples)
signal1 = (scipy.stats.norm.pdf(time, .40, .025) - scipy.stats.norm.pdf(time, .45, .025)) / 10
signal1 = signal1[:, np.newaxis] * np.ones((len(time), n_trials))
signal2 = (scipy.stats.norm.pdf(time, .43, .025) - scipy.stats.norm.pdf(time, .48, .025)) / 10
signal2 = signal2[:, np.newaxis] * np.ones((len(time), n_trials))
noise1 = np.random.normal(0, .2, size=(len(time), n_trials))
noise2 = np.random.normal(0, .1, size=(len(time), n_trials))
data1 = signal1 + noise1
data2 = signal2 + noise2
signals = np.stack((signal1, signal2), axis=-1)
data = np.stack((data1, data2), axis=-1)
fig, axis_handles = plt.subplots(5, 2, figsize=(12, 9))
axis_handles[0, 0].plot(time, signal1, color='blue');
axis_handles[0, 0].plot(time, signal2, color='green');
axis_handles[0, 0].set_xlabel('Time')
axis_handles[0, 1].plot(time, data1, color='blue');
axis_handles[0, 1].plot(time, data2, color='green');
axis_handles[0, 1].set_xlabel('Time')
m = Multitaper(signals,
sampling_frequency=sampling_frequency,
time_halfbandwidth_product=time_halfbandwidth_product,
start_time=time[0])
c = Connectivity(fourier_coefficients=m.fft(),
frequencies=m.frequencies,
time=m.time)
axis_handles[1, 0].plot(c.frequencies, c.power()[..., 0].squeeze())
axis_handles[1, 0].plot(c.frequencies, c.power()[..., 1].squeeze())
axis_handles[2, 0].plot(c.frequencies, c.coherence_magnitude()[..., 0, 1].squeeze())
axis_handles[3, 0].plot(c.frequencies, c.coherence_phase()[..., 0, 1].squeeze(), linestyle='None', marker='8')
psi = c.phase_slope_index(frequencies_of_interest=[2, 10], frequency_resolution=m.frequency_resolution)
axis_handles[4, 0].bar([1, 2], [psi[..., 0, 1].squeeze(), psi[..., 1, 0].squeeze()], color=['b', 'g']);
axis_handles[4, 0].set_xlim((0.5, 2.5))
axis_handles[4, 0].axhline(0, color='black');
m = Multitaper(data,
sampling_frequency=sampling_frequency,
time_halfbandwidth_product=time_halfbandwidth_product,
start_time=time[0])
c = Connectivity(fourier_coefficients=m.fft(),
frequencies=m.frequencies,
time=m.time)
axis_handles[1, 1].plot(c.frequencies, c.power()[..., 0].squeeze())
axis_handles[1, 1].plot(c.frequencies, c.power()[..., 1].squeeze())
axis_handles[2, 1].plot(c.frequencies, c.coherence_magnitude()[..., 0, 1].squeeze())
axis_handles[3, 1].plot(c.frequencies, c.coherence_phase()[..., 0, 1].squeeze(), linestyle='None', marker='8')
psi = c.phase_slope_index(frequencies_of_interest=[2, 10], frequency_resolution=m.frequency_resolution)
axis_handles[4, 1].bar([1, 2], [psi[..., 0, 1].squeeze(), psi[..., 1, 0].squeeze()], color=['b', 'g']);
axis_handles[4, 1].set_xlim((0.5, 2.5));
axis_handles[4, 1].axhline(0, color='black');
sampling_frequency = 1000
time_extent = (0, 1)
n_trials = 500
time_halfbandwidth_product = 1
n_time_samples = ((time_extent[1] - time_extent[0]) * sampling_frequency) + 1
time = np.linspace(time_extent[0], time_extent[1], num=n_time_samples)
signal1 = (scipy.stats.norm.pdf(time, .43, .025) - scipy.stats.norm.pdf(time, .48, .025)) / 10
signal1 = signal1[:, np.newaxis] * np.ones((len(time), n_trials))
signal2 = (scipy.stats.norm.pdf(time, .40, .025) - scipy.stats.norm.pdf(time, .45, .025)) / 10
signal2 = signal2[:, np.newaxis] * np.ones((len(time), n_trials))
noise1 = np.random.normal(0, .2, size=(len(time), n_trials))
noise2 = np.random.normal(0, .1, size=(len(time), n_trials))
data1 = signal1 + noise1
data2 = signal2 + noise2
signals = np.stack((signal1, signal2), axis=-1)
data = np.stack((data1, data2), axis=-1)
fig, axis_handles = plt.subplots(5, 2, figsize=(12, 9))
axis_handles[0, 0].plot(time, signal1, color='blue');
axis_handles[0, 0].plot(time, signal2, color='green');
axis_handles[0, 0].set_xlabel('Time')
axis_handles[0, 1].plot(time, data1, color='blue');
axis_handles[0, 1].plot(time, data2, color='green');
axis_handles[0, 1].set_xlabel('Time')
m = Multitaper(signals,
sampling_frequency=sampling_frequency,
time_halfbandwidth_product=time_halfbandwidth_product,
start_time=time[0])
c = Connectivity(fourier_coefficients=m.fft(),
frequencies=m.frequencies,
time=m.time)
axis_handles[1, 0].plot(c.frequencies, c.power()[..., 0].squeeze())
axis_handles[1, 0].plot(c.frequencies, c.power()[..., 1].squeeze())
axis_handles[2, 0].plot(c.frequencies, c.coherence_magnitude()[..., 0, 1].squeeze())
axis_handles[3, 0].plot(c.frequencies, c.coherence_phase()[..., 0, 1].squeeze(), linestyle='None', marker='8')
psi = c.phase_slope_index(frequencies_of_interest=[2, 10], frequency_resolution=m.frequency_resolution)
axis_handles[4, 0].bar([1, 2], [psi[..., 0, 1].squeeze(), psi[..., 1, 0].squeeze()], color=['b', 'g']);
axis_handles[4, 0].set_xlim((0.5, 2.5))
axis_handles[4, 0].axhline(0, color='black');
m = Multitaper(data,
sampling_frequency=sampling_frequency,
time_halfbandwidth_product=time_halfbandwidth_product,
start_time=time[0])
c = Connectivity(fourier_coefficients=m.fft(),
frequencies=m.frequencies,
time=m.time)
axis_handles[1, 1].plot(c.frequencies, c.power()[..., 0].squeeze())
axis_handles[1, 1].plot(c.frequencies, c.power()[..., 1].squeeze())
axis_handles[2, 1].plot(c.frequencies, c.coherence_magnitude()[..., 0, 1].squeeze())
axis_handles[3, 1].plot(c.frequencies, c.coherence_phase()[..., 0, 1].squeeze(), linestyle='None', marker='8')
psi = c.phase_slope_index(frequencies_of_interest=[2, 10], frequency_resolution=m.frequency_resolution)
axis_handles[4, 1].bar([1, 2], [psi[..., 0, 1].squeeze(), psi[..., 1, 0].squeeze()], color=['b', 'g']);
axis_handles[4, 1].set_xlim((0.5, 2.5));
axis_handles[4, 1].axhline(0, color='black');
sampling_frequency = 1000
time_extent = (0, 2)
n_trials = 500
time_halfbandwidth_product = 1
n_time_samples = ((time_extent[1] - time_extent[0]) * sampling_frequency) + 1
time = np.linspace(time_extent[0], time_extent[1], num=n_time_samples)
signal1 = (scipy.stats.norm.pdf(time, .43, .025) - scipy.stats.norm.pdf(time, .48, .025)) / 10
signal1 = signal1[:, np.newaxis] * np.ones((len(time), n_trials))
signal2 = (scipy.stats.norm.pdf(time, .40, .025) - scipy.stats.norm.pdf(time, .45, .025)) / 10
signal2 = signal2[:, np.newaxis] * np.ones((len(time), n_trials))
noise1 = np.random.normal(0, .2, size=(len(time), n_trials))
noise2 = np.random.normal(0, .1, size=(len(time), n_trials))
data1 = signal1 + noise1
data2 = signal2 + noise2
signals = np.stack((signal1, signal2), axis=-1)
data = np.stack((data1, data2), axis=-1)
fig, axis_handles = plt.subplots(2, 2, figsize=(12, 9), constrained_layout=True, sharex=True)
axis_handles[0, 0].plot(time, signal1, color='blue');
axis_handles[0, 0].plot(time, signal2, color='green');
axis_handles[0, 0].set_title('Signals')
axis_handles[0, 0].set_xlim(time_extent);
axis_handles[0, 1].plot(time, data1, color='blue');
axis_handles[0, 1].plot(time, data2, color='green');
axis_handles[0, 1].set_title('Signals')
axis_handles[0, 1].set_xlim(time_extent);
m = Multitaper(signals,
sampling_frequency=sampling_frequency,
time_halfbandwidth_product=time_halfbandwidth_product,
time_window_duration=0.500,
time_window_step=0.100,
start_time=time[0])
c = Connectivity(fourier_coefficients=m.fft(),
frequencies=m.frequencies,
time=m.time)
psi = c.phase_slope_index(frequencies_of_interest=[0, 15], frequency_resolution=m.frequency_resolution)
axis_handles[1, 0].plot(c.time + m.time_window_duration/2, psi[..., 0, 1], c.time + m.time_window_duration/2, psi[..., 1, 0])
axis_handles[1, 0].set_xlim(time_extent)
axis_handles[1, 0].set_xlabel('Time [s]')
axis_handles[1, 0].set_ylabel('Phase Slope Index')
m = Multitaper(data,
sampling_frequency=sampling_frequency,
time_halfbandwidth_product=time_halfbandwidth_product,
time_window_duration=0.500,
time_window_step=0.100,
start_time=time[0])
c = Connectivity(fourier_coefficients=m.fft(),
frequencies=m.frequencies,
time=m.time)
psi = c.phase_slope_index(frequencies_of_interest=[0, 15], frequency_resolution=m.frequency_resolution)
axis_handles[1, 1].plot(c.time + m.time_window_duration/2, psi[..., 0, 1], c.time + m.time_window_duration/2, psi[..., 1, 0]);
axis_handles[1, 1].set_xlim(time_extent);
axis_handles[1, 1].set_xlabel('Time [s]')
axis_handles[1, 1].set_ylabel('Phase Slope Index')
Text(0, 0.5, 'Phase Slope Index')
The advantage of canonical coherence is that it can be more statistically powerful than coherence because it is combining coherence from groups.
from itertools import product
time_halfbandwidth_product = 2
frequency_of_interest = 200
sampling_frequency = 1500
time_extent = (0, 2.400)
n_trials = 100
n_signals = 4
n_time_samples = int(((time_extent[1] - time_extent[0]) * sampling_frequency) + 1)
time = np.linspace(time_extent[0], time_extent[1], num=n_time_samples, endpoint=True)
signal = np.zeros((n_time_samples, n_trials, n_signals))
signal[:, :, 0:2] = (
np.sin(2 * np.pi * time * frequency_of_interest)[:, np.newaxis, np.newaxis] *
np.ones((1, n_trials, 2)))
other_signals = (n_signals + 1) // 2
n_other_signals = n_signals - other_signals
phase_offset = np.random.uniform(-np.pi, np.pi, size=(n_time_samples, n_trials, n_other_signals))
phase_offset[np.where(time > 1.5), :] = np.pi / 2
signal[:, :, other_signals:] = np.sin(
(2 * np.pi * time[:, np.newaxis, np.newaxis] * frequency_of_interest) + phase_offset)
noise = np.random.normal(0, 4, signal.shape)
m = Multitaper(signal + noise,
sampling_frequency=sampling_frequency,
time_halfbandwidth_product=time_halfbandwidth_product,
time_window_duration=0.080,
time_window_step=0.080,
start_time=time[0])
c = Connectivity(fourier_coefficients=m.fft(),
frequencies=m.frequencies,
time=m.time)
time_grid, freq_grid = np.meshgrid(
np.append(c.time, time_extent[-1]),
np.append(c.frequencies, m.nyquist_frequency))
fig, axes = plt.subplots(nrows=n_signals, ncols=n_signals, figsize=(15, 9))
meshes = list()
for ind1, ind2 in product(range(n_signals), range(n_signals)):
if ind1 == ind2:
vmin, vmax = c.power().min(), c.power().max()
else:
vmin, vmax = 0, 0.5
mesh = axes[ind1, ind2].pcolormesh(
time_grid, freq_grid, c.coherence_magnitude()[..., ind1, ind2].squeeze().T,
vmin=vmin, vmax=vmax, cmap='viridis')
meshes.append(mesh)
axes[ind1, ind2].set_ylim((0, 300))
axes[ind1, ind2].set_xlim(time_extent)
axes[ind1, ind2].axvline(1.5, color='black')
plt.suptitle('Coherence', y=1.02, fontsize=30)
plt.tight_layout()
cb = fig.colorbar(meshes[-2], ax=axes.ravel().tolist(), orientation='horizontal',
shrink=.5, aspect=15, pad=0.1, label='Coherence')
cb.outline.set_linewidth(0)
print('frequency resolution: {}'.format(m.frequency_resolution))
group_labels = (['a'] * (n_signals - n_other_signals)) + (['b'] * n_other_signals)
canonical_coherence, pair_labels = c.canonical_coherence(group_labels)
fig = plt.figure()
mesh = plt.pcolormesh(time_grid, freq_grid, canonical_coherence[..., 0, 1].squeeze().T,
vmin=0, vmax=0.5, cmap='viridis')
plt.suptitle('Canonical Coherence', y=1.02, fontsize=30)
cb = fig.colorbar(mesh, ax=plt.gca(), orientation='horizontal',
shrink=.5, aspect=15, pad=0.1, label='Coherence')
cb.outline.set_linewidth(0);
frequency resolution: 25.0
from itertools import product
time_halfbandwidth_product = 2
frequency_of_interest = 200
sampling_frequency = 1500
time_extent = (0, 2.400)
n_trials = 100
n_signals = 6
n_time_samples = int(((time_extent[1] - time_extent[0]) * sampling_frequency) + 1)
time = np.linspace(time_extent[0], time_extent[1], num=n_time_samples, endpoint=True)
signal = np.zeros((n_time_samples, n_trials, n_signals))
signal[:, :, 0:2] = (
np.sin(2 * np.pi * time * frequency_of_interest)[:, np.newaxis, np.newaxis] *
np.ones((1, n_trials, 2)))
other_signals = (n_signals + 1) // 2
n_other_signals = n_signals - other_signals
phase_offset = np.random.uniform(-np.pi, np.pi, size=(n_time_samples, n_trials, n_other_signals))
phase_offset[np.where(time > 1.5), :] = np.pi / 2
signal[:, :, other_signals:] = np.sin(
(2 * np.pi * time[:, np.newaxis, np.newaxis] * frequency_of_interest) + phase_offset)
noise = np.random.normal(10, 7, signal.shape)
m = Multitaper(signal + noise,
sampling_frequency=sampling_frequency,
time_halfbandwidth_product=time_halfbandwidth_product,
time_window_duration=0.080,
time_window_step=0.080,
start_time=time[0])
c = Connectivity(fourier_coefficients=m.fft(),
frequencies=m.frequencies,
time=m.time)
time_grid, freq_grid = np.meshgrid(
np.append(c.time, time_extent[-1]),
np.append(c.frequencies, m.nyquist_frequency))
fig, axes = plt.subplots(nrows=n_signals, ncols=n_signals, figsize=(15, 9))
meshes = list()
for ind1, ind2 in product(range(n_signals), range(n_signals)):
if ind1 == ind2:
vmin, vmax = c.power().min(), c.power().max()
else:
vmin, vmax = 0, 0.5
mesh = axes[ind1, ind2].pcolormesh(
time_grid, freq_grid, c.coherence_magnitude()[..., ind1, ind2].squeeze().T,
vmin=vmin, vmax=vmax, cmap='viridis')
meshes.append(mesh)
axes[ind1, ind2].set_ylim((0, 300))
axes[ind1, ind2].set_xlim(time_extent)
axes[ind1, ind2].axvline(1.5, color='black')
plt.suptitle('Coherence', y=1.02, fontsize=30)
plt.tight_layout()
cb = fig.colorbar(meshes[-2], ax=axes.ravel().tolist(), orientation='horizontal',
shrink=.5, aspect=15, pad=0.1, label='Coherence')
cb.outline.set_linewidth(0)
print('frequency resolution: {}'.format(m.frequency_resolution))
group_labels = (['a'] * (n_signals - n_other_signals)) + (['b'] * n_other_signals)
canonical_coherence, pair_labels = c.canonical_coherence(group_labels)
fig = plt.figure()
mesh = plt.pcolormesh(time_grid, freq_grid, canonical_coherence[..., 0, 1].squeeze().T,
vmin=0, vmax=0.5, cmap='viridis')
plt.xlabel('Time [s]')
plt.ylabel('Frequency [Hz]')
plt.suptitle('Canonical Coherence', y=1.02, fontsize=30)
cb = fig.colorbar(mesh, ax=plt.gca(), orientation='horizontal',
shrink=.5, aspect=15, pad=0.1, label='Coherence')
cb.outline.set_linewidth(0);
frequency resolution: 25.0
Global coherence finds the linear combinations of signals that maximizes the power at a given frequency.
from itertools import product
time_halfbandwidth_product = 2
frequency_of_interest = 200
sampling_frequency = 1500
time_extent = (0, 2.400)
n_trials = 100
n_signals = 6
n_time_samples = int(((time_extent[1] - time_extent[0]) * sampling_frequency) + 1)
time = np.linspace(time_extent[0], time_extent[1], num=n_time_samples, endpoint=True)
signal = np.zeros((n_time_samples, n_trials, n_signals))
signal[:, :, 0:2] = (
np.sin(2 * np.pi * time * frequency_of_interest)[:, np.newaxis, np.newaxis] *
np.ones((1, n_trials, 2)))
other_signals = (n_signals + 1) // 2
n_other_signals = n_signals - other_signals
phase_offset = np.random.uniform(-np.pi, np.pi, size=(n_time_samples, n_trials, n_other_signals))
phase_offset[np.where(time > 1.5), :] = np.pi / 2
signal[:, :, other_signals:] = np.sin(
(2 * np.pi * time[:, np.newaxis, np.newaxis] * frequency_of_interest) + phase_offset)
noise = np.random.normal(10, 7, signal.shape)
m = Multitaper(signal + noise,
sampling_frequency=sampling_frequency,
time_halfbandwidth_product=time_halfbandwidth_product,
time_window_duration=0.080,
time_window_step=0.080,
start_time=time[0])
c = Connectivity(fourier_coefficients=m.fft(),
frequencies=m.frequencies,
time=m.time)
global_coherence, unnormalized_global_coherence = c.global_coherence()
global_coherence.shape # n_time, n_frequencies, n_components
time_grid, freq_grid = np.meshgrid(
np.append(c.time, time_extent[-1]),
np.append(c.frequencies, m.nyquist_frequency))
plt.figure()
plt.pcolormesh(time_grid, freq_grid, global_coherence[:, m.frequencies >= 0, 0].T, shading='auto')
plt.title('Global Coherence (1st component)')
plt.xlabel('Time [s]')
plt.ylabel('Frequency [Hz]')
Text(0, 0.5, 'Frequency [Hz]')
The xarray interface provides three things:
coherence_magnitude = multitaper_connectivity(signal + noise,
sampling_frequency=sampling_frequency,
time_halfbandwidth_product=time_halfbandwidth_product,
time_window_duration=0.080,
time_window_step=0.080,
method='coherence_magnitude')
coherence_magnitude
<xarray.DataArray 'coherence_magnitude' (Time: 30, Frequency: 60, Source: 6, Target: 6)> array([[[[ nan, 2.08445759e-03, 6.38772392e-04, 7.25645676e-03, 2.18606451e-03, 3.57603936e-03], [2.08445759e-03, nan, 6.73329640e-03, 1.55010738e-03, 3.17180695e-04, 1.23702707e-03], [6.38772392e-04, 6.73329640e-03, nan, 3.43372943e-03, 3.67021000e-04, 6.25549185e-03], [7.25645676e-03, 1.55010738e-03, 3.43372943e-03, nan, 4.33401517e-03, 9.88034432e-03], [2.18606451e-03, 3.17180695e-04, 3.67021000e-04, 4.33401517e-03, nan, 4.74304847e-06], [3.57603936e-03, 1.23702707e-03, 6.25549185e-03, 9.88034432e-03, 4.74304847e-06, nan]], [[ nan, 1.95792847e-03, 4.78595112e-03, 6.32804088e-04, 9.75904756e-04, 2.22267229e-03], [1.95792847e-03, nan, 3.27527073e-03, 3.04599547e-03, 4.52339920e-03, 2.81317846e-03], [4.78595112e-03, 3.27527073e-03, nan, 9.95100471e-04, 2.22317542e-03, 2.50317006e-04], [6.32804088e-04, 3.04599547e-03, 9.95100471e-04, ... 5.48440707e-03, 1.76854804e-03, 6.92467748e-03], [6.43278006e-03, 7.71744653e-03, 5.48440707e-03, nan, 6.09645947e-04, 8.70175573e-04], [4.78705515e-03, 1.66259436e-03, 1.76854804e-03, 6.09645947e-04, nan, 1.76529904e-03], [2.40235883e-04, 1.12021124e-02, 6.92467748e-03, 8.70175573e-04, 1.76529904e-03, nan]], [[ nan, 1.01926980e-02, 5.71676520e-04, 3.32848656e-03, 3.38782168e-03, 1.10414148e-03], [1.01926980e-02, nan, 3.67052326e-03, 2.83449004e-03, 1.28761927e-03, 2.13895709e-03], [5.71676520e-04, 3.67052326e-03, nan, 3.77439290e-03, 5.20035000e-03, 2.51114305e-03], [3.32848656e-03, 2.83449004e-03, 3.77439290e-03, nan, 4.20354885e-04, 1.13101329e-03], [3.38782168e-03, 1.28761927e-03, 5.20035000e-03, 4.20354885e-04, nan, 2.23638303e-05], [1.10414148e-03, 2.13895709e-03, 2.51114305e-03, 1.13101329e-03, 2.23638303e-05, nan]]]]) Coordinates: * Time (Time) float64 0.0 0.08 0.16 0.24 0.32 ... 2.08 2.16 2.24 2.32 * Frequency (Frequency) float64 0.0 12.5 25.0 37.5 ... 712.5 725.0 737.5 * Source (Source) int64 0 1 2 3 4 5 * Target (Target) int64 0 1 2 3 4 5 Attributes: (12/15) mt_detrend_type: constant mt_frequency_resolution: 25.0 mt_is_low_bias: True mt_n_fft_samples: 120 mt_n_signals: 6 mt_n_tapers: 3 ... ... mt_nyquist_frequency: 750.0 mt_sampling_frequency: 1500 mt_start_time: 0 mt_time_halfbandwidth_product: 2 mt_time_window_duration: 0.08 mt_time_window_step: 0.08
array([[[[ nan, 2.08445759e-03, 6.38772392e-04, 7.25645676e-03, 2.18606451e-03, 3.57603936e-03], [2.08445759e-03, nan, 6.73329640e-03, 1.55010738e-03, 3.17180695e-04, 1.23702707e-03], [6.38772392e-04, 6.73329640e-03, nan, 3.43372943e-03, 3.67021000e-04, 6.25549185e-03], [7.25645676e-03, 1.55010738e-03, 3.43372943e-03, nan, 4.33401517e-03, 9.88034432e-03], [2.18606451e-03, 3.17180695e-04, 3.67021000e-04, 4.33401517e-03, nan, 4.74304847e-06], [3.57603936e-03, 1.23702707e-03, 6.25549185e-03, 9.88034432e-03, 4.74304847e-06, nan]], [[ nan, 1.95792847e-03, 4.78595112e-03, 6.32804088e-04, 9.75904756e-04, 2.22267229e-03], [1.95792847e-03, nan, 3.27527073e-03, 3.04599547e-03, 4.52339920e-03, 2.81317846e-03], [4.78595112e-03, 3.27527073e-03, nan, 9.95100471e-04, 2.22317542e-03, 2.50317006e-04], [6.32804088e-04, 3.04599547e-03, 9.95100471e-04, ... 5.48440707e-03, 1.76854804e-03, 6.92467748e-03], [6.43278006e-03, 7.71744653e-03, 5.48440707e-03, nan, 6.09645947e-04, 8.70175573e-04], [4.78705515e-03, 1.66259436e-03, 1.76854804e-03, 6.09645947e-04, nan, 1.76529904e-03], [2.40235883e-04, 1.12021124e-02, 6.92467748e-03, 8.70175573e-04, 1.76529904e-03, nan]], [[ nan, 1.01926980e-02, 5.71676520e-04, 3.32848656e-03, 3.38782168e-03, 1.10414148e-03], [1.01926980e-02, nan, 3.67052326e-03, 2.83449004e-03, 1.28761927e-03, 2.13895709e-03], [5.71676520e-04, 3.67052326e-03, nan, 3.77439290e-03, 5.20035000e-03, 2.51114305e-03], [3.32848656e-03, 2.83449004e-03, 3.77439290e-03, nan, 4.20354885e-04, 1.13101329e-03], [3.38782168e-03, 1.28761927e-03, 5.20035000e-03, 4.20354885e-04, nan, 2.23638303e-05], [1.10414148e-03, 2.13895709e-03, 2.51114305e-03, 1.13101329e-03, 2.23638303e-05, nan]]]])
array([0. , 0.08, 0.16, 0.24, 0.32, 0.4 , 0.48, 0.56, 0.64, 0.72, 0.8 , 0.88, 0.96, 1.04, 1.12, 1.2 , 1.28, 1.36, 1.44, 1.52, 1.6 , 1.68, 1.76, 1.84, 1.92, 2. , 2.08, 2.16, 2.24, 2.32])
array([ 0. , 12.5, 25. , 37.5, 50. , 62.5, 75. , 87.5, 100. , 112.5, 125. , 137.5, 150. , 162.5, 175. , 187.5, 200. , 212.5, 225. , 237.5, 250. , 262.5, 275. , 287.5, 300. , 312.5, 325. , 337.5, 350. , 362.5, 375. , 387.5, 400. , 412.5, 425. , 437.5, 450. , 462.5, 475. , 487.5, 500. , 512.5, 525. , 537.5, 550. , 562.5, 575. , 587.5, 600. , 612.5, 625. , 637.5, 650. , 662.5, 675. , 687.5, 700. , 712.5, 725. , 737.5])
array([0, 1, 2, 3, 4, 5])
array([0, 1, 2, 3, 4, 5])
coherence_magnitude.plot(col='Source', row='Target', x='Time');