Implmentation of R2C FFT

Cargo.toml

@@ -23,6 +23,7 @@ complex-nums = ["dep:num-complex", "dep:bytemuck"]
 [dev-dependencies]
 criterion = "0.5.1"
 fftw = "0.8.0"
+realfft = "3.3.0"
 rand = "0.8.5"
 utilities = { path = "utilities" }
 

benches/README.md

@@ -78,6 +78,30 @@ The generated images will be saved in your working directory.
 | Active Toolchain          | nightly-x86_64-unknown-linux-gnu (default)                                                      |
 | Rustc Version             | rustc 1.79.0-nightly (7f2fc33da 2024-04-22)                                                     |
 
+## Measuruing throughput with `criterion`
+
+0. Install gnuplot if you want `criterion` to use gnuplot as the plotting backend.
+   On macOS, just run:
+```bash
+brew install gnuplot 
+```
+
+1. Run benchmarks
+```bash
+cargo bench
+```
+
+2. Open the report using a browser. The report will be located in the `target` directory. 
+   For example, from the root of the `PhastFT` repository (on macOS) run:
+```bash
+open -a firefox target/criterion/r2c_versus_c2c/report/index.html
+```
+
+### Results
+
+![Alt text](../assets/lines.png)
+
+
 ## Profiling
 
 Navigate to the cloned repo:

benches/bench.rs

@@ -1,6 +1,11 @@
 use criterion::{black_box, criterion_group, criterion_main, BenchmarkId, Criterion, Throughput};
+use realfft::num_complex::Complex;
+use realfft::RealFftPlanner;
+use utilities::gen_random_signal;
+
 use num_traits::Float;
 use phastft::{
+    fft::r2c_fft_f64,
     fft_32_with_opts_and_plan, fft_64_with_opts_and_plan,
     options::Options,
     planner::{Direction, Planner32, Planner64},
@@ -9,9 +14,67 @@ use rand::{
     distributions::{Distribution, Standard},
     thread_rng, Rng,
 };
-use utilities::rustfft::num_complex::Complex;
 use utilities::rustfft::FftPlanner;
 
+fn benchmark_r2c_vs_c2c(c: &mut Criterion) {
+    let sizes = vec![1 << 10, 1 << 12, 1 << 14, 1 << 16, 1 << 18, 1 << 20];
+
+    let mut group = c.benchmark_group("r2c_versus_c2c");
+    for &size in &sizes {
+        group.throughput(Throughput::Elements(size as u64));
+
+        group.bench_with_input(BenchmarkId::new("r2c_fft", size), &size, |b, &size| {
+            let mut s_re = vec![0.0; size];
+            let mut s_im = vec![0.0; size];
+            gen_random_signal(&mut s_re, &mut s_im);
+            let mut output_re = vec![0.0; size];
+            let mut output_im = vec![0.0; size];
+
+            b.iter(|| {
+                r2c_fft_f64(
+                    black_box(&s_re),
+                    black_box(&mut output_re),
+                    black_box(&mut output_im),
+                );
+            });
+        });
+
+        group.bench_with_input(BenchmarkId::new("c2c_fft", size), &size, |b, &size| {
+            let mut s_re = vec![0.0; size];
+            let mut s_im = vec![0.0; size];
+            gen_random_signal(&mut s_re, &mut s_im);
+            s_im = vec![0.0; size];
+
+            let options = Options::guess_options(s_re.len());
+            let planner = Planner64::new(s_re.len(), Direction::Reverse);
+
+            b.iter(|| {
+                fft_64_with_opts_and_plan(
+                    black_box(&mut s_re),
+                    black_box(&mut s_im),
+                    black_box(&options),
+                    black_box(&planner),
+                );
+            });
+        });
+
+        group.bench_with_input(BenchmarkId::new("real_fft", size), &size, |b, &size| {
+            let mut s_re = vec![0.0; size];
+            let mut s_im = vec![0.0; size];
+            gen_random_signal(&mut s_re, &mut s_im);
+            let mut output = vec![Complex::default(); s_re.len() / 2 + 1];
+            let mut planner = RealFftPlanner::<f64>::new();
+
+            b.iter(|| {
+                planner
+                    .plan_fft_forward(s_re.len())
+                    .process(&mut s_re, &mut output)
+                    .expect("fft.process() failed!");
+            });
+        });
+    }
+}
+
 const LENGTHS: &[usize] = &[
     6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24,
 ];
@@ -170,11 +233,13 @@ fn benchmark_inverse_f64(c: &mut Criterion) {
     }
 }
 
