Repository: sappelhoff/pyprep
Pull request #66 (2021-04-16 18:19:20)

Add window-wise RANSAC to reduce RAM requirements
<!--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

(Updated) This PR does a few things:

- Adds support for window-wise RANSAC and makes it the default method
- Reworks channel-wise RANSAC to use the same pre-generated interpolation matrices as window-wise, which speeds things up and reduces the amount of divergent/redundant code between methods
- Reworks/refactors the channel-wise RANSAC functions a fair bit to better match the function names, argument names, variable names, and docstrings of their window-wise equivalence, for easier comparison and maintainability
- Adds a progress-printing function for window-wise RANSAC
- Overhauls the RANSAC unit tests considerably, including the addition of explicit comparison of the window-wise/channel-wise correlation matrices for both regular and MATLAB-strict RANSAC to make sure we get identical results for the same seed between methods (and that MATLAB strict values differ from non-strict ones as expected).

### Currently Missing
- [x] Progress indicators for RANSAC
- [x] Proper API for `matlab_strict` mode (handled in earlier PR)
- [ ] Any sort of optimization/parallelization for additional speed (I have some ideas, but will save for a future PR since this one's already pretty big)
- [x] Any sort of memory checking or channel chunking support (added memory check at start)
- [x] A clear idea of what to do about MATLAB PREP's weird correlation function (handled in earlier PR)

### Other Notes
- For my test file (`sub_002/ses-01` from that open meditation dataset), window-wise RANSAC is about 36% slower than the existing RANSAC approach (192 vs 141 seconds). However, RAM usage during RANSAC is now much lower (maxed out at 2.85 GB), which is nice but also suggests that window-wise RANSAC could be improved to take better advantage of RAM on beefier systems.

# 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


a-hurst: Whoops, accidentally submitted the PR before writing a description! Fixing that now.
sappelhoff: gonna do a proper review later, but regarding the weird matprep correlation --> is that maybe because the data are already mean-centered? ... i.e., mean(x) is already around 0, so x - mean(x) is unnecessary and you can just square x directly? Just a hunch ... could be completely off.
a-hurst: > gonna do a proper review later, but regarding the weird matprep correlation --> is that maybe because the data are already mean-centered? ... i.e., mean(x) is already around 0, so x - mean(x) is unnecessary and you can just square x directly? Just a hunch ... could be completely off.

That sounds like a good explanation. That said, it is generally going to be true that the real or predicted signal for each channel during any given 5-second window is centered on zero? I should probably compare the results directly with normal Pearson correlation to see how much of a difference it makes. Also, the "bad-by-correlation" code uses the normal numpy `corrcoef` correlation function and matches up identically with MATLAB PREP, so it's odd that they'd use this unique correlation method here but not there.

Regardless, if we ultimately decide that Pearson (or some other method of correlation) would make a better default then we'll probably need to adjust the RANSAC bad-correlation threshold as well, since the current one is tuned for the current calculation method. That makes me think this question might be better suited for a future issue/PR.
sappelhoff: > That sounds like a good explanation. That said, it is generally going to be true that the real or predicted signal for each channel during any given 5-second window is centered on zero? 

I don't know, it was a shot in the dark.

> Regardless, if we ultimately decide that Pearson (or some other method of correlation) would make a better default then we'll probably need to adjust the RANSAC bad-correlation threshold as well, since the current one is tuned for the current calculation method. That makes me think this question might be better suited for a future issue/PR.

yes, good point :)

> Just tested directly, and it looks like the the results are pretty similar (although still somewhat divergent) between the correlation method MatPREP uses and the one we were using before (np.corrcoef).

that *does* suggest that mean(x) is indeed close to zero, and it makes little difference whether it's subtracted from x or not. 

I like `np.corrcoef` better though 🤔 
a-hurst: > I like `np.corrcoef` better though 🤔

Me too. I'm going to make the MatPREP correlation method a `matlab_strict` option in my new PR and leave `np.corrcoef` as the default.
a-hurst: > @a-hurst are you taking a break from this PR to come back later? Or are there any blocking items that you need help with?

Since I've broken this down into a bunch of separate PRs now, the only new thing this PR adds is support for is window-wise RANSAC, which is slower than the current method but has much lower RAM demands for longer recordings. To finish this up, I'm going to rebase on the current master branch and see if I can improve performance at all. I'll try and do the former today!
sappelhoff: sounds great! --> but even if it's a bit slower, it'd still be a win for users with low RAM machines :+1: 
a-hurst: > sounds great! --> but even if it's a bit slower, it'd still be a win for users with low RAM machines 👍

