Repository: sappelhoff/pyprep
Pull request #72 (2021-04-29 01:48:32)

Skip RANSAC predictions for channels already confirmed as bad
<!--Thanks for contributing-->
<!--If this is your first time, please make sure to read the contributing guideline-->
<!--https://github.com/sappelhoff/pyprep/blob/master/.github/CONTRIBUTING.md-->

# PR Description

Started work on finishing up #66 and ended up noticing this in the process. Essentially, MATLAB PREP doesn't bother to do RANSAC predictions/correlations for channels that have already been flagged as unsuitable for use as RANSAC predictors (i.e., bad by correlation, deviation, or dropout). Instead, those rows in the correlation matrix get filled with ones, so those channels are never "bad-by-RANSAC" (since they're already confirmed bad by something else).

Anyway, this reduces the RAM requirements for channel-wise RANSAC and also ends up simplifying the existing RANSAC code quite a bit (since we no longer need to pass the lists of complete/good channels all the way down to `_get_ransac_pred` to transform good channel indices to complete channel indices).

# Merge Checklist

<!--To merge your PR we need to first take the following points into account.-->
<!--Please just leave this checklist untouched-->

- [ ] the PR has been reviewed and all comments are resolved
- [ ] all [CI][what-is-ci] checks pass
- [ ] (if applicable): the PR description includes the phrase `closes #<issue-number>` to [automatically close an issue][auto-close-documentation]
- [ ] (if applicable): bug fixes, new features, or [API][what-is-api] changes are documented in [whats_new.rst][whats-new-file]


[what-is-ci]: https://help.github.com/en/actions/building-and-testing-code-with-continuous-integration/about-continuous-integration
[auto-close-documentation]: https://help.github.com/en/github/managing-your-work-on-github/linking-a-pull-request-to-an-issue#linking-a-pull-request-to-an-issue-using-a-keyword
[what-is-api]: https://en.wikipedia.org/wiki/Application_programming_interface
[whats-new-file]: https://github.com/sappelhoff/pyprep/blob/master/docs/whats_new.rst


