fft_convolver/
utilities.rs

1use realfft::{FftNum, num_complex::Complex};
2
3pub fn next_power_of_2(value: usize) -> usize {
4    let mut new_value = 1;
5
6    while new_value < value {
7        new_value *= 2
8    }
9
10    new_value
11}
12
13pub fn complex_size(size: usize) -> usize {
14    (size / 2) + 1
15}
16
17pub fn copy_and_pad<F: FftNum>(dst: &mut [F], src: &[F], src_size: usize) {
18    assert!(dst.len() >= src_size);
19    dst[0..src_size].clone_from_slice(&src[0..src_size]);
20    dst[src_size..]
21        .iter_mut()
22        .for_each(|value| *value = F::zero());
23}
24
25pub fn complex_multiply_accumulate<F: FftNum>(
26    result: &mut [Complex<F>],
27    a: &[Complex<F>],
28    b: &[Complex<F>],
29) {
30    assert_eq!(result.len(), a.len());
31    assert_eq!(result.len(), b.len());
32    let len = result.len();
33    let end4 = 4 * (len / 4);
34    #[allow(clippy::identity_op)]
35    for i in (0..end4).step_by(4) {
36        result[i + 0].re =
37            result[i + 0].re + (a[i + 0].re * b[i + 0].re - a[i + 0].im * b[i + 0].im);
38        result[i + 1].re =
39            result[i + 1].re + (a[i + 1].re * b[i + 1].re - a[i + 1].im * b[i + 1].im);
40        result[i + 2].re =
41            result[i + 2].re + (a[i + 2].re * b[i + 2].re - a[i + 2].im * b[i + 2].im);
42        result[i + 3].re =
43            result[i + 3].re + (a[i + 3].re * b[i + 3].re - a[i + 3].im * b[i + 3].im);
44        result[i + 0].im =
45            result[i + 0].im + (a[i + 0].re * b[i + 0].im + a[i + 0].im * b[i + 0].re);
46        result[i + 1].im =
47            result[i + 1].im + (a[i + 1].re * b[i + 1].im + a[i + 1].im * b[i + 1].re);
48        result[i + 2].im =
49            result[i + 2].im + (a[i + 2].re * b[i + 2].im + a[i + 2].im * b[i + 2].re);
50        result[i + 3].im =
51            result[i + 3].im + (a[i + 3].re * b[i + 3].im + a[i + 3].im * b[i + 3].re);
52    }
53    for i in end4..len {
54        result[i].re = result[i].re + (a[i].re * b[i].re - a[i].im * b[i].im);
55        result[i].im = result[i].im + (a[i].re * b[i].im + a[i].im * b[i].re);
56    }
57}
58
59/// Computes the optimal tail block size for a two-stage convolver
60/// using García's formula (AES 2002).
61///
62/// For a given head block size N and impulse response length T, the optimal
63/// tail block size B minimizes the total computational cost of the double-FDL
64/// convolution scheme.
65pub fn compute_tail_block_size(head_block_size: usize, ir_len: usize) -> usize {
66    const K: f64 = 1.5;
67    let n = head_block_size as f64;
68    let t = ir_len as f64;
69    let kn = K * n / (2.0 * 2.0_f64.ln());
70    let b = -kn + (kn * kn + t * n).sqrt();
71    let b = next_power_of_2(b.ceil().max(1.0) as usize);
72    b.max(2 * head_block_size)
73}
74
75pub fn sum<F: FftNum>(result: &mut [F], a: &[F], b: &[F]) {
76    assert_eq!(result.len(), a.len());
77    assert_eq!(result.len(), b.len());
78    let len = result.len();
79    let end4 = 4 * (len / 4);
80    #[allow(clippy::identity_op)]
81    for i in (0..end4).step_by(4) {
82        result[i + 0] = a[i + 0] + b[i + 0];
83        result[i + 1] = a[i + 1] + b[i + 1];
84        result[i + 2] = a[i + 2] + b[i + 2];
85        result[i + 3] = a[i + 3] + b[i + 3];
86    }
87    for i in end4..len {
88        result[i] = a[i] + b[i];
89    }
90}
91
92#[cfg(test)]
93mod tests {
94    use realfft::num_complex::Complex;
95
96    use crate::utilities::complex_multiply_accumulate;
97    use crate::utilities::compute_tail_block_size;
98    use crate::utilities::copy_and_pad;
99    use crate::utilities::next_power_of_2;
100    use crate::utilities::sum;
101
102    #[test]
103    fn next_power_of_2_test() {
104        assert_eq!(128, next_power_of_2(122));
105        assert_eq!(1024, next_power_of_2(1000));
106        assert_eq!(1024, next_power_of_2(1024));
107        assert_eq!(1, next_power_of_2(1));
108    }
109
110    #[test]
111    fn copy_and_pad_test() {
112        let mut dst: Vec<f32> = vec![1.; 10];
113        let src: Vec<f32> = vec![2., 3., 4., 5., 6.];
114        copy_and_pad(&mut dst, &src, src.len());
115
116        assert_eq!(dst[0], 2.);
117        assert_eq!(dst[1], 3.);
118        assert_eq!(dst[2], 4.);
119        assert_eq!(dst[3], 5.);
120        assert_eq!(dst[4], 6.);
121        for num in &dst[5..] {
122            assert_eq!(*num, 0.);
123        }
124    }
125
126    #[test]
127    fn complex_mulitply_accumulate_test() {
128        let mut result: Vec<Complex<f32>> = vec![Complex::new(0., 0.); 10];
129
130        let a: Vec<Complex<f32>> = vec![
131            Complex::new(0., 9.),
132            Complex::new(1., 8.),
133            Complex::new(2., 7.),
134            Complex::new(3., 6.),
135            Complex::new(4., 5.),
136            Complex::new(5., 4.),
137            Complex::new(6., 3.),
138            Complex::new(7., 2.),
139            Complex::new(8., 1.),
140            Complex::new(9., 0.),
141        ];
142
143        let b: Vec<Complex<f32>> = vec![
144            Complex::new(9., 0.),
145            Complex::new(8., 1.),
146            Complex::new(7., 2.),
147            Complex::new(6., 3.),
148            Complex::new(5., 4.),
149            Complex::new(4., 5.),
150            Complex::new(3., 6.),
151            Complex::new(2., 7.),
152            Complex::new(1., 8.),
153            Complex::new(0., 9.),
154        ];
155        complex_multiply_accumulate(&mut result, &a, &b);
156
157        for num in &result {
158            assert_eq!(num.re, 0.);
159        }
160
161        assert_eq!(result[0].im, 81.);
162        assert_eq!(result[1].im, 65.);
163        assert_eq!(result[2].im, 53.);
164        assert_eq!(result[3].im, 45.);
165        assert_eq!(result[4].im, 41.);
166        assert_eq!(result[5].im, 41.);
167        assert_eq!(result[6].im, 45.);
168        assert_eq!(result[7].im, 53.);
169        assert_eq!(result[8].im, 65.);
170        assert_eq!(result[9].im, 81.);
171    }
172
173    #[test]
174    fn compute_tail_block_size_test() {
175        // For head=64, ir=131072 (3s reverb at 44.1kHz), should get a reasonable power of 2
176        let b = compute_tail_block_size(64, 131072);
177        assert!(b.is_power_of_two());
178        assert!(b >= 128); // at least 2 * head
179        assert!(b <= 8192); // reasonable upper bound
180
181        // Must always be at least 2 * head_block_size
182        let b = compute_tail_block_size(64, 100);
183        assert!(b >= 128);
184
185        // Small head, long IR
186        let b = compute_tail_block_size(32, 100000);
187        assert!(b.is_power_of_two());
188        assert!(b >= 64);
189    }
190
191    #[test]
192    fn sum_test() {
193        let mut result = vec![0.; 10];
194        let a = vec![0., 1., 2., 3., 4., 5., 6., 7., 8., 9.];
195        let b = vec![0., 6., 3., 1., 5., 3., 5., 1., 4., 0.];
196
197        sum(&mut result, &a, &b);
198
199        assert_eq!(result[0], 0.);
200        assert_eq!(result[1], 7.);
201        assert_eq!(result[2], 5.);
202        assert_eq!(result[3], 4.);
203        assert_eq!(result[4], 9.);
204        assert_eq!(result[5], 8.);
205        assert_eq!(result[6], 11.);
206        assert_eq!(result[7], 8.);
207        assert_eq!(result[8], 12.);
208        assert_eq!(result[9], 9.);
209    }
210}
fft_convolver/utilities.rs

fft_convolver/
utilities.rs