When I pull this together, should I make window-wise the new default then, with channel-wise still as an option for people with high-RAM workstations who want to process as quickly as possible?
sappelhoff: Sounds good!
codecov-commenter: # [Codecov](https://codecov.io/gh/sappelhoff/pyprep/pull/66?src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=Stefan+Appelhoff) Report
> Merging [#66](https://codecov.io/gh/sappelhoff/pyprep/pull/66?src=pr&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=Stefan+Appelhoff) (b9ce469) into [master](https://codecov.io/gh/sappelhoff/pyprep/commit/145a5a16d372c6317096871877e211fe01c51d54?el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=Stefan+Appelhoff) (145a5a1) will **decrease** coverage by `0.28%`.
> The diff coverage is `100.00%`.

> :exclamation: Current head b9ce469 differs from pull request most recent head 629950b. Consider uploading reports for the commit 629950b to get more accurate results
[![Impacted file tree graph](https://codecov.io/gh/sappelhoff/pyprep/pull/66/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/66?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      #66      +/-   ##
==========================================
- Coverage   97.31%   97.02%   -0.29%     
==========================================
  Files           7        7              
  Lines         634      606      -28     
==========================================
- Hits          617      588      -29     
- Misses         17       18       +1     
```


| [Impacted Files](https://codecov.io/gh/sappelhoff/pyprep/pull/66?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=Stefan+Appelhoff) | Coverage Δ | |
|---|---|---|
| [pyprep/find\_noisy\_channels.py](https://codecov.io/gh/sappelhoff/pyprep/pull/66/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=Stefan+Appelhoff#diff-cHlwcmVwL2ZpbmRfbm9pc3lfY2hhbm5lbHMucHk=) | `96.36% <ø> (-0.17%)` | :arrow_down: |
| [pyprep/ransac.py](https://codecov.io/gh/sappelhoff/pyprep/pull/66/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.14% <100.00%> (-0.13%)` | :arrow_down: |
| [pyprep/utils.py](https://codecov.io/gh/sappelhoff/pyprep/pull/66/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=Stefan+Appelhoff#diff-cHlwcmVwL3V0aWxzLnB5) | `98.55% <100.00%> (-0.52%)` | :arrow_down: |
| [pyprep/removeTrend.py](https://codecov.io/gh/sappelhoff/pyprep/pull/66/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=Stefan+Appelhoff#diff-cHlwcmVwL3JlbW92ZVRyZW5kLnB5) | `94.64% <0.00%> (-0.44%)` | :arrow_down: |
| [pyprep/reference.py](https://codecov.io/gh/sappelhoff/pyprep/pull/66/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=Stefan+Appelhoff#diff-cHlwcmVwL3JlZmVyZW5jZS5weQ==) | `96.39% <0.00%> (-0.16%)` | :arrow_down: |
| [pyprep/\_\_init\_\_.py](https://codecov.io/gh/sappelhoff/pyprep/pull/66/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=Stefan+Appelhoff#diff-cHlwcmVwL19faW5pdF9fLnB5) | `100.00% <0.00%> (ø)` | |
| [pyprep/prep\_pipeline.py](https://codecov.io/gh/sappelhoff/pyprep/pull/66/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=Stefan+Appelhoff#diff-cHlwcmVwL3ByZXBfcGlwZWxpbmUucHk=) | `100.00% <0.00%> (ø)` | |

------

[Continue to review full report at Codecov](https://codecov.io/gh/sappelhoff/pyprep/pull/66?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/66?src=pr&el=footer&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=Stefan+Appelhoff). Last update [145a5a1...629950b](https://codecov.io/gh/sappelhoff/pyprep/pull/66?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).

a-hurst: @sappelhoff Okay, so I think this is ready for review now: I haven't updated the `whats_new` yet and still need to finish exposing the public API for channel-wise RANSAC, but all the key parts are in place now and there are some things I wanted to get your thoughts on before finishing up.

The main thing I wanted to get feedback on is the repurposing of the `channel_wise` argument to `find_bad_by_ransac`: basically, I realized it would be confusing to have an argument for `window_wise` if that's going to be the default method, so I made `channel_wise` the argument to specify the channel-chunking method of RANSAC. I also added a new argument `max_chunk_size` to provide the same functionality as the old `channel_wise` argument (i.e. lets you run 1 channel per chunk) but with a bit more flexibility.

*However*, I realize that `channel_wise` is technically currently part of the public API and it's generally very bad practice to change the meaning of an argument between releases. My rationale here is that the current `channel_wise` functionality was never fully exposed (i.e., you couldn't select one-channel-at-a-time RANSAC from the `PrepPipeline` class, you could only use it if using `NoisyChannels` manually), and it was never formally introduced in the What's New, meaning it's unlikely to be used beyond the small circle of contributors here (@yjmantilla does this break anything for you?). Additionally, if I change the default value of `max_chunk_size` to 1 instead of none, `channel_wise` will behave exactly as it did before (though I think `None` is the more useful default).

Apart from that, my main changes are:

- Reworked channel-wise RANSAC to use the same pre-generated interpolation matrices as window-wise RANSAC, which speeds things up and reduces the amount of divergent/redundant code between methods
- Renamed the existing private channel-wise RANSAC functions to be consistent with the window-wise function names. Also rearranged their argument names/orders, variable names, and docstrings to match up more closely with their window-wise counterparts (for easier comparison/readability).
- Added a progress-printing function for window-wise RANSAC
- Overhauled the RANSAC unit tests considerably, including the addition of explicit comparison of the correlation matrices for window-wise/channel-wise for both regular and MATLAB-strict RANSAC to make sure window-wise and channel-wise produce identical results for the same seed.

Let me know what you think!
sappelhoff: >  I realize that channel_wise is technically currently part of the public API and it's generally very bad practice to change the meaning of an argument between releases.

generally yes, but with the immense changes that are coming in version `0.4.0`, it's totally fine to break a couple of things if afterwards everything is better. Strict semantic versioning would make us increment a "Major" version even (so 0.3.1 -> 1.0.0, instead of 0.3.1 -> 0.4.0), but a 1.0.0 release would be too much of a statement at this point, I think.

> and it was never formally introduced in the What's New, meaning it's unlikely to be used beyond the small circle of contributors here

another good argument why breaking would be fine here

> though I think None is the more useful default

+1 to None :+1: 


would be interested in your opinion as well here @yjmantilla 
a-hurst: @sappelhoff Just thinking of how to implement the public API for this now: do you think these kinds of settings should be part of the `prep_params` dict, or should they be separate arguments on their own?
sappelhoff: Probably separate parameters :thinking: and yes, the find_all_bads method should have these params as well.
yjmantilla: > @yjmantilla does this break anything for you?

It does not, and even if it would I would do a fork.
a-hurst: @sappelhoff I've now added a proper public API for `channel_wise` and `max_chunk_size` and updated the `whats_new.rst`. I think this is ready for final review!

I'm still open to adding chunking support for window-wise as you suggested, but I think that might be better to do in a future PR given the amount of change in this one as-is: from a cursory look at what's involved I'd have to refactor both channel-wise and window-wise RANSAC a fair bit to make that work (see: all the wrapping try/catching re: RAM-checking that's currently not used at all for window-wise after the initial check further up). While working on this PR I had some thoughts on how to clean that up too, but again I didn't want to overload this with too many peripheral changes.
a-hurst: > Cool! I have only one major comment for the docstr, the rest looks good to go!

Fixed in the latest commit, should be good to go now 🙂

Patch:
diff --git a/docs/whats_new.rst b/docs/whats_new.rst
index f28ca96..0b7cbce 100644
--- a/docs/whats_new.rst
+++ b/docs/whats_new.rst
@@ -41,6 +41,8 @@ Changelog
 - 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`)
+- Added a window-wise implementaion of RANSAC and made it the default method, reducing the typical RAM demands of robust re-referencing considerably, by `Austin Hurst`_ (:gh:`66`)
+- 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`)
 
 Bug
 ~~~
@@ -55,6 +57,7 @@ API
 - The permissible parameters for the following methods were removed and/or reordered: `ransac._ransac_correlations`, `ransac._run_ransac`, and `ransac._get_ransac_pred` methods, by `Yorguin Mantilla`_ (:gh:`51`)
 - The following methods have been moved to a new module named :mod:`~pyprep.ransac` and are now private: `NoisyChannels.ransac_correlations`, `NoisyChannels.run_ransac`, and `NoisyChannels.get_ransac_pred` methods, by `Yorguin Mantilla`_ (:gh:`51`)
 - The permissible parameters for the following methods were removed and/or reordered: `NoisyChannels.ransac_correlations`, `NoisyChannels.run_ransac`, and `NoisyChannels.get_ransac_pred` methods, by `Austin Hurst`_ and `Yorguin Mantilla`_ (:gh:`43`)
+- Changed the meaning of the argument `channel_wise` in :meth:`~pyprep.NoisyChannels.find_bad_by_ransac` to mean 'perform RANSAC across chunks of channels instead of window-wise', from its original meaning of 'perform channel-wise RANSAC one channel at a time', by `Austin Hurst`_ (:gh:`66`)
 
 
 .. _changes_0_3_1:
diff --git a/examples/run_ransac.py b/examples/run_ransac.py
index bfff7b0..93f2705 100644
--- a/examples/run_ransac.py
+++ b/examples/run_ransac.py
@@ -74,14 +74,15 @@ nd = NoisyChannels(raw)
 nd2 = NoisyChannels(raw)
 
 ###############################################################################
-# Find all bad channels and print a summary
+# Find all bad channels using channel-wise RANSAC and print a summary
 start_time = perf_counter()
-nd.find_bad_by_ransac()
+nd.find_bad_by_ransac(channel_wise=True)
 print("--- %s seconds ---" % (perf_counter() - start_time))
 
-# Repeat RANSAC in a channel wise manner. This is slower but needs less memory.
+# Repeat channel-wise RANSAC using a single channel at a time. This is slower
+# but needs less memory.
 start_time = perf_counter()
-nd2.find_bad_by_ransac(channel_wise=True)
+nd2.find_bad_by_ransac(channel_wise=True, max_chunk_size=1)
 print("--- %s seconds ---" % (perf_counter() - start_time))
 
 ###############################################################################
diff --git a/pyprep/find_noisy_channels.py b/pyprep/find_noisy_channels.py
index 03d02ef..6fbf9d9 100644
--- a/pyprep/find_noisy_channels.py
+++ b/pyprep/find_noisy_channels.py
@@ -144,22 +144,45 @@ class NoisyChannels:
             )
         return bads
 
-    def find_all_bads(self, ransac=True):
+    def find_all_bads(self, ransac=True, channel_wise=False, max_chunk_size=None):
         """Call all the functions to detect bad channels.
 
         This function calls all the bad-channel detecting functions.
 
         Parameters
         ----------
-        ransac : bool
-            To detect channels by ransac or not.
+        ransac : bool, optional
+            Whether RANSAC should be used for bad channel detection, in addition
+            to the other methods. RANSAC can detect bad channels that other
+            methods are unable to catch, but also slows down noisy channel
+            detection considerably. Defaults to ``True``.
+        channel_wise : bool, optional
+            Whether RANSAC should predict signals for chunks of channels over the
+            entire signal length ("channel-wise RANSAC", see `max_chunk_size`
+            parameter). If ``False``, RANSAC will instead predict signals for all
+            channels at once but over a number of smaller time windows instead of
+            over the entire signal length ("window-wise RANSAC"). Channel-wise
+            RANSAC generally has higher RAM demands than window-wise RANSAC
+            (especially if `max_chunk_size` is ``None``), but can be faster on
+            systems with lots of RAM to spare. Has no effect if not using RANSAC.
+            Defaults to ``False``.
+        max_chunk_size : {int, None}, optional
+            The maximum number of channels to predict at once during
+            channel-wise RANSAC. If ``None``, RANSAC will use the largest chunk
+            size that will fit into the available RAM, which may slow down
+            other programs on the host system. If using window-wise RANSAC
+            (the default) or not using RANSAC at all, this parameter has no
+            effect. Defaults to ``None``.
 
         """
         self.find_bad_by_nan_flat()
         self.find_bad_by_deviation()
         self.find_bad_by_SNR()
         if ransac:
-            self.find_bad_by_ransac()
+            self.find_bad_by_ransac(
+                channel_wise=channel_wise,
+                max_chunk_size=max_chunk_size
+            )
 
     def find_bad_by_nan_flat(self):
         """Detect channels that appear flat or have NaN values."""
@@ -409,6 +432,7 @@ class NoisyChannels:
         fraction_bad=0.4,
         corr_window_secs=5.0,
         channel_wise=False,
+        max_chunk_size=None,
     ):
         """Detect channels that are not predicted well by other channels.
 
@@ -447,10 +471,20 @@ class NoisyChannels:
             The duration (in seconds) of each RANSAC correlation window. Defaults
             to 5 seconds.
         channel_wise : bool, optional
-            Whether RANSAC should be performed one channel at a time (lower RAM
-            demands) or in chunks of as many channels as can fit into the
-            currently available RAM (faster). Defaults to ``False`` (i.e., using
-            the faster method).
+            Whether RANSAC should predict signals for chunks of channels over the
+            entire signal length ("channel-wise RANSAC", see `max_chunk_size`
+            parameter). If ``False``, RANSAC will instead predict signals for all
+            channels at once but over a number of smaller time windows instead of
+            over the entire signal length ("window-wise RANSAC"). Channel-wise
+            RANSAC generally has higher RAM demands than window-wise RANSAC
+            (especially if `max_chunk_size` is ``None``), but can be faster on
+            systems with lots of RAM to spare. Defaults to ``False``.
+        max_chunk_size : {int, None}, optional
+            The maximum number of channels to predict at once during
+            channel-wise RANSAC. If ``None``, RANSAC will use the largest chunk
+            size that will fit into the available RAM, which may slow down
+            other programs on the host system. If using window-wise RANSAC
+            (the default), this parameter has no effect. Defaults to ``None``.
 
         References
         ----------
@@ -479,6 +513,7 @@ class NoisyChannels:
             fraction_bad,
             corr_window_secs,
             channel_wise,
+            max_chunk_size,
             self.random_state,
             self.matlab_strict,
         )
diff --git a/pyprep/prep_pipeline.py b/pyprep/prep_pipeline.py
index 2c4a38b..6781b09 100644
--- a/pyprep/prep_pipeline.py
+++ b/pyprep/prep_pipeline.py
@@ -38,6 +38,22 @@ class PrepPipeline:
     ransac : bool, optional
         Whether or not to use RANSAC for noisy channel detection in addition to
         the other methods in :class:`~pyprep.NoisyChannels`. Defaults to True.
+    channel_wise : bool, optional
+        Whether RANSAC should predict signals for chunks of channels over the
+        entire signal length ("channel-wise RANSAC", see `max_chunk_size`
+        parameter). If ``False``, RANSAC will instead predict signals for all
+        channels at once but over a number of smaller time windows instead of
+        over the entire signal length ("window-wise RANSAC"). Channel-wise
+        RANSAC generally has higher RAM demands than window-wise RANSAC
+        (especially if `max_chunk_size` is ``None``), but can be faster on
+        systems with lots of RAM to spare. Has no effect if not using RANSAC.
+        Defaults to ``False``.
+    max_chunk_size : {int, None}, optional
+        The maximum number of channels to predict at once during channel-wise
+        RANSAC. If ``None``, RANSAC will use the largest chunk size that will
+        fit into the available RAM, which may slow down other programs on the
+        host system. If using window-wise RANSAC (the default) or not using
+        RANSAC at all, this parameter has no effect. Defaults to ``None``.
     random_state : {int, None, np.random.RandomState}, optional
         The random seed at which to initialize the class. If random_state is
         an int, it will be used as a seed for RandomState.
@@ -99,6 +115,8 @@ class PrepPipeline:
         prep_params,
         montage,
         ransac=True,
+        channel_wise=False,
+        max_chunk_size=None,
         random_state=None,
         filter_kwargs=None,
         matlab_strict=False,
@@ -133,7 +151,11 @@ class PrepPipeline:
         if self.prep_params["reref_chs"] == "eeg":
             self.prep_params["reref_chs"] = self.ch_names_eeg
         self.sfreq = self.raw_eeg.info["sfreq"]
-        self.ransac = ransac
+        self.ransac_settings = {
+            'ransac': ransac,
+            'channel_wise': channel_wise,
+            'max_chunk_size': max_chunk_size
+        }
         self.random_state = check_random_state(random_state)
         self.filter_kwargs = filter_kwargs
         self.matlab_strict = matlab_strict
@@ -189,9 +211,9 @@ class PrepPipeline:
         reference = Reference(
             self.raw_eeg,
             self.prep_params,
-            ransac=self.ransac,
             random_state=self.random_state,
-            matlab_strict=self.matlab_strict
+            matlab_strict=self.matlab_strict,
+            **self.ransac_settings
         )
         reference.perform_reference()
         self.raw_eeg = reference.raw
diff --git a/pyprep/ransac.py b/pyprep/ransac.py
index 5bc284f..3e1b86a 100644
--- a/pyprep/ransac.py
+++ b/pyprep/ransac.py
@@ -5,7 +5,8 @@ from mne.channels.interpolation import _make_interpolation_matrix
 from mne.utils import check_random_state
 
 from pyprep.utils import (
-    split_list, verify_free_ram, _get_random_subset, _mat_round, _correlate_arrays
+    _get_random_subset, _mat_round, _correlate_arrays, print_progress,
+    split_list, verify_free_ram
 )
 
 
@@ -21,6 +22,7 @@ def find_bad_by_ransac(
     fraction_bad=0.4,
     corr_window_secs=5.0,
     channel_wise=False,
+    max_chunk_size=None,
     random_state=None,
     matlab_strict=False,
 ):
@@ -70,10 +72,20 @@ def find_bad_by_ransac(
         The duration (in seconds) of each RANSAC correlation window. Defaults to
         5 seconds.
     channel_wise : bool, optional
-        Whether RANSAC should be performed one channel at a time (lower RAM
-        demands) or in chunks of as many channels as can fit into the currently
-        available RAM (faster). Defaults to ``False`` (i.e., using the faster
-        method).
+        Whether RANSAC should predict signals for chunks of channels over the
+        entire signal length ("channel-wise RANSAC", see `max_chunk_size`
+        parameter). If ``False``, RANSAC will instead predict signals for all
+        channels at once but over a number of smaller time windows instead of
+        over the entire signal length ("window-wise RANSAC"). Channel-wise
+        RANSAC generally has higher RAM demands than window-wise RANSAC
+        (especially if `max_chunk_size` is ``None``), but can be faster on
+        systems with lots of RAM to spare. Defaults to ``False``.
+    max_chunk_size : {int, None}, optional
+        The maximum number of channels to predict at once during channel-wise
+        RANSAC. If ``None``, RANSAC will use the largest chunk size that will
+        fit into the available RAM, which may slow down other programs on the
+        host system. If using window-wise RANSAC (the default), this parameter
+        has no effect. Defaults to ``None``.
     random_state : {int, None, np.random.RandomState}, optional
         The random seed with which to generate random samples of channels during
         RANSAC. If random_state is an int, it will be used as a seed for RandomState.
@@ -131,7 +143,11 @@ def find_bad_by_ransac(
 
     # Before running, make sure we have enough memory when using the
     # smallest possible chunk size
-    verify_free_ram(data, n_samples, 1)
+    if channel_wise:
+        verify_free_ram(data, n_samples, 1)
+    else:
+        window_size = int(sample_rate * corr_window_secs)
+        verify_free_ram(data[:, :window_size], n_samples, n_chans_good)
 
     # Generate random channel picks for each RANSAC sample
     random_ch_picks = []
@@ -142,25 +158,34 @@ def find_bad_by_ransac(
         picks = _get_random_subset(good_chans, n_pred_chns, rng)
         random_ch_picks.append(picks)
 
-    # Correlation windows setup
+    # Generate interpolation matrix for each RANSAC sample
+    interp_mats = _make_interpolation_matrices(random_ch_picks, chn_pos_good)
+
+    # Calculate the size (in frames) and count of correlation windows
     correlation_frames = corr_window_secs * sample_rate
-    correlation_window = np.arange(correlation_frames)
-    n = correlation_window.shape[0]
     signal_frames = data.shape[1]
     correlation_offsets = np.arange(
         0, (signal_frames - correlation_frames), correlation_frames
     )
-    w_correlation = correlation_offsets.shape[0]
+    win_size = int(correlation_frames)
+    win_count = correlation_offsets.shape[0]
 
-    # Preallocate
+    # Preallocate RANSAC correlation matrix
     n_chans_complete = len(complete_chn_labs)
-    channel_correlations = np.ones((w_correlation, n_chans_complete))
+    channel_correlations = np.ones((win_count, n_chans_complete))
     # Notice self.EEGData.shape[0] = self.n_chans_new
     # Is now data.shape[0] = n_chans_complete
     # They came from the same drop of channels
 
     print("Executing RANSAC\nThis may take a while, so be patient...")
 
+    # If enabled, run window-wise RANSAC
+    if not channel_wise:
+        # Get correlations between actual vs predicted signals for each RANSAC window
+        channel_correlations[:, good_idx] = _ransac_by_window(
+            data[good_idx, :], interp_mats, win_size, win_count, matlab_strict
+        )
+
     # Calculate smallest chunk size for each possible chunk count
     chunk_sizes = []
     chunk_count = 0
@@ -168,28 +193,28 @@ def find_bad_by_ransac(
         n_chunks = int(np.ceil(n_chans_good / i))
         if n_chunks != chunk_count:
             chunk_count = n_chunks
-            chunk_sizes.append(i)
+            if not max_chunk_size or i <= max_chunk_size:
+                chunk_sizes.append(i)
 
     chunk_size = chunk_sizes.pop()
     mem_error = True
     job = list(range(n_chans_good))
 
-    if channel_wise:
-        chunk_size = 1
-    while mem_error:
+    # If not using window-wise RANSAC, do channel-wise RANSAC
+    while mem_error and channel_wise:
         try:
             channel_chunks = split_list(job, chunk_size)
             total_chunks = len(channel_chunks)
             current = 1
             for chunk in channel_chunks:
-                channel_correlations[:, good_idx[chunk]] = _ransac_correlations(
+                interp_mats_for_chunk = [mat[chunk, :] for mat in interp_mats]
+                channel_correlations[:, good_idx[chunk]] = _ransac_by_channel(
+                    data[good_idx, :],
+                    interp_mats_for_chunk,
+                    win_size,
+                    win_count,
                     chunk,
                     random_ch_picks,
-                    chn_pos_good,
-                    data[good_idx, :],
-                    n_samples,
-                    n,
-                    w_correlation,
                     matlab_strict,
                 )
                 if chunk == channel_chunks[0]:
@@ -226,171 +251,277 @@ def find_bad_by_ransac(
     return bad_by_ransac, channel_correlations
 
 
-def _ransac_correlations(
+def _make_interpolation_matrices(random_ch_picks, chn_pos_good):
+    """Create an interpolation matrix for each RANSAC sample of channels.
+
+    This function takes the spatial coordinates of random subsets of currently-good
+    channels and uses them to predict what the signal will be at the spatial
+    coordinates of all other currently-good channels. The results of this process are
+    returned as matrices that can be multiplied with EEG data to generate predicted
+    signals.
+
+    Parameters
+    ----------
+    random_ch_picks : list of list of int
+        A list containing multiple random subsets of currently-good channels.
+    chn_pos_good : np.ndarray
+        3-D spatial coordinates of all currently-good channels.
+
+    Returns
+    -------
+    interpolation_mats : list of np.ndarray
+        A list of interpolation matrices, one for each random subset of channels.
+        Each matrix has the shape `[num_good_channels, num_good_channels]`, with the
+        number of good channels being inferred from the size of `ch_pos_good`.
+
+    Notes
+    -----
+    This function currently makes use of a private MNE function,
+    ``mne.channels.interpolation._make_interpolation_matrix``, to generate matrices.
+
+    """
+    n_chans_good = chn_pos_good.shape[0]
+
+    interpolation_mats = []
+    for sample in random_ch_picks:
+        mat = np.zeros((n_chans_good, n_chans_good))
+        subset_pos = chn_pos_good[sample, :]
+        mat[:, sample] = _make_interpolation_matrix(subset_pos, chn_pos_good)
+        interpolation_mats.append(mat)
+
+    return interpolation_mats
+
+
+def _ransac_by_window(data, interpolation_mats, win_size, win_count, matlab_strict):
+    """Calculate correlations of channels with their RANSAC-predicted values.
+
+    This function calculates RANSAC correlations for each RANSAC window
+    individually, requiring RAM equivalent to [channels * sample rate * seconds
+    per RANSAC window] to run. Generally, this method will use less RAM than
+    :func:`_ransac_by_channel`, with the exception of short recordings with high
+    electrode counts.
+
+    Parameters
+    ----------
+    data : np.ndarray
+        A 2-D array containing the EEG signals from all currently-good channels.
+    interpolation_mats : list of np.ndarray
+        A list of interpolation matrices, one for each RANSAC sample of channels.
+    win_size : int
+        Number of frames/samples of EEG data in each RANSAC correlation window.
+    win_count: int
+        Number of RANSAC correlation windows.
+    matlab_strict : bool
+        Whether or not RANSAC should strictly follow MATLAB PREP's internal
+        math, ignoring any improvements made in PyPREP over the original code.
+
+    Returns
+    -------
+    correlations : np.ndarray
+        Correlations of the given channels to their predicted values within each
+        RANSAC window.
+
+    """
+    ch_count = data.shape[0]
+    correlations = np.ones((win_count, ch_count))
+
+    for window in range(win_count):
+
+        # Print RANSAC progress in 10% increments
+        print_progress(window + 1, win_count, every=0.1)
+
+        # Get the current window of EEG data
+        start = window * win_size
+        end = (window + 1) * win_size
+        actual = data[:, start:end]
+
+        # Get the median RANSAC-predicted signal for each channel
+        predicted = _predict_median_signals(actual, interpolation_mats, matlab_strict)
+
+        # Calculate the actual vs predicted signal correlation for each channel
+        correlations[window, :] = _correlate_arrays(actual, predicted, matlab_strict)
+
+    return correlations
+
+
+def _predict_median_signals(window, interpolation_mats, matlab_strict=False):
+    """Calculate the median RANSAC-predicted signal for a given window of data.
+
+    Parameters
+    ----------
+    window : np.ndarray
+        A 2-D window of EEG data with the shape `[channels, samples]`.
+    interpolation_mats : list of np.ndarray
+        A set of channel interpolation matrices, one for each RANSAC sample of
+        channels.
+    matlab_strict : bool
+        Whether or not RANSAC should strictly follow MATLAB PREP's internal
+        math, ignoring any improvements made in PyPREP over the original code.
+
+    Returns
+    -------
+    predicted : np.ndarray
+        The median RANSAC-predicted EEG signal for the given window of data.
+
+    Notes
+    -----
+    In MATLAB PREP, the median signal is calculated by sorting the different
+    predictions for each EEG sample/channel from low to high and then taking the value
+    at the middle index (as calculated by ``int(n_ransac_samples / 2.0)``) for each.
+    Because this logic only returns the correct result for odd numbers of samples, the
+    current function will instead return the true median signal across predictions
+    unless strict MATLAB equivalence is requested.
+
+    """
+    ransac_samples = len(interpolation_mats)
+    merged_mats = np.concatenate(interpolation_mats, axis=0)
+
+    predictions_per_sample = np.reshape(
+        np.matmul(merged_mats, window),
+        (ransac_samples, window.shape[0], window.shape[1])
+    )
+
+    if matlab_strict:
+        # Match MATLAB's rounding logic (.5 always rounded up)
+        median_idx = int(_mat_round(ransac_samples / 2.0) - 1)
+        predictions_per_sample.sort(axis=0)
+        return predictions_per_sample[median_idx, :, :]
+    else:
+        return np.median(predictions_per_sample, axis=0)
+
+
+def _ransac_by_channel(
+    data,
+    interpolation_mats,
+    win_size,
+    win_count,
     chans_to_predict,
     random_ch_picks,
-    chn_pos_good,
-    data,
-    n_samples,
-    n,
-    w_correlation,
     matlab_strict,
 ):
-    """Get correlations of channels to their RANSAC-predicted values.
+    """Calculate correlations of channels with their RANSAC-predicted values.
+
+    This function calculates RANSAC correlations on one (or more) full channels
+    at once, requiring RAM equivalent to [channels per chunk * sample rate *
+    length of recording in seconds] to run. Generally, this method will use
+    more RAM than :func:`_ransac_by_window`, but may be faster for systems with
+    large amounts of RAM.
 
     Parameters
     ----------
-    chans_to_predict : list of int
-        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 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
     data : np.ndarray
-        2-D EEG data
-    n_samples : int
-        Number of samples used for computation of RANSAC.
-    n : int
-        Number of frames/samples of each window.
-    w_correlation: int
-        Number of windows.
+        A 2-D array containing the EEG signals from all currently-good channels.
+    interpolation_mats : list of np.ndarray
+        A set of channel interpolation matrices, one for each RANSAC sample of
+        channels.
+    win_size : int
+        Number of frames/samples of EEG data in each RANSAC correlation window.
+    win_count: int
+        Number of RANSAC correlation windows.
+    chans_to_predict : list of int
+        Indices of the channels to predict (as they appear in `data`) within the
+        current chunk.
+    random_ch_picks : list of list of int
+        A list containing multiple random subsets of currently-good channels.
     matlab_strict : bool
         Whether or not RANSAC should strictly follow MATLAB PREP's internal
         math, ignoring any improvements made in PyPREP over the original code.
 
     Returns
     -------
-    channel_correlations : np.ndarray
-        correlations of the given channels to their RANSAC-predicted values.
+    correlations : np.ndarray
+        Correlations of the given channels to their predicted values within each
+        RANSAC window.
 
     """
-    # Preallocate
-    channel_correlations = np.ones((w_correlation, len(chans_to_predict)))
+    # Preallocate RANSAC correlation matrix for current chunk
+    chunk_size = len(chans_to_predict)
+    correlations = np.ones((win_count, chunk_size))
 
-    # Make the ransac predictions
-    ransac_eeg = _run_ransac(
-        n_samples=n_samples,
-        random_ch_picks=random_ch_picks,
-        chn_pos=chn_pos_good[chans_to_predict, :],
-        chn_pos_good=chn_pos_good,
+    # Get median RANSAC predictions for each channel in the current chunk
+    predicted_chans = _predict_median_signals_channelwise(
         data=data,
+        interpolation_mats=interpolation_mats,
+        random_ch_picks=random_ch_picks,
+        chunk_size=len(chans_to_predict),
         matlab_strict=matlab_strict,
     )
 
     # Correlate ransac prediction and eeg data
 
     # For the actual data
-    data_window = data[chans_to_predict, : n * w_correlation]
-    data_window = data_window.reshape(len(chans_to_predict), w_correlation, n)
+    data_window = data[chans_to_predict, : win_size * win_count]
+    data_window = data_window.reshape(chunk_size, win_count, win_size)
     data_window = data_window.swapaxes(1, 0)
 
     # For the ransac predicted eeg
-    pred_window = ransac_eeg[: len(chans_to_predict), : n * w_correlation]
-    pred_window = pred_window.reshape(len(chans_to_predict), w_correlation, n)
+    pred_window = predicted_chans[: chunk_size, : win_size * win_count]
+    pred_window = pred_window.reshape(chunk_size, win_count, win_size)
     pred_window = pred_window.swapaxes(1, 0)
 
     # Perform correlations
-    for k in range(w_correlation):
+    for k in range(win_count):
         data_portion = data_window[k, :, :]
         pred_portion = pred_window[k, :, :]
         R = _correlate_arrays(data_portion, pred_portion, matlab_strict)
-        channel_correlations[k, :] = R
+        correlations[k, :] = R
 
-    return channel_correlations
+    return correlations
 
 
-def _run_ransac(
-    n_samples,
-    random_ch_picks,
-    chn_pos,
-    chn_pos_good,
+def _predict_median_signals_channelwise(
     data,
+    interpolation_mats,
+    random_ch_picks,
+    chunk_size,
     matlab_strict,
 ):
-    """Detect noisy channels apart from the ones described previously.
-
-    It creates a random subset of the so-far good channels
-    and predicts the values of the channels not in the subset.
+    """Calculate the median RANSAC-predicted signal for a given chunk of channels.
 
     Parameters
     ----------
-    n_samples : int
-        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 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
     data : np.ndarray
-        2-D EEG data
+        A 2-D array containing the EEG signals from all currently-good channels.
+    interpolation_mats : list of np.ndarray
+        A set of channel interpolation matrices, one for each RANSAC sample of
+        channels.
+    random_ch_picks : list of list of int
+        A list containing multiple random subsets of currently-good channels.
+    chunk_size : int
+        The number of channels to predict in the current chunk.
     matlab_strict : bool
         Whether or not RANSAC should strictly follow MATLAB PREP's internal
         math, ignoring any improvements made in PyPREP over the original code.
 
     Returns
     -------
-    ransac_eeg : np.ndarray
-        The EEG data predicted by RANSAC
+    predicted_chans : np.ndarray
+        The median RANSAC-predicted EEG signals for the given chunk of channels.
 
     """
     # n_chns, n_timepts = data.shape
     # 2 next lines should be equivalent but support single channel processing
+    ransac_samples = len(interpolation_mats)
     n_timepts = data.shape[1]
-    n_chns = chn_pos.shape[0]
 
     # Before running, make sure we have enough memory
-    verify_free_ram(data, n_samples, n_chns)
+    verify_free_ram(data, ransac_samples, chunk_size)
 
     # Memory seems to be fine ...
     # Make the predictions
-    eeg_predictions = np.zeros((n_chns, n_timepts, n_samples))
-    for sample in range(n_samples):
-        # Get the random channel selection for the current sample
+    eeg_predictions = np.zeros((chunk_size, n_timepts, ransac_samples))
+    for sample in range(ransac_samples):
+        # Get the random channels & interpolation matrix for the current sample
         reconstr_idx = random_ch_picks[sample]
-        eeg_predictions[..., sample] = _get_ransac_pred(
-            chn_pos, chn_pos_good, reconstr_idx, data
-        )
+        interp_mat = interpolation_mats[sample][:, reconstr_idx]
+        # Predict the EEG signals for the current RANSAC sample / channel chunk
+        eeg_predictions[..., sample] = np.matmul(interp_mat, data[reconstr_idx, :])
 
     # Form median from all predictions
     if matlab_strict:
         # Match MATLAB's rounding logic (.5 always rounded up)
-        median_idx = int(_mat_round(n_samples / 2.0) - 1)
+        median_idx = int(_mat_round(ransac_samples / 2.0) - 1)
         eeg_predictions.sort(axis=-1)
-        ransac_eeg = eeg_predictions[:, :, median_idx]
+        return 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, reconstr_idx, data):
-    """Perform RANSAC prediction.
-
-    Parameters
-    ----------
-    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
-    reconstr_idx : array_like
-        indexes of the channels in chn_pos_good to use for reconstruction
-    data : np.ndarray
-        2-D EEG data
-
-    Returns
-    -------
-    ransac_pred : np.ndarray
-        Single RANSAC prediction
-
-    """
-    # Get positions
-    reconstr_pos = chn_pos_good[reconstr_idx, :]
-
-    # Interpolate
-    interpol_mat = _make_interpolation_matrix(reconstr_pos, chn_pos)
-    ransac_pred = np.matmul(interpol_mat, data[reconstr_idx, :])
-
-    return ransac_pred
+        return np.median(eeg_predictions, axis=-1, overwrite_input=True)
diff --git a/pyprep/reference.py b/pyprep/reference.py
index 0cdddf3..5910867 100644
--- a/pyprep/reference.py
+++ b/pyprep/reference.py
@@ -32,6 +32,22 @@ class Reference:
     ransac : bool, optional
         Whether or not to use RANSAC for noisy channel detection in addition to
         the other methods in :class:`~pyprep.NoisyChannels`. Defaults to True.
+    channel_wise : bool, optional
+        Whether RANSAC should predict signals for chunks of channels over the
+        entire signal length ("channel-wise RANSAC", see `max_chunk_size`
+        parameter). If ``False``, RANSAC will instead predict signals for all
+        channels at once but over a number of smaller time windows instead of
+        over the entire signal length ("window-wise RANSAC"). Channel-wise
+        RANSAC generally has higher RAM demands than window-wise RANSAC
+        (especially if `max_chunk_size` is ``None``), but can be faster on
+        systems with lots of RAM to spare. Has no effect if not using RANSAC.
+        Defaults to ``False``.
+    max_chunk_size : {int, None}, optional
+        The maximum number of channels to predict at once during channel-wise
+        RANSAC. If ``None``, RANSAC will use the largest chunk size that will
+        fit into the available RAM, which may slow down other programs on the
+        host system. If using window-wise RANSAC (the default) or not using
+        RANSAC at all, this parameter has no effect. Defaults to ``None``.
     random_state : {int, None, np.random.RandomState}, optional
         The random seed at which to initialize the class. If random_state is
         an int, it will be used as a seed for RandomState.
@@ -51,7 +67,14 @@ class Reference:
     """
 
     def __init__(
-        self, raw, params, ransac=True, random_state=None, matlab_strict=False
+        self,
+        raw,
+        params,
+        ransac=True,
+        channel_wise=False,
+        max_chunk_size=None,
+        random_state=None,
+        matlab_strict=False
     ):
         """Initialize the class."""
         self.raw = raw.copy()
@@ -62,7 +85,11 @@ class Reference:
         self.reference_channels = params["ref_chs"]
         self.rereferenced_channels = params["reref_chs"]
         self.sfreq = self.raw.info["sfreq"]
-        self.ransac = ransac
+        self.ransac_settings = {
+            'ransac': ransac,
+            'channel_wise': channel_wise,
+            'max_chunk_size': max_chunk_size
+        }
         self.random_state = check_random_state(random_state)
         self._extra_info = {}
         self.matlab_strict = matlab_strict
@@ -103,7 +130,7 @@ class Reference:
         noisy_detector = NoisyChannels(
             self.raw, random_state=self.random_state, matlab_strict=self.matlab_strict
         )
-        noisy_detector.find_all_bads(ransac=self.ransac)
+        noisy_detector.find_all_bads(**self.ransac_settings)
 
         # Record Noisy channels and EEG before interpolation
         self.bad_before_interpolation = noisy_detector.get_bads(verbose=True)
@@ -141,7 +168,7 @@ class Reference:
         noisy_detector = NoisyChannels(
             self.raw, random_state=self.random_state, matlab_strict=self.matlab_strict
         )
-        noisy_detector.find_all_bads(ransac=self.ransac)
+        noisy_detector.find_all_bads(**self.ransac_settings)
         self.still_noisy_channels = noisy_detector.get_bads()
         self.raw.info["bads"] = self.still_noisy_channels
         self.noisy_channels_after_interpolation = {
@@ -186,7 +213,7 @@ class Reference:
             random_state=self.random_state,
             matlab_strict=self.matlab_strict
         )
-        noisy_detector.find_all_bads(ransac=self.ransac)
+        noisy_detector.find_all_bads(**self.ransac_settings)
         self.noisy_channels_original = {
             "bad_by_nan": noisy_detector.bad_by_nan,
             "bad_by_flat": noisy_detector.bad_by_flat,
@@ -243,7 +270,7 @@ class Reference:
                 matlab_strict=self.matlab_strict
             )
             # Detrend applied at the beginning of the function.
-            noisy_detector.find_all_bads(ransac=self.ransac)
+            noisy_detector.find_all_bads(**self.ransac_settings)
             self.noisy_channels["bad_by_nan"] = _union(
                 self.noisy_channels["bad_by_nan"], noisy_detector.bad_by_nan
             )
diff --git a/pyprep/utils.py b/pyprep/utils.py
index 51c9518..47aa8e9 100644
--- a/pyprep/utils.py
+++ b/pyprep/utils.py
@@ -362,6 +362,44 @@ def split_list(mylist, chunk_size):
     ]
 
 
+def print_progress(current, end, start=None, stepsize=1, every=0.1):
+    """Print the current progress in a loop.
+
+    Parameters
+    ----------
+    current: {int, float}
+        The index or numeric value of the current position in the loop.
+    end: {int, float}
+        The final index or numeric value in the loop.
+    start: {int, float, None}, optional
+        The first index or numeric value in the loop. If ``None``, the start
+        index will assumed to be `stepsize` (i.e., 3 if `stepsize` is 3).
+        Defaults to ``None``.
+    stepsize: {int, float}, optional
+        The fixed amount by which `current` increases every iteration of the
+        loop. Defaults to ``1``.
+    every: float, optional
+        The frequency with which to print progress updates during the loop,
+        as a proportion between 0 and 1, exclusive. Defaults to ``0.1``, which
+        prints a progress update after every 10%.
+
+    """
+    start = stepsize if not start else start
+    end = end - start + 1
+    current = current - start + 1
+
+    if current == 1:
+        print("Progress:", end=" ", flush=True)
+    elif current == end:
+        print("100%")
+    elif current > 0:
+        progress = float(current) / end
+        last = float(current - stepsize) / end
+        if int(progress / every) > int(last / every):
+            pct = int(progress / every) * every * 100
+            print("{0}%...".format(int(pct)), end=" ", flush=True)
+
+
 def make_random_mne_object(
     ch_names,
     ch_types,


Test patch:
diff --git a/tests/test_find_noisy_channels.py b/tests/test_find_noisy_channels.py
index 4e79558..c54d4ad 100644
--- a/tests/test_find_noisy_channels.py
+++ b/tests/test_find_noisy_channels.py
@@ -3,6 +3,7 @@ import numpy as np
 import pytest
 
 from pyprep.find_noisy_channels import NoisyChannels
+from pyprep.removeTrend import removeTrend
 
 
 @pytest.mark.usefixtures("raw", "montage")
@@ -111,48 +112,68 @@ def test_findnoisychannels(raw, montage):
     nd.find_bad_by_SNR()
     assert rand_chn_lab in nd.bad_by_SNR
 
-    # Test for finding bad channels by RANSAC
-    raw_tmp = raw.copy()
-    # Ransac identifies channels that go bad together and are highly correlated.
-    # Inserting highly correlated signal in channels 0 through 3 at 30 Hz
-    raw_tmp._data[0:6, :] = np.cos(2 * np.pi * raw.times * 30) * 1e-6
-    nd = NoisyChannels(raw_tmp, random_state=rng)
-    nd.find_bad_by_ransac()
-    bads = nd.bad_by_ransac
-    assert bads == raw_tmp.ch_names[0:6]
 
-    # Test for finding bad channels by matlab_strict RANSAC
-    raw_tmp = raw.copy()
-    # Ransac identifies channels that go bad together and are highly correlated.
-    # Inserting highly correlated signal in channels 0 through 3 at 30 Hz
-    raw_tmp._data[0:6, :] = np.cos(2 * np.pi * raw.times * 30) * 1e-6
-    nd = NoisyChannels(raw_tmp, random_state=rng, matlab_strict=True)
-    nd.find_bad_by_ransac()
-    bads = nd.bad_by_ransac
-    assert bads == raw_tmp.ch_names[0:6]
+@pytest.mark.usefixtures("raw", "montage")
+def test_find_bad_by_ransac(raw, montage):
+    """Test the RANSAC component of NoisyChannels."""
+    # Set a fixed random seed and a montage for the tests
+    rng = 435656
+    raw.set_montage(montage)
 
-    # Test for finding bad channels by channel-wise RANSAC
-    raw_tmp = raw.copy()
-    # Ransac identifies channels that go bad together and are highly correlated.
+    # RANSAC identifies channels that go bad together and are highly correlated.
     # Inserting highly correlated signal in channels 0 through 3 at 30 Hz
-    raw_tmp._data[0:6, :] = np.cos(2 * np.pi * raw.times * 30) * 1e-6
-    nd = NoisyChannels(raw_tmp, random_state=rng)
-    nd.find_bad_by_ransac(channel_wise=True)
-    bads = nd.bad_by_ransac
-    assert bads == raw_tmp.ch_names[0:6]
-
-    # Test not-enough-memory and n_samples type exceptions
     raw_tmp = raw.copy()
     raw_tmp._data[0:6, :] = np.cos(2 * np.pi * raw.times * 30) * 1e-6
-    nd = NoisyChannels(raw_tmp, random_state=rng)
+
+    # Pre-detrend data to save time during NoisyChannels initialization
+    raw_tmp._data = removeTrend(raw_tmp.get_data(), raw.info["sfreq"])
+
+    # Run different variations of RANSAC on the same data
+    test_matrix = {
+        # List items represent [matlab_strict, channel_wise, max_chunk_size]
+        'by_window': [False, False, None],
+        'by_channel': [False, True, None],
+        'by_channel_maxchunk': [False, True, 2],
+        'by_window_strict': [True, False, None],
+        'by_channel_strict': [True, True, None]
+    }
+    bads = {}
+    corr = {}
+    for name, args in test_matrix.items():
+        nd = NoisyChannels(
+            raw_tmp, do_detrend=False, random_state=rng, matlab_strict=args[0]
+        )
+        nd.find_bad_by_ransac(channel_wise=args[1], max_chunk_size=args[2])
+        # Save bad channels and RANSAC correlation matrix for later comparison
+        bads[name] = nd.bad_by_ransac
+        corr[name] = nd._extra_info['bad_by_ransac']['ransac_correlations']
+
+    # Test whether all methods detected bad channels properly
+    assert bads['by_window'] == raw_tmp.ch_names[0:6]
+    assert bads['by_channel'] == raw_tmp.ch_names[0:6]
+    assert bads['by_channel_maxchunk'] == raw_tmp.ch_names[0:6]
+    assert bads['by_window_strict'] == raw_tmp.ch_names[0:6]
+    assert bads['by_channel_strict'] == raw_tmp.ch_names[0:6]
+
+    # Make sure non-strict correlation matrices all match
+    assert np.allclose(corr['by_window'], corr['by_channel'])
+    assert np.allclose(corr['by_window'], corr['by_channel_maxchunk'])
+
+    # Make sure MATLAB-strict correlation matrices match
+    assert np.allclose(corr['by_window_strict'], corr['by_channel_strict'])
+
+    # Make sure strict and non-strict matrices differ
+    assert not np.allclose(corr['by_window'], corr['by_window_strict'])
 
     # Set n_samples very very high to trigger a memory error
     n_samples = int(1e100)
+    nd = NoisyChannels(raw_tmp, do_detrend=False, random_state=rng)
     with pytest.raises(MemoryError):
         nd.find_bad_by_ransac(n_samples=n_samples)
 
     # Set n_samples to a float to trigger a type error
     n_samples = 35.5
+    nd = NoisyChannels(raw_tmp, do_detrend=False, random_state=rng)
     with pytest.raises(TypeError):
         nd.find_bad_by_ransac(n_samples=n_samples)
 
diff --git a/tests/test_utils.py b/tests/test_utils.py
index 8cf3547..624fdfd 100644
--- a/tests/test_utils.py
+++ b/tests/test_utils.py
@@ -3,7 +3,7 @@ import numpy as np
 
 from pyprep.utils import (
     _mat_round, _mat_quantile, _mat_iqr, _get_random_subset, _correlate_arrays,
-    _eeglab_create_highpass
+    _eeglab_create_highpass, print_progress
 )
 
 
@@ -114,3 +114,43 @@ def test_eeglab_create_highpass():
     expected_val = 0.9961
     actual_val = vals[len(vals) // 2]
     assert np.isclose(expected_val, actual_val, atol=0.001)
+
+
+def test_print_progress(capsys):
+    """Test the function for printing progress updates within a loop."""
+    # Test printing start value
+    print_progress(1, 20)
+    captured = capsys.readouterr()
+    assert captured.out == "Progress: "
+
+    # Test printing end values
+    iterations = 27
+    for i in range(iterations):
+        print_progress(i + 1, iterations, every=0.2)
+    captured = capsys.readouterr()
+    assert captured.out == "Progress: 20%... 40%... 60%... 80%... 100%\n"
+
+    # Test printing of updates at right times
+    iterations = 176
+    for i in range(iterations):
+        print_progress(i + 1, iterations)
+        if (i + 1) == 17:
+            captured = capsys.readouterr()
+            assert captured.out == "Progress: "
+        elif (i + 1) == 18:
+            captured = capsys.readouterr()
+            assert captured.out == "10%... "
+            break
+
+    # Test shifted start value
+    iterations = 25
+    start = 5
+    for i in range(start, iterations + 1):
+        print_progress(i, iterations, start=start)
+        if i == 6:
+            captured = capsys.readouterr()
+            assert captured.out == "Progress: "
+        elif i == 7:
+            captured = capsys.readouterr()
+            assert captured.out == "10%... "
+            break