Merge pull request #7 from HEnquist/odds

Odds

.github/workflows/run_tests.yml

@@ -0,0 +1,55 @@
+on: [pull_request]
+
+name: CI
+
+jobs:
+  check:
+    name: Check+Test
+    runs-on: ubuntu-latest
+    strategy:
+      matrix:
+        rust:
+          - stable
+          - beta
+          - nightly
+          - 1.37
+    steps:
+      - name: Checkout sources
+        uses: actions/checkout@v2
+
+      - name: Install toolchain
+        uses: actions-rs/toolchain@v1
+        with:
+          profile: minimal
+          toolchain: ${{ matrix.rust }}
+          override: true
+
+      - name: Run cargo check
+        uses: actions-rs/cargo@v1
+        with:
+          command: check
+
+      - name: Run cargo test
+        uses: actions-rs/cargo@v1
+        with:
+          command: test
+  fmt:
+    name: Rustfmt
+    runs-on: ubuntu-latest
+    steps:
+      - name: Checkout sources
+        uses: actions/checkout@v2
+
+      - name: Install toolchain
+        uses: actions-rs/toolchain@v1
+        with:
+          profile: minimal
+          toolchain: stable
+          override: true
+          components: rustfmt
+
+      - name: Run cargo fmt
+        uses: actions-rs/cargo@v1
+        with:
+          command: fmt
+          args: -- --check

Cargo.toml

@@ -1,6 +1,6 @@
 [package]
 name = "realfft"
-version = "0.4.0"
+version = "1.0.0"
 authors = ["HEnquist <henrik.enquist@gmail.com>"]
 edition = "2018"
 description = "Real-to-complex FFT and complex-to-real iFFT for Rust"
@@ -17,6 +17,7 @@ rustfft = "5.0.0"
 
 [dev-dependencies]
 criterion = "0.3"
+rand = "0.8.1"
 
 [[bench]]
 name = "realfft"

README.md

@@ -1,94 +1,128 @@
-# RealFFT: Real-to-complex FFT and complex-to-real iFFT based on RustFFT
+# realfft
 
-This library is a wrapper for RustFFT that enables faster computations when the input data is real.
-It packs a 2N long real vector into an N long complex vector, which is transformed using a standard FFT.
+## RealFFT: Real-to-complex FFT and complex-to-real iFFT based on RustFFT
+
+This library is a wrapper for RustFFT that enables performing FFT of real-valued data.
+The API is designed to be as similar as possible to RustFFT.
+
+Using this library instead of RustFFT directly avoids the need of converting real-valued data to complex before performing a FFT.
+If the length is even, it also enables faster computations by using a complex FFT of half the length.
+It then packs a 2N long real vector into an N long complex vector, which is transformed using a standard FFT.
 It then post-processes the result to give only the first half of the complex spectrum, as an N+1 long complex vector.
 
 The iFFT goes through the same steps backwards, to transform an N+1 long complex spectrum to a 2N long real result.
 
 The speed increase compared to just converting the input to a 2N long complex vector
 and using a 2N long FFT depends on the length f the input data.
 The largest improvements are for long FFTs and for lengths over around 1000 elements there is an improvement of about a factor 2.
-The difference shrinks for shorter lengths, and around 100 elements there is no longer any difference.  
+The difference shrinks for shorter lengths, and around 30 elements there is no longer any difference.
 
-## Why use real-to-complex fft?
-### Using a complex-to-complex fft
-A simple way to get the fft of a rea values vector is to convert it to complex, and using a complex-to-complex fft.
+### Why use real-to-complex FFT?
+#### Using a complex-to-complex FFT
+A simple way to get the FFT of a rea values vector is to convert it to complex, and using a complex-to-complex FFT.
 
-Let's assume `x` is a 6 element long real vector: 
-```text
+Let's assume `x` is a 6 element long real vector:
+```
 x = [x0r, x1r, x2r, x3r, x4r, x5r]
 ```
 
-Converted to complex, using the notation `(xNr, xNi)` for the complex value `xN`, this becomes: 
-```text
+We now convert `x` to complex by adding an imaginary part with value zero. Using the notation `(xNr, xNi)` for the complex value `xN`, this becomes:
+```
 x_c = [(x0r, 0), (x1r, 0), (x2r, 0), (x3r, 0), (x4r, 0, (x5r, 0)]
 ```
 
-
-The general result of `X = FFT(x)` is:
-```text
-X = [(X0r, X0i), (X1r, X1i), (X2r, X2i), (X3r, X3i), (X4r, X4i), (X5r, X5i)]
+Performing a normal complex FFT, the result of `FFT(x_c)` is:
+```
+FFT(x_c) = [(X0r, X0i), (X1r, X1i), (X2r, X2i), (X3r, X3i), (X4r, X4i), (X5r, X5i)]
 ```
 
-However, because our `x` was real-valued, some of this is redundant:
-```text
+But because our `x_c` is real-valued (all imaginary parts are zero), some of this becomes redundant:
+```
 FFT(x_c) = [(X0r, 0), (X1r, X1i), (X2r, X2i), (X3r, 0), (X2r, -X2i), (X1r, -X1i)]
 ```
 