codecov-commenter: # [Codecov](https://codecov.io/gh/sappelhoff/pyprep/pull/72?src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=Stefan+Appelhoff) Report
> Merging [#72](https://codecov.io/gh/sappelhoff/pyprep/pull/72?src=pr&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=Stefan+Appelhoff) (1bf68ef) into [master](https://codecov.io/gh/sappelhoff/pyprep/commit/0fa6ac41acbe0e31228a54e3ec022505d2826137?el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=Stefan+Appelhoff) (0fa6ac4) will **decrease** coverage by `0.01%`.
> The diff coverage is `100.00%`.

[![Impacted file tree graph](https://codecov.io/gh/sappelhoff/pyprep/pull/72/graphs/tree.svg?width=650&height=150&src=pr&token=WniFfAHs8z&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=Stefan+Appelhoff)](https://codecov.io/gh/sappelhoff/pyprep/pull/72?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=Stefan+Appelhoff)

```diff
@@            Coverage Diff             @@
##           master      #72      +/-   ##
==========================================
- Coverage   97.15%   97.13%   -0.02%     
==========================================
  Files           7        7              
  Lines         597      594       -3     
==========================================
- Hits          580      577       -3     
  Misses         17       17              
```


| [Impacted Files](https://codecov.io/gh/sappelhoff/pyprep/pull/72?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=Stefan+Appelhoff) | Coverage Δ | |
|---|---|---|
| [pyprep/ransac.py](https://codecov.io/gh/sappelhoff/pyprep/pull/72/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=Stefan+Appelhoff#diff-cHlwcmVwL3JhbnNhYy5weQ==) | `97.27% <100.00%> (-0.08%)` | :arrow_down: |

------

[Continue to review full report at Codecov](https://codecov.io/gh/sappelhoff/pyprep/pull/72?src=pr&el=continue&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=Stefan+Appelhoff).
> **Legend** - [Click here to learn more](https://docs.codecov.io/docs/codecov-delta?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=Stefan+Appelhoff)
> `Δ = absolute <relative> (impact)`, `ø = not affected`, `? = missing data`
> Powered by [Codecov](https://codecov.io/gh/sappelhoff/pyprep/pull/72?src=pr&el=footer&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=Stefan+Appelhoff). Last update [0fa6ac4...1bf68ef](https://codecov.io/gh/sappelhoff/pyprep/pull/72?src=pr&el=lastupdated&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=Stefan+Appelhoff). Read the [comment docs](https://docs.codecov.io/docs/pull-request-comments?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=Stefan+Appelhoff).


Patch:
diff --git a/docs/matlab_differences.rst b/docs/matlab_differences.rst
index 17af766..85bf08a 100644
--- a/docs/matlab_differences.rst
+++ b/docs/matlab_differences.rst
@@ -18,6 +18,31 @@ internal math.
     :depth: 3
 
 
+Differences in Signal Detrending
+--------------------------------
+
+In the PREP pipeline, trends (i.e., slow drifts in EEG baseline signal) are
+temporarily removed from the data prior to adaptive line-noise removal
+as well as prior to bad channel detection via :class:`~pyprep.NoisyChannels`,
+which occurs at multiple points during robust re-referencing. This is done to
+improve the accuracy of both of these processes, which are sensitive to
+influence from trends in the signal.
+
+In MATLAB PREP, the default method of trend removal is to use EEGLAB's
+``pop_eegfiltnew``, which creates and applies an FIR high-pass filter to the
+data. MNE's :func:`mne.filter.filter_data` offers similar functionality, but
+uses slightly different filter creation math and a different filtering
+algorithm such that its results and subsequent :class:`~pyprep.NoisyChannels`
+values also differ slightly (on the order of ~0.002) for RANSAC correlations.
+
+Because the practical differences are small and MNE's filtering is fast and
+well-tested, PyPREP defaults to using :func:`mne.filter.filter_data` for
+high-pass trend removal. However, for exact numerical compatibility, PyPREP
+has a basic re-implementaion of EEGLAB's ``pop_eegfiltnew`` in Python that
+produces identical results to MATLAB PREP's ``removeTrend`` when
+``matlab_strict`` is set to ``True``.
+
+
 Differences in RANSAC
 ---------------------
 
diff --git a/docs/whats_new.rst b/docs/whats_new.rst
index b6d80c6..f28ca96 100644
--- a/docs/whats_new.rst
+++ b/docs/whats_new.rst
@@ -40,6 +40,7 @@ Changelog
 - Changed RANSAC's random channel sampling code to produce the same results as MATLAB PREP for the same random seed, additionally changing the default RANSAC sample size from 25% of all *good* channels (e.g. 15 for a 64-channel dataset with 4 bad channels) to 25% of *all* channels (e.g. 16 for the same dataset), by `Austin Hurst`_ (:gh:`62`)
 - Changed RANSAC so that "bad by high-frequency noise" channels are retained when making channel predictions (provided they aren't flagged as bad by any other metric), matching MATLAB PREP behaviour, by `Austin Hurst`_ (:gh:`64`)
 - Added a new flag ``matlab_strict`` to :class:`~pyprep.PrepPipeline`, :class:`~pyprep.Reference`, :class:`~pyprep.NoisyChannels`, and :func:`~pyprep.ransac.find_bad_by_ransac` for optionally matching MATLAB PREP's internal math as closely as possible, overriding areas where PyPREP attempts to improve on the original, by `Austin Hurst`_ (:gh:`70`)
+- Added a ``matlab_strict`` method for high-pass trend removal, exactly matching MATLAB PREP's values if ``matlab_strict`` is enabled, by `Austin Hurst`_ (:gh:`71`)
 
 Bug
 ~~~
@@ -47,6 +48,7 @@ Bug
 - Fixed "bad channel by flat" threshold in :meth:`NoisyChannels.find_bad_by_nan_flat` to be consistent with MATLAB PREP, by `Austin Hurst`_ (:gh:`60`)
 - Changed "bad channel by deviation" and "bad channel by correlation" detection code in :class:`NoisyChannels` to compute IQR and quantiles in the same manner as MATLAB, thus producing identical results to MATLAB PREP, by `Austin Hurst`_ (:gh:`57`)
 - Fixed a bug where EEG data was getting reshaped into RANSAC windows incorrectly (channel samples were not sequential), which was causing considerable variability and noise in RANSAC results, by `Austin Hurst`_ (:gh:`67`)
+- Fixed RANSAC to avoid making unnecessary signal predictions for known-bad channels, matching MATLAB behaviour and reducing RAM requirements, by `Austin Hurst`_ (:gh:`72`)
 
 API
 ~~~
diff --git a/pyprep/find_noisy_channels.py b/pyprep/find_noisy_channels.py
index d99aa38..03d02ef 100644
--- a/pyprep/find_noisy_channels.py
+++ b/pyprep/find_noisy_channels.py
@@ -52,7 +52,7 @@ class NoisyChannels:
         self.sample_rate = raw.info["sfreq"]
         if do_detrend:
             self.raw_mne._data = removeTrend(
-                self.raw_mne.get_data(), sample_rate=self.sample_rate
+                self.raw_mne.get_data(), self.sample_rate, matlab_strict=matlab_strict
             )
         self.matlab_strict = matlab_strict
 
diff --git a/pyprep/prep_pipeline.py b/pyprep/prep_pipeline.py
index ed0661b..2c4a38b 100644
--- a/pyprep/prep_pipeline.py
+++ b/pyprep/prep_pipeline.py
@@ -157,7 +157,9 @@ class PrepPipeline:
         # reference_channels = _set_diff(self.prep_params["ref_chs"], unusable_channels)
         # Step 1: 1Hz high pass filtering
         if len(self.prep_params["line_freqs"]) != 0:
-            self.EEG_new = removeTrend(self.EEG_raw, sample_rate=self.sfreq)
+            self.EEG_new = removeTrend(
+                self.EEG_raw, self.sfreq, matlab_strict=self.matlab_strict
+            )
 
             # Step 2: Removing line noise
             linenoise = self.prep_params["line_freqs"]
diff --git a/pyprep/ransac.py b/pyprep/ransac.py
index b6e0af5..5bc284f 100644
--- a/pyprep/ransac.py
+++ b/pyprep/ransac.py
@@ -111,7 +111,6 @@ def find_bad_by_ransac(
     # Exclude should be the bad channels from other methods
     # That is, identify all bad channels by other means
     good_idx = mne.pick_channels(list(complete_chn_labs), include=[], exclude=exclude)
-    good_chn_labs = complete_chn_labs[good_idx]
     n_chans_good = good_idx.shape[0]
     chn_pos_good = chn_pos[good_idx, :]
 
@@ -165,15 +164,15 @@ def find_bad_by_ransac(
     # Calculate smallest chunk size for each possible chunk count
     chunk_sizes = []
     chunk_count = 0
-    for i in range(1, n_chans_complete + 1):
-        n_chunks = int(np.ceil(n_chans_complete / i))
+    for i in range(1, n_chans_good + 1):
+        n_chunks = int(np.ceil(n_chans_good / i))
         if n_chunks != chunk_count:
             chunk_count = n_chunks
             chunk_sizes.append(i)
 
     chunk_size = chunk_sizes.pop()
     mem_error = True
-    job = list(range(n_chans_complete))
+    job = list(range(n_chans_good))
 
     if channel_wise:
         chunk_size = 1
@@ -183,14 +182,11 @@ def find_bad_by_ransac(
             total_chunks = len(channel_chunks)
             current = 1
             for chunk in channel_chunks:
-                channel_correlations[:, chunk] = _ransac_correlations(
+                channel_correlations[:, good_idx[chunk]] = _ransac_correlations(
                     chunk,
                     random_ch_picks,
-                    chn_pos,
                     chn_pos_good,
-                    good_chn_labs,
-                    complete_chn_labs,
-                    data,
+                    data[good_idx, :],
                     n_samples,
                     n,
                     w_correlation,
@@ -217,7 +213,7 @@ def find_bad_by_ransac(
                     "ested samples."
                 )
 
-    # Thresholding
+    # Calculate fractions of bad RANSAC windows for each channel
     thresholded_correlations = channel_correlations < corr_thresh
     frac_bad_corr_windows = np.mean(thresholded_correlations, axis=0)
 
@@ -233,10 +229,7 @@ def find_bad_by_ransac(
 def _ransac_correlations(
     chans_to_predict,
     random_ch_picks,
-    chn_pos,
     chn_pos_good,
-    good_chn_labs,
-    complete_chn_labs,
     data,
     n_samples,
     n,
@@ -251,15 +244,9 @@ def _ransac_correlations(
         Indexes of the channels to predict as they appear in chn_pos.
     random_ch_picks : list
         each element is a list of indexes of the channels (as they appear
-        in good_chn_labs) to use for reconstruction in each of the samples.
-    chn_pos : np.ndarray
-        3-D coordinates of the electrode positions to predict
+        in chn_pos_good) to use for reconstruction in each of the samples.
     chn_pos_good : np.ndarray
         3-D coordinates of all the channels not detected noisy so far
-    good_chn_labs : array_like
-        channel labels for the ch_pos_good channels
-    complete_chn_labs : array_like
-        labels of the channels in data in the same order
     data : np.ndarray
         2-D EEG data
     n_samples : int
@@ -285,10 +272,8 @@ def _ransac_correlations(
     ransac_eeg = _run_ransac(
         n_samples=n_samples,
         random_ch_picks=random_ch_picks,
-        chn_pos=chn_pos[chans_to_predict, :],
+        chn_pos=chn_pos_good[chans_to_predict, :],
         chn_pos_good=chn_pos_good,
-        good_chn_labs=good_chn_labs,
-        complete_chn_labs=complete_chn_labs,
         data=data,
         matlab_strict=matlab_strict,
     )
@@ -320,8 +305,6 @@ def _run_ransac(
     random_ch_picks,
     chn_pos,
     chn_pos_good,
-    good_chn_labs,
-    complete_chn_labs,
     data,
     matlab_strict,
 ):
@@ -336,15 +319,11 @@ def _run_ransac(
         number of interpolations from which a median will be computed
     random_ch_picks : list
         each element is a list of indexes of the channels (as they appear
-        in good_chn_labs) to use for reconstruction in each of the samples.
+        in chn_pos_good) to use for reconstruction in each of the samples.
     chn_pos : np.ndarray
         3-D coordinates of the electrode position
     chn_pos_good : np.ndarray
         3-D coordinates of all the channels not detected noisy so far
-    good_chn_labs : array_like
-        channel labels for the ch_pos_good channels
-    complete_chn_labs : array_like
-        labels of the channels in data in the same order
     data : np.ndarray
         2-D EEG data
     matlab_strict : bool
@@ -372,7 +351,7 @@ def _run_ransac(
         # Get the random channel selection for the current sample
         reconstr_idx = random_ch_picks[sample]
         eeg_predictions[..., sample] = _get_ransac_pred(
-            chn_pos, chn_pos_good, good_chn_labs, complete_chn_labs, reconstr_idx, data
+            chn_pos, chn_pos_good, reconstr_idx, data
         )
 
     # Form median from all predictions
@@ -383,12 +362,11 @@ def _run_ransac(
         ransac_eeg = eeg_predictions[:, :, median_idx]
     else:
         ransac_eeg = np.median(eeg_predictions, axis=-1, overwrite_input=True)
+
     return ransac_eeg
 
 
-def _get_ransac_pred(
-    chn_pos, chn_pos_good, good_chn_labs, complete_chn_labs, reconstr_idx, data
-):
+def _get_ransac_pred(chn_pos, chn_pos_good, reconstr_idx, data):
     """Perform RANSAC prediction.
 
     Parameters
@@ -397,12 +375,8 @@ def _get_ransac_pred(
         3-D coordinates of the electrode position
     chn_pos_good : np.ndarray
         3-D coordinates of all the channels not detected noisy so far
-    good_chn_labs : array_like
-        channel labels for the ch_pos_good channels
-    complete_chn_labs : array_like
-        labels of the channels in data in the same order
     reconstr_idx : array_like
-        indexes of the channels in good_chn_labs to use for reconstruction
+        indexes of the channels in chn_pos_good to use for reconstruction
     data : np.ndarray
         2-D EEG data
 
@@ -412,17 +386,11 @@ def _get_ransac_pred(
         Single RANSAC prediction
 
     """
-    # Get positions and according labels
-    reconstr_labels = good_chn_labs[reconstr_idx]
+    # Get positions
     reconstr_pos = chn_pos_good[reconstr_idx, :]
 
-    # Map the labels to their indices within the complete data
-    # Do not use mne.pick_channels, because it will return a sorted list.
-    reconstr_picks = [
-        list(complete_chn_labs).index(chn_lab) for chn_lab in reconstr_labels
-    ]
-
     # Interpolate
     interpol_mat = _make_interpolation_matrix(reconstr_pos, chn_pos)
-    ransac_pred = np.matmul(interpol_mat, data[reconstr_picks, :])
+    ransac_pred = np.matmul(interpol_mat, data[reconstr_idx, :])
+
     return ransac_pred
diff --git a/pyprep/reference.py b/pyprep/reference.py
index 92df288..0cdddf3 100644
--- a/pyprep/reference.py
+++ b/pyprep/reference.py
@@ -175,7 +175,9 @@ class Reference:
 
         """
         raw = self.raw.copy()
-        raw._data = removeTrend(raw.get_data(), sample_rate=self.sfreq)
+        raw._data = removeTrend(
+            raw.get_data(), self.sfreq, matlab_strict=self.matlab_strict
+        )
 
         # Determine unusable channels and remove them from the reference channels
         noisy_detector = NoisyChannels(
diff --git a/pyprep/removeTrend.py b/pyprep/removeTrend.py
index ced97e8..da3dba1 100644
--- a/pyprep/removeTrend.py
+++ b/pyprep/removeTrend.py
@@ -4,51 +4,76 @@ import logging
 import mne
 import numpy as np
 
+from pyprep.utils import _eeglab_create_highpass, _eeglab_fir_filter
+
 
 def removeTrend(
     EEG,
-    detrendType="High pass",
-    sample_rate=160.0,
+    sample_rate,
+    detrendType="high pass",
     detrendCutoff=1.0,
     detrendChannels=None,
+    matlab_strict=False,
 ):
-    """Perform high pass filtering or detrending.
+    """Remove trends (i.e., slow drifts in baseline) from an array of EEG data.
 
     Parameters
     ----------
     EEG : np.ndarray
-        The input EEG data.
-    detrendType : str
-        Type of detrending to be performed: high pass, high pass sinc, or local
-        detrending.
+        A 2-D array of EEG data to detrend.
     sample_rate : float
-        Rate at which the EEG data was sampled.
-    detrendCutoff : float
-        High pass cut-off frequency.
+        The sample rate (in Hz) of the input EEG data.
+    detrendType : str, optional
+        Type of detrending to be performed: must be one of 'high pass',
+        'high pass sinc, or 'local detrend'. Defaults to 'high pass'.
+    detrendCutoff : float, optional
+        The high-pass cutoff frequency (in Hz) to use for detrending. Defaults
+        to 1.0 Hz.
     detrendChannels : {list, None}, optional
-        List of all the channels that require detrending/filtering. If None,
-        all channels are used (default).
+        List of the indices of all channels that require detrending/filtering.
+        If ``None``, all channels are used (default).
+    matlab_strict : bool, optional
+        Whether or not detrending should strictly follow MATLAB PREP's internal
+        math, ignoring any improvements made in PyPREP over the original code
+        (see :ref:`matlab-diffs` for more details). Defaults to ``False``.
 
     Returns
     -------
     EEG : np.ndarray
-        Filtered/detrended EEG data.
+        A 2-D array containing the filtered/detrended EEG data.
 
     Notes
     -----
-    Filtering is implemented using the MNE filter function mne.filter.filter_data.
-    Local detrending is the python implementation of the chronux_2 runline command.
+    High-pass filtering is implemented using the MNE filter function
+    :func:``mne.filter.filter_data`` unless `matlab_strict` is ``True``, in
+    which case it is performed using a minimal re-implementation of EEGLAB's
+    ``pop_eegfiltnew``. Local detrending is performed using a Python
+    re-implementation of the ``runline`` function from the Chronux package for
+    MATLAB [1]_.
+
+    References
+    ----------
+    .. [1] http://chronux.org/
 
     """
     if len(EEG.shape) == 1:
         EEG = np.reshape(EEG, (1, EEG.shape[0]))
 
-    if detrendType == "High pass":
-        EEG = mne.filter.filter_data(
-            EEG, sfreq=sample_rate, l_freq=1, h_freq=None, picks=detrendChannels
-        )
+    if detrendType.lower() == "high pass":
+        if matlab_strict:
+            picks = detrendChannels if detrendChannels else range(EEG.shape[0])
+            filt = _eeglab_create_highpass(detrendCutoff, sample_rate)
+            EEG[picks, :] = _eeglab_fir_filter(EEG[picks, :], filt)
+        else:
+            EEG = mne.filter.filter_data(
+                EEG,
+                sfreq=sample_rate,
+                l_freq=detrendCutoff,
+                h_freq=None,
+                picks=detrendChannels
+            )
 
-    elif detrendType == "High pass sinc":
+    elif detrendType.lower() == "high pass sinc":
         fOrder = np.round(14080 * sample_rate / 512)
         fOrder = np.int(fOrder + fOrder % 2)
         EEG = mne.filter.filter_data(
@@ -60,7 +85,8 @@ def removeTrend(
             filter_length=fOrder,
             fir_window="blackman",
         )
-    elif detrendType == "Local detrend":
+
+    elif detrendType.lower() == "local detrend":
         if detrendChannels is None:
             detrendChannels = np.arange(0, EEG.shape[0])
         windowSize = 1.5 / detrendCutoff
@@ -82,29 +108,38 @@ def removeTrend(
             for ch in detrendChannels:
                 EEG[:, ch] = runline(EEG[:, ch], np.int(n), np.int(dn))
         EEG = np.transpose(EEG)
+
     else:
         logging.warning(
             "No filtering/detreding performed since the detrend type did not match"
         )
+
     return EEG
 
 
 def runline(y, n, dn):
-    """Implement chronux_2 runline command for performing local linear regression.
+    """Perform local linear regression on a channel of EEG data.
+
+    A re-implementation of the ``runline`` function from the Chronux package
+    for MATLAB [1]_.
 
     Parameters
     ----------
     y : np.ndarray
-        Input from one EEG channel.
+        A 1-D array of data from a single EEG channel.
     n : int
-        length of the detrending window.
+        Length of the detrending window.
     dn : int
-        length of the window step size.
+        Length of the window step size.
 
     Returns
     -------
     y: np.ndarray
-       Detrended EEG signal for one channel.
+       The detrended signal for the given EEG channel.
+
+    References
+    ----------
+    .. [1] http://chronux.org/
 
     """
     nt = y.shape[0]
diff --git a/pyprep/utils.py b/pyprep/utils.py
index 987d9ed..51c9518 100644
--- a/pyprep/utils.py
+++ b/pyprep/utils.py
@@ -6,6 +6,7 @@ import mne
 import numpy as np
 import scipy.interpolate
 from scipy.stats import iqr
+from scipy.signal import firwin, lfilter, lfilter_zi
 from psutil import virtual_memory
 
 
@@ -108,6 +109,104 @@ def _mat_iqr(arr, axis=None):
     return iqr(arr, rng=np.clip(iqr_adj, 0, 100), axis=axis)
 
 
+def _eeglab_create_highpass(cutoff, srate):
+    """Create a high-pass FIR filter using Hamming windows.
+
+    Parameters
+    ----------
+    cutoff : float
+        The lower pass-band edge of the filter, in Hz.
+    srate : float
+        The sampling rate of the EEG signal, in Hz.
+
+    Returns
+    -------
+    filter : np.ndarray
+        A 1-dimensional array of FIR filter coefficients.
+
+    Notes
+    -----
+    In MATLAB PREP, the internal ``removeTrend`` function uses EEGLAB's
+    ``pop_eegfiltnew`` to high-pass the EEG data to remove slow drifts.
+    Because MNE's ``mne.filter.filter_data`` and EEGLAB's ``pop_eegfiltnew``
+    calculate filter parameters slightly differently, this function is
+    used to precisely match EEGLAB & MATLAB PREP's filtering method.
+
+    """
+    TRANSITION_WIDTH_RATIO = 0.25
+    HAMMING_CONSTANT = 3.3  # note: not entirely clear what this represents
+
+    # Calculate parameters for constructing filter
+    trans_bandwidth = cutoff if cutoff < 2 else cutoff * TRANSITION_WIDTH_RATIO
+    order = HAMMING_CONSTANT / (trans_bandwidth / srate)
+    order = int(np.ceil(order / 2) * 2)  # ensure order is even
+    stop = cutoff - trans_bandwidth
+    transition = (stop + cutoff) / srate
+
+    # Generate highpass filter
+    N = order + 1
+    filt = np.zeros(N)
+    filt[N // 2] = 1
+    filt -= firwin(N, transition, window='hamming', nyq=1)
+    return filt
+
+
+def _eeglab_fir_filter(data, filt):
+    """Apply an FIR filter to a 2-D array of EEG data.
+
+    Parameters
+    ----------
+    data : np.ndarray
+        A 2-D array of EEG data to filter.
+    filt : np.ndarray
+        A 1-D array of FIR filter coefficients.
+
+    Returns
+    -------
+    filtered : np.ndarray
+        A 2-D array of FIR-filtered EEG data.
+
+    Notes
+    -----
+    Produces identical output to EEGLAB's ``firfilt`` function (for non-epoched
+    data). For internal use within :mod:`pyprep.removeTrend`.
+
+    """
+    # Initialize parameters for FIR filtering
+    frames_per_window = 2000
+    group_delay = int((len(filt) - 1) / 2)
+    n_samples = data.shape[1]
+    n_windows = int(np.ceil((n_samples - group_delay) / frames_per_window))
+    pad_len = min(group_delay, n_samples)
+
+    # Prepare initial state of filter, using padding at start of data
+    start_pad_idx = np.zeros(pad_len, dtype=np.uint8)
+    start_padded = np.concatenate(
+        (data[:, start_pad_idx], data[:, :pad_len]),
+        axis=1
+    )
+    zi_init = lfilter_zi(filt, 1) * np.take(start_padded, [0], axis=0)
+    _, zi = lfilter(filt, 1, start_padded, axis=1, zi=zi_init)
+
+    # Iterate over windows of signal, filtering in chunks
+    out = np.zeros_like(data)
+    for w in range(n_windows):
+        start = group_delay + w * frames_per_window
+        end = min(start + frames_per_window, n_samples)
+        start_out = start - group_delay
+        end_out = end - group_delay
+        out[:, start_out:end_out], zi = lfilter(
+            filt, 1, data[:, start:end], axis=1, zi=zi
+        )
+
+    # Finish filtering data, using padding at end to calculate final values
+    end_pad_idx = np.zeros(pad_len, dtype=np.uint8) + (n_samples - 1)
+    end, _ = lfilter(filt, 1, data[:, end_pad_idx], axis=1, zi=zi)
+    out[:, (n_samples - pad_len):] = end[:, (group_delay - pad_len):]
+
+    return out
+
+
 def _get_random_subset(x, size, rand_state):
     """Get a random subset of items from a list or array, without replacement.
 


Test patch:
diff --git a/tests/test_removeTrend.py b/tests/test_removeTrend.py
index 2bc88c0..9d06fef 100644
--- a/tests/test_removeTrend.py
+++ b/tests/test_removeTrend.py
@@ -17,10 +17,19 @@ def test_highpass():
     lowpass_filt2 = removeTrend.removeTrend(
         signal, detrendType="High pass", sample_rate=srate, detrendCutoff=1
     )
+    lowpass_filt3 = removeTrend.removeTrend(
+        signal,
+        detrendType="High pass",
+        sample_rate=srate,
+        detrendCutoff=1,
+        matlab_strict=True
+    )
     error1 = lowpass_filt1 - highfreq_signal
     error2 = lowpass_filt2 - highfreq_signal
+    error3 = lowpass_filt3 - highfreq_signal
     assert np.sqrt(np.mean(error1 ** 2)) < 0.1
     assert np.sqrt(np.mean(error2 ** 2)) < 0.1
+    assert np.sqrt(np.mean(error3 ** 2)) < 0.1
 
 
 def test_detrend():
diff --git a/tests/test_utils.py b/tests/test_utils.py
index a7711b8..8cf3547 100644
--- a/tests/test_utils.py
+++ b/tests/test_utils.py
@@ -2,7 +2,8 @@
 import numpy as np
 
 from pyprep.utils import (
-    _mat_round, _mat_quantile, _mat_iqr, _get_random_subset, _correlate_arrays
+    _mat_round, _mat_quantile, _mat_iqr, _get_random_subset, _correlate_arrays,
+    _eeglab_create_highpass
 )
 
 
@@ -92,3 +93,24 @@ def test_correlate_arrays():
     corr_expected = np.asarray([-0.0898, 0.0327, -0.1140])
     corr_actual = _correlate_arrays(a, b, matlab_strict=True)
     assert all(np.isclose(corr_expected, corr_actual, atol=0.001))
+
+
+def test_eeglab_create_highpass():
+    """Test EEGLAB-equivalent high-pass filter creation.
+
+    NOTE: EEGLAB values were obtained using breakpoints in ``pop_eegfiltnew``,
+    since filter creation and data filtering are both done in the same function.
+    Values here are first 4 values of the array ``b`` which contains the FIR
+    filter coefficents used by the function.
+
+    """
+    # Compare initial FIR filter coefficents with EEGLAB
+    expected_vals = [5.3691e-5, 5.4165e-5, 5.4651e-5, 5.5149e-5]
+    actual_vals = _eeglab_create_highpass(cutoff=1.0, srate=256)[:4]
+    assert all(np.isclose(expected_vals, actual_vals, atol=0.001))
+
+    # Compare middle FIR filter coefficent with EEGLAB
+    vals = _eeglab_create_highpass(cutoff=1.0, srate=256)
+    expected_val = 0.9961
+    actual_val = vals[len(vals) // 2]
+    assert np.isclose(expected_val, actual_val, atol=0.001)