Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

add Welch method #269

Merged
merged 1 commit into from
Dec 11, 2023
Merged

add Welch method #269

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
6 changes: 4 additions & 2 deletions py_neuromodulation/nm_features.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -11,7 +11,7 @@
nm_oscillatory,
nm_bursts,
nm_linelength,
nm_bispectra
nm_bispectra,
)


Expand Down Expand Up @@ -56,6 +56,8 @@ def __init__(
FeatureClass = nm_oscillatory.STFT
case "fft":
FeatureClass = nm_oscillatory.FFT
case "welch":
FeatureClass = nm_oscillatory.Welch
case "sharpwave_analysis":
FeatureClass = nm_sharpwaves.SharpwaveAnalyzer
case "fooof":
Expand Down Expand Up @@ -111,4 +113,4 @@ def estimate_features(self, data: np.ndarray) -> dict:
features_compute,
)

return features_compute
return features_compute
70 changes: 64 additions & 6 deletions py_neuromodulation/nm_oscillatory.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -137,7 +137,7 @@ def calc_feature(self, data: np.ndarray, features_compute: dict) -> dict:
feature_calc = np.mean(Z_ch)

if self.log_transform:
feature_calc = np.log(feature_calc)
feature_calc = np.log10(feature_calc)

if self.KF_dict:
feature_calc = self.update_KF(feature_calc, feature_name)
Expand All @@ -146,6 +146,57 @@ def calc_feature(self, data: np.ndarray, features_compute: dict) -> dict:
return features_compute


class Welch(OscillatoryFeature):
def __init__(
self,
settings: dict,
ch_names: Iterable[str],
sfreq: float,
) -> None:
super().__init__(settings, ch_names, sfreq)
if self.s["welch_settings"]["kalman_filter"]:
self.init_KF("welch")

self.log_transform = self.s["welch_settings"]["log_transform"]

self.feature_params = []
for ch_idx, ch_name in enumerate(self.ch_names):
for fband, f_range in self.f_ranges_dict.items():
feature_name = "_".join([ch_name, "welch", fband])
self.feature_params.append((ch_idx, feature_name, f_range))

@staticmethod
def test_settings(s: dict, ch_names: Iterable[str], sfreq: int | float):
OscillatoryFeature.test_settings_osc(
s, ch_names, sfreq, "welch_settings"
)

def calc_feature(self, data: np.ndarray, features_compute: dict) -> dict:
f, Z = signal.welch(
data,
fs=self.sfreq,
window="hann",
nperseg=self.sfreq,
noverlap=None,
)

for ch_idx, feature_name, f_range in self.feature_params:
Z_ch = Z[ch_idx]

if self.log_transform:
Z_ch = np.log10(Z_ch)

idx_range = np.where((f >= f_range[0]) & (f <= f_range[1]))[0]
feature_calc = np.mean(Z_ch[idx_range])

if self.KF_dict:
feature_calc = self.update_KF(feature_calc, feature_name)

features_compute[feature_name] = feature_calc

return features_compute


class STFT(OscillatoryFeature):
def __init__(
self,
Expand All @@ -158,6 +209,7 @@ def __init__(
self.init_KF("stft")

self.nperseg = int(self.s["stft_settings"]["windowlength_ms"])
self.log_transform = self.s["stft_settings"]["log_transform"]

self.feature_params = []
for ch_idx, ch_name in enumerate(self.ch_names):
Expand Down Expand Up @@ -188,6 +240,9 @@ def calc_feature(self, data: np.ndarray, features_compute: dict) -> dict:
if self.KF_dict:
feature_calc = self.update_KF(feature_calc, feature_name)

if self.log_transform:
feature_calc = np.log10(feature_calc)

features_compute[feature_name] = feature_calc

return features_compute
Expand Down Expand Up @@ -265,7 +320,9 @@ def test_settings(s: dict, ch_names: Iterable[str], sfreq: int | float):
].values()
), "Set at least one bandpower_feature to True."

for fband_name, seg_length_fband in s["bandpass_filter_settings"]["segment_lengths_ms"].items():
for fband_name, seg_length_fband in s["bandpass_filter_settings"][
"segment_lengths_ms"
].items():
assert isinstance(seg_length_fband, int), (
f"bandpass segment_lengths_ms for {fband_name} "
f"needs to be of type int, got {seg_length_fband}"
Expand All @@ -275,15 +332,16 @@ def test_settings(s: dict, ch_names: Iterable[str], sfreq: int | float):
f"segment length {seg_length_fband} needs to be smaller than "
f" s['segment_length_features_ms'] = {s['segment_length_features_ms']}"
)

for fband_name in list(s["frequency_ranges_hz"].keys()):
assert fband_name in list(s["bandpass_filter_settings"]["segment_lengths_ms"].keys()), (
assert fband_name in list(
s["bandpass_filter_settings"]["segment_lengths_ms"].keys()
), (
f"frequency range {fband_name} "
"needs to be defined in s['bandpass_filter_settings']['segment_lengths_ms']"
)

def calc_feature(self, data: np.ndarray, features_compute: dict) -> dict:

data = self.bandpass_filter.filter_data(data)

for (
Expand All @@ -297,7 +355,7 @@ def calc_feature(self, data: np.ndarray, features_compute: dict) -> dict:
) in self.feature_params:
if bp_feature == "activity":
if self.log_transform:
feature_calc = np.log(
feature_calc = np.log10(
np.var(data[ch_idx, f_band_idx, -seglen:])
)
else:
Expand Down
6 changes: 6 additions & 0 deletions py_neuromodulation/nm_settings.json
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -18,6 +18,7 @@
"bandpass_filter": false,
"stft": false,
"fft": true,
"welch": true,
"sharpwave_analysis": true,
"fooof": false,
"bursts": true,
Expand Down Expand Up @@ -90,6 +91,11 @@
"log_transform": true,
"kalman_filter": false
},
"welch_settings": {
"windowlength_ms": 1000,
"log_transform": true,
"kalman_filter": false
},
"stft_settings": {
"windowlength_ms": 500,
"log_transform": true,
Expand Down
Loading