-criterion_group!(
-    benches,
-    benchmark_forward_f32,
-    benchmark_inverse_f32,
-    benchmark_forward_f64,
-    benchmark_inverse_f64
-);
+fn critertion_benchmark(c: &mut Criterion) {
+    benchmark_forward_f32(c);
+    benchmark_inverse_f32(c);
+    benchmark_forward_f64(c);
+    benchmark_inverse_f64(c);
+    benchmark_r2c_vs_c2c(c);
+}
+
+criterion_group!(benches, critertion_benchmark,);
 criterion_main!(benches);

examples/profile.rs

@@ -1,29 +1,38 @@
 use std::env;
 use std::str::FromStr;
 
-use utilities::gen_random_signal;
-
-use phastft::fft_64_with_opts_and_plan;
-use phastft::options::Options;
-use phastft::planner::{Direction, Planner64};
-
-fn benchmark_fft_64(n: usize) {
+use phastft::fft::r2c_fft_f64;
+// use phastft::fft_64;
+// use phastft::fft_64_with_opts_and_plan;
+// use phastft::options::Options;
+// use phastft::planner::{Direction, Planner64};
+// use utilities::gen_random_signal;
+
+// fn benchmark_fft_64(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 planner = Planner64::new(reals.len(), Direction::Forward);
+//     let opts = Options::guess_options(reals.len());
+//
+//     fft_64_with_opts_and_plan(&mut reals, &mut imags, &opts, &planner);
+// }
+
+fn benchmark_r2c_fft(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 planner = Planner64::new(reals.len(), Direction::Forward);
-    let opts = Options::guess_options(reals.len());
-
-    fft_64_with_opts_and_plan(&mut reals, &mut imags, &opts, &planner);
+    let reals: Vec<f64> = (1..=big_n).map(|i| i as f64).collect();
+    let mut output_re = vec![0.0; big_n];
+    let mut output_im = vec![0.0; big_n];
+    r2c_fft_f64(&reals, &mut output_re, &mut output_im);
 }
 
 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_fft_64(n);
+    benchmark_r2c_fft(n);
+    // benchmark_fft_64(n);
 }

pyphastft/audio_visual.py

