scirs2_fft/
real_planner.rs

1//! Real FFT planner with trait object support
2//!
3//! This module provides trait object interfaces for real-to-complex and complex-to-real
4//! FFT operations, matching the API patterns used by realfft crate for easier migration.
5//!
6//! # Features
7//!
8//! - `RealToComplex` trait for forward real-to-complex FFT operations
9//! - `ComplexToReal` trait for inverse complex-to-real FFT operations
10//! - `RealFftPlanner` for creating and caching FFT plans
11//! - Support for both f32 and f64 precision
12//! - Thread-safe plan caching with Arc<dyn Trait>
13//!
14//! # Examples
15//!
16//! ```
17//! use scirs2_fft::real_planner::{RealFftPlanner, RealToComplex, ComplexToReal};
18//! use std::sync::Arc;
19//!
20//! // Create a planner
21//! let mut planner = RealFftPlanner::<f64>::new();
22//!
23//! // Plan forward FFT
24//! let forward_fft = planner.plan_fft_forward(1024);
25//!
26//! // Plan inverse FFT
27//! let inverse_fft = planner.plan_fft_inverse(1024);
28//!
29//! // Use in struct (common VoiRS pattern)
30//! struct AudioProcessor {
31//!     forward: Arc<dyn RealToComplex<f64>>,
32//!     backward: Arc<dyn ComplexToReal<f64>>,
33//! }
34//! ```
35
36use crate::error::{FFTError, FFTResult};
37use scirs2_core::numeric::Complex;
38use scirs2_core::numeric::Float;
39use std::fmt::Debug;
40use std::sync::{Arc, Mutex};
41
42/// Trait for real-to-complex FFT operations
43///
44/// This trait defines the interface for forward FFT transforms that convert
45/// real-valued input data to complex-valued frequency domain output.
46pub trait RealToComplex<T: Float>: Send + Sync {
47    /// Process a real-valued input and produce complex-valued output
48    ///
49    /// # Arguments
50    ///
51    /// * `input` - Real-valued input samples
52    /// * `output` - Complex-valued frequency domain output (length = input.len()/2 + 1)
53    fn process(&self, input: &[T], output: &mut [Complex<T>]);
54
55    /// Get the length of the input this FFT is configured for
56    fn len(&self) -> usize;
57
58    /// Check if this FFT is empty (length 0)
59    fn is_empty(&self) -> bool {
60        self.len() == 0
61    }
62
63    /// Get the length of the output this FFT produces
64    fn output_len(&self) -> usize {
65        self.len() / 2 + 1
66    }
67}
68
69/// Trait for complex-to-real FFT operations
70///
71/// This trait defines the interface for inverse FFT transforms that convert
72/// complex-valued frequency domain data back to real-valued time domain output.
73pub trait ComplexToReal<T: Float>: Send + Sync {
74    /// Process a complex-valued input and produce real-valued output
75    ///
76    /// # Arguments
77    ///
78    /// * `input` - Complex-valued frequency domain samples (length = output.len()/2 + 1)
79    /// * `output` - Real-valued time domain output
80    fn process(&self, input: &[Complex<T>], output: &mut [T]);
81
82    /// Get the length of the output this IFFT is configured for
83    fn len(&self) -> usize;
84
85    /// Check if this IFFT is empty (length 0)
86    fn is_empty(&self) -> bool {
87        self.len() == 0
88    }
89
90    /// Get the length of the input this IFFT expects
91    fn input_len(&self) -> usize {
92        self.len() / 2 + 1
93    }
94}
95
96/// Real FFT plan implementation for f64
97struct RealFftPlanF64 {
98    length: usize,
99    planner: Arc<Mutex<rustfft::FftPlanner<f64>>>,
100}
101
102impl RealFftPlanF64 {
103    fn new(length: usize, planner: Arc<Mutex<rustfft::FftPlanner<f64>>>) -> Self {
104        Self { length, planner }
105    }
106}
107
108impl RealToComplex<f64> for RealFftPlanF64 {
109    fn process(&self, input: &[f64], output: &mut [Complex<f64>]) {
110        // Validate input/output lengths
111        assert_eq!(
112            input.len(),
113            self.length,
114            "Input length {} doesn't match plan length {}",
115            input.len(),
116            self.length
117        );
118        assert_eq!(
119            output.len(),
120            self.output_len(),
121            "Output length {} doesn't match expected length {}",
122            output.len(),
123            self.output_len()
124        );
125
126        // Convert input to complex for full FFT
127        let mut buffer: Vec<Complex<f64>> = input.iter().map(|&x| Complex::new(x, 0.0)).collect();
128
129        // Get FFT plan and process
130        let mut planner = self.planner.lock().unwrap();
131        let fft = planner.plan_fft_forward(self.length);
132        drop(planner); // Release lock before processing
133
134        fft.process(&mut buffer);
135
136        // Copy first n/2 + 1 elements to output (real FFT property)
137        output.copy_from_slice(&buffer[..self.output_len()]);
138    }
139
140    fn len(&self) -> usize {
141        self.length
142    }
143}
144
145/// Inverse real FFT plan implementation for f64
146struct InverseRealFftPlanF64 {
147    length: usize,
148    planner: Arc<Mutex<rustfft::FftPlanner<f64>>>,
149}
150
151impl InverseRealFftPlanF64 {
152    fn new(length: usize, planner: Arc<Mutex<rustfft::FftPlanner<f64>>>) -> Self {
153        Self { length, planner }
154    }
155}
156
157impl ComplexToReal<f64> for InverseRealFftPlanF64 {
158    fn process(&self, input: &[Complex<f64>], output: &mut [f64]) {
159        // Validate input/output lengths
160        assert_eq!(
161            input.len(),
162            self.input_len(),
163            "Input length {} doesn't match expected length {}",
164            input.len(),
165            self.input_len()
166        );
167        assert_eq!(
168            output.len(),
169            self.length,
170            "Output length {} doesn't match plan length {}",
171            output.len(),
172            self.length
173        );
174
175        // Reconstruct full spectrum using Hermitian symmetry
176        let mut buffer: Vec<Complex<f64>> = Vec::with_capacity(self.length);
177        buffer.extend_from_slice(input);
178
179        // Add conjugate symmetric part
180        let start_idx = if self.length.is_multiple_of(2) {
181            input.len() - 1
182        } else {
183            input.len()
184        };
185
186        for i in (1..start_idx).rev() {
187            buffer.push(input[i].conj());
188        }
189
190        // Pad to full length if needed
191        while buffer.len() < self.length {
192            buffer.push(Complex::new(0.0, 0.0));
193        }
194
195        // Get inverse FFT plan and process
196        let mut planner = self.planner.lock().unwrap();
197        let ifft = planner.plan_fft_inverse(self.length);
198        drop(planner); // Release lock before processing
199
200        ifft.process(&mut buffer);
201
202        // Extract real parts and normalize
203        let scale = 1.0 / self.length as f64;
204        for (i, &val) in buffer.iter().enumerate() {
205            output[i] = val.re * scale;
206        }
207    }
208
209    fn len(&self) -> usize {
210        self.length
211    }
212}
213
214/// Real FFT plan implementation for f32
215struct RealFftPlanF32 {
216    length: usize,
217    planner: Arc<Mutex<rustfft::FftPlanner<f32>>>,
218}
219
220impl RealFftPlanF32 {
221    fn new(length: usize, planner: Arc<Mutex<rustfft::FftPlanner<f32>>>) -> Self {
222        Self { length, planner }
223    }
224}
225
226impl RealToComplex<f32> for RealFftPlanF32 {
227    fn process(&self, input: &[f32], output: &mut [Complex<f32>]) {
228        // Validate input/output lengths
229        assert_eq!(
230            input.len(),
231            self.length,
232            "Input length {} doesn't match plan length {}",
233            input.len(),
234            self.length
235        );
236        assert_eq!(
237            output.len(),
238            self.output_len(),
239            "Output length {} doesn't match expected length {}",
240            output.len(),
241            self.output_len()
242        );
243
244        // Convert input to complex for full FFT
245        let mut buffer: Vec<Complex<f32>> = input.iter().map(|&x| Complex::new(x, 0.0)).collect();
246
247        // Get FFT plan and process
248        let mut planner = self.planner.lock().unwrap();
249        let fft = planner.plan_fft_forward(self.length);
250        drop(planner); // Release lock before processing
251
252        fft.process(&mut buffer);
253
254        // Copy first n/2 + 1 elements to output (real FFT property)
255        output.copy_from_slice(&buffer[..self.output_len()]);
256    }
257
258    fn len(&self) -> usize {
259        self.length
260    }
261}
262
263/// Inverse real FFT plan implementation for f32
264struct InverseRealFftPlanF32 {
265    length: usize,
266    planner: Arc<Mutex<rustfft::FftPlanner<f32>>>,
267}
268
269impl InverseRealFftPlanF32 {
270    fn new(length: usize, planner: Arc<Mutex<rustfft::FftPlanner<f32>>>) -> Self {
271        Self { length, planner }
272    }
273}
274
275impl ComplexToReal<f32> for InverseRealFftPlanF32 {
276    fn process(&self, input: &[Complex<f32>], output: &mut [f32]) {
277        // Validate input/output lengths
278        assert_eq!(
279            input.len(),
280            self.input_len(),
281            "Input length {} doesn't match expected length {}",
282            input.len(),
283            self.input_len()
284        );
285        assert_eq!(
286            output.len(),
287            self.length,
288            "Output length {} doesn't match plan length {}",
289            output.len(),
290            self.length
291        );
292
293        // Reconstruct full spectrum using Hermitian symmetry
294        let mut buffer: Vec<Complex<f32>> = Vec::with_capacity(self.length);
295        buffer.extend_from_slice(input);
296
297        // Add conjugate symmetric part
298        let start_idx = if self.length.is_multiple_of(2) {
299            input.len() - 1
300        } else {
301            input.len()
302        };
303
304        for i in (1..start_idx).rev() {
305            buffer.push(input[i].conj());
306        }
307
308        // Pad to full length if needed
309        while buffer.len() < self.length {
310            buffer.push(Complex::new(0.0, 0.0));
311        }
312
313        // Get inverse FFT plan and process
314        let mut planner = self.planner.lock().unwrap();
315        let ifft = planner.plan_fft_inverse(self.length);
316        drop(planner); // Release lock before processing
317
318        ifft.process(&mut buffer);
319
320        // Extract real parts and normalize
321        let scale = 1.0 / self.length as f32;
322        for (i, &val) in buffer.iter().enumerate() {
323            output[i] = val.re * scale;
324        }
325    }
326
327    fn len(&self) -> usize {
328        self.length
329    }
330}
331
332/// Real FFT planner for creating and managing FFT plans
333///
334/// This planner creates reusable FFT plans optimized for real-valued input/output.
335/// Plans are thread-safe and can be shared across threads using Arc.
336///
337/// # Type Parameters
338///
339/// * `T` - Float type (f32 or f64)
340///
341/// # Examples
342///
343/// ```
344/// use scirs2_fft::real_planner::RealFftPlanner;
345///
346/// let mut planner = RealFftPlanner::<f64>::new();
347/// let forward = planner.plan_fft_forward(1024);
348/// let inverse = planner.plan_fft_inverse(1024);
349/// ```
350pub struct RealFftPlanner<T: Float> {
351    _phantom: std::marker::PhantomData<T>,
352}
353
354impl RealFftPlanner<f64> {
355    /// Create a new planner for f64 precision
356    pub fn new() -> Self {
357        Self {
358            _phantom: std::marker::PhantomData,
359        }
360    }
361
362    /// Create a forward FFT plan for real-to-complex transformation
363    ///
364    /// # Arguments
365    ///
366    /// * `length` - Length of the input real-valued array
367    ///
368    /// # Returns
369    ///
370    /// Arc-wrapped trait object implementing RealToComplex
371    pub fn plan_fft_forward(&mut self, length: usize) -> Arc<dyn RealToComplex<f64>> {
372        let planner = Arc::new(Mutex::new(rustfft::FftPlanner::<f64>::new()));
373        Arc::new(RealFftPlanF64::new(length, planner))
374    }
375
376    /// Create an inverse FFT plan for complex-to-real transformation
377    ///
378    /// # Arguments
379    ///
380    /// * `length` - Length of the output real-valued array
381    ///
382    /// # Returns
383    ///
384    /// Arc-wrapped trait object implementing ComplexToReal
385    pub fn plan_fft_inverse(&mut self, length: usize) -> Arc<dyn ComplexToReal<f64>> {
386        let planner = Arc::new(Mutex::new(rustfft::FftPlanner::<f64>::new()));
387        Arc::new(InverseRealFftPlanF64::new(length, planner))
388    }
389}
390
391impl Default for RealFftPlanner<f64> {
392    fn default() -> Self {
393        Self::new()
394    }
395}
396
397impl RealFftPlanner<f32> {
398    /// Create a new planner for f32 precision
399    pub fn new() -> Self {
400        Self {
401            _phantom: std::marker::PhantomData,
402        }
403    }
404
405    /// Create a forward FFT plan for real-to-complex transformation
406    ///
407    /// # Arguments
408    ///
409    /// * `length` - Length of the input real-valued array
410    ///
411    /// # Returns
412    ///
413    /// Arc-wrapped trait object implementing RealToComplex
414    pub fn plan_fft_forward(&mut self, length: usize) -> Arc<dyn RealToComplex<f32>> {
415        let planner = Arc::new(Mutex::new(rustfft::FftPlanner::<f32>::new()));
416        Arc::new(RealFftPlanF32::new(length, planner))
417    }
418
419    /// Create an inverse FFT plan for complex-to-real transformation
420    ///
421    /// # Arguments
422    ///
423    /// * `length` - Length of the output real-valued array
424    ///
425    /// # Returns
426    ///
427    /// Arc-wrapped trait object implementing ComplexToReal
428    pub fn plan_fft_inverse(&mut self, length: usize) -> Arc<dyn ComplexToReal<f32>> {
429        let planner = Arc::new(Mutex::new(rustfft::FftPlanner::<f32>::new()));
430        Arc::new(InverseRealFftPlanF32::new(length, planner))
431    }
432}
433
434impl Default for RealFftPlanner<f32> {
435    fn default() -> Self {
436        Self::new()
437    }
438}
439
440#[cfg(test)]
441mod tests {
442    use super::*;
443    use scirs2_core::numeric::Complex64;
444    use std::f64::consts::PI;
445
446    #[test]
447    fn test_real_fft_planner_f64() {
448        let mut planner = RealFftPlanner::<f64>::new();
449        let forward = planner.plan_fft_forward(8);
450        let inverse = planner.plan_fft_inverse(8);
451
452        // Test input
453        let input = vec![1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0];
454        let mut spectrum = vec![Complex64::new(0.0, 0.0); 5]; // 8/2 + 1 = 5
455
456        // Forward transform
457        forward.process(&input, &mut spectrum);
458
459        // Check DC component
460        let sum: f64 = input.iter().sum();
461        assert!((spectrum[0].re - sum).abs() < 1e-10);
462        assert!(spectrum[0].im.abs() < 1e-10);
463
464        // Inverse transform
465        let mut recovered = vec![0.0; 8];
466        inverse.process(&spectrum, &mut recovered);
467
468        // Check round-trip accuracy
469        for (i, (&orig, &recov)) in input.iter().zip(recovered.iter()).enumerate() {
470            assert!(
471                (orig - recov).abs() < 1e-10,
472                "Mismatch at index {}: {} vs {}",
473                i,
474                orig,
475                recov
476            );
477        }
478    }
479
480    #[test]
481    fn test_real_fft_planner_f32() {
482        let mut planner = RealFftPlanner::<f32>::new();
483        let forward = planner.plan_fft_forward(8);
484        let inverse = planner.plan_fft_inverse(8);
485
486        // Test input
487        let input = vec![1.0f32, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0];
488        let mut spectrum = vec![Complex::new(0.0f32, 0.0); 5]; // 8/2 + 1 = 5
489
490        // Forward transform
491        forward.process(&input, &mut spectrum);
492
493        // Inverse transform
494        let mut recovered = vec![0.0f32; 8];
495        inverse.process(&spectrum, &mut recovered);
496
497        // Check round-trip accuracy (lower precision for f32)
498        for (i, (&orig, &recov)) in input.iter().zip(recovered.iter()).enumerate() {
499            assert!(
500                (orig - recov).abs() < 1e-5,
501                "Mismatch at index {}: {} vs {}",
502                i,
503                orig,
504                recov
505            );
506        }
507    }
508
509    #[test]
510    fn test_sine_wave_fft() {
511        let mut planner = RealFftPlanner::<f64>::new();
512        let length = 128;
513        let forward = planner.plan_fft_forward(length);
514
515        // Generate sine wave at frequency bin 5
516        let freq_index = 5;
517        let input: Vec<f64> = (0..length)
518            .map(|i| (2.0 * PI * freq_index as f64 * i as f64 / length as f64).sin())
519            .collect();
520
521        let mut spectrum = vec![Complex64::new(0.0, 0.0); length / 2 + 1];
522        forward.process(&input, &mut spectrum);
523
524        // Check that energy is concentrated at the expected frequency
525        let magnitudes: Vec<f64> = spectrum.iter().map(|c| c.norm()).collect();
526
527        // Find peak
528        let (peak_idx, &peak_mag) = magnitudes
529            .iter()
530            .enumerate()
531            .max_by(|(_, a), (_, b)| a.partial_cmp(b).unwrap())
532            .unwrap();
533
534        assert_eq!(
535            peak_idx, freq_index,
536            "Peak should be at frequency index {}",
537            freq_index
538        );
539        assert!(peak_mag > length as f64 / 4.0, "Peak magnitude too small");
540    }
541
542    #[test]
543    fn test_plan_properties() {
544        let mut planner = RealFftPlanner::<f64>::new();
545        let forward = planner.plan_fft_forward(1024);
546
547        assert_eq!(forward.len(), 1024);
548        assert_eq!(forward.output_len(), 513); // 1024/2 + 1
549        assert!(!forward.is_empty());
550    }
551
552    #[test]
553    fn test_voirs_usage_pattern() {
554        // This test demonstrates the VoiRS usage pattern with Arc<dyn Trait>
555        struct AudioProcessor {
556            forward: Arc<dyn RealToComplex<f64>>,
557            backward: Arc<dyn ComplexToReal<f64>>,
558        }
559
560        impl AudioProcessor {
561            fn new(size: usize) -> Self {
562                let mut planner = RealFftPlanner::<f64>::new();
563                Self {
564                    forward: planner.plan_fft_forward(size),
565                    backward: planner.plan_fft_inverse(size),
566                }
567            }
568
569            fn process(&self, input: &[f64]) -> Vec<f64> {
570                let mut spectrum = vec![Complex64::new(0.0, 0.0); self.forward.output_len()];
571                self.forward.process(input, &mut spectrum);
572
573                let mut output = vec![0.0; self.backward.len()];
574                self.backward.process(&spectrum, &mut output);
575
576                output
577            }
578        }
579
580        let processor = AudioProcessor::new(16);
581        let input = vec![
582            1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 11.0, 12.0, 13.0, 14.0, 15.0, 16.0,
583        ];
584        let output = processor.process(&input);
585
586        // Verify round-trip
587        for (i, (&orig, &recov)) in input.iter().zip(output.iter()).enumerate() {
588            assert!((orig - recov).abs() < 1e-10, "Mismatch at {}", i);
589        }
590    }
591}
scirs2_fft/real_planner.rs

scirs2_fft/
real_planner.rs
