Repository: sappelhoff/pyprep (language: python)

Issue:
Testing Matlab PREP in our CI? Directly comparing with pyprep
This is not something to do right now, but as far as I heard it's now possible to test actual MATLAB in some CI services, see for reference: https://github.com/bids-standard/bids-matlab/issues/68#issuecomment-717187963

That could mean we can in our CI

1. download some lightweight data
1. run PREP on it (matlab), save outputs as build artifacts
1. run pyprep on it (python), save outputs as build artifacts
1. run another python script that compares the two build artifacts
Numeric equivalence with MATLAB PREP
This is a general tracking issue for the numeric equivalence of PyPREP vs Matlab PREP at each stage of the pipeline. I figure this will be useful for keeping a bird's eye view of what's been tested and what hasn't, as well as what parts have been confirmed to differ.

Guide: items with a checkmark have been tested and confirmed as producing identical results to MATLAB. Items with no checkmark and a corresponding issue # have been tested and confirmed as diverging. Items with no checkmark and no issue number have yet to be formally tested.

### NoisyChannels Methods

 - [x] Bad by NaN
 - [x] Bad by flat (fixed by #60)
 - [x] Bad by deviation (fixed by #57)
 - [x] Bad by HF noise
 - [x] Bad by correlation (fixed by #57)
 - [x] Bad by low SNR (fixed by #57)
 - [x] Bad by dropout (fixed by #81)
 - [x] Bad by RANSAC
      - [x] 'Bad' channel exclusion (fixed by #64)
      - [x] Channel picks (with fixed seed) (fixed by #62)
      - [x] Spherical interpolation algorithm (with identical channel locations, produces identical matrix 🎉)
      - [x] Channel interpolation (fixed by #70)
      - [x] Channel correlation (fixed by #67 and #70)
      - [x] Channel thresholding

### Other Functions

- [x] removeTrend (fixed by #71)
- [ ] CleanLine (`mne.filter.notch_filter`) (#82: PyPREP's line noise removal is different, seems to be more aggressive)
- [x] Robust re-referencing
     - [x] Initial average referencing (fixed by #78)
     - [x] Reference removal
     - [x] Bad channel identification (fixed by #78, #89, and #93)
     - [x] Bad channel interpolation (fixed by #96)
     - [x] Final re-referencing (fixed by #92)
- [x] Final bad channel interpolation (fixed by #92 and #96)





Discussion / hints:
Direct comparison during CI like this is a good idea, but I think we can probably do it without actually adding any MATLAB code to the CI: as long as the MATLAB PREP results are consistent with a random seed, we could run MATLAB PREP on a file locally, upload the resulting .set to Google Drive or a similar service that allows downloading files with curl/wget, and then have the CI just download it prior to running PyPREP tests.

Also, I remember @yjmantilla mentioning that PyPREP intentionally differs in some minor ways from MATLAB PREP because they do some weird things with medians or rounding or something. To make tests like this work, we'd need to add a special flag of some kind that would make PyPREP to do the same weird things for comparison purposes so that we can be sure the only functional differences are in those areas.
> but I think we can probably do it without actually adding any MATLAB code to the CI: as long as the MATLAB PREP results are consistent with a random seed, we could run MATLAB PREP on a file locally, upload the resulting .set to Google Drive or a similar service that allows downloading files with curl/wget, and then have the CI just download it prior to running PyPREP tests.

agreed, that'd be a good first step

Having it properly implemented in CI however serves the additional purposes of:
- having the MATLAB pipeline documented transparently and easily accessible
- having the opportunity to be notified about upstream changes in MATLAB prep

> because they do some weird things with medians or rounding or something.

yes, I think they didn't take the proper median for an even number of channels ... not sure if that has been fixed by now. We should check it out if we find the time. I did a quick search of autoreject issues and prep (matlab) issues and PRs, but couldn't find anything. @jasmainak originally detected that problem, maybe he knows where to find it again.

>  To make tests like this work, we'd need to add a special flag of some kind that would make PyPREP to do the same weird things for comparison purposes so that we can be sure the only functional differences are in those areas.

yes, that'd need some clever fixes ... such as flags. or preferably, selecting a test dataset where that specific median-difference doesn't matter between PREP and PYPREP.


>  I did a quick search of autoreject issues and prep (matlab) issues and PRs, but couldn't find anything

there are none. Autoreject used to be a private repo at the time the RANSAC in PREP was implemented and I didn't want to go into a rabbit hole fixing PREP. Had some ideas to improve RANSAC in PREP but never got around to them :)
@sappelhoff 

> yes, that'd need some clever fixes ... such as flags. or preferably, selecting a test dataset where that specific median-difference doesn't matter between PREP and PYPREP.

In this one I would go for the flag, that is , if someone wants to go with the matlab's prep way they can do it through it. I guess we should really open the issue there to see whats going on with that. Maybe they will correct it and there won't be any need for the flag.

Also I don't know how easy it would be to find such a dataset... maybe is super easy but the flag would suffice in case it is hard to find.

With #45 the CI comparison should become easier.
 

yes, would be fine with me. I guess the first job would be to go back to MATLAB PREP and find the ominous false-median-computation.

... and then all other points in code where pyprep and prep deviate. THEN we can think about how best to allow users to change pyprep from "best practice in our opinion" to original PREP.


jasmainak reported [these lines](https://github.com/VisLab/EEG-Clean-Tools/blob/48c2e670cf47ba998edae45a9b49c74595bddb3a/PrepPipeline/utilities/findNoisyChannels.m#L352-L353) in #21 


@sappelhoff I've never used it, but there's a promising-looking package called [oct2py](https://pypi.org/project/oct2py/) that lets you run MATLAB functions via Octave through Python, passing and receiving data as numpy arrays.

If we can get it working with the basics of PREP, it could be handy for both CI and for our validation efforts since it would  let us substitute individual PyPREP functions with their MATLAB PREP equivalents to compare things like interpolation directly between packages.
@sappelhoff @yjmantilla A quick update I thought you'd be interested in: during my recent burst of PyPREP debugging I spent a few hours hacking away at things and got EEGLAB + PREP (somewhat) running via oct2py:

![Screen Shot 2021-03-22 at 12 22 25 AM](https://user-images.githubusercontent.com/18648066/111937032-b2dd0c00-8aa5-11eb-8688-00112681d044.png)
Most of the work was figuring out how to convert MNE Raw objects into the EEGLAB format, which now seems to be working well enough for PREP to not error out immediately. A quick rundown of some problems to eventually overcome:

- PREP is not natively Octave-compatible: it relies on a couple odd MATLAB-only functions that I needed to substitute or comment out. Fine for testing locally, but would require workarounds for CI usage.
- The findNoisyChannels bit that does the channel correlations before RANSAC is *slooooooooooooowwwwwwww*, like 4-5 hours per iteration slow. Guessing this is an Octave issue as well? At its current pace, it's going to take a full day to do robust referencing on a file that takes maybe 6-10 minutes for the same via PyPREP.

In other words, it seems to work, but there seem to be some changes needed to the MATLAB code itself before CI comparison is doable.
nice, thanks! :sunglasses:

Now did you by any chance do your changes by
1. forking PREP
2. making a new branch on your fork and adjusting for octave

?

Because then you could push that branch up to GitHub and open a PR to your own fork --> that'll let us inspect the git diff, so we can better see the extent of the adjustments :-)

it's surprising that it's so ultra slow now. Perhaps profiling the function would be a good idea, to see exactly where this slowing comes from. :thinking: 
@sappelhoff I'll get to doing a proper fork eventually: right now I'm still ironing out what's Octave-compatible and what isn't (which is made a lot more difficult by every change/fix taking half a day chugging along to test). This morning I woke up to this error after two successful iterations of findNoisyChannels so I guess that's the next thing to figure out:

```
oct2py.utils.Oct2PyError: Octave evaluation error:
error: structure has no member 'etc'
error: called from:
    pop_select at line 196, column 1
    eeg_interp at line 158, column 19
    interpolateChannels at line 49, column 12
    robustReference at line 80, column 19
    performReference>doRobustPost at line 118, column 22
    performReference at line 76, column 5
```

Because I'm lazy and don't know MATLAB too well, I'm going to try "print-by-line" profiling here to get a better sense of what's going on with the correlation speed. That said, if you know of a better/easier way of profiling MATLAB code that's Octave-compatible, please let me know!
I see - and I think there's a very high chance to get lost by **not** following the workflow of:

fork, then clone, then checkout a new branch, and then adding your changes/adjustments incrementally with new commits that you then push go GitHub (on your fork), using meaningful commit messages

... so just a recommendation to start with that better early than late :)

apart from that I don't have good advice unfortunately :-(
So far the changes are minimal, but of course you're right: I'll start tracking things more formally as I go.

Also, found the cause of the slowness: the calculation of the correlation coefficients for each 1-second window of the file takes 2-to-3 seconds per window for some reason? 

```
windowCorrelation = corrcoef(eegPortion);
```

Since the file is about 2700 seconds long, that adds up to several hours. No clue why it's so slow here (I'm assuming it's faster in regular MATLAB, since it's massively faster in PyPREP), but I'll look into it!
@sappelhoff Just a quick status update: I've given up on oct2py for the time being given how horribly slow Octave seems to be at findNoisyChannels, but I realized I've got a site licence for MATLAB from my institution so I installed it and got MatPREP working smoothly using @yjmantilla's script they shared a few months back.

In doing so, I noticed that MatPREP retains a lot of debug info regarding the findNoisyChannels process that PyPREP doesn't: e.g. the full matrices of regular correlations and RANSAC correlations for each channel/window, the HF noise z-scores for each channel, the fractions of "bad" RANSAC windows for each channel, etc. What I did yesterday was:

1. Write a Python script to import and parse the PREP debug info .mat files into Python/csv format
2. Adapt PyPREP to collect most of the same debug information
3. Modify the MatPREP script to only perform CleanLine and save the resulting data for each participant to a .set/.fdt
4. Modify my PyPREP script to import the .set/.fdt files, apply the EEGLAB montage, and perform robust re-referencing, saving the debug info to pickle files for later parsing. This theoretically allows for the closest comparison between PyPREP's NoisyChannels and MatPREP's possible, and also saves a ton of time waiting for cleanline during testing.

With the detailed debug information for both implementations saved such that they can be directly compared, hopefully I can get to the bottom of the NoisyChannels differences!

**EDIT**: From what I can tell, it looks like PyPREP and MatPREP start to diverge as soon as removeTrend is done: the values in the input data are both identical immediately before (I set breakpoints in both and checked), but immediately after the results diverge quite a bit:

 First sample of first 10 channels in MATLAB PREP before trend removal:
```
signal.data(:, 1)

ans =

   1.0e+04 *

   -0.3335
    1.0037
   -0.3581
   -0.1795
    0.9913
   -0.7414
   -0.7198
   -0.5665
    0.5557
   -1.2237
```

 First sample of first 10 channels in MATLAB PREP after trend removal:
```
signal2.data(:, 1)

ans =

   -3.7068
   -4.8012
   -0.8008
   -2.5809
   -1.5205
    0.1093
    1.9928
   -4.6906
  -10.2859
    1.6814
```

First sample of first 10 channels in PyPREP before trend removal:
```
raw_copy._data[:10, 0]
array([-0.00333541,  0.01003745, -0.00358109, -0.00179516,  0.00991251,
       -0.00741443, -0.00719762, -0.00566537,  0.00555661, -0.01223686])
```

First sample of first 10 channels in PyPREP after trend removal:
```
raw_filtered._data[:10, 0]
array([-4.11658012e-19, -1.79973325e-18,  6.98802181e-19,  2.89261751e-19,
        2.50128829e-18, -9.28348110e-19,  7.57247455e-19,  3.55753838e-19,
        3.82858892e-19, -7.01343280e-19])
```

As you can see, there's a huge order-of-magnitude shift with PyPREP's removeTrend that's much smaller in MatPREP, and the values also don't line up. Looks like our filter is doing something different from theirs, although I'm not quite sure what that is yet.
That's great that you got your hands on a MATLAB license! Looking forward to what you may find :-)
@sappelhoff Okay, so once #70 is merged and I dig a bit more into #55 to see if it's an easy fix, I think we'll be at the point where we can start adding MATLAB comparison for `NoisyChannels`, `removeTrend`, and potentially `Reference` into the CI! I guess the question now is what's the best way to do it.

Since Oct2Py seems to have been a dead-end in terms of PREP compatibility, here's my current rough proposal:

1. Create a separate Git repository with a MATLAB script for running all stages of PREP on a given file and saving the resulting `.set`/`.fdt`s following each stage (plus a `.mat` containing the extra info described in #68).
2. Set up CI for the above repository to run the script and save all the output files as GitHub Release artifacts.
3. Set up PyPREP to (optionally?) automatically grab the latest release artifacts from that repository during pytest and use the data for comparison

The advantages of this approach (in my mind) would be that we wouldn't need to re-run MATLAB PREP every time we run a CI job, and we also wouldn't need to depend on PyPREP's CI service (or anyone running unit tests locally) as having a properly-configured MATLAB environment. Separate `.set/.fdt` sets after each major step (post-cleanline, post-detrend, post-rereference, post-interpolation) would mean that a difference in an earlier part (e.g. CleanLine) won't automatically result in comparisons of later parts (e.g. NoisyChannels) failing, since we could load the relevant `.set/.fdt` before each to make sure PyPREP's starting with the same data for each stage.

I've already got a Python script that reads MATLAB PREP's extra info (saved into a .mat) and transforms it into a useful Python dict, so that part's already done, and my local MATLAB script for testing could be quickly adapted for the purposes described above. I think the main headache will be automating the setup of the MATLAB environment within a CI context: that said, apparently [GitHub Actions](https://github.com/matlab-actions/overview) and [Travis CI](https://docs.travis-ci.com/user/languages/matlab/) both have some form of MATLAB support, so I guess we'll see.



- +1 for using Matlab in GitHub actions ... let's resort to Travis only when GH actions turns out to be difficult
- re: separate repos: I think I would prefer everything to live in one (this) repo
- why not setup a matlab CI flow that installs EEGLAB and a recent version of PREP and then runs on a dataset of a dataset of our choice, and then the results are passed on to the next step in the CI where we run pyprep and then compare results?

Let me explain what I mean with a brief circleci example:

- in this CI step, we build the docs, and "store" the results: https://github.com/mne-tools/mne-bids/blob/7a4e65f9847ba6e94d1cd469e8472f024f6329e8/.circleci/config.yml#L35-L37
- in a later CI step, we get back our stored results and continue working on them: https://github.com/mne-tools/mne-bids/blob/7a4e65f9847ba6e94d1cd469e8472f024f6329e8/.circleci/config.yml#L68-L69
- I don't know what these are called in GH actions from the top of my head, but this concept exists for all serious CIs

>  Separate .set/.fdt sets after each major step 

I think that's a nice idea, my only concern is that we also start testing IO roundtrip accuracies there and I hope that won't trip us up. --> sidenote: Let's only use `.set` and not the set+fdt combo (I personally find the one file approach nicer, and with the new EEGLAB set style of storing the matlab EEGLAB data at the "root"  of the set file (which is a mat file) instead of an EEG "struct", one can also access metadata without loading the full binary data. FYI: https://github.com/mne-tools/mne-python/issues/8864#issuecomment-779375469
>  why not setup a matlab CI flow that installs EEGLAB and a recent version of PREP and then runs on a dataset of a dataset of our choice, and then the results are passed on to the next step in the CI where we run pyprep and then compare results?

Two questions on this:

1. Would the caching of build artifacts like this mean that it'd only re-run the MATLAB script any time it was changed, instead of every time the CI tests are run?
2. Would the artifacts be accessible from all the different CI bots, (Windows/macOS/Linux) or just the one it was run in? I ask because MATLAB via both GitHub Actions and Travis CI are both Linux-only. I guess the MATLAB PREP comparison could be its own separate runner in the CI matrix.

My main rationale behind the separate-repo-artifacts approach was that it would get around potential issues 1) and 2) easily, and also allow for running of local unit tests of MatPREP vs. PyPREP without a full MATLAB/EEGLAB/PREP install (using `urllib` or similar to just download the artifacts to a temporary folder for local tests). If 1 & 2 aren't issues, then doing everything in a single repository definitely makes sense.

> my only concern is that we also start testing IO roundtrip accuracies there and I hope that won't trip us up.

From my local testing, the EEG float values seem to be identical in an MNE-imported .set vs the same .set in EEGLAB (other than the volts/microVolts conversion) so this thankfully shouldn't be an issue!

> Let's only use .set and not the set+fdt combo

Ah cool, didn't even know that was possible! Single-file is definitely much cleaner to work with. The MATLAB/EEGLAB world is still very foreign to me...
> Would the caching of build artifacts like this mean that it'd only re-run the MATLAB script any time it was changed, instead of every time the CI tests are run?

no, that'd be an "every time" approach instead of computing once per release. Why would that be a problem? The results should stay the same if all versions are properly pinned, and a seed is set, right?

The only issue I could potentially see is (a) waiting time, (b) computation is expensive and pollutes the environment

> Would the artifacts be accessible from all the different CI bots, (Windows/macOS/Linux) or just the one it was run in?

no, pretty sure that the "artifacts" (in the way we talk about it currently) are only available to the current CI "bot" (matrix).

But I think comparing PREP and PyPREP for equivalence would be fine enough to do only under one specific set of parameters, as long as we continue to run all other unit tests under more diverse conditions.

> and also allow for running of local unit tests of MatPREP vs. PyPREP without a full MATLAB/EEGLAB/PREP install

yes, that is a fair point. Needs to be balanced against the increased amount of infrastructure (a whole separate repo). How would you want to do the new repo? under your user account? mine? a new pyprep GH account that we then also transfer sappelhoff/pyprep to?
> I think comparing PREP and PyPREP for equivalence would be fine enough to do only under one specific set of parameters

Agreed, I definitely think it would be overkill to multiple test files/configurations: setting up just one will be hard enough! I only mentioned artifacts across CIs to wonder if comparison tests could be on macOS/Windows runners too given that the MATLAB CI only works on Linux instances. I guess there'd be no reason to expect differences between OSes in terms of MATLAB differences, though.

> that'd be an "every time" approach instead of computing once per release. Why would that be a problem?

From my reading into the MATLAB GitHub Actions, [it looks like](https://github.com/matlab-actions/setup-matlab/blob/main/src/install.ts) it needs to download & install MATLAB on the instance every single time it's run, so I'd imagine it would slow down CI considerably if the results can't be cached for future runs.

Based on your thoughts/feedback, here's my current proposal: for now, I'll try and set up a MATLAB PREP script with CI and artifact saving in a new repository of my own, just to see what it's like in terms of speed/functionality. If it's pretty quick/painless and wouldn't add much time to the PyPREP CI, I'll take my work and make it into a PyPREP PR so everything's in one repository. If it's excessively slow and/or messy, we'll keep it as a separate repository that rarely needs to be run. If that happens, I can transfer the repository to you or a separate PyPREP organization (if you'd prefer control), but I'm also happy to keep it as my own if you'd prefer not to be responsible for its maintenance. Does that sound like a reasonable plan?

Sounds great! Good luck and let me know if I can help.
@sappelhoff Spent most of the day hacking away at this, and after much hacking and debugging I have a basic proof-of-concept working: https://github.com/a-hurst/matprep_artifacts/runs/2478319171?check_suite_focus=true

Right now it's just generating the .mat with the internal NoisyChannels info and it's still in a hacked-until-it-worked state, but definitely a relief to see the basics up and running!

EDIT: Got a chance to tidy it up, now it's got a README and generates single-file `.set` files at each main stage across the pipleine in addition the noisy_info `.mat`. All that's left is to set up a "release-on-tag" mechanism and finalize the choice of the test file (using S004R04.edf from the `eegbci` dataset as pyprep's current S004R01.edf it uses for tests). I opted for the other file because it's twice as long and thus there's a better chance of RANSAC having enough windows to detect bads properly, but I'm also open to S004R01.edf for speed/consistency.
I just took a look --> that looks awesome! Super clean / clear! I really like that you separated config/settings/setup and the .m script.

Choice of the data is up to you, I haven't worked with either of the datasets.

You could add a GH actions "badge" to the README that shows whether CI is passing and that acts as a link to the most recent run.
> I just took a look --> that looks awesome! Super clean / clear! I really like that you separated config/settings/setup and the .m script.

Thanks! One of my biggest long-standing annoyances with scripting in the sciences is that people like to bury all the configurable setting switches throughout their scripts (or don't make key things configurable at all) instead of having it all in one clear and obvious place, so I always try to make a  separate `settings.R/.m/.py` script whenever it makes sense.

If you've got no strong preference on the test file I'll move along with S004R04.edf and see how things go on the Python end! The nice thing about the script setup is that the test file could be changed at any time and it wouldn't affect the Python side at all: since it's using `1_matprep_raw.set` as its starting point, it'll always be starting with the same data as MATLAB PREP.
regarding the organization of repos etc., we could also:

- leave pyprep under sappelhoff
- leave matprep_artifacts under a-hurst

... but include `matprep_artifacts` in `pyprep` as a [git submodule](https://git-scm.com/book/en/v2/Git-Tools-Submodules) ... to keep a record of it within pyprep.

I thinking about this also in the direction of longevity of the project: For release 0.4 I'd want to switch on the [Zenodo tracking](https://guides.github.com/activities/citable-code/)  of the repo, i.e., to permanently store/archive a copy of the full code on an independent site. And it'd be cool if the matprep artifacts could be part of that archive, which I think they would via the git submodules solution. :thinking: :thought_balloon:  (but have never tested)

WDYT?


EDIT (even before sending off): Actually git submodules seem to not be included in the zenodo archive :confused: it's an open issue: https://github.com/zenodo/zenodo/issues/1036

So at least you know my thought process now. If anyone can come up with a good solution for proper archival and cross-referencing of "all things pyprep", that'd be great.
Hey @a-hurst this is amazing, you are doing a lot of great work for pyprep, thanks :-) I sent you a "collaborator" invitation for this GitHub repo, which will give you a couple of more permissions that will hopefully be helpful for organizing labels, milestones, reviews, etc.
@sappelhoff Glad I can help! I've definitely learned a lot through the process so far. Thanks for the collaborator invite, and for your helpful feedback on my posts and PRs!
@sappelhoff @yjmantilla do either of you have any data with channel drop-outs, just so I can test that and (hopefully) check it off? Or would that be really easy to simulate for a given recording? As an EEG novice (my background is almost entirely behavioural psych), I'm assuming that dropouts are what happens when the amplifier loses connection with an electrode for some reason, but I'm not sure how that ends up looking in the data.
> do either of you have any data with channel drop-outs, just so I can test that and (hopefully) check it off?

Mh I don't have such data at hand right now (nothing comes to mind), but for the data that I worked with (BrainVision ".vhdr+.vmrk+.eeg" [format](https://github.com/bids-standard/pybv)), I never encountered NaN values during recording or analysis.

I did have some recordings where an electrode got loose and then connection was bad ... and other cases where an electrode itself was broken.

That usually resulted in very flat (very low variance and amplitude) signals.
> but for the data that I worked with (BrainVision ".vhdr+.vmrk+.eeg" format), I never encountered NaN values during recording or analysis.

I don't think EDF/BDF files have any concept of NaN either, since all data is stored exclusively as integers. Maybe EEGLAB generates NaNs under certain circumstances that MNE doesn't? If that's the case, maybe for PyPREP dropouts should be based on windows where the signal is totally flat instead of windows where the signal is NaN. The original PREP publication reports a dataset they used where they had considerable dropouts, so maybe I should fire them an email...
Once numerical equivalence with matlab PREP is reached, we can consider removing the "ALPHA SOFTWARE" warning from the README:

> ALPHA SOFTWARE. This package is currently in its early stages of iteration. It may change both its internals or its user-facing API in the near future. Any feedback and ideas on how to improve either of these is more than welcome! Use this software at your own risk.
Seems I unpinned this issue without me noticing (maybe I miss-clicked somewhere (?)). Pinning it back.

Fix patch:
diff --git a/docs/whats_new.rst b/docs/whats_new.rst
index 7984feae..22107be7 100644
--- a/docs/whats_new.rst
+++ b/docs/whats_new.rst
@@ -45,6 +45,9 @@ Changelog
 - Added `max_chunk_size` parameter for specifying the maximum chunk size to use for channel-wise RANSAC, allowing more control over PyPREP RAM usage, by `Austin Hurst`_ (:gh:`66`)
 - Changed :class:`~pyprep.Reference` to exclude "bad-by-SNR" channels from initial average referencing, matching MATLAB PREP behaviour, by `Austin Hurst`_ (:gh:`78`)
 - Changed :class:`~pyprep.Reference` to only flag "unusable" channels (bad by flat, NaNs, or low SNR) from the first pass of noisy detection for permanent exclusion from the reference signal, matching MATLAB PREP behaviour, by `Austin Hurst`_ (:gh:`78`)
+- Added a framework for automated testing of PyPREP's components against their MATLAB PREP counterparts (using ``.mat`` and ``.set`` files generated with the `matprep_artifacts`_ script), helping verify that the two PREP implementations are numerically equivalent when `matlab_strict` is ``True``, by `Austin Hurst`_ (:gh:`79`)
+
+.. matprep_artifacts: https://github.com/a-hurst/matprep_artifacts
 
 Bug
 ~~~
@@ -54,6 +57,7 @@ Bug
 - 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`)
 - Fixed a bug in :meth:`NoisyChannels.find_bad_by_correlation` that prevented it from being able to handle channels with dropouts (intermittent flat regions), by `Austin Hurst`_ (:gh:`81`).
+- Fixed :class:`~pyprep.NoisyChannels` so that it always runs "bad channel by NaN" and "bad channel by flat" detection, preventing these channels from causing problems with other :class:`~pyprep.NoisyChannels` methods, by `Austin Hurst`_ (:gh:`79`)
 
 API
 ~~~
diff --git a/pyprep/find_noisy_channels.py b/pyprep/find_noisy_channels.py
index 045bbca2..83757676 100644
--- a/pyprep/find_noisy_channels.py
+++ b/pyprep/find_noisy_channels.py
@@ -49,6 +49,7 @@ def __init__(self, raw, do_detrend=True, random_state=None, matlab_strict=False)
         assert isinstance(raw, mne.io.BaseRaw)
 
         self.raw_mne = raw.copy()
+        self.raw_mne.pick_types(eeg=True)
         self.sample_rate = raw.info["sfreq"]
         if do_detrend:
             self.raw_mne._data = removeTrend(
@@ -56,18 +57,6 @@ def __init__(self, raw, do_detrend=True, random_state=None, matlab_strict=False)
             )
         self.matlab_strict = matlab_strict
 
-        self.EEGData = self.raw_mne.get_data(picks="eeg")
-        self.EEGData_beforeFilt = self.EEGData
-        self.ch_names_original = np.asarray(raw.info["ch_names"])
-        self.n_chans_original = len(self.ch_names_original)
-        self.n_chans_new = self.n_chans_original
-        self.original_dimensions = np.shape(self.EEGData)
-        self.new_dimensions = self.original_dimensions
-        self.original_channels = np.arange(self.original_dimensions[0])
-        self.new_channels = self.original_channels
-        self.ch_names_new = self.ch_names_original
-        self.channels_interpolate = self.original_channels
-
         # Extra data for debugging
         self._extra_info = {
             'bad_by_deviation': {},
@@ -90,6 +79,25 @@ def __init__(self, raw, do_detrend=True, random_state=None, matlab_strict=False)
         self.bad_by_dropout = []
         self.bad_by_ransac = []
 
+        # Get original EEG channel names, channel count & samples
+        ch_names = np.asarray(self.raw_mne.info["ch_names"])
+        self.ch_names_original = ch_names
+        self.n_chans_original = len(ch_names)
+        self.n_samples = raw._data.shape[1]
+
+        # Before anything else, flag bad-by-NaNs and bad-by-flats
+        self.find_bad_by_nan_flat()
+        bads_by_nan_flat = self.bad_by_nan + self.bad_by_flat
+
+        # Make a subset of the data containing only usable EEG channels
+        self.usable_idx = np.isin(ch_names, bads_by_nan_flat, invert=True)
+        self.EEGData = self.raw_mne.get_data(picks=ch_names[self.usable_idx])
+        self.EEGData_beforeFilt = self.EEGData
+
+        # Get usable EEG channel names & channel counts
+        self.ch_names_new = np.asarray(ch_names[self.usable_idx])
+        self.n_chans_new = len(self.ch_names_new)
+
     def get_bads(self, verbose=False):
         """Get a list of all bad channels.
 
@@ -186,35 +194,23 @@ def find_all_bads(self, ransac=True, channel_wise=False, max_chunk_size=None):
 
     def find_bad_by_nan_flat(self):
         """Detect channels that appear flat or have NaN values."""
-        nan_channel_mask = [False] * self.original_dimensions[0]
-        no_signal_channel_mask = [False] * self.original_dimensions[0]
+        # Get all EEG channels from original copy of data
+        EEGData = self.raw_mne.get_data()
 
-        FLAT_THRESHOLD = 1e-15  # corresponds to 10e-10 µV in MATLAB PREP
-        for i in range(0, self.original_dimensions[0]):
-            nan_channel_mask[i] = np.sum(np.isnan(self.EEGData[i, :])) > 0
-        for i in range(0, self.original_dimensions[0]):
-            no_signal_channel_mask[i] = (
-                robust.mad(self.EEGData[i, :], c=1) < FLAT_THRESHOLD
-                or np.std(self.EEGData[i, :]) < FLAT_THRESHOLD
-            )
-        nan_channels = self.channels_interpolate[nan_channel_mask]
-        flat_channels = self.channels_interpolate[no_signal_channel_mask]
+        # Detect channels containing any NaN values
+        nan_channel_mask = np.isnan(np.sum(EEGData, axis=1))
+        nan_channels = self.ch_names_original[nan_channel_mask]
 
-        nans_no_data_channels = np.union1d(nan_channels, flat_channels)
-        self.channels_interpolate = np.setdiff1d(
-            self.channels_interpolate, nans_no_data_channels
-        )
+        # Detect channels with flat or extremely weak signals
+        FLAT_THRESHOLD = 1e-15  # corresponds to 10e-10 µV in MATLAB PREP
+        flat_by_mad = robust.mad(EEGData, axis=1, c=1) < FLAT_THRESHOLD
+        flat_by_stdev = np.std(EEGData, axis=1) < FLAT_THRESHOLD
+        flat_channel_mask = flat_by_mad | flat_by_stdev
+        flat_channels = self.ch_names_original[flat_channel_mask]
 
-        for i in range(0, len(nan_channels)):
-            self.bad_by_nan.append(self.ch_names_original[nan_channels[i]])
-        for i in range(0, len(flat_channels)):
-            self.bad_by_flat.append(self.ch_names_original[flat_channels[i]])
-
-        self.raw_mne.drop_channels(list(set(self.bad_by_nan + self.bad_by_flat)))
-        self.EEGData = self.raw_mne.get_data(picks="eeg")
-        self.ch_names_new = np.asarray(self.raw_mne.info["ch_names"])
-        self.n_chans_new = len(self.ch_names_new)
-        self.new_dimensions = np.shape(self.EEGData)
+        # Update names of bad channels by NaN or flat signal
+        self.bad_by_nan = nan_channels.tolist()
+        self.bad_by_flat = flat_channels.tolist()
 
     def find_bad_by_deviation(self, deviation_threshold=5.0):
         """Robust z-score of the robust standard deviation for each channel is calculated.
@@ -227,22 +223,22 @@ def find_bad_by_deviation(self, deviation_threshold=5.0):
             z-score threshold above which channels will be labelled bad.
 
         """
-        deviation_channel_mask = [False] * (self.new_dimensions[0])
-        channel_deviation = np.zeros(self.new_dimensions[0])
-        for i in range(0, self.new_dimensions[0]):
+        deviation_channel_mask = [False] * (self.n_chans_original)
+        robust_channel_deviation = np.zeros(self.n_chans_original, dtype=np.float64)
+        channel_deviation = np.zeros(self.n_chans_new)
+        for i in range(0, self.n_chans_new):
             channel_deviation[i] = 0.7413 * _mat_iqr(self.EEGData[i, :])
         channel_deviationSD = 0.7413 * _mat_iqr(channel_deviation)
         channel_deviationMedian = np.nanmedian(channel_deviation)
-        robust_channel_deviation = np.divide(
+        robust_channel_deviation[self.usable_idx] = np.divide(
             np.subtract(channel_deviation, channel_deviationMedian), channel_deviationSD
         )
-        for i in range(0, self.new_dimensions[0]):
+        for i in range(0, self.n_chans_original):
             deviation_channel_mask[i] = abs(
                 robust_channel_deviation[i]
             ) > deviation_threshold or np.isnan(robust_channel_deviation[i])
-        deviation_channels = self.channels_interpolate[deviation_channel_mask]
-        for i in range(0, len(deviation_channels)):
-            self.bad_by_deviation.append(self.ch_names_original[deviation_channels[i]])
+        deviation_channels = self.ch_names_original[deviation_channel_mask]
+        self.bad_by_deviation = deviation_channels.tolist()
         self._extra_info['bad_by_deviation'].update({
             'median_channel_deviation': channel_deviationMedian,
             'channel_deviation_sd': channel_deviationSD,
@@ -266,6 +262,7 @@ def find_bad_by_hfnoise(self, HF_zscore_threshold=5.0):
         """
         data_tmp = np.transpose(self.EEGData)
         dimension = np.shape(data_tmp)
+        zscore_HFNoise = np.zeros(self.n_chans_original, dtype=np.float64)
         if self.sample_rate > 100:
             EEG_filt = np.zeros((dimension[0], dimension[1]))
             bandpass_filter = filter_design(
@@ -284,25 +281,19 @@ def find_bad_by_hfnoise(self, HF_zscore_threshold=5.0):
                 np.median(np.absolute(np.subtract(noisiness, np.median(noisiness))))
                 * 1.4826
             )
-            zscore_HFNoise = np.divide(
+            zscore_HFNoise[self.usable_idx] = np.divide(
                 np.subtract(noisiness, noisiness_median), noiseSD
             )
-            HFnoise_channel_mask = [False] * self.new_dimensions[0]
-            for i in range(0, self.new_dimensions[0]):
-                HFnoise_channel_mask[i] = zscore_HFNoise[
-                    i
-                ] > HF_zscore_threshold or np.isnan(zscore_HFNoise[i])
-            HFNoise_channels = self.channels_interpolate[HFnoise_channel_mask]
+            hf_mask = np.isnan(zscore_HFNoise) | (zscore_HFNoise > HF_zscore_threshold)
+            HFNoise_channels = self.ch_names_original[hf_mask]
         else:
             EEG_filt = data_tmp
             noisiness_median = 0
-            # noisinessSD = 1
-            zscore_HFNoise = np.zeros((self.new_dimensions[0], 1))
-            HFNoise_channels = []
+            noiseSD = 1
+            HFNoise_channels = np.asarray([])
         self.EEGData_beforeFilt = data_tmp
         self.EEGData = np.transpose(EEG_filt)
-        for i in range(0, len(HFNoise_channels)):
-            self.bad_by_hf_noise.append(self.ch_names_original[HFNoise_channels[i]])
+        self.bad_by_hf_noise = HFNoise_channels.tolist()
         self._extra_info['bad_by_hf_noise'].update({
             'median_channel_noisiness': noisiness_median,
             'channel_noisiness_sd': noiseSD,
@@ -337,30 +328,29 @@ def find_bad_by_correlation(
 
         """
         self.find_bad_by_hfnoise()  # since filtering is performed there
+        usable_idx = self.usable_idx
         correlation_frames = correlation_secs * self.sample_rate
         correlation_window = np.arange(correlation_frames)
         correlation_offsets = np.arange(
-            1, (self.new_dimensions[1] - correlation_frames), correlation_frames
+            1, (self.n_samples - correlation_frames), correlation_frames
         )
         w_correlation = len(correlation_offsets)
-        maximum_correlations = np.ones((self.original_dimensions[0], w_correlation))
-        drop_out = np.zeros((self.new_dimensions[0], w_correlation))
-        channel_correlation = np.ones((w_correlation, self.new_dimensions[0]))
-        noiselevels = np.zeros((w_correlation, self.new_dimensions[0]))
-        channel_deviations = np.zeros((w_correlation, self.new_dimensions[0]))
-        drop = np.zeros((w_correlation, self.new_dimensions[0]))
+        channel_correlation = np.ones((w_correlation, self.n_chans_original))
+        noiselevels = np.zeros((w_correlation, self.n_chans_original))
+        channel_deviations = np.zeros((w_correlation, self.n_chans_original))
+        drop = np.zeros((w_correlation, self.n_chans_original))
         len_correlation_window = len(correlation_window)
         EEGData = np.transpose(self.EEGData)
         EEG_new_win = np.reshape(
             np.transpose(EEGData[0 : len_correlation_window * w_correlation, :]),
-            (self.new_dimensions[0], len_correlation_window, w_correlation),
+            (self.n_chans_new, len_correlation_window, w_correlation),
             order="F",
         )
         data_win = np.reshape(
             np.transpose(
                 self.EEGData_beforeFilt[0 : len_correlation_window * w_correlation, :]
             ),
-            (self.new_dimensions[0], len_correlation_window, w_correlation),
+            (self.n_chans_new, len_correlation_window, w_correlation),
             order="F",
         )
         for k in range(0, w_correlation):
@@ -371,41 +361,38 @@ def find_bad_by_correlation(
             abs_corr = np.abs(
                 np.subtract(window_correlation, np.diag(np.diag(window_correlation)))
             )
-            channel_correlation[k, :] = _mat_quantile(abs_corr, 0.98, axis=0)
+            channel_correlation[k, usable_idx] = _mat_quantile(abs_corr, 0.98, axis=0)
             with np.errstate(invalid='ignore'):  # suppress divide-by-zero warnings
-                noiselevels[k, :] = np.divide(
+                noiselevels[k, usable_idx] = np.divide(
                     robust.mad(np.subtract(data_portion, eeg_portion), c=1),
                     robust.mad(eeg_portion, c=1),
                 )
-            channel_deviations[k, :] = 0.7413 * _mat_iqr(data_portion, axis=0)
+            channel_deviations[k, usable_idx] = 0.7413 * _mat_iqr(data_portion, axis=0)
         for i in range(0, w_correlation):
-            for j in range(0, self.new_dimensions[0]):
+            for j in range(0, self.n_chans_original):
                 drop[i, j] = np.int(
                     np.isnan(channel_correlation[i, j]) or np.isnan(noiselevels[i, j])
                 )
                 if drop[i, j] == 1:
                     channel_correlation[i, j] = 0
                     noiselevels[i, j] = 0
-        maximum_correlations[self.channels_interpolate, :] = np.transpose(
-            channel_correlation
-        )
-        drop_out[:] = np.transpose(drop)
+
+        # Flag channels with above-threshold fractions of bad correlation windows
+        maximum_correlations = np.transpose(channel_correlation)
         thresholded_correlations = maximum_correlations < correlation_threshold
         thresholded_correlations = thresholded_correlations.astype(int)
         fraction_BadCorrelationWindows = np.mean(thresholded_correlations, axis=1)
+        bad_correlation_mask = fraction_BadCorrelationWindows > frac_bad
+        bad_correlation_channels = self.ch_names_original[bad_correlation_mask]
+        self.bad_by_correlation = bad_correlation_channels.tolist()
+
+        # Flag channels with above-threshold fractions of drop-out windows
+        drop_out = np.transpose(drop)
         fraction_BadDropOutWindows = np.mean(drop_out, axis=1)
+        dropout_mask = fraction_BadDropOutWindows > frac_bad
+        dropout_channels = self.ch_names_original[dropout_mask]
+        self.bad_by_dropout = dropout_channels.tolist()
 
-        bad_correlation_channels_idx = np.argwhere(
-            fraction_BadCorrelationWindows > frac_bad
-        )
-        bad_correlation_channels_name = self.ch_names_original[
-            bad_correlation_channels_idx.astype(int)
-        ]
-        self.bad_by_correlation = [i[0] for i in bad_correlation_channels_name]
-
-        dropout_channels_idx = np.argwhere(fraction_BadDropOutWindows > frac_bad)
-        dropout_channels_name = self.ch_names_original[dropout_channels_idx.astype(int)]
-        self.bad_by_dropout = [i[0] for i in dropout_channels_name]
         self._extra_info['bad_by_correlation'] = {
             'max_correlations': maximum_correlations,
             'median_max_correlations': np.median(maximum_correlations, axis=1),
@@ -503,11 +490,11 @@ def find_bad_by_ransac(
             self.bad_by_deviation +
             self.bad_by_dropout
         )
-        self.bad_by_ransac, ch_correlations = find_bad_by_ransac(
+        self.bad_by_ransac, ch_correlations_usable = find_bad_by_ransac(
             self.EEGData,
             self.sample_rate,
             self.ch_names_new,
-            self.raw_mne._get_channel_positions(),
+            self.raw_mne._get_channel_positions()[self.usable_idx, :],
             exclude_from_ransac,
             n_samples,
             fraction_good,
@@ -519,6 +506,12 @@ def find_bad_by_ransac(
             self.random_state,
             self.matlab_strict,
         )
+
+        # Reshape correlation matrix to match original channel count
+        n_ransac_windows = ch_correlations_usable.shape[0]
+        ch_correlations = np.ones((n_ransac_windows, self.n_chans_original))
+        ch_correlations[:, self.usable_idx] = ch_correlations_usable
+
         self._extra_info['bad_by_ransac'] = {
             'ransac_correlations': ch_correlations,
             'bad_window_fractions': np.mean(ch_correlations < corr_thresh, axis=0)


Test patch:
diff --git a/tests/test_find_noisy_channels.py b/tests/test_find_noisy_channels.py
index 3246fad1..8d45ae3a 100644
--- a/tests/test_find_noisy_channels.py
+++ b/tests/test_find_noisy_channels.py
@@ -118,6 +118,12 @@ def test_findnoisychannels(raw, montage):
     nd.find_bad_by_hfnoise()
     assert rand_chn_lab in nd.bad_by_hf_noise
 
+    # Test for high freq noise detection when sample rate < 100 Hz
+    raw_tmp.resample(80)  # downsample to 80 Hz
+    nd = NoisyChannels(raw_tmp, random_state=rng)
+    nd.find_bad_by_hfnoise()
+    assert len(nd.bad_by_hf_noise) == 0
+
     # Test for signal to noise ratio in EEG data
     raw_tmp = raw.copy()
     m, n = raw_tmp._data.shape
diff --git a/tests/test_matprep_compare.py b/tests/test_matprep_compare.py
new file mode 100644
index 00000000..9f191eaf
--- /dev/null
+++ b/tests/test_matprep_compare.py
@@ -0,0 +1,305 @@
+"""Compare PyPREP results to MATLAB PREP."""
+from urllib.request import urlopen
+import numpy as np
+import scipy
+import mne
+
+import pytest
+
+from pyprep.find_noisy_channels import NoisyChannels
+from pyprep.removeTrend import removeTrend
+
+
+# Define some fixtures for things that will be used across multiple tests
+
+@pytest.fixture(scope='session')
+def matprep_artifacts(tmpdir_factory):
+    """Fixture for downloading & using MATLAB PREP artifacts from CI.
+
+    Downloads the latest set of CI-generated MATLAB PREP artifacts and saves
+    them to a temporary folder for use with various tests, returning the paths
+    to each artifact in a {name: path} dict. The temporary folder will be
+    automatically cleaned up by pytest once all tests have completed.
+
+    """
+    base_url = "https://github.com/a-hurst/matprep_artifacts/releases/latest/download/"
+    artifacts = [
+        "1_matprep_raw.set",
+        "2_matprep_removetrend.set",
+        "3_matprep_cleanline.set",
+        "4_matprep_pre_reference.set",
+        "5_matprep_post_reference.set",
+        "matprep_info.mat"
+    ]
+    dirpath = tmpdir_factory.mktemp("artifacts")
+
+    artifact_paths = {}
+    for f in artifacts:
+        artifact_name = f.split(".")[0]
+        outfile = str(dirpath.join(f))
+        dl = urlopen(base_url + f)
+        with open(outfile, 'wb') as out:
+            out.write(dl.read())
+        artifact_paths[artifact_name] = outfile
+
+    return artifact_paths
+
+
+@pytest.fixture(scope='session')
+def matprep_noisy(matprep_artifacts):
+    """Import and preprocess artifact containing MATLAB PREP runtime info.
+
+    This fixture only parses and retains data from the first pass of noisy
+    channel detection during re-referencing, since it's easiest to compare with
+    PyPREP. It also adds a new key to the imported struct, 'bads', which
+    contains the names of the channels flagged as bad by each detection
+    method (as opposed to just channel indices).
+
+    """
+    # Read in and parse noisy channel info artifact from MATLAB PREP
+    info_path = matprep_artifacts['matprep_info']
+    matprep_info = scipy.io.loadmat(info_path, simplify_cells=True)['prep_info']
+    matprep_noisy_all = matprep_info['reference']
+    matprep_noisy = matprep_noisy_all['noisyStatisticsOriginal']
+
+    # Gather bad channel names from MatPREP, converting numbers to labels
+    ch_names = matprep_info['originalChannelLabels']
+    bad_types = {
+        'badChannelsFromNaNs': 'by_nan',
+        'badChannelsFromNoData': 'by_flat',
+        'badChannelsFromDeviation': 'by_deviation',
+        'badChannelsFromHFNoise': 'by_hf_noise',
+        'badChannelsFromCorrelation': 'by_correlation',
+        'badChannelsFromLowSNR': 'by_SNR',
+        'badChannelsFromDropOuts': 'by_dropout',
+        'badChannelsFromRansac': 'by_ransac',
+        'all': 'all',
+    }
+    matprep_bads = {}
+    for bad_type, name in bad_types.items():
+        bads_idx = matprep_noisy['noisyChannels'][bad_type]
+        bads_idx = [bads_idx] if isinstance(bads_idx, int) else bads_idx
+        matprep_bads[name] = [ch_names[i-1] for i in bads_idx]
+    matprep_noisy['bads'] = matprep_bads
+
+    return matprep_noisy
+
+
+@pytest.fixture(scope='session')
+def pyprep_noisy(matprep_artifacts):
+    """Get original NoisyChannels results for comparison with MATLAB PREP.
+
+    This fixture uses an artifact from MATLAB PREP of the CleanLined and
+    detrended EEG signal right before MATLAB PREP runs its first iteration of
+    NoisyChannels during re-referncing. As such, any differences in test results
+    will be due to actual differences in the noisy channel detection code rather
+    than differences at an earlier stage of the pipeline.
+
+    """
+    # Import pre-reference MATLAB PREP data
+    preref_path = matprep_artifacts['4_matprep_pre_reference']
+    matprep_preref = mne.io.read_raw_eeglab(preref_path, preload=True)
+
+    # Run NoisyChannels on MATLAB data and extract internal noisy info
+    matprep_seed = 435656
+    pyprep_noisy = NoisyChannels(
+        matprep_preref,
+        do_detrend=False,
+        random_state=matprep_seed,
+        matlab_strict=True
+    )
+    pyprep_noisy.find_all_bads()
+
+    return pyprep_noisy
+
+
+# Define MATLAB comparison tests for each main component of PyPREP
+
+@pytest.mark.usefixtures('matprep_artifacts')
+def test_compare_removeTrend(matprep_artifacts):
+    """Test the numeric equivalence of removeTrend to MATLAB PREP."""
+    # Get paths of MATLAB .set files
+    raw_path = matprep_artifacts['1_matprep_raw']
+    detrend_path = matprep_artifacts['2_matprep_removetrend']
+
+    # Load relevant MATLAB data
+    matprep_raw = mne.io.read_raw_eeglab(raw_path, preload=True)
+    matprep_detrended = mne.io.read_raw_eeglab(detrend_path, preload=True)
+    sample_rate = matprep_raw.info["sfreq"]
+
+    # Apply removeTrend to raw artifact to get expected and actual signals
+    expected = matprep_detrended._data
+    actual = removeTrend(
+        matprep_raw._data, sample_rate, detrendType="high pass", matlab_strict=True
+    )
+
+    # Check MATLAB equivalence at start of recording
+    win_size = 500  # window of samples to check
+    assert np.allclose(
+        actual[:, :win_size], expected[:, :win_size], equal_nan=True
+    )
+
+    # Check MATLAB equivalence in middle of recording
+    win_start = int(actual.shape[1] / 2)
+    win_end = win_start + win_size
+    assert np.allclose(
+        actual[:, win_start:win_end], expected[:, win_start:win_end], equal_nan=True
+    )
+
+    # Check MATLAB equivalence at end of recording
+    assert np.allclose(
+        actual[:, -win_size:], expected[:, -win_size:], equal_nan=True
+    )
+
+
+class TestCompareNoisyChannels(object):
+    """Compare the results of NoisyChannels to the equivalent MatPREP code.
+
+    These comparisons use input data that's already had adaptive line noise
+    removal and high-pass trend removal done to the signal, so any differences
+    found will be due to differences in the noisy channel detection code rather
+    than any differences at an earlier stage of the pipeline.
+
+    """
+
+    def test_bad_by_nan(self, pyprep_noisy, matprep_noisy):
+        """Compare bad-by-NaN results between PyPREP and MatPREP."""
+        # NOTE: The current test artifacts contain no channels with NaN values
+        # (when does that ever happen?), meaning that this may not be testing
+        # anything useful
+        assert pyprep_noisy.bad_by_nan == matprep_noisy['bads']['by_nan']
+
+    def test_bad_by_flat(self, pyprep_noisy, matprep_noisy):
+        """Compare bad-by-flat results between PyPREP and MatPREP."""
+        # NOTE: The current test artifacts contain no flat channels, meaning
+        # that this may not be testing anything useful
+        assert pyprep_noisy.bad_by_flat == matprep_noisy['bads']['by_flat']
+
+    def test_bad_by_deviation(self, pyprep_noisy, matprep_noisy):
+        """Compare bad-by-deviation results between PyPREP and MatPREP."""
+        # Gather PyPREP deviation info and MATLAB equivalents
+        mat = matprep_noisy
+        matprep_info = {
+            'median_channel_deviation': mat['channelDeviationMedian'],
+            'channel_deviation_sd': mat['channelDeviationSD'],
+            'robust_channel_deviations': mat['robustChannelDeviation'],
+            'channel_deviations': mat['channelDeviations']
+        }
+        pyprep_info = pyprep_noisy._extra_info['bad_by_deviation']
+
+        # Compare overall medians and SDs for channel deviations
+        median_dev = 'median_channel_deviation'
+        dev_sd = 'channel_deviation_sd'
+        assert np.isclose(pyprep_info[median_dev] * 1e6, matprep_info[median_dev])
+        assert np.isclose(pyprep_info[dev_sd] * 1e6, matprep_info[dev_sd])
+
+        # Compare robust deviations across channels
+        dev_by_chan = 'robust_channel_deviations'
+        assert np.allclose(pyprep_info[dev_by_chan], matprep_info[dev_by_chan])
+
+        # Compare windows of channel deviations across recording
+        chan_devs = 'channel_deviations'
+        assert np.allclose(pyprep_info[chan_devs] * 1e6, matprep_info[chan_devs].T)
+
+        # Compare names of bad-by-deviation channels
+        assert pyprep_noisy.bad_by_deviation == matprep_noisy['bads']['by_deviation']
+
+    def test_bad_by_hf_noise(self, pyprep_noisy, matprep_noisy):
+        """Compare bad-by-HF-noise results between PyPREP and MatPREP."""
+        # Gather PyPREP high-frequency noise info and MATLAB equivalents
+        mat = matprep_noisy
+        matprep_info = {
+            'median_channel_noisiness': mat['noisinessMedian'],
+            'channel_noisiness_sd': mat['noisinessSD'],
+            'hf_noise_zscores': mat['zscoreHFNoise'],
+            'noise_levels': mat['noiseLevels']
+        }
+        pyprep_info = pyprep_noisy._extra_info['bad_by_hf_noise']
+
+        # Compare overall medians and SDs for channel noisiness
+        median_noise = 'median_channel_noisiness'
+        noise_sd = 'channel_noisiness_sd'
+        assert np.isclose(pyprep_info[median_noise], matprep_info[median_noise])
+        assert np.isclose(pyprep_info[noise_sd], matprep_info[noise_sd])
+
+        # Compare noisiness z-scores across channels
+        TOL = 1e-5  # NOTE: Some diffs > 1e-6 (default), maybe slight diff in math?
+        noise_z = 'hf_noise_zscores'
+        assert np.allclose(pyprep_info[noise_z], matprep_info[noise_z], atol=TOL)
+
+        # Compare windows of noisiness per channel across recording
+        noise_win = 'noise_levels'
+        assert np.allclose(pyprep_info[noise_win], matprep_info[noise_win].T)
+
+        # Compare names of bad-by-HF-noise channels
+        assert pyprep_noisy.bad_by_hf_noise == matprep_noisy['bads']['by_hf_noise']
+
+    def test_bad_by_correlation(self, pyprep_noisy, matprep_noisy):
+        """Compare bad-by-correlation results between PyPREP and MatPREP."""
+        # Gather PyPREP correlation info and MATLAB equivalents
+        mat = matprep_noisy
+        matprep_info = {
+            'max_correlations': mat['maximumCorrelations'],
+            'median_max_correlations': mat['medianMaxCorrelation']
+        }
+        pyprep_info = pyprep_noisy._extra_info['bad_by_correlation']
+
+        # Compare median maximum correlations across channels
+        med_max_corr = 'median_max_correlations'
+        assert np.allclose(pyprep_info[med_max_corr], matprep_info[med_max_corr])
+
+        # Compare max correlations per channel across recording
+        max_corr = 'max_correlations'
+        assert np.allclose(pyprep_info[max_corr], matprep_info[max_corr])
+
+        # Compare names of bad-by-correlation channels
+        matprep_bad_by_corr = matprep_noisy['bads']['by_correlation']
+        assert pyprep_noisy.bad_by_correlation == matprep_bad_by_corr
+
+    def test_bad_by_SNR(self, pyprep_noisy, matprep_noisy):
+        """Compare bad-by-SNR results between PyPREP and MatPREP."""
+        pyprep_bads_snr = sorted(pyprep_noisy.bad_by_SNR)
+        matprep_bads_snr = sorted(matprep_noisy['bads']['by_SNR'])
+        assert pyprep_bads_snr == matprep_bads_snr
+
+    def test_bad_by_dropout(self, pyprep_noisy, matprep_noisy):
+        """Compare bad-by-dropout results between PyPREP and MatPREP."""
+        # NOTE: The current test artifacts contain no channels with dropouts,
+        # meaning that this may not be testing anything useful
+
+        # Gather PyPREP and MATLAB PREP dropout info
+        matprep_dropouts = matprep_noisy['dropOuts']
+        pyprep_dropouts = pyprep_noisy._extra_info['bad_by_dropout']['dropouts']
+
+        # Compare dropout windows per channel across recording
+        assert np.allclose(pyprep_dropouts, matprep_dropouts)
+
+        # Compare names of bad-by-dropout channels
+        assert pyprep_noisy.bad_by_dropout == matprep_noisy['bads']['by_dropout']
+
+    def test_bad_by_ransac(self, pyprep_noisy, matprep_noisy):
+        """Compare bad-by-RANSAC results between PyPREP and MatPREP."""
+        # Gather PyPREP RANSAC correlation info and MATLAB equivalents
+        mat = matprep_noisy
+        matprep_info = {
+            'ransac_correlations': mat['ransacCorrelations'],
+            'bad_window_fractions': mat['ransacBadWindowFraction']
+        }
+        pyprep_info = pyprep_noisy._extra_info['bad_by_ransac']
+
+        # Compare fractions of bad RANSAC windows across channels
+        bad_fracs = 'bad_window_fractions'
+        assert np.allclose(pyprep_info[bad_fracs], matprep_info[bad_fracs])
+
+        # Compare RANSAC window correlations per channel across recording
+        ransac_corr = 'ransac_correlations'
+        assert np.allclose(pyprep_info[ransac_corr], matprep_info[ransac_corr].T)
+
+        # Compare names of bad-by-RANSAC channels
+        assert pyprep_noisy.bad_by_ransac == matprep_noisy['bads']['by_ransac']
+
+    def test_all_bads(self, pyprep_noisy, matprep_noisy):
+        """Compare names of all bad channels between PyPREP and MatPREP."""
+        pyprep_bads_all = sorted(pyprep_noisy.get_bads())
+        matprep_bads_all = sorted(matprep_noisy['bads']['all'])
+        assert pyprep_bads_all == matprep_bads_all