rustfft/algorithm/
bluesteins_algorithm.rs

1use std::sync::Arc;
2
3use num_complex::Complex;
4use num_traits::Zero;
5
6use crate::array_utils;
7use crate::common::{fft_error_inplace, fft_error_outofplace};
8use crate::{common::FftNum, twiddles, FftDirection};
9use crate::{Direction, Fft, Length};
10
11/// Implementation of Bluestein's Algorithm
12///
13/// This algorithm computes an arbitrary-sized FFT in O(nlogn) time. It does this by converting this size-N FFT into a
14/// size-M FFT where M >= 2N - 1.
15///
16/// The choice of M is very important for the performance of Bluestein's Algorithm. The most obvious choice is the next-largest
17/// power of two -- but if there's a smaller/faster FFT size that satisfies the `>= 2N - 1` requirement, that will significantly
18/// improve this algorithm's overall performance.
19///
20/// ~~~
21/// // Computes a forward FFT of size 1201, using Bluestein's Algorithm
22/// use rustfft::algorithm::BluesteinsAlgorithm;
23/// use rustfft::{Fft, FftPlanner};
24/// use rustfft::num_complex::Complex;
25///
26/// let mut buffer = vec![Complex{ re: 0.0f32, im: 0.0f32 }; 1201];
27///
28/// // We need to find an inner FFT whose size is greater than 1201*2 - 1.
29/// // The size 2401 (7^4) satisfies this requirement, while also being relatively fast.
30/// let mut planner = FftPlanner::new();
31/// let inner_fft = planner.plan_fft_forward(2401);
32///
33/// let fft = BluesteinsAlgorithm::new(1201, inner_fft);
34/// fft.process(&mut buffer);
35/// ~~~
36///
37/// Bluesteins's Algorithm is relatively expensive compared to other FFT algorithms. Benchmarking shows that it is up to
38/// an order of magnitude slower than similar composite sizes. In the example size above of 1201, benchmarking shows
39/// that it takes 5x more time to compute than computing a FFT of size 1200 via a step of MixedRadix.
40
41pub struct BluesteinsAlgorithm<T> {
42    inner_fft: Arc<dyn Fft<T>>,
43
44    inner_fft_multiplier: Box<[Complex<T>]>,
45    twiddles: Box<[Complex<T>]>,
46
47    len: usize,
48    direction: FftDirection,
49}
50
51impl<T: FftNum> BluesteinsAlgorithm<T> {
52    /// Creates a FFT instance which will process inputs/outputs of size `len`. `inner_fft.len()` must be >= `len * 2 - 1`
53    ///
54    /// Note that this constructor is quite expensive to run; This algorithm must compute a FFT using `inner_fft` within the
55    /// constructor. This further underlines the fact that Bluesteins Algorithm is more expensive to run than other
56    /// FFT algorithms
57    ///
58    /// # Panics
59    /// Panics if `inner_fft.len() < len * 2 - 1`.
60    pub fn new(len: usize, inner_fft: Arc<dyn Fft<T>>) -> Self {
61        let inner_fft_len = inner_fft.len();
62        assert!(len * 2 - 1 <= inner_fft_len, "Bluestein's algorithm requires inner_fft.len() >= self.len() * 2 - 1. Expected >= {}, got {}", len * 2 - 1, inner_fft_len);
63
64        // when computing FFTs, we're going to run our inner multiply pairise by some precomputed data, then run an inverse inner FFT. We need to precompute that inner data here
65        let inner_fft_scale = T::one() / T::from_usize(inner_fft_len).unwrap();
66        let direction = inner_fft.fft_direction();
67
68        // Compute twiddle factors that we'll run our inner FFT on
69        let mut inner_fft_input = vec![Complex::zero(); inner_fft_len];
70        twiddles::fill_bluesteins_twiddles(
71            &mut inner_fft_input[..len],
72            direction.opposite_direction(),
73        );
74
75        // Scale the computed twiddles and copy them to the end of the array
76        inner_fft_input[0] = inner_fft_input[0] * inner_fft_scale;
77        for i in 1..len {
78            let twiddle = inner_fft_input[i] * inner_fft_scale;
79            inner_fft_input[i] = twiddle;
80            inner_fft_input[inner_fft_len - i] = twiddle;
81        }
82
83        //Compute the inner fft
84        let mut inner_fft_scratch = vec![Complex::zero(); inner_fft.get_inplace_scratch_len()];
85        inner_fft.process_with_scratch(&mut inner_fft_input, &mut inner_fft_scratch);
86
87        // also compute some more mundane twiddle factors to start and end with
88        let mut twiddles = vec![Complex::zero(); len];
89        twiddles::fill_bluesteins_twiddles(&mut twiddles, direction);
90
91        Self {
92            inner_fft: inner_fft,
93
94            inner_fft_multiplier: inner_fft_input.into_boxed_slice(),
95            twiddles: twiddles.into_boxed_slice(),
96
97            len,
98            direction,
99        }
100    }
101
102    fn perform_fft_inplace(&self, input: &mut [Complex<T>], scratch: &mut [Complex<T>]) {
103        let (inner_input, inner_scratch) = scratch.split_at_mut(self.inner_fft_multiplier.len());
104
105        // Copy the buffer into our inner FFT input. the buffer will only fill part of the FFT input, so zero fill the rest
106        for ((buffer_entry, inner_entry), twiddle) in input
107            .iter()
108            .zip(inner_input.iter_mut())
109            .zip(self.twiddles.iter())
110        {
111            *inner_entry = *buffer_entry * *twiddle;
112        }
113        for inner in (&mut inner_input[input.len()..]).iter_mut() {
114            *inner = Complex::zero();
115        }
116
117        // run our inner forward FFT
118        self.inner_fft
119            .process_with_scratch(inner_input, inner_scratch);
120
121        // Multiply our inner FFT output by our precomputed data. Then, conjugate the result to set up for an inverse FFT
122        for (inner, multiplier) in inner_input.iter_mut().zip(self.inner_fft_multiplier.iter()) {
123            *inner = (*inner * *multiplier).conj();
124        }
125
126        // inverse FFT. we're computing a forward but we're massaging it into an inverse by conjugating the inputs and outputs
127        self.inner_fft
128            .process_with_scratch(inner_input, inner_scratch);
129
130        // copy our data back to the buffer, applying twiddle factors again as we go. Also conjugate inner_input to complete the inverse FFT
131        for ((buffer_entry, inner_entry), twiddle) in input
132            .iter_mut()
133            .zip(inner_input.iter())
134            .zip(self.twiddles.iter())
135        {
136            *buffer_entry = inner_entry.conj() * twiddle;
137        }
138    }
139
140    fn perform_fft_out_of_place(
141        &self,
142        input: &mut [Complex<T>],
143        output: &mut [Complex<T>],
144        scratch: &mut [Complex<T>],
145    ) {
146        let (inner_input, inner_scratch) = scratch.split_at_mut(self.inner_fft_multiplier.len());
147
148        // Copy the buffer into our inner FFT input. the buffer will only fill part of the FFT input, so zero fill the rest
149        for ((buffer_entry, inner_entry), twiddle) in input
150            .iter()
151            .zip(inner_input.iter_mut())
152            .zip(self.twiddles.iter())
153        {
154            *inner_entry = *buffer_entry * *twiddle;
155        }
156        for inner in inner_input.iter_mut().skip(input.len()) {
157            *inner = Complex::zero();
158        }
159
160        // run our inner forward FFT
161        self.inner_fft
162            .process_with_scratch(inner_input, inner_scratch);
163
164        // Multiply our inner FFT output by our precomputed data. Then, conjugate the result to set up for an inverse FFT
165        for (inner, multiplier) in inner_input.iter_mut().zip(self.inner_fft_multiplier.iter()) {
166            *inner = (*inner * *multiplier).conj();
167        }
168
169        // inverse FFT. we're computing a forward but we're massaging it into an inverse by conjugating the inputs and outputs
170        self.inner_fft
171            .process_with_scratch(inner_input, inner_scratch);
172
173        // copy our data back to the buffer, applying twiddle factors again as we go. Also conjugate inner_input to complete the inverse FFT
174        for ((buffer_entry, inner_entry), twiddle) in output
175            .iter_mut()
176            .zip(inner_input.iter())
177            .zip(self.twiddles.iter())
178        {
179            *buffer_entry = inner_entry.conj() * twiddle;
180        }
181    }
182}
183boilerplate_fft!(
184    BluesteinsAlgorithm,
185    |this: &BluesteinsAlgorithm<_>| this.len, // FFT len
186    |this: &BluesteinsAlgorithm<_>| this.inner_fft_multiplier.len()
187        + this.inner_fft.get_inplace_scratch_len(), // in-place scratch len
188    |this: &BluesteinsAlgorithm<_>| this.inner_fft_multiplier.len()
189        + this.inner_fft.get_inplace_scratch_len()  // out of place scratch len
190);
191
192#[cfg(test)]
193mod unit_tests {
194    use super::*;
195    use crate::algorithm::Dft;
196    use crate::test_utils::check_fft_algorithm;
197    use std::sync::Arc;
198
199    #[test]
200    fn test_bluesteins_scalar() {
201        for &len in &[3, 5, 7, 11, 13] {
202            test_bluesteins_with_length(len, FftDirection::Forward);
203            test_bluesteins_with_length(len, FftDirection::Inverse);
204        }
205    }
206
207    fn test_bluesteins_with_length(len: usize, direction: FftDirection) {
208        let inner_fft = Arc::new(Dft::new(
209            (len * 2 - 1).checked_next_power_of_two().unwrap(),
210            direction,
211        ));
212        let fft = BluesteinsAlgorithm::new(len, inner_fft);
213
214        check_fft_algorithm::<f32>(&fft, len, direction);
215    }
216}
rustfft/algorithm/bluesteins_algorithm.rs

rustfft/algorithm/
bluesteins_algorithm.rs
