rill_fft/
real_fft.rs

1// rill-fft/src/real_fft.rs
2//! Real-valued FFT using a complex FFT with packing/unpacking.
3//!
4//! Transforms `N` real samples into `N/2 + 1` complex frequency bins.
5//! The inverse transform reconstructs `N` real samples from the complex bins.
6
7use num_complex::Complex;
8use rill_core::Transcendental;
9
10use crate::complex_fft::ComplexFft;
11
12/// Real-valued FFT.
13///
14/// Uses a half-size complex FFT internally via the two-for-one method.
15/// Transforms `N` real samples into `N/2 + 1` complex bins (only the
16/// non-redundant half of the spectrum). The Nyquist bin (`N/2`) and
17/// DC bin (`0`) are purely real.
18///
19/// # Panics
20///
21/// Panics if `size` is not a power of two or less than 4.
22pub struct RealFft<T: Transcendental> {
23    size: usize,
24    half_size: usize,
25    complex_fft: ComplexFft<T>,
26    scratch: Vec<Complex<T>>,
27}
28
29impl<T: Transcendental> RealFft<T> {
30    /// Create a new real FFT for the given size.
31    ///
32    /// # Panics
33    ///
34    /// Panics if `size` is not a power of two or is less than 4.
35    pub fn new(size: usize) -> Self {
36        assert!(
37            size.is_power_of_two(),
38            "FFT size must be a power of two, got {size}"
39        );
40        assert!(size >= 4, "FFT size must be at least 4, got {size}");
41
42        let half_size = size / 2;
43        let complex_fft = ComplexFft::new(half_size);
44        let scratch = vec![Complex::new(T::ZERO, T::ZERO); half_size];
45
46        Self {
47            size,
48            half_size,
49            complex_fft,
50            scratch,
51        }
52    }
53
54    /// Returns the FFT size (number of real input samples).
55    pub fn size(&self) -> usize {
56        self.size
57    }
58
59    /// Forward real FFT.
60    ///
61    /// Transforms `input` (N real samples) into `output` (N/2 + 1 complex bins).
62    ///
63    /// # Panics
64    ///
65    /// Panics if `input.len() != self.size()` or `output.len() != self.half_size + 1`.
66    pub fn forward(&mut self, input: &[T], output: &mut [Complex<T>]) {
67        assert_eq!(
68            input.len(),
69            self.size,
70            "input length ({}) must match FFT size ({})",
71            input.len(),
72            self.size
73        );
74        assert_eq!(
75            output.len(),
76            self.half_size + 1,
77            "output length ({}) must be half_size + 1 ({})",
78            output.len(),
79            self.half_size + 1
80        );
81
82        for i in 0..self.half_size {
83            self.scratch[i] = Complex::new(input[2 * i], input[2 * i + 1]);
84        }
85
86        self.complex_fft.forward(&mut self.scratch);
87
88        let z = &self.scratch;
89        let z0 = z[0];
90        output[0] = Complex::new(z0.re + z0.im, T::ZERO);
91        output[self.half_size] = Complex::new(z0.re - z0.im, T::ZERO);
92
93        let twopi = T::from_f64(2.0 * std::f64::consts::PI);
94        let size_t = T::from_usize(self.size);
95        for k in 1..self.half_size {
96            let theta = twopi * T::from_usize(k) / size_t;
97            let w_cos = theta.cos();
98            let w_sin = theta.sin();
99
100            let a = z[k];
101            let b = Complex::new(z[self.half_size - k].re, -z[self.half_size - k].im);
102
103            let c_re = a.re + b.re;
104            let c_im = a.im + b.im;
105            let d_re = a.re - b.re;
106            let d_im = a.im - b.im;
107
108            let half = T::from_f64(0.5);
109            output[k].re = half * (c_re - w_sin * d_re + w_cos * d_im);
110            output[k].im = half * (c_im - w_sin * d_im - w_cos * d_re);
111        }
112    }
113
114    /// Inverse real FFT.
115    ///
116    /// Reconstructs `output` (N real samples) from `input` (N/2 + 1 complex bins).
117    /// This is the exact inverse of `forward()`.
118    ///
119    /// # Panics
120    ///
121    /// Panics if `input.len() != self.half_size + 1` or `output.len() != self.size()`.
122    pub fn inverse(&mut self, input: &[Complex<T>], output: &mut [T]) {
123        assert_eq!(
124            input.len(),
125            self.half_size + 1,
126            "input length ({}) must be half_size + 1 ({})",
127            input.len(),
128            self.half_size + 1
129        );
130        assert_eq!(
131            output.len(),
132            self.size,
133            "output length ({}) must match FFT size ({})",
134            output.len(),
135            self.size
136        );
137
138        let x = input;
139
140        let half = T::from_f64(0.5);
141        self.scratch[0] = Complex::new(
142            half * (x[0].re + x[self.half_size].re),
143            half * (x[0].re - x[self.half_size].re),
144        );
145
146        let twopi = T::from_f64(2.0 * std::f64::consts::PI);
147        let size_t = T::from_usize(self.size);
148        for k in 1..self.half_size {
149            let theta = twopi * T::from_usize(k) / size_t;
150            let w_cos = theta.cos();
151            let w_sin = theta.sin();
152
153            let a = x[k];
154            let b = Complex::new(x[self.half_size - k].re, -x[self.half_size - k].im);
155
156            let c_re = a.re + b.re;
157            let c_im = a.im + b.im;
158            let d_re = a.re - b.re;
159            let d_im = a.im - b.im;
160
161            let half = T::from_f64(0.5);
162            self.scratch[k].re = half * (c_re - w_sin * d_re - w_cos * d_im);
163            self.scratch[k].im = half * (c_im - w_sin * d_im + w_cos * d_re);
164        }
165
166        self.complex_fft.inverse(&mut self.scratch);
167
168        for i in 0..self.half_size {
169            output[2 * i] = self.scratch[i].re;
170            output[2 * i + 1] = self.scratch[i].im;
171        }
172    }
173}
174
175#[cfg(test)]
176mod tests {
177    use super::*;
178
179    #[test]
180    fn test_real_fft_size_4_roundtrip() {
181        let mut fft = RealFft::<f32>::new(4);
182        let input = [1.0f32, 2.0, 3.0, 4.0];
183        let mut spectrum = vec![Complex::new(0.0, 0.0); 3];
184        let mut output = [0.0f32; 4];
185
186        fft.forward(&input, &mut spectrum);
187        fft.inverse(&spectrum, &mut output);
188
189        for (a, b) in input.iter().zip(output.iter()) {
190            assert!((a - b).abs() < 1e-3, "expected {a}, got {b}");
191        }
192    }
193
194    #[test]
195    fn test_real_fft_size_8_roundtrip() {
196        let mut fft = RealFft::<f32>::new(8);
197        let input: Vec<f32> = (0..8).map(|i| (i as f32 * 0.7).sin()).collect();
198        let mut spectrum = vec![Complex::new(0.0, 0.0); 5];
199        let mut output = vec![0.0f32; 8];
200
201        fft.forward(&input, &mut spectrum);
202        fft.inverse(&spectrum, &mut output);
203
204        for (a, b) in input.iter().zip(output.iter()) {
205            assert!((a - b).abs() < 1e-3, "expected {a}, got {b}");
206        }
207    }
208
209    #[test]
210    fn test_real_fft_size_1024_roundtrip() {
211        let mut fft = RealFft::<f32>::new(1024);
212        let input: Vec<f32> = (0..1024)
213            .map(|i| {
214                let x = i as f32 * 0.05;
215                x.sin() + 0.5 * (x * 2.3).cos()
216            })
217            .collect();
218        let mut spectrum = vec![Complex::new(0.0, 0.0); 513];
219        let mut output = vec![0.0f32; 1024];
220
221        fft.forward(&input, &mut spectrum);
222        fft.inverse(&spectrum, &mut output);
223
224        for (a, b) in input.iter().zip(output.iter()) {
225            assert!((a - b).abs() < 5e-4, "at index: expected {a}, got {b}");
226        }
227    }
228
229    #[test]
230    fn test_real_fft_size_16_roundtrip() {
231        let mut fft = RealFft::<f32>::new(16);
232        let input: Vec<f32> = (0..16).map(|i| i as f32).collect();
233        let mut spectrum = vec![Complex::new(0.0, 0.0); 9];
234        let mut output = vec![0.0f32; 16];
235
236        fft.forward(&input, &mut spectrum);
237        fft.inverse(&spectrum, &mut output);
238
239        for (a, b) in input.iter().zip(output.iter()) {
240            assert!((a - b).abs() < 5e-4, "expected {a}, got {b}");
241        }
242    }
243
244    #[test]
245    fn test_real_fft_f64_roundtrip() {
246        let mut fft = RealFft::<f64>::new(8);
247        let input: Vec<f64> = (0..8).map(|i| (i as f64 * 0.5).sin()).collect();
248        let mut spectrum = vec![Complex::new(0.0, 0.0); 5];
249        let mut output = vec![0.0f64; 8];
250
251        fft.forward(&input, &mut spectrum);
252        fft.inverse(&spectrum, &mut output);
253
254        for (a, b) in input.iter().zip(output.iter()) {
255            assert!((a - b).abs() < 1e-10, "expected {a}, got {b}");
256        }
257    }
258
259    #[test]
260    fn test_real_fft_dc_input() {
261        let mut fft = RealFft::<f32>::new(16);
262        let input = [2.5f32; 16];
263        let mut spectrum = vec![Complex::new(0.0, 0.0); 9];
264
265        fft.forward(&input, &mut spectrum);
266
267        assert!((spectrum[0].re - 40.0).abs() < 1e-3);
268        assert!(spectrum[0].im.abs() < 1e-3);
269    }
270
271    #[test]
272    fn test_real_fft_nyquist_bin() {
273        let mut fft = RealFft::<f32>::new(8);
274        let input = [1.0f32, -1.0, 1.0, -1.0, 1.0, -1.0, 1.0, -1.0];
275        let mut spectrum = vec![Complex::new(0.0, 0.0); 5];
276
277        fft.forward(&input, &mut spectrum);
278
279        assert!(spectrum[4].re.abs() > 0.1);
280    }
281}
rill_fft/real_fft.rs

rill_fft/
real_fft.rs