@@ -0,0 +1,111 @@
+import sys
+import numpy as np
+import pyaudio
+import pyqtgraph as pg
+from pyphastft import rfft
+from pyqtgraph.Qt import QtWidgets, QtCore
+
+class RealTimeAudioSpectrum(QtWidgets.QWidget):
+	def __init__(self, parent=None):
+		super(RealTimeAudioSpectrum, self).__init__(parent)
+		self.n_fft_bins = 1024
+		self.n_display_bins = 64
+		self.sample_rate = 44100
+		self.smoothing_factor = 0.1  # Fine-tuned smoothing factor
+		self.ema_fft_data = np.zeros(self.n_display_bins)
+		self.init_ui()
+		self.init_audio_stream()
+
+	def init_ui(self):
+		self.layout = QtWidgets.QVBoxLayout(self)
+		self.plot_widget = pg.PlotWidget()
+		self.layout.addWidget(self.plot_widget)
+
+		self.plot_widget.setBackground("k")
+		self.plot_item = self.plot_widget.getPlotItem()
+		self.plot_item.setTitle(
+				"Real-Time Audio Spectrum Visualizer powered by PhastFT",
+				color="w",
+				size="16pt",
+				)
+
+		self.plot_item.getAxis("left").hide()
+		self.plot_item.getAxis("bottom").hide()
+
+		self.plot_item.setXRange(0, self.sample_rate / 2, padding=0)
+		self.plot_item.setYRange(0, 1, padding=0)
+
+		self.bar_width = (self.sample_rate / 2) / self.n_display_bins * 0.8
+		self.freqs = np.linspace(
+				0, self.sample_rate / 2, self.n_display_bins, endpoint=False
+				)
+
+		self.bar_graph = pg.BarGraphItem(
+				x=self.freqs,
+				height=np.zeros(self.n_display_bins),
+				width=self.bar_width,
+				brush=pg.mkBrush("m"),
+				)
+		self.plot_item.addItem(self.bar_graph)
+
+		self.timer = QtCore.QTimer()
+		self.timer.timeout.connect(self.update)
+		self.timer.start(50)
+
+	def init_audio_stream(self):
+		self.p = pyaudio.PyAudio()
+		self.stream = self.p.open(
+				format=pyaudio.paFloat32,
+				channels=1,
+				rate=self.sample_rate,
+				input=True,
+				frames_per_buffer=self.n_fft_bins,
+				stream_callback=self.audio_callback,
+				)
+		self.stream.start_stream()
+
+	def audio_callback(self, in_data, frame_count, time_info, status):
+		audio_data = np.frombuffer(in_data, dtype=np.float32)
+		audio_data = np.ascontiguousarray(audio_data, dtype=np.float64)
+		reals, imags = rfft(audio_data, direction="f")
+		reals = np.ascontiguousarray(reals)
+		imags = np.ascontiguousarray(imags)
+		fft_magnitude = np.sqrt(reals**2 + imags**2)[: self.n_fft_bins // 2]
+
+		new_fft_data = np.interp(
+				np.linspace(0, len(fft_magnitude), self.n_display_bins),
+				np.arange(len(fft_magnitude)),
+				fft_magnitude,
+				)
+
+		new_fft_data = np.log1p(new_fft_data)  # Apply logarithmic scaling
+
+		self.ema_fft_data = self.ema_fft_data * self.smoothing_factor + new_fft_data * (
+				1.0 - self.smoothing_factor
+				)
+
+		return in_data, pyaudio.paContinue
+
+	def update(self):
+		# Normalize the FFT data to ensure it fits within the display range
+		max_value = np.max(self.ema_fft_data)
+		if max_value > 0:
+			normalized_fft_data = self.ema_fft_data / max_value
+		else:
+			normalized_fft_data = self.ema_fft_data
+
+		self.bar_graph.setOpts(height=normalized_fft_data, width=self.bar_width)
+
+	def closeEvent(self, event):
+		self.stream.stop_stream()
+		self.stream.close()
+		self.p.terminate()
+		event.accept()
+
+if __name__ == "__main__":
+	app = QtWidgets.QApplication(sys.argv)
+	window = RealTimeAudioSpectrum()
+	window.show()
+	sys.exit(app.exec_())
+
+

pyphastft/requirements.txt

@@ -0,0 +1,7 @@
+maturin==1.6.0
+pip==24.0
+PyAudio==0.2.14
+pylint==3.2.3
+pyphastft==0.2.1
+PyQt5==5.15.10
+pyqtgraph==0.13.7

pyphastft/src/lib.rs

@@ -1,5 +1,5 @@
-use numpy::PyReadwriteArray1;
-use phastft::{fft_64 as fft_64_rs, planner::Direction};
+use numpy::{PyReadonlyArray1, PyReadwriteArray1};
+use phastft::{fft::r2c_fft_f64, fft_64 as fft_64_rs, planner::Direction};
 use pyo3::prelude::*;
 
 #[pyfunction]
@@ -18,9 +18,27 @@ fn fft(mut reals: PyReadwriteArray1<f64>, mut imags: PyReadwriteArray1<f64>, dir
     );
 }
 
+#[pyfunction]
+fn rfft(reals: PyReadonlyArray1<f64>, direction: char) -> (Vec<f64>, Vec<f64>) {
+    assert!(direction == 'f' || direction == 'r');
+    let _dir = if direction == 'f' {
+        Direction::Forward
+    } else {
+        Direction::Reverse
+    };
+
+    let big_n = reals.as_slice().unwrap().len();
+
+    let mut output_re = vec![0.0; big_n];
+    let mut output_im = vec![0.0; big_n];
+    r2c_fft_f64(reals.as_slice().unwrap(), &mut output_re, &mut output_im);
+    (output_re, output_im)
+}
+
 /// A Python module implemented in Rust.
 #[pymodule]
 fn pyphastft(_py: Python, m: &PyModule) -> PyResult<()> {
     m.add_function(wrap_pyfunction!(fft, m)?)?;
+    m.add_function(wrap_pyfunction!(rfft, m)?)?;
     Ok(())
 }