-As we can see, the output contains a fair bit of redundant data. But it still takes time for the FFT to calculate these values. Converting the input data to complex also takes a little bit of time.
+The last two values are the complex conjugates of `X1` and `X2`. Additionally, `X0i` and `X3i` are zero.
+As we can see, the output contains 6 independent values, and the rest is redundant.
+But it still takes time for the FFT to calculate the redundant values.
+Converting the input data to complex also takes a little bit of time.
+
+If the length of `x` instead had been 7, result would have been:
+```
+FFT(x_c) = [(X0r, 0), (X1r, X1i), (X2r, X2i), (X3r, X3i), (X3r, -X3i), (X2r, -X2i), (X1r, -X1i)]
+```
+
+The result is similar, but this time there is no zero at `X3i`. Also in this case we got the same number of indendent values as we started with.
 
-### real-to-complex
-Using a real-to-complex fft removes the need for converting the input data to complex.
+#### Real-to-complex
+Using a real-to-complex FFT removes the need for converting the input data to complex.
 It also avoids caclulating the redundant output values.
 
-The result is: 
-```text
+The result for 6 elements is:
+```
 RealFFT(x) = [(X0r, 0), (X1r, X1i), (X2r, X2i), (X3r, 0)]
 ```
 
-This is the data layout output by the real-to-complex fft, and the one expected as input to the complex-to-real ifft.
+The result for 7 elements is:
+```
+RealFFT(x) = [(X0r, 0), (X1r, X1i), (X2r, X2i), (X3r, X3i)]
+```
+
+This is the data layout output by the real-to-complex FFT, and the one expected as input to the complex-to-real iFFT.
 
-## Scaling
+### Scaling
 RealFFT matches the behaviour of RustFFT and does not normalize the output of either FFT of iFFT. To get normalized results, each element must be scaled by `1/sqrt(length)`. If the processing involves both an FFT and an iFFT step, it is advisable to merge the two normalization steps to a single, by scaling by `1/length`.
 
-## Documentation
+### Documentation
 
 The full documentation can be generated by rustdoc. To generate and view it run:
-```text
+```
 cargo doc --open
 ```
 
-## Benchmarks
+### Benchmarks
 
 To run a set of benchmarks comparing real-to-complex FFT with standard complex-to-complex, type:
-```text
+```
 cargo bench
 ```
 The results are printed while running, and are compiled into an html report containing much more details.
 To view, open `target/criterion/report/index.html` in a browser.
 
-## Example
+### Example
 Transform a vector, and then inverse transform the result.
 ```rust
-use realfft::{ComplexToReal, RealToComplex};
+use realfft::RealFftPlanner;
 use rustfft::num_complex::Complex;
 use rustfft::num_traits::Zero;
 
-// make dummy input vector, spectrum and output vectors
-let mut indata = vec![0.0f64; 256];
-let mut spectrum: Vec<Complex<f64>> = vec![Complex::zero(); 129];
-let mut outdata: Vec<f64> = vec![0.0; 256];
+let length = 256;
 
-//create an FFT and forward transform the input data
-let mut r2c = RealToComplex::<f64>::new(256).unwrap();
+// make a planner
+let mut real_planner = RealFftPlanner::<f64>::new();
+
+// create a FFT
+let r2c = real_planner.plan_fft_forward(length);
+// make input and output vectors
+let mut indata = r2c.make_input_vec();
+let mut spectrum = r2c.make_output_vec();
+
+// Are they the length we expect?
+assert_eq!(indata.len(), length);
+assert_eq!(spectrum.len(), length/2+1);
+
+// Forward transform the input data
 r2c.process(&mut indata, &mut spectrum).unwrap();
 
