diff --git a/phys2bids/cli/run.py b/phys2bids/cli/run.py index a7770171e..5ad83915b 100644 --- a/phys2bids/cli/run.py +++ b/phys2bids/cli/run.py @@ -1,7 +1,6 @@ # -*- coding: utf-8 -*- """Parser for phys2bids.""" - import argparse from phys2bids import __version__ @@ -134,6 +133,27 @@ def _get_parser(): "or just one TR if it is consistent throughout the session.", default=None, ) + optional.add_argument( + "-freq", + "--frequency", + dest="freq", + type=float, + help="Frequency of input data. Now specifically for GE gep files. " + "Units are Hertz, Hz or samples/sec " + "GE defaults to PPG of 100Hz, Resp at 25Hz and ECG at 1000Hz ", + default=None, + ) + + optional.add_argument( + "-pt", + "--pre_time", + dest="pretime", + type=float, + help="Pre-run time offset of input data. Now specifically for GE gep files. " + "Time offset in seconds. Default for GE is 30 seconds. ", + default=30.0, + ) + optional.add_argument( "-thr", "--threshold", diff --git a/phys2bids/io.py b/phys2bids/io.py index d4a5e3903..2efe57836 100644 --- a/phys2bids/io.py +++ b/phys2bids/io.py @@ -3,6 +3,7 @@ """phys2bids interfaces for loading extension files.""" import logging +import pdb import warnings from copy import deepcopy from itertools import groupby @@ -200,19 +201,22 @@ def read_header_and_channels(filename): """ header = [] # Read in the header until it's numbers + with open(filename, "r") as f: for n, line in enumerate(f): line = line.rstrip("\n").split("\t") if line[-1] == "": line.remove("") try: - float(line[0]) + list(map(float, line[0].split(" "))) + # float(line[0]) break except ValueError: header.append(line) continue + # Read in the rest paying attention to possible differences - if "Interval=" in header[0]: + if "Interval=" in str(header[0]): # Not specifying delimiters will ignore comments channel_list = np.genfromtxt(filename, skip_header=n) elif "acq" in header[0][0]: @@ -261,29 +265,31 @@ def extract_header_items(header): If Labchart headers cannot be processed If files are not in acq or txt format """ + # check header is not empty and detect if it is in labchart or Acqknoledge format if len(header) == 0: raise NotImplementedError("Files without header are not supported yet") - elif "Interval=" in header[0]: + elif "Interval=" in header[0][0]: LGR.info("phys2bids detected that your file is in Labchart format") interval = None orig_names = None range_list = None for line in header: - if "Interval=" in line: - interval = line[1].split(" ") - if "ChannelTitle=" in line: - orig_names = line[1:] - if "Range=" in line: - range_list = line[1:] + ll = line[0] + if "Interval=" in ll: + interval = ll.split()[1:] + if "ChannelTitle=" in ll: + orig_names = ll.split()[1:] + if "Range=" in ll: + range_list = ll.split()[1:] if None in [interval, orig_names, range_list]: raise NotImplementedError(OPEN_ISSUE) orig_units = [] for item in range_list: - orig_units.append(item.split(" ")[1]) + orig_units.append(item.split()[0]) elif "acq" in header[0][0]: LGR.info("phys2bids detected that your file is in AcqKnowledge format") @@ -450,7 +456,7 @@ def load_mat(filename, chtrig=0): return BlueprintInput(timeseries, freq, names, units, chtrig) -def load_gep(filename): +def load_gep(filename, inifreq=None, pretime=30.0): """ Populate object phys_input from GE physiological files. @@ -491,11 +497,15 @@ def load_gep(filename): # Initiate lists of column names and units with time and trigger names = ["time", "trigger"] - units = ["s", "mV"] # Assuming recording units are mV... + units = ["s", "mV", "au"] # Assuming recording units are mV... + + # start debugging here + # pdb.set_trace() # Add column for file given by user if "PPGData" in filename: freq = [100, 100, 100] + # freq = [50, 50, 50] names.append("cardiac") elif "RESPData" in filename: freq = [25, 25, 25] @@ -504,13 +514,18 @@ def load_gep(filename): freq = [1000, 1000, 1000] names.append("cardiac") + # if frequency specified in call, use that instead + if inifreq: + freq = [inifreq, inifreq, inifreq] + # Load in user file data data = [np.loadtxt(filename)] # Calculate time in seconds for first input (starts from -30s) interval = 1 / freq[0] duration = data[0].shape[0] * interval - t_ch = np.ogrid[-30 : duration - 30 : interval] + # t_ch = np.ogrid[-30 : duration - 30 : interval] + t_ch = np.ogrid[0:duration:interval] # Find and add additional data files filename = Path(filename) @@ -519,20 +534,33 @@ def load_gep(filename): if not len(fnames) == 0: for fname in fnames: if "PPGData" in fname: - freq.append(100) + if inifreq: + freq.append(inifreq) + else: + freq.append(100) names.append("cardiac") data.append(np.loadtxt(fname)) elif "RESPData" in fname: - freq.append(25) + if inifreq: + freq.append(inifreq) + else: + freq.append(25) names.append("respiratory") data.append(np.loadtxt(fname)) elif "ECGData" in fname: - freq.append(1000) + if inifreq: + freq.append(inifreq) + else: + freq.append(1000) names.append("cardiac") data.append(np.loadtxt(fname)) # Create trigger channel - trigger = np.hstack((np.zeros(int(30 / interval)), np.ones(int((duration - 30) / interval)))) + trigger = np.hstack( + (np.zeros(int(pretime / interval)), np.ones(int((duration - pretime) / interval))) + ) + # trigger = np.hstack((np.zeros(int(30 / interval)), np.ones(int((duration - 30) / interval)))) + # trigger = np.hstack( np.ones(int((duration) / interval))) # Create final list of timeseries timeseries = [t_ch, trigger] diff --git a/phys2bids/phys2bids.py b/phys2bids/phys2bids.py index 5161792bf..e8b5686dc 100644 --- a/phys2bids/phys2bids.py +++ b/phys2bids/phys2bids.py @@ -147,6 +147,8 @@ def phys2bids( yml="", debug=False, quiet=False, + freq=None, + pretime=None, ): """ Run main workflow of phys2bids. @@ -259,7 +261,7 @@ def phys2bids( elif ftype == "gep": from phys2bids.io import load_gep - phys_in = load_gep(infile) + phys_in = load_gep(infile, inifreq=freq, pretime=pretime) LGR.info("Checking that units of measure are BIDS compatible") for index, unit in enumerate(phys_in.units): diff --git a/phys2bids/physio_obj.py b/phys2bids/physio_obj.py index 131cb7ce9..6056ebb4e 100644 --- a/phys2bids/physio_obj.py +++ b/phys2bids/physio_obj.py @@ -518,13 +518,16 @@ def check_trigger_amount(self, thr=None, num_timepoints_expected=0, tr=0): if thr is None: # If trigger channels are binary # (i.e., "on" is a higher value and "off" is a lower value) - # and each "on" and "off" are each always approzimately the same value + # and each "on" and "off" are each always approximately the same value # then any value above the mean is "on" and every value below the mean # is "off". thr = np.mean(trigger) flag = 1 - timepoints = trigger > thr - num_timepoints_found = np.count_nonzero(np.ediff1d(timepoints.astype(np.int8)) > 0) + # drg - modified to >= rather than > + timepoints = trigger >= thr + # drg - also removed threshold difference that looks at consecutive diffs + # num_timepoints_found = np.count_nonzero(np.ediff1d(timepoints.astype(np.int8)) > 0) + num_timepoints_found = np.count_nonzero(timepoints.astype(np.int8) > 0) if flag == 1: LGR.info( f"The number of timepoints according to the std_thr method "