From 753559e7d2d9710ffc816ed995fccadaf0755108 Mon Sep 17 00:00:00 2001 From: DRomBas Date: Wed, 10 Mar 2021 12:22:24 +0100 Subject: [PATCH 1/9] add time-domain trigger detection --- phys2bids/physio_obj.py | 26 ++++++++++++++++++-------- 1 file changed, 18 insertions(+), 8 deletions(-) diff --git a/phys2bids/physio_obj.py b/phys2bids/physio_obj.py index 0eea26765..c960e4766 100644 --- a/phys2bids/physio_obj.py +++ b/phys2bids/physio_obj.py @@ -245,10 +245,6 @@ def __init__(self, timeseries, freq, ch_name, units, trigger_idx, self.time_offset = deepcopy(time_offset) if trigger_idx == 0: self.auto_trigger_selection() - else: - if ch_name[trigger_idx] not in TRIGGER_NAMES: - LGR.info('Trigger channel name is not in our trigger channel name alias list. ' - 'Please make sure you choose the proper channel.') @property def ch_amount(self): @@ -594,11 +590,25 @@ def auto_trigger_selection(self): 'Please run phys2bids specifying the -chtrig argument.') else: self.trigger_idx = indexes[0] - LGR.info(f'{self.ch_name[self.trigger_idx]} selected as trigger channel') else: - raise Exception('No trigger channel automatically found. Please run phys2bids ' - 'specifying the -chtrig argument.') - + # Time-domain automatic trigger detection + # Initialize distance array + mean_d = [np.nan]*(self.ch_amount-1) + # Loop through channels + for n in range(1, self.ch_amount): + # Get timeseries + s = self.timeseries[n] + # Normalize to [0,1] + s = (s-min(s))/(max(s)-min(s)) + # Calculate distance to the closest signal limit (min or max) + d = np.minimum(abs(s-max(s)), abs(s-min(s))) + # Store the mean distance + mean_d[n-1] = np.mean(d) + + # Set the trigger as the channel with smaller distance + self.trigger_idx = np.argmin(mean_d) + 1 + + LGR.info(f'{self.ch_name[self.trigger_idx]} selected as trigger channel') class BlueprintOutput(): """ From 381bfe0d480df105c4e3b2deb2b07dfc01cfbb21 Mon Sep 17 00:00:00 2001 From: DRomBas Date: Wed, 10 Mar 2021 12:47:59 +0100 Subject: [PATCH 2/9] add basic channel classification --- phys2bids/phys2bids.py | 4 ++++ phys2bids/physio_obj.py | 32 ++++++++++++++++++++++++++++++++ 2 files changed, 36 insertions(+) diff --git a/phys2bids/phys2bids.py b/phys2bids/phys2bids.py index 3d5ab6495..c6303d0f1 100644 --- a/phys2bids/phys2bids.py +++ b/phys2bids/phys2bids.py @@ -252,6 +252,10 @@ def phys2bids(filename, info=False, indir='.', outdir='.', heur_file=None, if ch_name: LGR.info('Renaming channels with given names') phys_in.rename_channels(ch_name) + else: + LGR.info('Renaming channels automatically') + phys_in.auto_rename_channels() + # Checking acquisition type via user's input if tr is not None and num_timepoints_expected is not None: diff --git a/phys2bids/physio_obj.py b/phys2bids/physio_obj.py index c960e4766..934d277e9 100644 --- a/phys2bids/physio_obj.py +++ b/phys2bids/physio_obj.py @@ -610,6 +610,38 @@ def auto_trigger_selection(self): LGR.info(f'{self.ch_name[self.trigger_idx]} selected as trigger channel') + def auto_rename_channels(self): + LGR.info('Running automatic channel detection.') + + # Loop through channels excluding the trigger + for n in range(1,self.ch_amount): + # Skip trigger + if n == self.trigger_idx: + continue + + # Get timeseries and sampling frequency + s = self.timeseries[n] + T = 1/self.freq[n] + + # Remove DC component + s = s - np.mean(s) + + # Compute power spectrum + sf = np.power(np.abs(np.fft.fft(s)), 2) + f = np.array(np.fft.fftfreq(len(s), d=T)) + + # Get power in respiratory band + pow_resp = sum(sf[(f > 0) & (f < 0.5)]) + + # Get power in cardiac band + pow_cardiac = sum(sf[(f>0.5) & (f<4)]) + + # Classify as cardiac or respiratory + if pow_resp > pow_cardiac: + self.ch_name[n] = 'respiratory' + else: + self.ch_name[n] = 'cardiac' + class BlueprintOutput(): """ Main output object for phys2bids. From f327777fad7b51518c9ae39b9c8119c58772bac6 Mon Sep 17 00:00:00 2001 From: dombas Date: Wed, 7 Apr 2021 09:22:40 +0200 Subject: [PATCH 3/9] configure only trigger automatic detection --- phys2bids/phys2bids.py | 4 --- phys2bids/physio_obj.py | 45 ++++++------------------------ phys2bids/tests/test_physio_obj.py | 6 ---- 3 files changed, 8 insertions(+), 47 deletions(-) diff --git a/phys2bids/phys2bids.py b/phys2bids/phys2bids.py index c6303d0f1..3d5ab6495 100644 --- a/phys2bids/phys2bids.py +++ b/phys2bids/phys2bids.py @@ -252,10 +252,6 @@ def phys2bids(filename, info=False, indir='.', outdir='.', heur_file=None, if ch_name: LGR.info('Renaming channels with given names') phys_in.rename_channels(ch_name) - else: - LGR.info('Renaming channels automatically') - phys_in.auto_rename_channels() - # Checking acquisition type via user's input if tr is not None and num_timepoints_expected is not None: diff --git a/phys2bids/physio_obj.py b/phys2bids/physio_obj.py index 934d277e9..90d93ca5c 100644 --- a/phys2bids/physio_obj.py +++ b/phys2bids/physio_obj.py @@ -245,6 +245,10 @@ def __init__(self, timeseries, freq, ch_name, units, trigger_idx, self.time_offset = deepcopy(time_offset) if trigger_idx == 0: self.auto_trigger_selection() + else: + if ch_name[trigger_idx] not in TRIGGER_NAMES: + LGR.info('Trigger channel name is not in our trigger channel name alias list. ' + 'Please make sure you choose the proper channel.') @property def ch_amount(self): @@ -554,7 +558,8 @@ def auto_trigger_selection(self): Find a trigger index matching the channels with a regular expresion. It compares the channel name with the the regular expressions stored - in TRIGGER_NAMES. + in TRIGGER_NAMES. If that fails a time-domain recognition of the + trigger signal is performed. Parameters ---------- @@ -565,9 +570,6 @@ def auto_trigger_selection(self): Exception More than one possible trigger channel was automatically found. - Exception - No trigger channel automatically found - Notes ----- Outcome: @@ -599,9 +601,9 @@ def auto_trigger_selection(self): # Get timeseries s = self.timeseries[n] # Normalize to [0,1] - s = (s-min(s))/(max(s)-min(s)) + s = (s - min(s))/(max(s) - min(s)) # Calculate distance to the closest signal limit (min or max) - d = np.minimum(abs(s-max(s)), abs(s-min(s))) + d = np.minimum(abs(s - max(s)), abs(s - min(s))) # Store the mean distance mean_d[n-1] = np.mean(d) @@ -610,37 +612,6 @@ def auto_trigger_selection(self): LGR.info(f'{self.ch_name[self.trigger_idx]} selected as trigger channel') - def auto_rename_channels(self): - LGR.info('Running automatic channel detection.') - - # Loop through channels excluding the trigger - for n in range(1,self.ch_amount): - # Skip trigger - if n == self.trigger_idx: - continue - - # Get timeseries and sampling frequency - s = self.timeseries[n] - T = 1/self.freq[n] - - # Remove DC component - s = s - np.mean(s) - - # Compute power spectrum - sf = np.power(np.abs(np.fft.fft(s)), 2) - f = np.array(np.fft.fftfreq(len(s), d=T)) - - # Get power in respiratory band - pow_resp = sum(sf[(f > 0) & (f < 0.5)]) - - # Get power in cardiac band - pow_cardiac = sum(sf[(f>0.5) & (f<4)]) - - # Classify as cardiac or respiratory - if pow_resp > pow_cardiac: - self.ch_name[n] = 'respiratory' - else: - self.ch_name[n] = 'cardiac' class BlueprintOutput(): """ diff --git a/phys2bids/tests/test_physio_obj.py b/phys2bids/tests/test_physio_obj.py index fb6d61f53..22dbf15b3 100644 --- a/phys2bids/tests/test_physio_obj.py +++ b/phys2bids/tests/test_physio_obj.py @@ -284,12 +284,6 @@ def test_auto_trigger_selection(caplog): test_units, test_chtrig) assert phys_in.trigger_idx == 1 # test when no trigger is found - test_chn_name = ['time', 'TRIGGAH', 'half', 'CO2', 'CO 2', 'strigose'] - with raises(Exception) as errorinfo: - phys_in = po.BlueprintInput(test_timeseries, test_freq, test_chn_name, - test_units, test_chtrig) - assert 'No trigger channel automatically found' in str(errorinfo.value) - # test when no trigger is found test_chn_name = ['time', 'trigger', 'TRIGGER', 'CO2', 'CO 2', 'strigose'] with raises(Exception) as errorinfo: phys_in = po.BlueprintInput(test_timeseries, test_freq, test_chn_name, From 23c49d8348bbc402e8a1fa28f77cee171b3f5906 Mon Sep 17 00:00:00 2001 From: dombas Date: Wed, 7 Apr 2021 10:45:33 +0200 Subject: [PATCH 4/9] TEST:add test for time-domain recognition of the trigger --- phys2bids/physio_obj.py | 3 ++- phys2bids/tests/test_physio_obj.py | 30 +++++++++++++++++++++++++++++- 2 files changed, 31 insertions(+), 2 deletions(-) diff --git a/phys2bids/physio_obj.py b/phys2bids/physio_obj.py index 90d93ca5c..43b85d645 100644 --- a/phys2bids/physio_obj.py +++ b/phys2bids/physio_obj.py @@ -555,7 +555,8 @@ def print_info(self, filename): def auto_trigger_selection(self): """ - Find a trigger index matching the channels with a regular expresion. + Find a trigger index matching the channels with a regular expresion or + by using a time-domain recognition of the trigger. It compares the channel name with the the regular expressions stored in TRIGGER_NAMES. If that fails a time-domain recognition of the diff --git a/phys2bids/tests/test_physio_obj.py b/phys2bids/tests/test_physio_obj.py index 22dbf15b3..7f44e506c 100644 --- a/phys2bids/tests/test_physio_obj.py +++ b/phys2bids/tests/test_physio_obj.py @@ -261,7 +261,7 @@ def test_BlueprintOutput(): assert blueprint_out == blueprint_out -def test_auto_trigger_selection(caplog): +def test_auto_trigger_selection_text(caplog): """Test auto_trigger_selection.""" test_time = np.array([0, 1, 2, 3, 4]) test_trigger = np.array([0, 1, 2, 3, 4]) @@ -289,3 +289,31 @@ def test_auto_trigger_selection(caplog): phys_in = po.BlueprintInput(test_timeseries, test_freq, test_chn_name, test_units, test_chtrig) assert 'More than one possible trigger channel' in str(errorinfo.value) + + +def test_auto_trigger_selection_time(caplog): + """Test auto_trigger_selection in time domain.""" + # Simulate 10 s of a trigger, O2 and ECG + T = 10 + nSamp = 100 + fs = nSamp/T + test_time = np.linspace(0, T, nSamp) + test_freq = [fs, fs, fs, fs] + + # O2 as a sinusoidal of 0.5 Hz + test_O2 = np.sin(2*np.pi*0.5*test_time) + # ECG as a sinusoidal with 1.5 Hz + test_ecg = np.sin(2*np.pi*1.5*test_time) + # Trigger as a binary signal + test_trigger = np.zeros(nSamp) + test_trigger[1:nSamp:4] = 1 + + test_timeseries = [test_time, test_O2, test_ecg, test_trigger] + test_chn_name = ['time', 'O2', 'ecg', 'tiger'] + test_units = ['s', 'V', 'V', 'V'] + + # test when trigger is 0 and the trigger is not recognized by text matching: + test_chtrig = 0 + phys_in = po.BlueprintInput(test_timeseries, test_freq, test_chn_name, + test_units, test_chtrig) + assert phys_in.trigger_idx == 3 From b06eb3c0db7b7f6e395ddc2d25bf21c13d486673 Mon Sep 17 00:00:00 2001 From: dombas Date: Wed, 7 Apr 2021 11:21:55 +0200 Subject: [PATCH 5/9] STYLE: fix lintin issues --- phys2bids/physio_obj.py | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/phys2bids/physio_obj.py b/phys2bids/physio_obj.py index 43b85d645..4c4c60366 100644 --- a/phys2bids/physio_obj.py +++ b/phys2bids/physio_obj.py @@ -555,8 +555,7 @@ def print_info(self, filename): def auto_trigger_selection(self): """ - Find a trigger index matching the channels with a regular expresion or - by using a time-domain recognition of the trigger. + Find a trigger index automatically. It compares the channel name with the the regular expressions stored in TRIGGER_NAMES. If that fails a time-domain recognition of the @@ -596,17 +595,17 @@ def auto_trigger_selection(self): else: # Time-domain automatic trigger detection # Initialize distance array - mean_d = [np.nan]*(self.ch_amount-1) + mean_d = [np.nan] * (self.ch_amount - 1) # Loop through channels for n in range(1, self.ch_amount): # Get timeseries s = self.timeseries[n] # Normalize to [0,1] - s = (s - min(s))/(max(s) - min(s)) + s = (s - min(s)) / (max(s) - min(s)) # Calculate distance to the closest signal limit (min or max) d = np.minimum(abs(s - max(s)), abs(s - min(s))) # Store the mean distance - mean_d[n-1] = np.mean(d) + mean_d[n - 1] = np.mean(d) # Set the trigger as the channel with smaller distance self.trigger_idx = np.argmin(mean_d) + 1 From 00661247dd98a669dd2afb9d22efb25bedbc750a Mon Sep 17 00:00:00 2001 From: dombas Date: Thu, 15 Apr 2021 18:49:58 +0200 Subject: [PATCH 6/9] TEST: relaunch CI --- phys2bids/tests/test_physio_obj.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/phys2bids/tests/test_physio_obj.py b/phys2bids/tests/test_physio_obj.py index 7f44e506c..9e4104540 100644 --- a/phys2bids/tests/test_physio_obj.py +++ b/phys2bids/tests/test_physio_obj.py @@ -291,7 +291,7 @@ def test_auto_trigger_selection_text(caplog): assert 'More than one possible trigger channel' in str(errorinfo.value) -def test_auto_trigger_selection_time(caplog): +def test_auto_trigger_selection_time(): """Test auto_trigger_selection in time domain.""" # Simulate 10 s of a trigger, O2 and ECG T = 10 From 296da3d7d7c1d2a91b56b15ea07a4e640b4eab31 Mon Sep 17 00:00:00 2001 From: dombas Date: Fri, 16 Apr 2021 18:49:58 +0200 Subject: [PATCH 7/9] RF: replace for loop by numpy array computation --- phys2bids/physio_obj.py | 30 +++++++++++++++--------------- phys2bids/tests/test_physio_obj.py | 2 +- 2 files changed, 16 insertions(+), 16 deletions(-) diff --git a/phys2bids/physio_obj.py b/phys2bids/physio_obj.py index 4c4c60366..0b8481275 100644 --- a/phys2bids/physio_obj.py +++ b/phys2bids/physio_obj.py @@ -594,21 +594,21 @@ def auto_trigger_selection(self): self.trigger_idx = indexes[0] else: # Time-domain automatic trigger detection - # Initialize distance array - mean_d = [np.nan] * (self.ch_amount - 1) - # Loop through channels - for n in range(1, self.ch_amount): - # Get timeseries - s = self.timeseries[n] - # Normalize to [0,1] - s = (s - min(s)) / (max(s) - min(s)) - # Calculate distance to the closest signal limit (min or max) - d = np.minimum(abs(s - max(s)), abs(s - min(s))) - # Store the mean distance - mean_d[n - 1] = np.mean(d) - - # Set the trigger as the channel with smaller distance - self.trigger_idx = np.argmin(mean_d) + 1 + + # Create numpy array with all channels (excluding time) + channel_ts = np.array(self.timeseries[1:]) + + # Normalize each signal to [0,1] + min_ts = np.min(channel_ts, axis=1)[:, None] + max_ts = np.max(channel_ts, axis=1)[:, None] + channel_ts = (channel_ts - min_ts)/(max_ts - min_ts) + + # Compute distance to the closest signal limit (0 or 1) + distance = np.minimum(abs(channel_ts - 0), abs(channel_ts - 1)) + distance_mean = np.mean(distance, axis=1) + + # Set the trigger as the channel with the smallest distance + self.trigger_idx = np.argmin(distance_mean) + 1 LGR.info(f'{self.ch_name[self.trigger_idx]} selected as trigger channel') diff --git a/phys2bids/tests/test_physio_obj.py b/phys2bids/tests/test_physio_obj.py index 9e4104540..2c5f6dfcc 100644 --- a/phys2bids/tests/test_physio_obj.py +++ b/phys2bids/tests/test_physio_obj.py @@ -312,7 +312,7 @@ def test_auto_trigger_selection_time(): test_chn_name = ['time', 'O2', 'ecg', 'tiger'] test_units = ['s', 'V', 'V', 'V'] - # test when trigger is 0 and the trigger is not recognized by text matching: + # test when chtrig is 0 and the trigger is not recognized by text matching: test_chtrig = 0 phys_in = po.BlueprintInput(test_timeseries, test_freq, test_chn_name, test_units, test_chtrig) From 82e306c64349e79662175a5032384f1b3db78124 Mon Sep 17 00:00:00 2001 From: dombas Date: Fri, 16 Apr 2021 18:53:15 +0200 Subject: [PATCH 8/9] STYLE: add missing whitespace --- phys2bids/physio_obj.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/phys2bids/physio_obj.py b/phys2bids/physio_obj.py index 0b8481275..faeb190e3 100644 --- a/phys2bids/physio_obj.py +++ b/phys2bids/physio_obj.py @@ -601,7 +601,7 @@ def auto_trigger_selection(self): # Normalize each signal to [0,1] min_ts = np.min(channel_ts, axis=1)[:, None] max_ts = np.max(channel_ts, axis=1)[:, None] - channel_ts = (channel_ts - min_ts)/(max_ts - min_ts) + channel_ts = (channel_ts - min_ts) / (max_ts - min_ts) # Compute distance to the closest signal limit (0 or 1) distance = np.minimum(abs(channel_ts - 0), abs(channel_ts - 1)) From e860c5d26c9ac5f688be118d245b96c7c2294238 Mon Sep 17 00:00:00 2001 From: dombas Date: Tue, 11 May 2021 20:58:55 +0200 Subject: [PATCH 9/9] RF: make the method robust agains nan --- phys2bids/physio_obj.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/phys2bids/physio_obj.py b/phys2bids/physio_obj.py index faeb190e3..e6467fd06 100644 --- a/phys2bids/physio_obj.py +++ b/phys2bids/physio_obj.py @@ -608,7 +608,7 @@ def auto_trigger_selection(self): distance_mean = np.mean(distance, axis=1) # Set the trigger as the channel with the smallest distance - self.trigger_idx = np.argmin(distance_mean) + 1 + self.trigger_idx = np.nanargmin(distance_mean) + 1 LGR.info(f'{self.ch_name[self.trigger_idx]} selected as trigger channel')