Added the pyphotometry preprocessing steps. There is a function called 'import_ppd' that is needed to import .ppd data. I included it in the description of the processppdphotometry function. Is this correct?

Author: j2figuer

src/jdb_to_nwb/convert_photometry.py

@@ -4,7 +4,7 @@
 import numpy as np
 import warnings
 import scipy.io
-from scipy.signal import butter, lfilter, hilbert
+from scipy.signal import butter, lfilter, hilbert, filtfilt
 from scipy.sparse import diags, eye, csc_matrix
 from scipy.sparse.linalg import spsolve
 from sklearn.linear_model import Lasso
@@ -336,6 +336,141 @@ def airPLS(data, lambda_=1e8, max_iterations=50):
         weights[-1] = weights[0]
     return baseline
 
+def processppdphotometry(ppd_file_path):
+    """
+    Process ppd file from pyPhotometry into a dict
+    requires the import_ppd function from the data_import.py script
+    The code needed is below.
+
+        def import_ppd(file_path, low_pass=20, high_pass=15.9):
+        Function to import pyPhotometry binary data files into Python. The high_pass
+        and low_pass arguments determine the frequency in Hz of highpass and lowpass
+        filtering applied to the filtered analog signals. To disable highpass or lowpass
+        filtering set the respective argument to None.  Returns a dictionary with the
+        following items:
+            'filename'      - Data filename
+            'subject_ID'    - Subject ID
+            'date_time'     - Recording start date and time (ISO 8601 format string)
+            'end_time'      - Recording end date and time (ISO 8601 format string)
+            'mode'          - Acquisition mode
+            'sampling_rate' - Sampling rate (Hz)
+            'LED_current'   - Current for LEDs 1 and 2 (mA)
+            'version'       - Version number of pyPhotometry
+            'analog_1'      - Raw analog signal 1 (volts)
+            'analog_2'      - Raw analog signal 2 (volts)
+            'analog_3'      - Raw analog signal 3 (if present, volts)
+            'analog_1_filt' - Filtered analog signal 1 (volts)
+            'analog_2_filt' - Filtered analog signal 2 (volts)
+            'analog_3_filt' - Filtered analog signal 2 (if present, volts)
+            'digital_1'     - Digital signal 1
+            'digital_2'     - Digital signal 2 (if present)
+            'pulse_inds_1'  - Locations of rising edges on digital input 1 (samples).
+            'pulse_inds_2'  - Locations of rising edges on digital input 2 (samples).
+            'pulse_times_1' - Times of rising edges on digital input 1 (ms).
+            'pulse_times_2' - Times of rising edges on digital input 2 (ms).
+            'time'          - Time of each sample relative to start of recording (ms)
+
+    with open(file_path, "rb") as f:
+        header_size = int.from_bytes(f.read(2), "little")
+        data_header = f.read(header_size)
+        data = np.frombuffer(f.read(), dtype=np.dtype("<u2"))
+    # Extract header information
+    header_dict = json.loads(data_header)
+    volts_per_division = header_dict["volts_per_division"]
+    sampling_rate = header_dict["sampling_rate"]
+    # Extract signals.
+    analog = data >> 1  # Analog signal is most significant 15 bits.
+    digital = ((data & 1) == 1).astype(int)  # Digital signal is least significant bit.
+    # Alternating samples are different signals.
+    if "n_analog_signals" in header_dict.keys():
+        n_analog_signals = header_dict["n_analog_signals"]
+        n_digital_signals = header_dict["n_digital_signals"]
+    else:  # Pre version 1.0 data file.
+        n_analog_signals = 2
+        n_digital_signals = 2
+    analog_1 = analog[::n_analog_signals] * volts_per_division[0]
+    analog_2 = analog[1::n_analog_signals] * volts_per_division[1]
+    analog_3 = analog[2::n_analog_signals] * volts_per_division[0] if n_analog_signals == 3 else None
+    digital_1 = digital[::n_analog_signals]
+    digital_2 = digital[1::n_analog_signals] if n_digital_signals == 2 else None
+    time = np.arange(analog_1.shape[0]) * 1000 / sampling_rate  # Time relative to start of recording (ms).
+    # Filter signals with specified high and low pass frequencies (Hz).
+    if low_pass and high_pass:
+        b, a = butter(2, np.array([high_pass, low_pass]) / (0.5 * sampling_rate), "bandpass")
+    elif low_pass:
+        b, a = butter(2, low_pass / (0.5 * sampling_rate), "low")
+    elif high_pass:
+        b, a = butter(2, high_pass / (0.5 * sampling_rate), "high")
+    if low_pass or high_pass:
+        analog_1_filt = filtfilt(b, a, analog_1)
+        analog_2_filt = filtfilt(b, a, analog_2)
+        analog_3_filt = filtfilt(b, a, analog_3) if n_analog_signals == 3 else None
+    else:
+        analog_1_filt = analog_2_filt = analog_3_filt = None
+    # Extract rising edges for digital inputs.
+    pulse_inds_1 = 1 + np.where(np.diff(digital_1) == 1)[0]
+    pulse_inds_2 = 1 + np.where(np.diff(digital_2) == 1)[0] if n_digital_signals == 2 else None
+    pulse_times_1 = pulse_inds_1 * 1000 / sampling_rate
+    pulse_times_2 = pulse_inds_2 * 1000 / sampling_rate if n_digital_signals == 2 else None
+    # Return signals + header information as a dictionary.
+    data_dict = {
+        "filename": os.path.basename(file_path),
+        "analog_1": analog_1,
+        "analog_2": analog_2,
+        "analog_1_filt": analog_1_filt,
+        "analog_2_filt": analog_2_filt,
+        "digital_1": digital_1,
+        "digital_2": digital_2,
+        "pulse_inds_1": pulse_inds_1,
+        "pulse_inds_2": pulse_inds_2,
+        "pulse_times_1": pulse_times_1,
+        "pulse_times_2": pulse_times_2,
+        "time": time,
+    }
+    if n_analog_signals == 3:
+        data_dict.update(
+            {
+                "analog_3": analog_3,
+                "analog_3_filt": analog_3_filt,
+            }
+        )
+    data_dict.update(header_dict)
+    return data_dict
+    """
+    ppd_data = import_ppd(ppd_file_path)  
+
+    raw_green =  pd.Series(ppd_data['analog_1'])
+    raw_red = pd.Series(ppd_data ['analog_2'])
+    raw_405 = pd.Series(ppd_data['analog_3'])
+    
+    relative_raw_signal = raw_green / raw_405   
+
+    sampling_rate = ppd_data['sampling_rate']
+
+    b,a = butter(2, 10, btype='low', fs=sampling_rate)
+    GACh_denoised = filtfilt(b,a, raw_green)
+    rDA3m_denoised = filtfilt(b,a, raw_red)
+    ratio_denoised = filtfilt(b,a, relative_raw_signal)
+    denoised_405 = filtfilt(b,a, raw_405)
+
+    # high pass at 0.001Hz which removes the drift due to bleaching, but will also remove any physiological variation in the signal on very slow timescales.
+    b,a = butter(2, 0.001, btype='high', fs=sampling_rate)
+    GACh_highpass = filtfilt(b,a, GACh_denoised, padtype='even')
+    rDA3m_highpass = filtfilt(b,a, rDA3m_denoised, padtype='even')
+    ratio_highpass = filtfilt(b,a, ratio_denoised, padtype='even')
+    highpass_405 = filtfilt(b,a, denoised_405, padtype='even')
+
+    green_zscored = np.divide(np.subtract(GACh_highpass,GACh_highpass.mean()),GACh_highpass.std())
+
+    red_zscored = np.divide(np.subtract(rDA3m_highpass,rDA3m_highpass.mean()),rDA3m_highpass.std())
+
+    zscored_405 = np.divide(np.subtract(highpass_405,highpass_405.mean()),highpass_405.std())
+
+    ratio_zscored = np.divide(np.subtract(ratio_highpass,ratio_highpass.mean()),ratio_highpass.std())
+
+    return {"green_zscored": green_zscored, "red_zscored": red_zscored, "ratio_zscored": ratio_zscored, "zscored_405": zscored_405}
+    pass
+
 
 def add_photometry_metadata(nwbfile: NWBFile, metadata: dict):
     # TODO for Ryan - add photometry metadata to NWB :)
@@ -395,9 +530,11 @@ def add_photometry(nwbfile: NWBFile, metadata: dict):
     elif "ppd_file_path" in metadata["photometry"]:
         # Process ppd file from pyPhotometry
         print("Processing ppd file from pyPhotometry...")
-        
+        ppd_file_path = metadata["photometry"]["ppd_file_path"]
+        signals = processppdphotometry(ppd_file_path)
         # TODO for Jose - add pyPhotometry processing here!! 
         # Probably add the processing functions above and just call them here
+
         raise NotImplementedError("pyPhotometry processing is not yet implemented.")
 
     else: