You signed in with another tab or window. Reload to refresh your session.You signed out in another tab or window. Reload to refresh your session.You switched accounts on another tab or window. Reload to refresh your session.Dismiss alert
docs: promote Rust to first-class implementation in top-level README
Rust was previously framed as 'optional Python acceleration' but the
crate ships both a standalone CLI (rust/src/main.rs) and a PyO3 Python
extension, so it belongs alongside Matlab / Python / R in the
implementation list. Reorganizes the TOC + section ordering
accordingly and folds the Rust-vs-Python equivalence note into the
Numerical Differences section.
Copy file name to clipboardExpand all lines: README.md
+25-40Lines changed: 25 additions & 40 deletions
Display the source diff
Display the rich diff
Original file line number
Diff line number
Diff line change
@@ -6,8 +6,8 @@
6
6
*[General Information](#general-information)
7
7
*[Matlab Implementation](#matlab-implementation)
8
8
*[Python Implementation](#python-implementation)
9
+
*[Rust Implementation](#rust-implementation)
9
10
*[R Implementation](#r-implementation)
10
-
*[Rust Backend (optional acceleration for Python)](#rust-backend-optional-acceleration-for-python)
11
11
*[Parameters](#parameters)
12
12
*[Numerical Differences Between Implementations](#numerical-differences-between-implementations)
13
13
*[Status](#status)
@@ -16,7 +16,7 @@
16
16
*[Contact](#contact)
17
17
18
18
## General Information
19
-
This repository contains Matlab, Python, and R implementations of the multitaper spectrogram analysis described in the paper ["Sleep Neurophysiological Dynamics Through the Lens of Multitaper Spectral Analysis"](https://prerau.bwh.harvard.edu/publications/Physiology_Bethesda_2017_Prerau.pdf)<sup>1</sup>. Multitaper spectral estimation was developed in the early 1980s by David Thomson<sup>2</sup> and has been shown to have superior statistical properties compared with single-taper spectral estimates<sup>3,4</sup>. The multitaper method works by averaging together multiple independent spectra estimated from a single segment of data. The innovation of the multitaper method is that, instead of using a single-taper function to compute the spectrum, it uses multiple taper functions called discrete prolate spheroidal sequences (DPSS). Because DPSS tapers are uncorrelated with each other, they can be averaged together as if they were independent trials of the same condition, producing a spectrum with reduced variance compared to periodogram and single-taper estimation.
19
+
This repository contains Matlab, Python, Rust, and R implementations of the multitaper spectrogram analysis described in the paper ["Sleep Neurophysiological Dynamics Through the Lens of Multitaper Spectral Analysis"](https://prerau.bwh.harvard.edu/publications/Physiology_Bethesda_2017_Prerau.pdf)<sup>1</sup>. Multitaper spectral estimation was developed in the early 1980s by David Thomson<sup>2</sup> and has been shown to have superior statistical properties compared with single-taper spectral estimates<sup>3,4</sup>. The multitaper method works by averaging together multiple independent spectra estimated from a single segment of data. The innovation of the multitaper method is that, instead of using a single-taper function to compute the spectrum, it uses multiple taper functions called discrete prolate spheroidal sequences (DPSS). Because DPSS tapers are uncorrelated with each other, they can be averaged together as if they were independent trials of the same condition, producing a spectrum with reduced variance compared to periodogram and single-taper estimation.
20
20
21
21
Find videos describing the theory of spectral estimation and demonstrating how multitaper spectral estimation works at [http://sleepeeg.org/multitaper](http://sleepeeg.org/multitaper) on the Prerau Lab website.
22
22
@@ -33,52 +33,36 @@ Find videos describing the theory of spectral estimation and demonstrating how m
33
33
<br/>
34
34
35
35
## Python Implementation
36
-
***multitaper_spectrogram_python.py**: optimized implementation in Python with batched `scipy.fft.rfft` (multi-threaded via `workers=-1`), optional multiprocessing, and transparent acceleration via a Rust backend when installed.
36
+
***multitaper_spectrogram_python.py**: optimized implementation in Python with batched `scipy.fft.rfft` (multi-threaded via `workers=-1`), optional multiprocessing, and transparent acceleration via the Rust implementation when installed as a Python extension (`use_rust` kwarg).
37
37
***requirements.txt**: contains names and versions of non-standard library Python packages required to run multitaper_spectrogram_python.py
38
38
* See [the python implementation folder](./python/README.md) for usage information and the `use_rust` kwarg.
39
39
40
40
<br/>
41
41
42
-
## R Implementation
43
-
***multitaper_spectrogram_R.R**: baseline implementation in R with option for multiprocessing
44
-
* See [the R implementation folder](https://github.com/preraulab/multitaper_toolbox/tree/master/R) for usage information and other details of the Matlab implementation
42
+
## Rust Implementation
43
+
***multitaper_rs**: a Rust crate (`rust/src/lib.rs`) + CLI binary (`rust/src/main.rs`) implementing the core multitaper compute loop — per-window detrend, taper multiply, zero-padded rFFT (`realfft` / pocketfft), weighted power averaging across tapers, and one-sided PSD scaling. `rayon` parallelism over windows. `'unity'` and `'eigen'` taper weightings are implemented; `'adapt'` is not.
44
+
***Two use modes**, sharing the same compiled core:
45
+
1.**Standalone CLI** — `cargo build --release` then `target/release/multitaper_rs --data data.npy --tapers tapers.npy --fs 200 ...`. Reads / writes `.npy`, prints a `TIMING` line. DPSS tapers must be pre-computed (e.g. with `scipy.signal.windows.dpss`) and passed via `.npy`.
46
+
2.**Python extension** — the same crate exposes a PyO3 module. Install with `cd rust && maturin develop --release`; afterwards `multitaper_spectrogram(data, fs, use_rust=None)` transparently routes through Rust, falling back to pure Python when the extension is missing. Force with `use_rust=True`/`False`. A one-time banner to stderr announces the active backend.
47
+
***Speedup (x86_64 Python 3.9, numpy 1.26 under Rosetta):**
48
+
49
+
| Duration | Python (s) | Rust (s) | Speedup |
50
+
|---:|---:|---:|---:|
51
+
| 1 h | 0.60 | 0.03 | 19× |
52
+
| 4 h | 2.52 | 0.12 | 21× |
53
+
| 10 h | 6.60 | 0.68 | 10× |
54
+
| 24 h | 17.76 | 1.63 | 11× |
55
+
56
+
On native arm64 Python 3.11 (numpy 2.4 with Apple Accelerate) the pure-Python path is already much faster, so the Rust speedup is smaller (≈5–9× across the same durations) but still substantial.
57
+
***Equivalence:** across three parameter sets the Rust and pure-Python outputs agree to a max absolute difference of 1.67e-16 (bit-identical to f64 round-off).
58
+
***Limitations:**`weighting='adapt'` is not yet implemented (Python wrapper falls back automatically). DPSS tapers are not generated by the Rust crate — they must be pre-computed by the caller (Python: `scipy.signal.windows.dpss`; MATLAB: `dpss`).
59
+
* See [the rust implementation folder](./rust/README.md) for build instructions, CLI flags, and API details.
45
60
46
61
<br/>
47
62
48
-
## Rust Backend (optional acceleration for Python)
49
-
50
-
The Python module ships with an optional Rust backend that transparently accelerates the core compute loop (per-window detrend, taper multiply, zero-padded rFFT, weighted power sum, and one-sided PSD scaling). Callers do not need to change anything: if the compiled extension is importable, `multitaper_spectrogram(...)` uses it; if not, it falls back to the pure-Python path automatically. A one-time banner is printed to stderr announcing which backend is active.
51
-
52
-
**Install:**
53
-
```
54
-
cd multitaper/rust
55
-
pip install maturin
56
-
maturin develop --release
57
-
```
58
-
59
-
**Opt in / out:**
60
-
```python
61
-
multitaper_spectrogram(data, fs, use_rust=None) # auto (default)
62
-
multitaper_spectrogram(data, fs, use_rust=True) # force Rust (error if not installed)
63
-
multitaper_spectrogram(data, fs, use_rust=False) # force pure Python
64
-
```
65
-
66
-
**Speedup (x86_64 Python 3.9, numpy 1.26 under Rosetta):**
67
-
68
-
| Duration | Python (s) | Rust (s) | Speedup |
69
-
|---:|---:|---:|---:|
70
-
| 1 h | 0.60 | 0.03 | 19× |
71
-
| 4 h | 2.52 | 0.12 | 21× |
72
-
| 10 h | 6.60 | 0.68 | 10× |
73
-
| 24 h | 17.76 | 1.63 | 11× |
74
-
75
-
On native arm64 Python 3.11 (numpy 2.4 with Apple Accelerate) the pure-Python path is already much faster, so the Rust speedup is smaller (≈5–9× across the same durations) but still substantial.
76
-
77
-
**Equivalence:** across three parameter sets the Rust and pure-Python outputs agree to a max absolute difference of 1.67e-16 (bit-identical to f64 round-off).
78
-
79
-
**Limitations:**`weighting='adapt'` is not yet implemented in Rust; the wrapper falls back automatically. DPSS tapers are still computed by `scipy.signal.windows.dpss` and passed into Rust.
80
-
81
-
See [python/README.md](./python/README.md) and [rust/README.md](./rust/README.md) for details.
63
+
## R Implementation
64
+
***multitaper_spectrogram_R.R**: baseline implementation in R with option for multiprocessing
65
+
* See [the R implementation folder](https://github.com/preraulab/multitaper_toolbox/tree/master/R) for usage information and other details of the R implementation
82
66
83
67
<br/>
84
68
@@ -97,6 +81,7 @@ The spectral parameters used in all implementations of the multitaper spectrogra
97
81
98
82
## Numerical Differences Between Implementations
99
83
* In data point comparisons Matlab and Python implementation results tend to agree on average with precision on the order of at most 10^-13 with SD of at most 10^-10.
84
+
* Rust and Python implementations agree to 1.67e-16 max abs-diff (f64 round-off — effectively bit-identical) across the validated parameter sets.
100
85
* In data point comparisons R and Python implementation results tend to agree on average with precision on the order of at most 10^-11 with SD of at most 10^-9.
0 commit comments