Merge pull request #56 from dj-sciops/update-lfp

fix: major refactor and update methods in LFP processing

Author: ttngu207

CHANGELOG.md

@@ -3,6 +3,11 @@
 Observes [Semantic Versioning](https://semver.org/spec/v2.0.0.html) standard and
 [Keep a Changelog](https://keepachangelog.com/en/1.0.0/) convention.
 
+
+## [0.6.0] - 2025-05-31
+
+- Fix - Major refactor and fix methods in `ephys.LFP`  
+
 ## [0.5.1] - 2025-04-23
 
 - Add - `execution_duration` to `LFP` table

element_array_ephys/ephys_no_curation.py

@@ -239,7 +239,7 @@ class Trace(dj.Part):
         -> master
         -> probe.ElectrodeConfig.Electrode
         ---
-        lfp              : blob@datajoint-blob
+        lfp              : blob@datajoint-blob # uV
         """
 
     @property
@@ -252,143 +252,154 @@ def key_source(self):
         )
 
     def make(self, key):
+        """Compute broadband LFP signals for each electrode.
+
+        Args:
+            key (dict): EphysSession primary key.
+
+        Raises:
+            ValueError: If the trace duration is not within the expected range.
+            OSError: If there is an error when loading the file.
+
+        Logic:
+            - Fetch the probe information for the given ephys session.
+            - Fetch the electrode configuration for the given probe.
+            - Fetch the raw data files for the given ephys session.
+            - Check for missing files or short trace durations in min
+            - Design notch filter to remove powerline noise that contaminates the LFP
+            - Downsample the signal with `decimate` and apply an anti-aliasing FIR filter
+        """
         execution_time = datetime.now(timezone.utc)
 
+        # Define constants
         TARGET_SAMPLING_RATE = 2500  # Hz
         POWERLINE_NOISE_FREQ = 60  # Hz
-        LFP_DURATION = 30  # minutes
+        MAX_DURATION_MINUTES = 30  # Minutes
 
+        # Check if the trace duration is within the expected range
         duration = (key["end_time"] - key["start_time"]).total_seconds() / 60  # minutes
         assert (
-            duration <= LFP_DURATION
-        ), f"LFP sessions cannot exceeds {LFP_DURATION} minutes in duration."
+            duration <= MAX_DURATION_MINUTES
+        ), f"LFP session duration {duration} min > max session duration {MAX_DURATION_MINUTES} min"
 
+        # Fetch the raw data files for the given ephys session
         query = (
             EphysRawFile
             & f"file_time BETWEEN '{key['start_time']}' AND '{key['end_time']}'"
         )
         if not query:
-            logger.info(
-                f"No raw data file found. Skip populating ephys.LFP for <{key}>"
-            )
-        else:
-            logger.info(f"Populating ephys.LFP for <{key}>")
+            logger.info(f"No raw data file found. Skipping LFP for <{key}>")
+            return
 
-            # Get probe info
-            probe_info = (EphysSessionProbe & key).fetch1()
-            probe_type = (probe.Probe & {"probe": probe_info["probe"]}).fetch1(
-                "probe_type"
-            )
+        logger.info(f"Populating ephys.LFP for <{key}>")
 
-            electrode_query = probe.ElectrodeConfig.Electrode & (
-                probe.ElectrodeConfig & {"probe_type": probe_type}
-            )
-
-            # Filter for used electrodes. If probe_info["used_electrodes"] is None, it means all electrodes were used.
-            if probe_info["used_electrodes"]:
-                electrode_query &= (
-                    f'electrode IN {tuple(probe_info["used_electrodes"])}'
-                )
+        # Fetch the probe information for the given ephys session
+        probe_info = (EphysSessionProbe & key).fetch1()
+        probe_type = (probe.Probe & {"probe": probe_info["probe"]}).fetch1("probe_type")
+        electrode_query = probe.ElectrodeConfig.Electrode & (
+            probe.ElectrodeConfig & {"probe_type": probe_type}
+        )
 
-            header = {}
-            lfp_concat = np.array([], dtype=np.float64)
+        # Fetch the electrode configuration for the given probe
+        # Filter for used electrodes. If probe_info["used_electrodes"] is None, it means all electrodes were used.
+        if probe_info["used_electrodes"]:
+            electrode_query &= f"electrode IN {tuple(probe_info['used_electrodes'])}"
 
-            for file_relpath in query.fetch("file_path", order_by="file_time"):
-                file = find_full_path(get_ephys_root_data_dir(), file_relpath)
+        header = {}
+        lfp_concat = []
 
-                try:
-                    data = intanrhdreader.load_file(file)
-                except OSError:
-                    raise OSError(f"OS error occurred when loading file {file.name}")
+        # Iterate over the raw data files for the given ephys session to load the data
+        for file_relpath in query.fetch("file_path", order_by="file_time"):
+            file = find_full_path(get_ephys_root_data_dir(), file_relpath)
+            try:
+                data = intanrhdreader.load_file(file)
+            except OSError:
+                raise OSError(f"OS error occurred when loading file {file.name}")
+
+            if not header:
+                header = data.pop("header")
+                lfp_sampling_rate = header["sample_rate"]
+                powerline_noise_freq = (
+                    header["notch_filter_frequency"] or POWERLINE_NOISE_FREQ
+                )  # in Hz
+
+                # Calculate downsampling factor
+                true_ratio = lfp_sampling_rate / TARGET_SAMPLING_RATE
+                downsample_factor = int(np.round(true_ratio))
+
+                # Check if the ratio is within 1% of an integer (1% tolerance)
+                if not np.isclose(true_ratio, downsample_factor, rtol=0.01, atol=1e-8):
+                    raise ValueError(
+                        f"Downsampling factor {true_ratio} is too far from an integer. Check LFP sampling rates."
+                    )
 
-                if not header:
-                    header = data.pop("header")
-                    lfp_sampling_rate = header["sample_rate"]
-                    powerline_noise_freq = (
-                        header["notch_filter_frequency"] or POWERLINE_NOISE_FREQ
-                    )  # in Hz
-                    downsample_factor = int(lfp_sampling_rate / TARGET_SAMPLING_RATE)
+                # Get LFP indices (row index of the LFP matrix to be used)
+                lfp_indices = np.array(electrode_query.fetch("channel_idx"), dtype=int)
+                port_indices = np.array(
+                    [
+                        ind
+                        for ind, ch in enumerate(data["amplifier_channels"])
+                        if ch["port_prefix"] == probe_info["port_id"]
+                    ]
+                )
+                lfp_indices = np.sort(port_indices[lfp_indices])
 
-                    # Get LFP indices (row index of the LFP matrix to be used)
-                    lfp_indices = np.array(
-                        electrode_query.fetch("channel_idx"), dtype=int
-                    )
-                    port_indices = np.array(
-                        [
-                            ind
-                            for ind, ch in enumerate(data["amplifier_channels"])
-                            if ch["port_prefix"] == probe_info["port_id"]
-                        ]
-                    )
-                    lfp_indices = np.sort(port_indices[lfp_indices])
+                self.insert1({**key, "lfp_sampling_rate": TARGET_SAMPLING_RATE})
 
-                    self.insert1(
-                        {
-                            **key,
-                            "lfp_sampling_rate": TARGET_SAMPLING_RATE,
-                        }
-                    )
+                # Get LFP channels
+                channels = np.array(
+                    [
+                        ch["native_channel_name"]
+                        for ch in data["amplifier_channels"]
+                        if ch["port_prefix"]
+                    ]
+                )[lfp_indices]
 
-                    channels = np.array(
-                        [
-                            ch["native_channel_name"]
-                            for ch in data["amplifier_channels"]
-                            if ch["port_prefix"]
-                        ]
-                    )[lfp_indices]
+                # Get channel to electrode mapping
+                electrode_df = electrode_query.fetch(format="frame").reset_index()
+                channel_to_electrode_map = dict(
+                    zip(electrode_df["channel_idx"], electrode_df["electrode"])
+                )
 
-                    electrode_df = electrode_query.fetch(format="frame").reset_index()
+                channel_to_electrode_map = {
+                    f'{probe_info["port_id"]}-{int(channel):03d}': electrode
+                    for channel, electrode in channel_to_electrode_map.items()
+                }
 
-                    channel_to_electrode_map = dict(
-                        zip(electrode_df["channel_idx"], electrode_df["electrode"])
-                    )
+            lfps = data.pop("amplifier_data")[lfp_indices]
+            lfp_concat.append(lfps)
 
-                    channel_to_electrode_map = {
-                        f'{probe_info["port_id"]}-{int(channel):03d}': electrode
-                        for channel, electrode in channel_to_electrode_map.items()
-                    }
+        full_lfp = np.hstack(lfp_concat)
 
-                lfps = data.pop("amplifier_data")[lfp_indices]
-                lfp_concat = (
-                    lfps if lfp_concat.size == 0 else np.hstack((lfp_concat, lfps))
-                )
-                del data
+        # Check if the trace duration is within the expected range
+        trace_duration = full_lfp.shape[1] / lfp_sampling_rate / 60  # in min
+        if abs(trace_duration - duration) > 0.5:
+            raise ValueError(
+                f"Trace duration mismatch: expected {duration}, got {trace_duration} min"
+            )
 
-            # Check for missing files or short trace durations in min
-            trace_duration = lfp_concat.shape[1] / lfp_sampling_rate / 60  # in min
-            if trace_duration != duration:
-                raise ValueError(
-                    f"Trace length ({trace_duration} min) is less than session duration"
-                )
+        # Design notch filter
+        notch_b, notch_a = signal.iirnotch(
+            w0=powerline_noise_freq, Q=30, fs=lfp_sampling_rate
+        )
 
-            # Single insert in loop to mitigate potential memory issue.
-            for ch, lfp in zip(channels, lfp_concat):
-                # Powerline noise removal
-                b_notch, a_notch = signal.iirnotch(
-                    w0=powerline_noise_freq, Q=30, fs=TARGET_SAMPLING_RATE
-                )
-                lfp = signal.filtfilt(b_notch, a_notch, lfp)
+        for ch_idx, raw_lfp in enumerate(full_lfp):
 
-                # Lowpass filter
-                b_butter, a_butter = signal.butter(
-                    N=4, Wn=1000, btype="lowpass", fs=TARGET_SAMPLING_RATE
-                )
-                lfp = signal.filtfilt(b_butter, a_butter, lfp)
+            # Apply notch filter
+            lfp = signal.filtfilt(notch_b, notch_a, raw_lfp)
 
-                # Downsample the signal
-                lfp = lfp[::downsample_factor]
+            # Downsample the signal with `decimate`
+            lfp = signal.decimate(lfp, downsample_factor, ftype="fir", zero_phase=True)
 
-                self.Trace.insert1(
-                    {
-                        **key,
-                        "electrode_config_hash": electrode_df["electrode_config_hash"][
-                            0
-                        ],
-                        "probe_type": electrode_df["probe_type"][0],
-                        "electrode": channel_to_electrode_map[ch],
-                        "lfp": lfp,
-                    }
-                )
+            self.Trace.insert1(
+                {
+                    **key,
+                    "electrode_config_hash": electrode_df["electrode_config_hash"][0],
+                    "probe_type": electrode_df["probe_type"][0],
+                    "electrode": channel_to_electrode_map[ch_idx],
+                    "lfp": lfp,
+                }
+            )
 
         self.update1(
             {

element_array_ephys/version.py

@@ -1,3 +1,3 @@
 """Package metadata."""
 
-__version__ = "0.5.1"
+__version__ = "0.6.0"