-// create an iFFT and inverse transform the spectum
-let mut c2r = ComplexToReal::<f64>::new(256).unwrap();
-c2r.process(&spectrum, &mut outdata).unwrap();
+// create an iFFT and an output vector
+let c2r = real_planner.plan_fft_inverse(length);
+let mut outdata = c2r.make_output_vec();
+assert_eq!(outdata.len(), length);
+
+c2r.process(&mut spectrum, &mut outdata).unwrap();
 ```
 
-## Compatibility
+### Compatibility
+
+The `realfft` crate requires rustc version 1.37 or newer.
 
-The `realfft` crate requires rustc version 1.37 or newer.
+License: MIT

benches/realfft.rs

@@ -2,7 +2,7 @@ use criterion::{criterion_group, criterion_main, Bencher, BenchmarkId, Criterion
 extern crate realfft;
 extern crate rustfft;
 
-use realfft::RealToComplex;
+use realfft::RealFftPlanner;
 use rustfft::num_complex::Complex;
 
 /// Times just the FFT execution (not allocation and pre-calculation)
@@ -24,7 +24,8 @@ fn bench_fft(b: &mut Bencher, len: usize) {
 }
 
 fn bench_realfft(b: &mut Bencher, len: usize) {
-    let mut fft = RealToComplex::<f64>::new(len).unwrap();
+    let mut planner = RealFftPlanner::<f64>::new();
+    let fft = planner.plan_fft_forward(len);
 
     let mut signal = vec![0_f64; len];
     let mut spectrum = vec![
@@ -34,27 +35,97 @@ fn bench_realfft(b: &mut Bencher, len: usize) {
         };
         len / 2 + 1
     ];
-    b.iter(|| fft.process(&mut signal, &mut spectrum));
+    let mut scratch = vec![Complex::from(0.0); fft.get_scratch_len()];
+    b.iter(|| fft.process_with_scratch(&mut signal, &mut spectrum, &mut scratch));
+}
+
+/// Times just the FFT execution (not allocation and pre-calculation)
+/// for a given length
+fn bench_ifft(b: &mut Bencher, len: usize) {
+    let mut planner = rustfft::FftPlanner::new();
+    let fft = planner.plan_fft_inverse(len);
+    let mut scratch = vec![Complex::from(0.0); fft.get_outofplace_scratch_len()];
+
+    let mut signal = vec![
+        Complex {
+            re: 0_f64,
+            im: 0_f64
+        };
+        len
+    ];
+    let mut spectrum = signal.clone();
+    b.iter(|| fft.process_outofplace_with_scratch(&mut signal, &mut spectrum, &mut scratch));
 }
 
-fn bench_pow2(c: &mut Criterion) {
-    let mut group = c.benchmark_group("Powers of 2");
-    for i in [64, 128, 256, 512, 4096, 65536].iter() {
+fn bench_realifft(b: &mut Bencher, len: usize) {
+    let mut planner = RealFftPlanner::<f64>::new();
+    let fft = planner.plan_fft_inverse(len);
+
+    let mut signal = vec![0_f64; len];
+    let mut spectrum = vec![
+        Complex {
+            re: 0_f64,
+            im: 0_f64
+        };
+        len / 2 + 1
+    ];
+    let mut scratch = vec![Complex::from(0.0); fft.get_scratch_len()];
+    b.iter(|| fft.process_with_scratch(&mut spectrum, &mut signal, &mut scratch));
+}
+
+fn bench_pow2_fw(c: &mut Criterion) {
+    let mut group = c.benchmark_group("Fw Powers of 2");
+    for i in [8, 16, 32, 64, 128, 256, 1024, 4096, 65536].iter() {
         group.bench_with_input(BenchmarkId::new("Complex", i), i, |b, i| bench_fft(b, *i));
         group.bench_with_input(BenchmarkId::new("Real", i), i, |b, i| bench_realfft(b, *i));
     }
     group.finish();
 }
 
-fn bench_pow7(c: &mut Criterion) {
-    let mut group = c.benchmark_group("Powers of 7");
-    for i in [2 * 343, 2 * 2401, 2 * 16807].iter() {
-        group.bench_with_input(BenchmarkId::new("Complex", i), i, |b, i| bench_fft(b, *i));
-        group.bench_with_input(BenchmarkId::new("Real", i), i, |b, i| bench_realfft(b, *i));
+fn bench_pow2_inv(c: &mut Criterion) {
+    let mut group = c.benchmark_group("Inv Powers of 2");
+    for i in [8, 16, 32, 64, 128, 256, 1024, 4096, 65536].iter() {
+        group.bench_with_input(BenchmarkId::new("Complex", i), i, |b, i| bench_ifft(b, *i));
+        group.bench_with_input(BenchmarkId::new("Real", i), i, |b, i| bench_realifft(b, *i));
+    }
+    group.finish();
+}
+
+//fn bench_pow7(c: &mut Criterion) {
+//    let mut group = c.benchmark_group("Powers of 7");
+//    for i in [2 * 343, 2 * 2401, 2 * 16807].iter() {
+//        group.bench_with_input(BenchmarkId::new("Complex", i), i, |b, i| bench_fft(b, *i));
+//        group.bench_with_input(BenchmarkId::new("Real", i), i, |b, i| bench_realfft(b, *i));
+//    }
+//    group.finish();
+//}
+
+fn bench_range_fw(c: &mut Criterion) {
+    let mut group = c.benchmark_group("Fw Range 1022-1025");
+    for i in 1022..1026 {
+        group.bench_with_input(BenchmarkId::new("Complex", i), &i, |b, i| bench_fft(b, *i));
+        group.bench_with_input(BenchmarkId::new("Real", i), &i, |b, i| bench_realfft(b, *i));
+    }
+    group.finish();
+}
+
+fn bench_range_inv(c: &mut Criterion) {
+    let mut group = c.benchmark_group("Inv Range 1022-1025");
+    for i in 1022..1026 {
+        group.bench_with_input(BenchmarkId::new("Complex", i), &i, |b, i| bench_ifft(b, *i));
+        group.bench_with_input(BenchmarkId::new("Real", i), &i, |b, i| {
+            bench_realifft(b, *i)
+        });
     }
     group.finish();
 }
 
-criterion_group!(benches, bench_pow2, bench_pow7);
+criterion_group!(
+    benches,
+    bench_pow2_fw,
+    bench_range_fw,
+    bench_pow2_inv,
+    bench_range_inv
+);
 
 criterion_main!(benches);