davrot/pytutorial
1
0
Fork 0
mirror of https://github.com/davrot/pytutorial.git synced 2025-04-19 05:36:42 +02:00
Code Issues Projects Releases Packages Wiki Activity Actions

Update README.md

Browse source 
Signed-off-by: David Rotermund <54365609+davrot@users.noreply.github.com>
This commit is contained in:
David Rotermund 2024-02-15 15:41:14 +01:00 • committed by GitHub
parent 391fc432b4
commit 21e33d248d
No known key found for this signature in database
GPG key ID: B5690EEEBB952194
1 changed files with 69 additions and 64 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

117
data_analysis/spectral_coherence/README.md
View file

@ -193,7 +193,6 @@ plt.show()
import numpy as np
import matplotlib.pyplot as plt
import pywt
from tqdm import trange
# Calculate the wavelet scales we requested
@ -301,44 +300,17 @@ def calculate_wavelet_tf_complex_coeffs(
return (complex_spectrum, frequency_axis, t)
# Parameters for the wavelet transform
number_of_frequences: int = 25 # frequency bands
frequency_range_min: float = 5 # Hz
frequency_range_max: float = 200 # Hz
dt: float = 1.0 / 1000.0
def calculate_spectral_coherence(
n_trials: int,
y_a: np.ndarray,
y_b: np.ndarray,
number_of_frequences: int,
frequency_range_min: float,
frequency_range_max: float,
dt: float,
) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
# I want more trials
f_base: float = 50
f_delta: float = 1
delay_enable: bool = False
# Test data ->
rng = np.random.default_rng(1)
n_t: int = 1000
n_trials: int = 100
t: np.ndarray = np.arange(0, n_t) * dt
amplitude: float = 0.75
x = dt * 2.0 * np.pi * (f_base + f_delta * 2 * (rng.random((n_t, n_trials)) - 0.5))
x = np.cumsum(x, axis=0)
y_a: np.ndarray = np.sin(x)
y_a[:400, :] = 0.0
y_a[-400:, :] = 0.0
y_a = y_a + amplitude * rng.random((n_t, n_trials))
y_a -= y_a.mean(axis=0, keepdims=True)
y_a /= y_a.std(axis=0, keepdims=True)
# <- Test data
if delay_enable:
y_b: np.ndarray = np.roll(y_a, shift=250, axis=0)
else:
y_b = y_a.copy()
for trial_id in trange(0, n_trials):
for trial_id in range(0, n_trials):
wave_data_a, frequency_axis, t = calculate_wavelet_tf_complex_coeffs(
data=y_a[..., trial_id],
number_of_frequences=number_of_frequences,
@ -379,31 +351,64 @@ for trial_id in trange(0, n_trials):
norm_data_a += np.abs(wave_data_a) ** 2
norm_data_b += np.abs(wave_data_b) ** 2
calculation /= float(n_trials)
norm_data_a /= float(n_trials)
norm_data_b /= float(n_trials)
calculation /= float(n_trials)
norm_data_a /= float(n_trials)
norm_data_b /= float(n_trials)
coherence = np.abs(calculation) ** 2 / (norm_data_a * norm_data_b)
coherence = np.abs(calculation) ** 2 / (norm_data_a * norm_data_b)
y_reduction_to_ticks: int = 10
x_reduction_to_ticks: int = 10
y_round: int = 1
x_round: int = 1
return np.nanmean(coherence, axis=-1), frequency_axis, t
freq_ticks, freq_values = get_y_ticks(
reduction_to_ticks=y_reduction_to_ticks,
frequency_axis=frequency_axis,
round=y_round,
)
time_ticks, time_values = get_x_ticks(
reduction_to_ticks=x_reduction_to_ticks,
# Parameters for the wavelet transform
number_of_frequences: int = 25 # frequency bands
frequency_range_min: float = 5 # Hz
frequency_range_max: float = 200 # Hz
dt: float = 1.0 / 1000.0
# I want more trials
f_base: float = 50
f_delta: float = 1
delay_enable: bool = False
# Test data ->
rng = np.random.default_rng(1)
n_t: int = 1000
n_trials: int = 100
t: np.ndarray = np.arange(0, n_t) * dt
amplitude: float = 0.75
x = dt * 2.0 * np.pi * (f_base + f_delta * 2 * (rng.random((n_t, n_trials)) - 0.5))
x = np.cumsum(x, axis=0)
y_a: np.ndarray = np.sin(x)
y_a[:400, :] = 0.0
y_a[-400:, :] = 0.0
y_a = y_a + amplitude * rng.random((n_t, n_trials))
y_a -= y_a.mean(axis=0, keepdims=True)
y_a /= y_a.std(axis=0, keepdims=True)
# <- Test data
if delay_enable:
y_b: np.ndarray = np.roll(y_a, shift=250, axis=0)
else:
y_b = y_a.copy()
coherence, frequency_axis, t = calculate_spectral_coherence(
n_trials=n_trials,
y_a=y_a,
y_b=y_b,
number_of_frequences=number_of_frequences,
frequency_range_min=frequency_range_min,
frequency_range_max=frequency_range_max,
dt=dt,
number_of_timesteps=t.shape[0],
round=x_round,
)
plt.plot(frequency_axis, np.nanmean(coherence, axis=-1))
plt.plot(frequency_axis, coherence)
plt.ylabel("Spectral Coherence")
plt.xlabel("Frequency [Hz]")
plt.show()