Add more examples and fix typos in readme/plots

- Added rust bindings of FFTW as an example, which will be used for
  benchmarks

- Add fftw (rust bindings) crate as a dev-dependency

- Add an example of using pyphastft to reproduce an example use case of
  FFT from the FFT wikipedia page

- Fix typos in the README and distinguish pyphastft from phastft in the
  python benchmarks plots

Cargo.toml

@@ -13,6 +13,7 @@ exclude = ["assets", "scripts", "benches"]
 
 [dev-dependencies]
 utilities = { path = "utilities" }
+fftw = "0.8.0"
 
 [profile.release]
 codegen-units = 1

README.md

@@ -36,16 +36,16 @@ Transform (FFT) library written in pure Rust.
 PhastFT is designed around the capabilities and limitations of modern hardware (that is, anything made in the last 10
 years or so).
 
-The two major bottlenecks in FFT are the **CPU cycles** and **memory accesses.**
+The two major bottlenecks in FFT are the **CPU cycles** and **memory accesses**.
 
 We picked an efficient, general-purpose FFT algorithm. Our implementation can make use of latest CPU features such as
 AVX-512, but performs well even without them.
 
 Our key insight for speeding up memory accesses is that FFT is equivalent to applying gates to all qubits in `[0, n)`.
-This creates to oppurtunity to leverage the same memory access patterns as
+This creates the opportunity to leverage the same memory access patterns as
 a [high-performance quantum state simulator](https://github.com/QuState/spinoza).
 
-We also use the Cache-Optimal Bit Reveral
+We also use the Cache-Optimal Bit Reversal
 Algorithm ([COBRA](https://csaws.cs.technion.ac.il/~itai/Courses/Cache/bit.pdf))
 on large datasets and optionally run it on 2 parallel threads, accelerating it even further.
 

benches/benchmark_plots.py

@@ -22,8 +22,8 @@ def read_file(filepath: str) -> list[float]:
     return y
 
 
-def get_figure_of_interest(vals: list[int]) -> float:
-    return np.median(vals)
+def get_figure_of_interest(vals: list[float]) -> float:
+    return np.min(vals)
 
 
 def build_and_clean_data(

benches/py_benchmarks.py

@@ -118,8 +118,8 @@ def plot_elapsed_times(data: dict) -> None:
     phastft_timings = np.asarray(data["phastft_time"])
 
     plt.plot(index, np_fft_timings, label="NumPy FFT", lw=0.8)
-    plt.plot(index, pyfftw_timings, label="PyFFTW FFT", lw=0.8)
-    plt.plot(index, phastft_timings, label="PhastFT", lw=0.8)
+    plt.plot(index, pyfftw_timings, label="pyFFTW", lw=0.8)
+    plt.plot(index, phastft_timings, label="pyPhastFT", lw=0.8)
 
     plt.title("pyPhastFT vs. pyFFTW vs. NumPy FFT")
     plt.xticks(fontsize=8, rotation=-45)
@@ -143,7 +143,7 @@ def grouped_bar_plot(data: dict, start=0, end=1):
         {
             "NumPy fft": np.ones(len(index)),
             "pyFFTW": pyfftw_timings / np_fft_timings,
-            "PhastFT": phastft_timings / np_fft_timings,
+            "pyPhastFT": phastft_timings / np_fft_timings,
         },
         index=index,
     )

examples/fftwrb.rs

@@ -0,0 +1,47 @@
+use std::{env, ptr::slice_from_raw_parts_mut, str::FromStr};
+
+use fftw::{
+    array::AlignedVec,
+    plan::{C2CPlan, C2CPlan64},
+    types::{Flag, Sign},
+};
+use utilities::{gen_random_signal, rustfft::num_complex::Complex};
+
+fn benchmark_fftw(n: usize) {
+    let big_n = 1 << n;
+
+    let mut reals = vec![0.0; big_n];
+    let mut imags = vec![0.0; big_n];
+
+    gen_random_signal(&mut reals, &mut imags);
+    let mut nums = AlignedVec::new(big_n);
+    reals
+        .drain(..)
+        .zip(imags.drain(..))
+        .zip(nums.iter_mut())
+        .for_each(|((re, im), z)| *z = Complex::new(re, im));
+
+    let now = std::time::Instant::now();
+    C2CPlan64::aligned(
+        &[big_n],
+        Sign::Backward,
+        Flag::DESTROYINPUT | Flag::ESTIMATE,
+    )
+    .unwrap()
+    .c2c(
+        // SAFETY: See above comment.
+        unsafe { &mut *slice_from_raw_parts_mut(nums.as_mut_ptr(), big_n) },
+        &mut nums,
+    )
+    .unwrap();
+    let elapsed = now.elapsed().as_micros();
+    println!("{elapsed}");
+}
+
+fn main() {
+    let args: Vec<String> = env::args().collect();
+    assert_eq!(args.len(), 2, "Usage {} <n>", args[0]);
+
+    let n = usize::from_str(&args[1]).unwrap();
+    benchmark_fftw(n);
+}

pyphastft/example.py

@@ -0,0 +1,52 @@
+import matplotlib.pyplot as plt
+import numpy as np
+
+from pyphastft import fft
+
+
+def main():
+    fs = 100  # Sampling frequency (100 samples/second for this synthetic example)
+    t_max = 6  # maximum time in "seconds"
+
+    # Find the next lower power of 2 for the number of samples
+    n_samples = 2 ** int(np.log2(t_max * fs))
+
+    t = np.linspace(0, n_samples / fs, n_samples, endpoint=False)  # Adjusted time vector
+
+    # Generate the signal
+    s_re = 2 * np.sin(2 * np.pi * t) + np.sin(2 * np.pi * 10 * t)
+    s_im = np.ascontiguousarray([0.0] * len(s_re), dtype=np.float64)
+
+    # Plot the original signal
+    plt.figure(figsize=(10, 7))
+
+    plt.subplot(2, 1, 1)
+    plt.plot(t, s_re, label="f(x) = 2sin(x) + sin(10x)")
+    plt.title("signal: f(x) = 2sin(x) + sin(10x)")
+    plt.xlabel("time [seconds]")
+    plt.ylabel("f(x)")
+    plt.legend()
+
+    # Perform FFT
+    fft(s_re, s_im, direction="f")
+
+    # Plot the magnitude spectrum of the FFT result
+    plt.subplot(2, 1, 2)
+    plt.plot(
+        np.abs(s_re),
+        label="frequency spectrum",
+    )
+    plt.title("Signal after FFT")
+    plt.xlabel("frequency (in Hz)")
+    plt.ylabel("|FFT(f(x))|")
+
+    # only show up to 11 Hz as in the wiki example
+    plt.xlim(0, 11)
+
+    plt.legend()
+    plt.tight_layout()
+    plt.savefig("wiki_fft_example.png", dpi=600)
+
+
+if __name__ == "__main__":
+    main()