From d740d55ea5162dd4c600fd3d02c70f1f58df259e Mon Sep 17 00:00:00 2001 From: Luni-4 Date: Sun, 19 Apr 2020 01:12:00 +0200 Subject: [PATCH 1/2] Make fft library safe --- fft/src/dradb.rs | 646 +++++++++++++ fft/src/dradf.rs | 558 +++++++++++ fft/src/fftwrap.rs | 100 +- fft/src/lib.rs | 2 + fft/src/smallft.rs | 2157 +++++------------------------------------- fft/tests/fftwrap.rs | 76 +- fft/tests/smallft.rs | 9 +- 7 files changed, 1493 insertions(+), 2055 deletions(-) create mode 100644 fft/src/dradb.rs create mode 100644 fft/src/dradf.rs diff --git a/fft/src/dradb.rs b/fft/src/dradb.rs new file mode 100644 index 0000000..f2c4cb3 --- /dev/null +++ b/fft/src/dradb.rs @@ -0,0 +1,646 @@ +pub(crate) fn dradb2( + ido: i32, + l1: i32, + cc: &[f32], + ch: &mut [f32], + wa1: &[f32], +) { + let t0 = l1 * ido; + let mut t1 = 0 as i32; + let mut t2 = 0 as i32; + let mut t3 = (ido << 1 as i32) - 1 as i32; + + for _ in 0..l1 { + ch[t1 as usize] = cc[t2 as usize] + cc[t3 as usize + t2 as usize]; + ch[t1 as usize + t0 as usize] = + cc[t2 as usize] - cc[t3 as usize + t2 as usize]; + t1 += ido; + t2 = t1 << 1 as i32; + } + + if ido < 2 { + return; + } + + if ido != 2 { + t1 = 0 as i32; + t2 = 0 as i32; + for _ in 0..l1 { + t3 = t1; + let mut t4 = t2; + let mut t5 = t4 + (ido << 1 as i32); + let mut t6 = t0 + t1; + for i in (2..ido as usize).step_by(2) { + t3 += 2 as i32; + t4 += 2 as i32; + t5 -= 2 as i32; + t6 += 2 as i32; + ch[t3 as usize - 1] = + cc[t4 as usize - 1] + cc[t5 as usize - 1]; + let tr2 = cc[t4 as usize - 1] - cc[t5 as usize - 1]; + ch[t3 as usize] = cc[t4 as usize] - cc[t5 as usize]; + let ti2 = cc[t4 as usize] + cc[t5 as usize]; + ch[t6 as usize - 1] = wa1[i - 2] * tr2 - wa1[i - 1] * ti2; + ch[t6 as usize] = wa1[i - 2] * ti2 + wa1[i - 1] * tr2; + } + t1 += ido; + t2 = t1 << 1 as i32; + } + + if ido % 2 == 1 { + return; + } + } + + t1 = ido - 1 as i32; + t2 = ido - 1 as i32; + for _ in 0..l1 { + ch[t1 as usize] = cc[t2 as usize] + cc[t2 as usize]; + ch[t1 as usize + t0 as usize] = + -(cc[t2 as usize + 1] + cc[t2 as usize + 1]); + t1 += ido; + t2 += ido << 1 as i32; + } +} + +pub(crate) fn dradb3( + ido: i32, + l1: i32, + cc: &[f32], + ch: &mut [f32], + wa1: &[f32], + wa2: &[f32], +) { + const TAUR: f32 = -0.5; + const TAUI: f32 = 0.866_025_403_784_438_6; + + let t0 = l1 * ido; + let mut t1 = 0 as i32; + let t2 = t0 << 1 as i32; + let mut t3 = ido << 1 as i32; + let t4 = ido + (ido << 1 as i32); + let mut t5 = 0 as i32; + + for _ in 0..l1 { + let tr2 = cc[t3 as usize - 1] + cc[t3 as usize - 1]; + let cr2 = cc[t5 as usize] + TAUR * tr2; + ch[t1 as usize] = cc[t5 as usize] + tr2; + let ci3 = TAUI * (cc[t3 as usize] + cc[t3 as usize]); + ch[t1 as usize + t0 as usize] = cr2 - ci3; + ch[t1 as usize + t2 as usize] = cr2 + ci3; + t1 += ido; + t3 += t4; + t5 += t4; + } + + if ido == 1 { + return; + } + + t1 = 0 as i32; + t3 = ido << 1 as i32; + for _ in 0..l1 { + let mut t7 = t1 + (t1 << 1 as i32); + t5 = t7 + t3; + let mut t6 = t5; + let mut t8 = t1; + let mut t9 = t1 + t0; + let mut t10 = t9 + t0; + for i in (2..ido as usize).step_by(2) { + t5 += 2 as i32; + t6 -= 2 as i32; + t7 += 2 as i32; + t8 += 2 as i32; + t9 += 2 as i32; + t10 += 2 as i32; + let tr2 = cc[t5 as usize - 1] + cc[t6 as usize - 1]; + let cr2 = cc[t7 as usize - 1] + TAUR * tr2; + ch[t8 as usize - 1] = cc[t7 as usize - 1] + tr2; + let ti2 = cc[t5 as usize] - cc[t6 as usize]; + let ci2 = cc[t7 as usize] + TAUR * ti2; + ch[t8 as usize] = cc[t7 as usize] + ti2; + let cr3 = TAUI * (cc[t5 as usize - 1] - cc[t6 as usize - 1]); + let ci3 = TAUI * (cc[t5 as usize] + cc[t6 as usize]); + let dr2 = cr2 - ci3; + let dr3 = cr2 + ci3; + let di2 = ci2 + cr3; + let di3 = ci2 - cr3; + ch[t9 as usize - 1] = wa1[i - 2] * dr2 - wa1[i - 1] * di2; + ch[t9 as usize] = wa1[i - 2] * di2 + wa1[i - 1] * dr2; + ch[t10 as usize - 1] = wa2[i - 2] * dr3 - wa2[i - 1] * di3; + ch[t10 as usize] = wa2[i - 2] * di3 + wa2[i - 1] * dr3; + } + t1 += ido; + } +} + +pub(crate) fn dradb4( + ido: i32, + l1: i32, + cc: &[f32], + ch: &mut [f32], + wa1: &[f32], + wa2: &[f32], + wa3: &[f32], +) { + const SQRT2: f32 = std::f32::consts::SQRT_2; + + let t0 = l1 * ido; + let mut t1 = 0 as i32; + let mut t2 = ido << 2 as i32; + let mut t3 = 0 as i32; + let t6 = ido << 1 as i32; + + for _ in 0..l1 { + let mut t4 = t3 + t6; + let mut t5 = t1; + let tr3 = cc[t4 as usize - 1] + cc[t4 as usize - 1]; + let tr4 = cc[t4 as usize] + cc[t4 as usize]; + t4 += t6; + let tr1 = cc[t3 as usize] - cc[t4 as usize - 1]; + let tr2 = cc[t3 as usize] + cc[t4 as usize - 1]; + ch[t5 as usize] = tr2 + tr3; + t5 += t0; + ch[t5 as usize] = tr1 - tr4; + t5 += t0; + ch[t5 as usize] = tr2 - tr3; + t5 += t0; + ch[t5 as usize] = tr1 + tr4; + t1 += ido; + t3 += t2; + } + + if ido < 2 { + return; + } + + if ido != 2 { + t1 = 0 as i32; + for _ in 0..l1 { + t2 = t1 << 2 as i32; + t3 = t2 + t6; + let mut t4 = t3; + let mut t5 = t4 + t6; + let mut t7 = t1; + for i in (2..ido as usize).step_by(2) { + t2 += 2 as i32; + t3 += 2 as i32; + t4 -= 2 as i32; + t5 -= 2 as i32; + t7 += 2 as i32; + let ti1 = cc[t2 as usize] + cc[t5 as usize]; + let ti2 = cc[t2 as usize] - cc[t5 as usize]; + let ti3 = cc[t3 as usize] - cc[t4 as usize]; + let tr4 = cc[t3 as usize] + cc[t4 as usize]; + let tr1 = cc[t2 as usize - 1] - cc[t5 as usize - 1]; + let tr2 = cc[t2 as usize - 1] + cc[t5 as usize - 1]; + let ti4 = cc[t3 as usize - 1] - cc[t4 as usize - 1]; + let tr3 = cc[t3 as usize - 1] + cc[t4 as usize - 1]; + ch[t7 as usize - 1] = tr2 + tr3; + let cr3 = tr2 - tr3; + ch[t7 as usize] = ti2 + ti3; + let ci3 = ti2 - ti3; + let cr2 = tr1 - tr4; + let cr4 = tr1 + tr4; + let ci2 = ti1 + ti4; + let ci4 = ti1 - ti4; + let mut t8 = t7 + t0; + ch[t8 as usize - 1] = wa1[i - 2] * cr2 - wa1[i - 1] * ci2; + ch[t8 as usize] = wa1[i - 2] * ci2 + wa1[i - 1] * cr2; + t8 += t0; + ch[t8 as usize - 1] = wa2[i - 2] * cr3 - wa2[i - 1] * ci3; + ch[t8 as usize] = wa2[i - 2] * ci3 + wa2[i - 1] * cr3; + t8 += t0; + ch[t8 as usize - 1] = wa3[i - 2] * cr4 - wa3[i - 1] * ci4; + ch[t8 as usize] = wa3[i - 2] * ci4 + wa3[i - 1] * cr4; + } + t1 += ido; + } + + if ido % 2 == 1 { + return; + } + } + + t1 = ido; + t2 = ido << 2 as i32; + let mut t3 = ido - 1 as i32; + let mut t4 = ido + (ido << 1 as i32); + for _ in 0..l1 { + let mut t5 = t3; + let ti1 = cc[t1 as usize] + cc[t4 as usize]; + let ti2 = cc[t4 as usize] - cc[t1 as usize]; + let tr1 = cc[t1 as usize - 1] - cc[t4 as usize - 1]; + let tr2 = cc[t1 as usize - 1] + cc[t4 as usize - 1]; + ch[t5 as usize] = tr2 + tr2; + t5 += t0; + ch[t5 as usize] = SQRT2 * (tr1 - ti1); + t5 += t0; + ch[t5 as usize] = ti2 + ti2; + t5 += t0; + ch[t5 as usize] = -SQRT2 * (tr1 + ti1); + t3 += ido; + t1 += t2; + t4 += t2; + } +} + +#[inline(always)] +fn dradbg_l102( + ido: i32, + nbd: i32, + ipp2: i32, + ipph: i32, + l1: i32, + t0: i32, + t10: i32, + cc: &[f32], + ch: &mut [f32], +) { + if ido != 1 { + if nbd < l1 { + let mut t1 = 0 as i32; + let mut t2 = ipp2 * t0; + let mut t7 = 0 as i32; + for _ in 1..ipph { + t1 += t0; + t2 -= t0; + let mut t3 = t1; + let mut t4 = t2; + t7 += ido << 1 as i32; + let mut t8 = t7; + let mut t9 = t7; + for _ in (2..ido as usize).step_by(2) { + t3 += 2 as i32; + t4 += 2 as i32; + t8 += 2 as i32; + t9 -= 2 as i32; + let mut t5 = t3; + let mut t6 = t4; + let mut t11 = t8; + let mut t12 = t9; + for _ in 0..l1 { + ch[t5 as usize - 1] = + cc[t11 as usize - 1] + cc[t12 as usize - 1]; + ch[t6 as usize - 1] = + cc[t11 as usize - 1] - cc[t12 as usize - 1]; + ch[t5 as usize] = cc[t11 as usize] - cc[t12 as usize]; + ch[t6 as usize] = cc[t11 as usize] + cc[t12 as usize]; + t5 += ido; + t6 += ido; + t11 += t10; + t12 += t10; + } + } + } + } else { + let mut t1 = 0 as i32; + let mut t2 = ipp2 * t0; + let mut t7 = 0 as i32; + for _ in 1..ipph { + t1 += t0; + t2 -= t0; + let mut t3 = t1; + let mut t4 = t2; + t7 += ido << 1 as i32; + let mut t8 = t7; + for _ in 0..l1 { + let mut t5 = t3; + let mut t6 = t4; + let mut t9 = t8; + let mut t11 = t8; + for _ in (2..ido as usize).step_by(2) { + t5 += 2 as i32; + t6 += 2 as i32; + t9 += 2 as i32; + t11 -= 2 as i32; + ch[t5 as usize - 1] = + cc[t9 as usize - 1] + cc[t11 as usize - 1]; + ch[t6 as usize - 1] = + cc[t9 as usize - 1] - cc[t11 as usize - 1]; + ch[t5 as usize] = cc[t9 as usize] - cc[t11 as usize]; + ch[t6 as usize] = cc[t9 as usize] + cc[t11 as usize]; + } + t3 += ido; + t4 += ido; + t8 += t10; + } + } + } + } +} + +#[inline(always)] +fn dradbg_l103( + ido: i32, + nbd: i32, + ipp2: i32, + ipph: i32, + l1: i32, + t0: i32, + cc: &[f32], + ch: &mut [f32], +) { + if ido != 1 { + if nbd < l1 { + let mut t1 = 0 as i32; + let mut t2 = ipp2 * t0; + for _ in 1..ipph { + t1 += t0; + t2 -= t0; + let mut t3 = t1; + let mut t4 = t2; + for _ in (2..ido as usize).step_by(2) { + t3 += 2 as i32; + t4 += 2 as i32; + let mut t5 = t3; + let mut t6 = t4; + for _ in 0..l1 { + ch[t5 as usize - 1] = + cc[t5 as usize - 1] - cc[t6 as usize]; + ch[t6 as usize - 1] = + cc[t5 as usize - 1] + cc[t6 as usize]; + ch[t5 as usize] = + cc[t5 as usize] + cc[t6 as usize - 1]; + ch[t6 as usize] = + cc[t5 as usize] - cc[t6 as usize - 1]; + t5 += ido; + t6 += ido; + } + } + } + } else { + let mut t1 = 0 as i32; + let mut t2 = ipp2 * t0; + for _ in 1..ipph { + t1 += t0; + t2 -= t0; + let mut t3 = t1; + let mut t4 = t2; + for _ in 0..l1 { + let mut t5 = t3; + let mut t6 = t4; + for _ in (2..ido as usize).step_by(2) { + t5 += 2 as i32; + t6 += 2 as i32; + ch[t5 as usize - 1] = + cc[t5 as usize - 1] - cc[t6 as usize]; + ch[t6 as usize - 1] = + cc[t5 as usize - 1] + cc[t6 as usize]; + ch[t5 as usize] = + cc[t5 as usize] + cc[t6 as usize - 1]; + ch[t6 as usize] = + cc[t5 as usize] - cc[t6 as usize - 1]; + } + t3 += ido; + t4 += ido; + } + } + } + } +} + +#[inline(always)] +fn dradbg_l104( + ido: i32, + nbd: i32, + ip: i32, + l1: i32, + t0: i32, + cc: &mut [f32], + ch: &[f32], + wa: &[f32], +) { + if nbd > l1 { + let mut is = -ido - 1 as i32; + let mut t1 = 0 as i32; + for _ in 1..ip { + is += ido; + t1 += t0; + let mut t2 = t1; + for _ in 0..l1 { + let mut idij = is; + let mut t3 = t2; + for _ in (2..ido as usize).step_by(2) { + idij += 2 as i32; + t3 += 2 as i32; + cc[t3 as usize - 1] = wa[idij as usize - 1] + * ch[t3 as usize - 1] + - wa[idij as usize] * ch[t3 as usize]; + cc[t3 as usize] = wa[idij as usize - 1] * ch[t3 as usize] + + wa[idij as usize] * ch[t3 as usize - 1]; + } + t2 += ido; + } + } + } else { + let mut is = -ido - 1 as i32; + let mut t1 = 0 as i32; + for _ in 1..ip { + is += ido; + t1 += t0; + let mut idij = is; + let mut t2 = t1; + for _ in (2..ido as usize).step_by(2) { + t2 += 2 as i32; + idij += 2 as i32; + let mut t3 = t2; + for _ in 0..l1 { + cc[t3 as usize - 1] = wa[idij as usize - 1] + * ch[t3 as usize - 1] + - wa[idij as usize] * ch[t3 as usize]; + cc[t3 as usize] = wa[idij as usize - 1] * ch[t3 as usize] + + wa[idij as usize] * ch[t3 as usize - 1]; + t3 += ido; + } + } + } + } +} + +pub(crate) fn dradbg( + ido: i32, + ip: i32, + l1: i32, + idl1: i32, + cc: &mut [f32], + ch: &mut [f32], + wa: &[f32], +) { + const TPI: f32 = 6.283_185_307_179_586; + + let t10 = ip * ido; + let t0 = l1 * ido; + let arg = TPI / ip as f32; + let dcp = f64::cos(arg as f64) as f32; + let dsp = f64::sin(arg as f64) as f32; + let nbd = ido - 1 as i32 >> 1 as i32; + let ipp2 = ip; + let ipph = ip + 1 as i32 >> 1 as i32; + + if ido < l1 { + let mut t1 = 0 as i32; + for _ in 0..ido { + let mut t2 = t1; + let mut t3 = t1; + for _ in 0..l1 { + ch[t2 as usize] = cc[t3 as usize]; + t2 += ido; + t3 += t10; + } + t1 += 1; + } + } else { + let mut t1 = 0 as i32; + let mut t2 = 0 as i32; + for _ in 0..l1 { + let mut t3 = t1; + let mut t4 = t2; + for _ in 0..ido { + ch[t3 as usize] = cc[t4 as usize]; + t3 += 1; + t4 += 1; + } + t1 += ido; + t2 += t10; + } + } + + let mut t1 = 0 as i32; + let mut t2 = ipp2 * t0; + let mut t5 = ido << 1 as i32; + let t7 = t5; + + for _ in 1..ipph { + t1 += t0; + t2 -= t0; + let mut t3 = t1; + let mut t4 = t2; + let mut t6 = t5; + for _ in 0..l1 { + ch[t3 as usize] = cc[t6 as usize - 1] + cc[t6 as usize - 1]; + ch[t4 as usize] = cc[t6 as usize] + cc[t6 as usize]; + t3 += ido; + t4 += ido; + t6 += t10; + } + t5 += t7; + } + + dradbg_l102(ido, nbd, ipp2, ipph, l1, t0, t10, cc, ch); + + let ar1: f32 = 1.0; + let ai1: f32 = 0.0; + let mut t1 = 0 as i32; + let mut t2 = ipp2 * idl1; + let t9 = t2; + let t3 = (ip - 1 as i32) * idl1; + + for _ in 1..ipph { + t1 += idl1; + t2 -= idl1; + let ar1h = dcp * ar1 - dsp * ai1; + let ai1 = dcp * ai1 + dsp * ar1; + let ar1 = ar1h; + let mut t4 = t1; + let mut t5 = t2; + let mut t6 = 0 as i32; + let mut t7 = idl1; + let mut t8 = t3; + + for _ in 0..idl1 { + let fresh13 = t6; + t6 = t6 + 1; + let fresh14 = t7; + t7 = t7 + 1; + let fresh15 = t4; + t4 = t4 + 1; + cc[fresh15 as usize] = + ch[fresh13 as usize] + ar1 * ch[fresh14 as usize]; + let fresh16 = t8; + t8 = t8 + 1; + let fresh17 = t5; + t5 = t5 + 1; + cc[fresh17 as usize] = ai1 * ch[fresh16 as usize]; + } + let dcc = ar1; + let ds2 = ai1; + let ar2 = ar1; + let ai2 = ai1; + let mut t6 = idl1; + let mut t7 = t9 - idl1; + + for _ in 2..ipph { + t6 += idl1; + t7 -= idl1; + let ar2h = dcc * ar2 - ds2 * ai2; + let ai2 = dcc * ai2 + ds2 * ar2; + let ar2 = ar2h; + let mut t4 = t1; + let mut t5 = t2; + let mut t11 = t6; + let mut t12 = t7; + + for _ in 0..idl1 { + let fresh18 = t11; + t11 = t11 + 1; + let fresh19 = t4; + t4 = t4 + 1; + cc[fresh19 as usize] += ar2 * ch[fresh18 as usize]; + let fresh20 = t12; + t12 = t12 + 1; + let fresh21 = t5; + t5 = t5 + 1; + cc[fresh21 as usize] += ai2 * ch[fresh20 as usize]; + } + } + } + + t1 = 0 as i32; + for _ in 1..ipph { + t1 += idl1; + let mut t2 = t1; + for ik in 0..idl1 { + let fresh22 = t2; + t2 = t2 + 1; + ch[ik as usize] += ch[fresh22 as usize]; + } + } + + t1 = 0 as i32; + t2 = ipp2 * t0; + for _ in 1..ipph { + t1 += t0; + t2 -= t0; + let mut t3 = t1; + let mut t4 = t2; + for _ in 0..l1 { + ch[t3 as usize] = cc[t3 as usize] - cc[t4 as usize]; + ch[t4 as usize] = cc[t3 as usize] + cc[t4 as usize]; + t3 += ido; + t4 += ido; + } + } + + dradbg_l103(ido, nbd, ipp2, ipph, l1, t0, cc, ch); + + if ido == 1 { + return; + } + + for ik in 0..idl1 { + cc[ik as usize] = ch[ik as usize]; + } + + t1 = 0 as i32; + for _ in 1..ip { + t1 += t0; + t2 = t1; + for _ in 0..l1 { + cc[t2 as usize] = ch[t2 as usize]; + t2 += ido; + } + } + + dradbg_l104(ido, nbd, ip, l1, t0, cc, ch, wa); +} diff --git a/fft/src/dradf.rs b/fft/src/dradf.rs new file mode 100644 index 0000000..249ce1d --- /dev/null +++ b/fft/src/dradf.rs @@ -0,0 +1,558 @@ +pub(crate) fn dradf2( + ido: i32, + l1: i32, + cc: &mut [f32], + ch: &mut [f32], + wa1: &mut [f32], +) { + let mut t1 = 0; + let mut t2 = l1 * ido; + let t0 = t2; + let mut t3 = ido << 1 as i32; + + for _ in 0..l1 { + ch[(t1 as usize) << 1] = cc[t1 as usize] + cc[t2 as usize]; + ch[(t1 as usize) << 1 + t3 as usize - 1] = + cc[t1 as usize] - cc[t2 as usize]; + t1 += ido; + t2 += ido; + } + + if ido < 2 { + return; + } + + if ido != 2 { + t1 = 0; + t2 = t0; + for _ in 0..l1 { + let mut t3 = t2; + let mut t4 = (t1 << 1) + (ido << 1); + let mut t5 = t1; + let mut t6 = t1 + t1; + for i in (2..ido as usize).step_by(2) { + t3 += 2; + t4 -= 2; + t5 += 2; + t6 += 2; + let tr2 = wa1[i - 2] * cc[t3 as usize - 1] + + wa1[i - 1] * cc[t3 as usize]; + let ti2 = wa1[i - 2] * cc[t3 as usize] + - wa1[i - 1] * cc[t3 as usize - 1]; + ch[t6 as usize] = cc[t5 as usize] + ti2; + ch[t4 as usize] = ti2 - cc[t5 as usize]; + ch[t6 as usize - 1] = cc[t5 as usize - 1] + tr2; + ch[t4 as usize - 1] = cc[t5 as usize - 1] - tr2; + } + t1 += ido; + t2 += ido; + } + + if ido % 2 == 1 { + return; + } + } + + t1 = ido; + t2 = t1 - 1; + t3 = t2; + t2 += t0; + + for _ in 0..l1 { + ch[t1 as usize] = -cc[t2 as usize]; + ch[t1 as usize - 1] = cc[t3 as usize]; + t1 += ido << 1; + t2 += ido; + t3 += ido; + } +} + +pub(crate) fn dradf4( + ido: i32, + l1: i32, + cc: &mut [f32], + ch: &mut [f32], + wa1: &[f32], + wa2: &[f32], + wa3: &[f32], +) { + const HSQT2: f32 = 0.707_106_781_186_547_52; + + let t0 = l1 * ido; + let mut t1 = t0; + let mut t4 = t1 << 1; + let mut t2 = t1 + (t1 << 1); + let mut t3 = 0; + + for _ in 0..l1 { + let tr1 = cc[t1 as usize] + cc[t2 as usize]; + let tr2 = cc[t3 as usize] + cc[t4 as usize]; + let mut t5 = t3 << 2; + ch[t5 as usize] = tr1 + tr2; + ch[(ido << 2) as usize + t5 as usize - 1] = tr2 - tr1; + t5 += ido << 1; + ch[t5 as usize - 1] = cc[t3 as usize] - cc[t4 as usize]; + ch[t5 as usize] = cc[t2 as usize] - cc[t1 as usize]; + t1 += ido; + t2 += ido; + t3 += ido; + t4 += ido; + } + + if ido < 2 { + return; + } + + if ido != 2 { + t1 = 0; + for _ in 0..l1 { + t2 = t1; + t4 = t1 << 2; + let t6 = ido << 1; + let mut t5 = t6 + t4; + for i in (2..ido as usize).step_by(2) { + t2 += 2; + t3 = t2; + t4 += 2; + t5 -= 2; + t3 += t0; + let cr2 = wa1[i - 2] * cc[t3 as usize - 1] + + wa1[i - 1] * cc[t3 as usize]; + let ci2 = wa1[i - 2] * cc[t3 as usize] + - wa1[i - 1] * cc[t3 as usize - 1]; + t3 += t0; + let cr3 = wa2[i - 2] * cc[t3 as usize - 1] + + wa2[i - 1] * cc[t3 as usize]; + let ci3 = wa2[i - 2] * cc[t3 as usize] + - wa2[i - 1] * cc[t3 as usize - 1]; + t3 += t0; + let cr4 = wa3[i - 2] * cc[t3 as usize - 1] + + wa3[i - 1] * cc[t3 as usize]; + let ci4 = wa3[i - 2] * cc[t3 as usize] + - wa3[i - 1] * cc[t3 as usize - 1]; + let tr1 = cr2 + cr4; + let tr4 = cr4 - cr2; + let ti1 = ci2 + ci4; + let ti4 = ci2 - ci4; + let ti2 = cc[t2 as usize] + ci3; + let ti3 = cc[t2 as usize] - ci3; + let tr2 = cc[t2 as usize - 1] + cr3; + let tr3 = cc[t2 as usize - 1] - cr3; + ch[t4 as usize - 1] = tr1 + tr2; + ch[t4 as usize] = ti1 + ti2; + ch[t5 as usize - 1] = tr3 - ti4; + ch[t5 as usize] = tr4 - ti3; + ch[t4 as usize + t6 as usize - 1] = ti4 + tr3; + ch[t4 as usize + t6 as usize] = tr4 + ti3; + ch[t5 as usize + t6 as usize - 1] = tr2 - tr1; + ch[t5 as usize + t6 as usize] = ti1 - ti2; + } + t1 += ido; + } + + if ido & 1 != 0 { + return; + } + } + + t1 = t0 + ido - 1; + t2 = t1 + (t0 << 1); + t3 = ido << 2; + t4 = ido; + let t5 = ido << 1; + let mut t6 = ido; + + for _ in 0..l1 { + let ti1 = -HSQT2 * (cc[t1 as usize] + cc[t2 as usize]); + let tr1 = HSQT2 * (cc[t1 as usize] - cc[t2 as usize]); + ch[t4 as usize - 1] = tr1 + cc[t6 as usize - 1]; + ch[t4 as usize + t5 as usize - 1] = cc[t6 as usize - 1] - tr1; + ch[t4 as usize] = ti1 - cc[t1 as usize + t0 as usize]; + ch[t4 as usize + t5 as usize] = ti1 + cc[t1 as usize + t0 as usize]; + t1 += ido; + t2 += ido; + t4 += t3; + t6 += ido; + } +} + +#[inline(always)] +fn dradfg_l102( + ido: i32, + nbd: i32, + ip: i32, + l1: i32, + t0: i32, + cc: &[f32], + wa: &[f32], + ch: &mut [f32], +) { + let mut is = -ido; + let mut t1 = 0 as i32; + if nbd > l1 { + for _ in 1..ip { + t1 += t0; + is += ido; + let mut t2 = -ido + t1; + for _ in 0..l1 { + let mut idij = is - 1 as i32; + t2 += ido; + let mut t3 = t2; + for _ in (2..ido as usize).step_by(2) { + idij += 2 as i32; + t3 += 2 as i32; + ch[t3 as usize - 1] = wa[idij as usize - 1] + * cc[t3 as usize - 1] + + wa[idij as usize] * cc[t3 as usize]; + ch[t3 as usize] = wa[idij as usize - 1] * cc[t3 as usize] + - wa[idij as usize] * cc[t3 as usize - 1]; + } + } + } + } else { + for _ in 0..ip { + is += ido; + let mut idij = is - 1 as i32; + t1 += t0; + let mut t2 = t1; + for _ in (2..ido as usize).step_by(2) { + idij += 2 as i32; + t2 += 2 as i32; + let mut t3 = t2; + for _ in 0..l1 { + ch[t3 as usize - 1] = wa[idij as usize - 1] + * cc[t3 as usize - 1] + + wa[idij as usize] * cc[t3 as usize]; + ch[t3 as usize] = wa[idij as usize - 1] * cc[t3 as usize] + - wa[idij as usize] * cc[t3 as usize - 1]; + t3 += ido; + } + } + } + } +} + +#[inline(always)] +fn dradfg_l103( + ido: i32, + nbd: i32, + ipp2: i32, + ipph: i32, + l1: i32, + t0: i32, + ch: &[f32], + cc: &mut [f32], +) { + let mut t1 = 0 as i32; + let mut t2 = ipp2 * t0; + if nbd < l1 { + for _ in 1..ipph { + t1 += t0; + t2 -= t0; + let mut t3 = t1; + let mut t4 = t2; + for _ in (2..ido as usize).step_by(2) { + t3 += 2 as i32; + t4 += 2 as i32; + let mut t5 = t3 - ido; + let mut t6 = t4 - ido; + for _ in 0..l1 { + t5 += ido; + t6 += ido; + cc[t5 as usize - 1] = + ch[t5 as usize - 1] + ch[t6 as usize - 1]; + cc[t6 as usize - 1] = ch[t5 as usize] - ch[t6 as usize]; + cc[t5 as usize] = ch[t5 as usize] + ch[t6 as usize]; + cc[t6 as usize] = + ch[t6 as usize - 1] - ch[t5 as usize - 1]; + } + } + } + } else { + for _ in 1..ipph { + t1 += t0; + t2 -= t0; + let mut t3 = t1; + let mut t4 = t2; + for _ in 0..l1 { + let mut t5 = t3; + let mut t6 = t4; + for _ in (2..ido as usize).step_by(2) { + t5 += 2 as i32; + t6 += 2 as i32; + cc[t5 as usize - 1] = + ch[t5 as usize - 1] + ch[t6 as usize - 1]; + cc[t6 as usize - 1] = ch[t5 as usize] - ch[t6 as usize]; + cc[t5 as usize] = ch[t5 as usize] + ch[t6 as usize]; + cc[t6 as usize] = + ch[t6 as usize - 1] - ch[t5 as usize - 1]; + } + t3 += ido; + t4 += ido; + } + } + } +} + +#[inline(always)] +fn dradfg_l104( + ido: i32, + nbd: i32, + idp2: i32, + ipp2: i32, + ipph: i32, + l1: i32, + t0: i32, + t2: i32, + t10: i32, + ch: &[f32], + cc: &mut [f32], +) { + if nbd < l1 { + let mut t1 = -ido; + let mut t3 = 0 as i32; + let mut t4 = 0 as i32; + let mut t5 = ipp2 * t0; + for _ in 1..ipph { + t1 += t2; + t3 += t2; + t4 += t0; + t5 -= t0; + for i in (2..ido as i32).step_by(2) { + let mut t6 = idp2 + t1 - i; + let mut t7 = i + t3; + let mut t8 = i + t4; + let mut t9 = i + t5; + for _ in 0..l1 { + cc[t7 as usize - 1] = + ch[t8 as usize - 1] + ch[t9 as usize - 1]; + cc[t6 as usize - 1] = + ch[t8 as usize - 1] - ch[t9 as usize - 1]; + cc[t7 as usize] = ch[t8 as usize] + ch[t9 as usize]; + cc[t6 as usize] = ch[t9 as usize] - ch[t8 as usize]; + t6 += t10; + t7 += t10; + t8 += ido; + t9 += ido; + } + } + } + } else { + let mut t1 = -ido; + let mut t3 = 0 as i32; + let mut t4 = 0 as i32; + let mut t5 = ipp2 * t0; + for _ in 1..ipph { + t1 += t2; + t3 += t2; + t4 += t0; + t5 -= t0; + let mut t6 = t1; + let mut t7 = t3; + let mut t8 = t4; + let mut t9 = t5; + for _ in 0..l1 { + for i in (2..ido as usize).step_by(2) { + let ic = idp2 - i as i32; + cc[i + t7 as usize - 1] = + ch[i + t8 as usize - 1] + ch[i + t9 as usize - 1]; + cc[ic as usize + t6 as usize - 1] = + ch[i + t8 as usize - 1] - ch[i + t9 as usize - 1]; + cc[i + t7 as usize] = + ch[i + t8 as usize] + ch[i + t9 as usize]; + cc[ic as usize + t6 as usize] = + ch[i + t9 as usize] - ch[i + t8 as usize]; + } + t6 += t10; + t7 += t10; + t8 += ido; + t9 += ido; + } + } + } +} + +pub(crate) fn dradfg( + ido: i32, + ip: i32, + l1: i32, + idl1: i32, + cc: &mut [f32], + ch: &mut [f32], + wa: &[f32], +) { + const TPI: f32 = 6.283_185_307_179_586; + + let arg = TPI / ip as f32; + let dcp = f64::cos(arg as f64) as f32; + let dsp = f64::sin(arg as f64) as f32; + let ipph = ip + 1 >> 1; + let ipp2 = ip; + let idp2 = ido; + let nbd = ido - 1 >> 1; + let t0 = l1 * ido; + let t10 = ip * ido; + + if ido != 1 { + for ik in 0..idl1 { + ch[ik as usize] = cc[ik as usize]; + } + + let mut t1 = 0; + for _ in 1..ip { + t1 += t0; + let mut t2 = t1; + for _ in 0..l1 { + ch[t2 as usize] = cc[t2 as usize]; + t2 += ido; + } + } + + dradfg_l102(ido, nbd, ip, l1, t0, cc, wa, ch); + + dradfg_l103(ido, nbd, ipp2, ipph, l1, t0, ch, cc); + } + + for ik in 0..idl1 { + cc[ik as usize] = ch[ik as usize]; + } + + let mut t1 = 0 as i32; + let mut t2 = ipp2 * idl1; + for _ in 1..ipph { + t1 += t0; + t2 -= t0; + let mut t3 = t1 - ido; + let mut t4 = t2 - ido; + for _ in 0..l1 { + t3 += ido; + t4 += ido; + cc[t3 as usize] = ch[t3 as usize] + ch[t4 as usize]; + cc[t4 as usize] = ch[t4 as usize] - ch[t3 as usize]; + } + } + + let mut ar1: f32 = 1.0; + let mut ai1: f32 = 0.0; + t1 = 0 as i32; + t2 = ipp2 * idl1; + let mut t3 = (ip - 1 as i32) * idl1; + + for _ in 1..ipph { + t1 += idl1; + t2 -= idl1; + let ar1h = dcp * ar1 - dsp * ai1; + ai1 = dcp * ai1 + dsp * ar1; + ar1 = ar1h; + let mut t4 = t1; + let mut t5 = t2; + let mut t6 = t3; + let mut t7 = idl1; + for ik in 0..idl1 { + let fresh2 = t7; + t7 = t7 + 1; + let fresh3 = t4; + t4 = t4 + 1; + ch[fresh3 as usize] = cc[ik as usize] + ar1 * cc[fresh2 as usize]; + let fresh4 = t6; + t6 = t6 + 1; + let fresh5 = t5; + t5 = t5 + 1; + ch[fresh5 as usize] = ai1 * cc[fresh4 as usize]; + } + let dc2 = ar1; + let ds2 = ai1; + let mut ar2 = ar1; + let mut ai2 = ai1; + t4 = idl1; + t5 = (ipp2 - 1 as i32) * idl1; + for _ in 2..ipph { + t4 += idl1; + t5 -= idl1; + let ar2h = dc2 * ar2 - ds2 * ai2; + ai2 = dc2 * ai2 + ds2 * ar2; + ar2 = ar2h; + let mut t6 = t1; + let mut t7 = t2; + let mut t8 = t4; + let mut t9 = t5; + for _ in 0..idl1 { + let fresh6 = t8; + t8 = t8 + 1; + let fresh7 = t6; + t6 = t6 + 1; + ch[fresh7 as usize] += ar2 * cc[fresh6 as usize]; + let fresh8 = t9; + t9 = t9 + 1; + let fresh9 = t7; + t7 = t7 + 1; + ch[fresh9 as usize] += ai2 * cc[fresh8 as usize]; + } + } + } + + t1 = 0 as i32; + for _ in 1..ipph { + t1 += idl1; + t2 = t1; + for ik in 0..idl1 { + let fresh10 = t2; + t2 = t2 + 1; + ch[ik as usize] += cc[fresh10 as usize]; + } + } + + if ido < l1 { + for i in 0..ido { + t1 = i; + t2 = i; + for _ in 0..l1 { + cc[t2 as usize] = ch[t1 as usize]; + t1 += ido; + t2 += t10; + } + } + } else { + t1 = 0 as i32; + t2 = 0 as i32; + for _ in 0..l1 { + t3 = t1; + let mut t4 = t2; + for _ in 0..ido { + let fresh11 = t3; + t3 = t3 + 1; + let fresh12 = t4; + t4 = t4 + 1; + cc[fresh12 as usize] = ch[fresh11 as usize]; + } + t1 += ido; + t2 += t10; + } + } + + t1 = 0 as i32; + t2 = ido << 1 as i32; + t3 = 0 as i32; + let mut t4 = ipp2 * t0; + for _ in 1..ipph { + t1 += t2; + t3 += t0; + t4 -= t0; + let mut t5 = t1; + let mut t6 = t3; + let mut t7 = t4; + for _ in 0..l1 { + cc[t5 as usize - 1] = ch[t6 as usize]; + cc[t5 as usize] = ch[t7 as usize]; + t5 += t10; + t6 += ido; + t7 += ido; + } + } + + if ido == 1 { + return; + } + + dradfg_l104(ido, nbd, idp2, ipp2, ipph, l1, t0, t2, t10, ch, cc); +} diff --git a/fft/src/fftwrap.rs b/fft/src/fftwrap.rs index 165d6b6..e8773d9 100644 --- a/fft/src/fftwrap.rs +++ b/fft/src/fftwrap.rs @@ -1,18 +1,3 @@ -#![allow( - dead_code, - mutable_transmutes, - non_camel_case_types, - non_snake_case, - non_upper_case_globals, - unused_assignments, - unused_mut -)] - -use std::{ - ffi::c_void, - os::raw::{c_double, c_float, c_int}, -}; - use crate::smallft::*; /* Copyright (C) 2005-2006 Jean-Marc Valin @@ -48,79 +33,46 @@ use crate::smallft::*; SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. */ -#[no_mangle] -pub unsafe extern "C" fn spx_fft_init(mut size: c_int) -> *mut c_void { - let mut table = Box::new(drft_lookup::new(size as usize)); - - return Box::into_raw(table) as *mut c_void; -} - -#[no_mangle] -pub unsafe extern "C" fn spx_fft_destroy(mut table: *mut c_void) { - let table = table as *mut drft_lookup; - - Box::from_raw(table); +pub fn spx_fft_init(size: usize) -> DrftLookup { + DrftLookup::new(size) } -#[no_mangle] -pub unsafe extern "C" fn spx_fft( - mut table: *mut c_void, - mut in_0: *mut c_float, - mut out: *mut c_float, -) { +pub fn spx_fft(table: &mut DrftLookup, in_0: &mut [f32], out: &mut [f32]) { + let scale = (1.0f64 / table.n as f64) as f32; if in_0 == out { - let mut i: c_int = 0; - let mut scale: c_float = - (1.0f64 / (*(table as *mut drft_lookup)).n as c_double) as c_float; eprintln!("FFT should not be done in-place"); - i = 0 as c_int; - while i < (*(table as *mut drft_lookup)).n { - *out.offset(i as isize) = scale * *in_0.offset(i as isize); - i += 1 - } - } else { - let mut i_0: c_int = 0; - let mut scale_0: c_float = - (1.0f64 / (*(table as *mut drft_lookup)).n as c_double) as c_float; - i_0 = 0 as c_int; - while i_0 < (*(table as *mut drft_lookup)).n { - *out.offset(i_0 as isize) = scale_0 * *in_0.offset(i_0 as isize); - i_0 += 1 - } } - spx_drft_forward(table as *mut drft_lookup, out); + + out.iter_mut() + .zip(in_0.iter()) + .take(table.n as usize) + .for_each(|(o, i)| *o = scale * *i); + + spx_drft_forward(table, out); } -#[no_mangle] -pub unsafe extern "C" fn spx_ifft( - mut table: *mut c_void, - mut in_0: *mut c_float, - mut out: *mut c_float, -) { + +pub fn spx_ifft(table: &mut DrftLookup, in_0: &mut [f32], out: &mut [f32]) { if in_0 == out { eprintln!("FFT should not be done in-place"); } else { - let mut i: c_int = 0; - i = 0 as c_int; - while i < (*(table as *mut drft_lookup)).n { - *out.offset(i as isize) = *in_0.offset(i as isize); - i += 1 - } + out.copy_from_slice(&in_0[..table.n as usize]); } - spx_drft_backward(table as *mut drft_lookup, out); + + spx_drft_backward(table, out); } -#[no_mangle] -pub unsafe extern "C" fn spx_fft_float( - mut table: *mut c_void, - mut in_0: *mut c_float, - mut out: *mut c_float, + +pub fn spx_fft_float( + table: &mut DrftLookup, + in_0: &mut [f32], + out: &mut [f32], ) { spx_fft(table, in_0, out); } -#[no_mangle] -pub unsafe extern "C" fn spx_ifft_float( - mut table: *mut c_void, - mut in_0: *mut c_float, - mut out: *mut c_float, + +pub fn spx_ifft_float( + table: &mut DrftLookup, + in_0: &mut [f32], + out: &mut [f32], ) { spx_ifft(table, in_0, out); } diff --git a/fft/src/lib.rs b/fft/src/lib.rs index 1e16137..74e6239 100644 --- a/fft/src/lib.rs +++ b/fft/src/lib.rs @@ -1,3 +1,5 @@ +mod dradb; +mod dradf; mod fftwrap; mod smallft; diff --git a/fft/src/smallft.rs b/fft/src/smallft.rs index 1ab1409..09e6744 100644 --- a/fft/src/smallft.rs +++ b/fft/src/smallft.rs @@ -1,1000 +1,192 @@ -#![allow( - dead_code, - mutable_transmutes, - non_camel_case_types, - non_snake_case, - non_upper_case_globals, - unused_assignments, - unused_mut -)] +use crate::dradb::*; +use crate::dradf::*; -use std::os::raw::{c_double, c_float, c_int}; - -/* ******************************************************************* -* * -* THIS FILE IS PART OF THE OggVorbis SOFTWARE CODEC SOURCE CODE. * -* USE, DISTRIBUTION AND REPRODUCTION OF THIS LIBRARY SOURCE IS * -* GOVERNED BY A BSD-STYLE SOURCE LICENSE INCLUDED WITH THIS SOURCE * -* IN 'COPYING'. PLEASE READ THESE TERMS BEFORE DISTRIBUTING. * -* * -* THE OggVorbis SOURCE CODE IS (C) COPYRIGHT 1994-2001 * -* by the XIPHOPHORUS Company http://www.xiph.org/ * -* * -******************************************************************** - -function: fft transform -last mod: $Id: smallft.h,v 1.3 2003/09/16 18:35:45 jm Exp $ - -********************************************************************/ -/* * - @file smallft.h - @brief Discrete Rotational Fourier Transform (DRFT) -*/ -/* * Discrete Rotational Fourier Transform lookup */ #[derive(Clone)] -#[repr(C)] -pub struct drft_lookup { - pub n: c_int, +pub struct DrftLookup { + pub n: usize, pub trigcache: Vec, pub splitcache: Vec, } -impl drft_lookup { +impl DrftLookup { pub fn new(n: usize) -> Self { - let mut l = Self { - n: n as c_int, - trigcache: vec![0.0; 3 * (n as usize)], + let mut drft = Self { + n: n, + trigcache: vec![0.0; 3 * n], splitcache: vec![0; 32], }; - unsafe { - fdrffti(n, &mut l.trigcache, &mut l.splitcache); - } - l + fdrffti(n, &mut drft.trigcache, &mut drft.splitcache); + + drft } } -/* ******************************************************************* -* * -* THIS FILE IS PART OF THE OggVorbis SOFTWARE CODEC SOURCE CODE. * -* USE, DISTRIBUTION AND REPRODUCTION OF THIS LIBRARY SOURCE IS * -* GOVERNED BY A BSD-STYLE SOURCE LICENSE INCLUDED WITH THIS SOURCE * -* IN 'COPYING'. PLEASE READ THESE TERMS BEFORE DISTRIBUTING. * -* * -* THE OggVorbis SOURCE CODE IS (C) COPYRIGHT 1994-2001 * -* by the XIPHOPHORUS Company http://www.xiph.org/ * -* * -******************************************************************** - -function: *unnormalized* fft transform -last mod: $Id: smallft.c,v 1.19 2003/10/08 05:12:37 jm Exp $ - -********************************************************************/ -/* FFT implementation from OggSquish, minus cosine transforms, - * minus all but radix 2/4 case. In Vorbis we only need this - * cut-down version. - * - * To do more than just power-of-two sized vectors, see the full - * version I wrote for NetLib. - * - * Note that the packing is a little strange; rather than the FFT r/i - * packing following R_0, I_n, R_1, I_1, R_2, I_2 ... R_n-1, I_n-1, - * it follows R_0, R_1, I_1, R_2, I_2 ... R_n-1, I_n-1, I_n like the - * FORTRAN version - */ -unsafe extern "C" fn drfti1(wa: &mut [f32], ifac: &mut [i32]) { - static mut ntryh: [c_int; 4] = - [4 as c_int, 2 as c_int, 3 as c_int, 5 as c_int]; - static mut tpi: c_float = 6.28318530717958648f32; +#[inline(always)] +fn drfti1_c_10244(ifac: &mut [i32], n: i32, nf: &mut i32) { + const NTRYH: [i32; 4] = [4, 2, 3, 5]; - let mut n = wa.len() as i32; - let mut wa = wa.as_mut_ptr(); - let mut ifac = ifac.as_mut_ptr(); + let mut j: i32 = -1; + let mut ntry = 0; + let mut nl = n; - let mut arg: c_float = 0.; - let mut argh: c_float = 0.; - let mut argld: c_float = 0.; - let mut fi: c_float = 0.; - let mut ntry: c_int = 0 as c_int; - let mut i: c_int = 0; - let mut j: c_int = -(1 as c_int); - let mut k1: c_int = 0; - let mut l1: c_int = 0; - let mut l2: c_int = 0; - let mut ib: c_int = 0; - let mut ld: c_int = 0; - let mut ii: c_int = 0; - let mut ip: c_int = 0; - let mut is: c_int = 0; - let mut nq: c_int = 0; - let mut nr: c_int = 0; - let mut ido: c_int = 0; - let mut ipm: c_int = 0; - let mut nfm1: c_int = 0; - let mut nl: c_int = n; - let mut nf: c_int = 0 as c_int; 'c_10244: loop { j += 1; - if j < 4 as c_int { - ntry = ntryh[j as usize] + if j < 4 { + ntry = NTRYH[j as usize] } else { - ntry += 2 as c_int + ntry += 2 } loop { - nq = nl / ntry; - nr = nl - ntry * nq; - if nr != 0 as c_int { + let nq = nl / ntry; + let nr = nl - ntry * nq; + if nr != 0 { break; } - nf += 1; - *ifac.offset((nf + 1 as c_int) as isize) = ntry; + *nf += 1; + ifac[*nf as usize + 1] = ntry; nl = nq; - if !(ntry != 2 as c_int) { - if !(nf == 1 as c_int) { - i = 1 as c_int; - while i < nf { - ib = nf - i + 1 as c_int; - *ifac.offset((ib + 1 as c_int) as isize) = - *ifac.offset(ib as isize); - i += 1 + if ntry == 2 { + if *nf != 1 { + for i in 1..*nf { + let ib = *nf - i + 1 as i32; + ifac[ib as usize + 1] = ifac[ib as usize]; } - *ifac.offset(2 as c_int as isize) = 2 as c_int + ifac[2] = 2; } } - if !(nl != 1 as c_int) { + if nl == 1 { break 'c_10244; } } } - *ifac.offset(0 as c_int as isize) = n; - *ifac.offset(1 as c_int as isize) = nf; - argh = tpi / n as c_float; - is = 0 as c_int; - nfm1 = nf - 1 as c_int; - l1 = 1 as c_int; - if nfm1 == 0 as c_int { +} + +fn drfti1(wa: &mut [f32], ifac: &mut [i32]) { + const TPI: f32 = 6.283_185_307_179_586_48; + + let n = wa.len() as i32; + let mut nf = 0; + + drfti1_c_10244(ifac, n, &mut nf); + + ifac[0] = n; + ifac[1] = nf; + let nfm1 = nf - 1; + let argh = TPI / n as f32; + let mut is = 0; + let mut l1 = 1; + + if nfm1 == 0 { return; } - k1 = 0 as c_int; - while k1 < nfm1 { - ip = *ifac.offset((k1 + 2 as c_int) as isize); - ld = 0 as c_int; - l2 = l1 * ip; - ido = n / l2; - ipm = ip - 1 as c_int; - j = 0 as c_int; + + for k1 in 0..nfm1 { + let ip = ifac[k1 as usize + 2]; + let l2 = l1 * ip; + let ido = n / l2; + let ipm = ip - 1; + let mut ld = 0; + let mut j = 0; while j < ipm { ld += l1; - i = is; - argld = ld as c_float * argh; - fi = 0.0f32; - ii = 2 as c_int; + let argld = ld as f32 * argh; + let mut i = is; + let mut ii = 2; + let mut fi = 0.0f32; while ii < ido { fi += 1.0f32; - arg = fi * argld; + let arg = fi * argld; let fresh0 = i; i = i + 1; - *wa.offset(fresh0 as isize) = - f64::cos(arg as c_double) as c_float; + wa[fresh0] = f64::cos(arg as f64) as f32; let fresh1 = i; i = i + 1; - *wa.offset(fresh1 as isize) = - f64::sin(arg as c_double) as c_float; - ii += 2 as c_int + wa[fresh1] = f64::sin(arg as f64) as f32; + ii += 2 as i32 } - is += ido; + is += ido as usize; j += 1 } l1 = l2; - k1 += 1 } } -unsafe extern "C" fn fdrffti( - n: usize, - wsave: &mut [f32], - mut ifac: &mut [i32], -) { + +fn fdrffti(n: usize, wsave: &mut [f32], ifac: &mut [i32]) { if n == 1 { return; } drfti1(&mut wsave[n..n * 2], ifac); } -unsafe extern "C" fn dradf2( - mut ido: c_int, - mut l1: c_int, - mut cc: *mut c_float, - mut ch: *mut c_float, - mut wa1: *mut c_float, -) { - let mut i: c_int = 0; - let mut k: c_int = 0; - let mut ti2: c_float = 0.; - let mut tr2: c_float = 0.; - let mut t0: c_int = 0; - let mut t1: c_int = 0; - let mut t2: c_int = 0; - let mut t3: c_int = 0; - let mut t4: c_int = 0; - let mut t5: c_int = 0; - let mut t6: c_int = 0; - t1 = 0 as c_int; - t2 = l1 * ido; - t0 = t2; - t3 = ido << 1 as c_int; - k = 0 as c_int; - while k < l1 { - *ch.offset((t1 << 1 as c_int) as isize) = - *cc.offset(t1 as isize) + *cc.offset(t2 as isize); - *ch.offset(((t1 << 1 as c_int) + t3 - 1 as c_int) as isize) = - *cc.offset(t1 as isize) - *cc.offset(t2 as isize); - t1 += ido; - t2 += ido; - k += 1 - } - if ido < 2 as c_int { - return; - } - if !(ido == 2 as c_int) { - t1 = 0 as c_int; - t2 = t0; - k = 0 as c_int; - while k < l1 { - t3 = t2; - t4 = (t1 << 1 as c_int) + (ido << 1 as c_int); - t5 = t1; - t6 = t1 + t1; - i = 2 as c_int; - while i < ido { - t3 += 2 as c_int; - t4 -= 2 as c_int; - t5 += 2 as c_int; - t6 += 2 as c_int; - tr2 = *wa1.offset((i - 2 as c_int) as isize) - * *cc.offset((t3 - 1 as c_int) as isize) - + *wa1.offset((i - 1 as c_int) as isize) - * *cc.offset(t3 as isize); - ti2 = *wa1.offset((i - 2 as c_int) as isize) - * *cc.offset(t3 as isize) - - *wa1.offset((i - 1 as c_int) as isize) - * *cc.offset((t3 - 1 as c_int) as isize); - *ch.offset(t6 as isize) = *cc.offset(t5 as isize) + ti2; - *ch.offset(t4 as isize) = ti2 - *cc.offset(t5 as isize); - *ch.offset((t6 - 1 as c_int) as isize) = - *cc.offset((t5 - 1 as c_int) as isize) + tr2; - *ch.offset((t4 - 1 as c_int) as isize) = - *cc.offset((t5 - 1 as c_int) as isize) - tr2; - i += 2 as c_int - } - t1 += ido; - t2 += ido; - k += 1 - } - if ido % 2 as c_int == 1 as c_int { - return; - } - } - t1 = ido; - t2 = t1 - 1 as c_int; - t3 = t2; - t2 += t0; - k = 0 as c_int; - while k < l1 { - *ch.offset(t1 as isize) = -*cc.offset(t2 as isize); - *ch.offset((t1 - 1 as c_int) as isize) = *cc.offset(t3 as isize); - t1 += ido << 1 as c_int; - t2 += ido; - t3 += ido; - k += 1 - } -} -unsafe extern "C" fn dradf4( - mut ido: c_int, - mut l1: c_int, - mut cc: *mut c_float, - mut ch: *mut c_float, - mut wa1: *mut c_float, - mut wa2: *mut c_float, - mut wa3: *mut c_float, -) { - static mut hsqt2: c_float = 0.70710678118654752f32; - let mut i: c_int = 0; - let mut k: c_int = 0; - let mut t0: c_int = 0; - let mut t1: c_int = 0; - let mut t2: c_int = 0; - let mut t3: c_int = 0; - let mut t4: c_int = 0; - let mut t5: c_int = 0; - let mut t6: c_int = 0; - let mut ci2: c_float = 0.; - let mut ci3: c_float = 0.; - let mut ci4: c_float = 0.; - let mut cr2: c_float = 0.; - let mut cr3: c_float = 0.; - let mut cr4: c_float = 0.; - let mut ti1: c_float = 0.; - let mut ti2: c_float = 0.; - let mut ti3: c_float = 0.; - let mut ti4: c_float = 0.; - let mut tr1: c_float = 0.; - let mut tr2: c_float = 0.; - let mut tr3: c_float = 0.; - let mut tr4: c_float = 0.; - t0 = l1 * ido; - t1 = t0; - t4 = t1 << 1 as c_int; - t2 = t1 + (t1 << 1 as c_int); - t3 = 0 as c_int; - k = 0 as c_int; - while k < l1 { - tr1 = *cc.offset(t1 as isize) + *cc.offset(t2 as isize); - tr2 = *cc.offset(t3 as isize) + *cc.offset(t4 as isize); - t5 = t3 << 2 as c_int; - *ch.offset(t5 as isize) = tr1 + tr2; - *ch.offset(((ido << 2 as c_int) + t5 - 1 as c_int) as isize) = - tr2 - tr1; - t5 += ido << 1 as c_int; - *ch.offset((t5 - 1 as c_int) as isize) = - *cc.offset(t3 as isize) - *cc.offset(t4 as isize); - *ch.offset(t5 as isize) = - *cc.offset(t2 as isize) - *cc.offset(t1 as isize); - t1 += ido; - t2 += ido; - t3 += ido; - t4 += ido; - k += 1 - } - if ido < 2 as c_int { - return; - } - if !(ido == 2 as c_int) { - t1 = 0 as c_int; - k = 0 as c_int; - while k < l1 { - t2 = t1; - t4 = t1 << 2 as c_int; - t6 = ido << 1 as c_int; - t5 = t6 + t4; - i = 2 as c_int; - while i < ido { - t2 += 2 as c_int; - t3 = t2; - t4 += 2 as c_int; - t5 -= 2 as c_int; - t3 += t0; - cr2 = *wa1.offset((i - 2 as c_int) as isize) - * *cc.offset((t3 - 1 as c_int) as isize) - + *wa1.offset((i - 1 as c_int) as isize) - * *cc.offset(t3 as isize); - ci2 = *wa1.offset((i - 2 as c_int) as isize) - * *cc.offset(t3 as isize) - - *wa1.offset((i - 1 as c_int) as isize) - * *cc.offset((t3 - 1 as c_int) as isize); - t3 += t0; - cr3 = *wa2.offset((i - 2 as c_int) as isize) - * *cc.offset((t3 - 1 as c_int) as isize) - + *wa2.offset((i - 1 as c_int) as isize) - * *cc.offset(t3 as isize); - ci3 = *wa2.offset((i - 2 as c_int) as isize) - * *cc.offset(t3 as isize) - - *wa2.offset((i - 1 as c_int) as isize) - * *cc.offset((t3 - 1 as c_int) as isize); - t3 += t0; - cr4 = *wa3.offset((i - 2 as c_int) as isize) - * *cc.offset((t3 - 1 as c_int) as isize) - + *wa3.offset((i - 1 as c_int) as isize) - * *cc.offset(t3 as isize); - ci4 = *wa3.offset((i - 2 as c_int) as isize) - * *cc.offset(t3 as isize) - - *wa3.offset((i - 1 as c_int) as isize) - * *cc.offset((t3 - 1 as c_int) as isize); - tr1 = cr2 + cr4; - tr4 = cr4 - cr2; - ti1 = ci2 + ci4; - ti4 = ci2 - ci4; - ti2 = *cc.offset(t2 as isize) + ci3; - ti3 = *cc.offset(t2 as isize) - ci3; - tr2 = *cc.offset((t2 - 1 as c_int) as isize) + cr3; - tr3 = *cc.offset((t2 - 1 as c_int) as isize) - cr3; - *ch.offset((t4 - 1 as c_int) as isize) = tr1 + tr2; - *ch.offset(t4 as isize) = ti1 + ti2; - *ch.offset((t5 - 1 as c_int) as isize) = tr3 - ti4; - *ch.offset(t5 as isize) = tr4 - ti3; - *ch.offset((t4 + t6 - 1 as c_int) as isize) = ti4 + tr3; - *ch.offset((t4 + t6) as isize) = tr4 + ti3; - *ch.offset((t5 + t6 - 1 as c_int) as isize) = tr2 - tr1; - *ch.offset((t5 + t6) as isize) = ti1 - ti2; - i += 2 as c_int - } - t1 += ido; - k += 1 - } - if ido & 1 as c_int != 0 { - return; - } - } - t1 = t0 + ido - 1 as c_int; - t2 = t1 + (t0 << 1 as c_int); - t3 = ido << 2 as c_int; - t4 = ido; - t5 = ido << 1 as c_int; - t6 = ido; - k = 0 as c_int; - while k < l1 { - ti1 = -hsqt2 * (*cc.offset(t1 as isize) + *cc.offset(t2 as isize)); - tr1 = hsqt2 * (*cc.offset(t1 as isize) - *cc.offset(t2 as isize)); - *ch.offset((t4 - 1 as c_int) as isize) = - tr1 + *cc.offset((t6 - 1 as c_int) as isize); - *ch.offset((t4 + t5 - 1 as c_int) as isize) = - *cc.offset((t6 - 1 as c_int) as isize) - tr1; - *ch.offset(t4 as isize) = ti1 - *cc.offset((t1 + t0) as isize); - *ch.offset((t4 + t5) as isize) = ti1 + *cc.offset((t1 + t0) as isize); - t1 += ido; - t2 += ido; - t4 += t3; - t6 += ido; - k += 1 - } -} -unsafe extern "C" fn dradfg( - mut ido: c_int, - mut ip: c_int, - mut l1: c_int, - mut idl1: c_int, - mut cc: *mut c_float, - mut c1: *mut c_float, - mut c2: *mut c_float, - mut ch: *mut c_float, - mut ch2: *mut c_float, - mut wa: *mut c_float, + +#[inline(always)] +fn drftf1_l102( + ip: i32, + ido: i32, + l1: i32, + idl1: i32, + iw: i32, + na: &mut i32, + c: &mut [f32], + ch: &mut [f32], + wa: &mut [f32], ) { - static mut tpi: c_float = 6.283185307179586f32; - let mut idij: c_int = 0; - let mut ipph: c_int = 0; - let mut i: c_int = 0; - let mut j: c_int = 0; - let mut k: c_int = 0; - let mut l: c_int = 0; - let mut ic: c_int = 0; - let mut ik: c_int = 0; - let mut is: c_int = 0; - let mut t0: c_int = 0; - let mut t1: c_int = 0; - let mut t2: c_int = 0; - let mut t3: c_int = 0; - let mut t4: c_int = 0; - let mut t5: c_int = 0; - let mut t6: c_int = 0; - let mut t7: c_int = 0; - let mut t8: c_int = 0; - let mut t9: c_int = 0; - let mut t10: c_int = 0; - let mut dc2: c_float = 0.; - let mut ai1: c_float = 0.; - let mut ai2: c_float = 0.; - let mut ar1: c_float = 0.; - let mut ar2: c_float = 0.; - let mut ds2: c_float = 0.; - let mut nbd: c_int = 0; - let mut dcp: c_float = 0.; - let mut arg: c_float = 0.; - let mut dsp: c_float = 0.; - let mut ar1h: c_float = 0.; - let mut ar2h: c_float = 0.; - let mut idp2: c_int = 0; - let mut ipp2: c_int = 0; - arg = tpi / ip as c_float; - dcp = f64::cos(arg as c_double) as c_float; - dsp = f64::sin(arg as c_double) as c_float; - ipph = ip + 1 as c_int >> 1 as c_int; - ipp2 = ip; - idp2 = ido; - nbd = ido - 1 as c_int >> 1 as c_int; - t0 = l1 * ido; - t10 = ip * ido; - if !(ido == 1 as c_int) { - ik = 0 as c_int; - while ik < idl1 { - *ch2.offset(ik as isize) = *c2.offset(ik as isize); - ik += 1 - } - t1 = 0 as c_int; - j = 1 as c_int; - while j < ip { - t1 += t0; - t2 = t1; - k = 0 as c_int; - while k < l1 { - *ch.offset(t2 as isize) = *c1.offset(t2 as isize); - t2 += ido; - k += 1 - } - j += 1 - } - is = -ido; - t1 = 0 as c_int; - if nbd > l1 { - j = 1 as c_int; - while j < ip { - t1 += t0; - is += ido; - t2 = -ido + t1; - k = 0 as c_int; - while k < l1 { - idij = is - 1 as c_int; - t2 += ido; - t3 = t2; - i = 2 as c_int; - while i < ido { - idij += 2 as c_int; - t3 += 2 as c_int; - *ch.offset((t3 - 1 as c_int) as isize) = *wa - .offset((idij - 1 as c_int) as isize) - * *c1.offset((t3 - 1 as c_int) as isize) - + *wa.offset(idij as isize) - * *c1.offset(t3 as isize); - *ch.offset(t3 as isize) = *wa - .offset((idij - 1 as c_int) as isize) - * *c1.offset(t3 as isize) - - *wa.offset(idij as isize) - * *c1.offset((t3 - 1 as c_int) as isize); - i += 2 as c_int - } - k += 1 - } - j += 1 + let iw_temp = iw as usize - 1; + if ip != 4 { + if ip != 2 { + if ido == 1 { + *na = 1 - *na; } - } else { - j = 1 as c_int; - while j < ip { - is += ido; - idij = is - 1 as c_int; - t1 += t0; - t2 = t1; - i = 2 as c_int; - while i < ido { - idij += 2 as c_int; - t2 += 2 as c_int; - t3 = t2; - k = 0 as c_int; - while k < l1 { - *ch.offset((t3 - 1 as c_int) as isize) = *wa - .offset((idij - 1 as c_int) as isize) - * *c1.offset((t3 - 1 as c_int) as isize) - + *wa.offset(idij as isize) - * *c1.offset(t3 as isize); - *ch.offset(t3 as isize) = *wa - .offset((idij - 1 as c_int) as isize) - * *c1.offset(t3 as isize) - - *wa.offset(idij as isize) - * *c1.offset((t3 - 1 as c_int) as isize); - t3 += ido; - k += 1 - } - i += 2 as c_int - } - j += 1 - } - } - t1 = 0 as c_int; - t2 = ipp2 * t0; - if nbd < l1 { - j = 1 as c_int; - while j < ipph { - t1 += t0; - t2 -= t0; - t3 = t1; - t4 = t2; - i = 2 as c_int; - while i < ido { - t3 += 2 as c_int; - t4 += 2 as c_int; - t5 = t3 - ido; - t6 = t4 - ido; - k = 0 as c_int; - while k < l1 { - t5 += ido; - t6 += ido; - *c1.offset((t5 - 1 as c_int) as isize) = *ch - .offset((t5 - 1 as c_int) as isize) - + *ch.offset((t6 - 1 as c_int) as isize); - *c1.offset((t6 - 1 as c_int) as isize) = - *ch.offset(t5 as isize) - *ch.offset(t6 as isize); - *c1.offset(t5 as isize) = - *ch.offset(t5 as isize) + *ch.offset(t6 as isize); - *c1.offset(t6 as isize) = *ch - .offset((t6 - 1 as c_int) as isize) - - *ch.offset((t5 - 1 as c_int) as isize); - k += 1 - } - i += 2 as c_int - } - j += 1 + if *na != 0 { + dradfg(ido, ip, l1, idl1, ch, c, &wa[iw_temp..]); + *na = 0; + } else { + dradfg(ido, ip, l1, idl1, c, ch, &wa[iw_temp..]); + *na = 1; } + } else if *na != 0 { + dradf2(ido, l1, ch, c, &mut wa[iw_temp..]); } else { - j = 1 as c_int; - while j < ipph { - t1 += t0; - t2 -= t0; - t3 = t1; - t4 = t2; - k = 0 as c_int; - while k < l1 { - t5 = t3; - t6 = t4; - i = 2 as c_int; - while i < ido { - t5 += 2 as c_int; - t6 += 2 as c_int; - *c1.offset((t5 - 1 as c_int) as isize) = *ch - .offset((t5 - 1 as c_int) as isize) - + *ch.offset((t6 - 1 as c_int) as isize); - *c1.offset((t6 - 1 as c_int) as isize) = - *ch.offset(t5 as isize) - *ch.offset(t6 as isize); - *c1.offset(t5 as isize) = - *ch.offset(t5 as isize) + *ch.offset(t6 as isize); - *c1.offset(t6 as isize) = *ch - .offset((t6 - 1 as c_int) as isize) - - *ch.offset((t5 - 1 as c_int) as isize); - i += 2 as c_int - } - t3 += ido; - t4 += ido; - k += 1 - } - j += 1 - } - } - } - ik = 0 as c_int; - while ik < idl1 { - *c2.offset(ik as isize) = *ch2.offset(ik as isize); - ik += 1 - } - t1 = 0 as c_int; - t2 = ipp2 * idl1; - j = 1 as c_int; - while j < ipph { - t1 += t0; - t2 -= t0; - t3 = t1 - ido; - t4 = t2 - ido; - k = 0 as c_int; - while k < l1 { - t3 += ido; - t4 += ido; - *c1.offset(t3 as isize) = - *ch.offset(t3 as isize) + *ch.offset(t4 as isize); - *c1.offset(t4 as isize) = - *ch.offset(t4 as isize) - *ch.offset(t3 as isize); - k += 1 - } - j += 1 - } - ar1 = 1.0f32; - ai1 = 0.0f32; - t1 = 0 as c_int; - t2 = ipp2 * idl1; - t3 = (ip - 1 as c_int) * idl1; - l = 1 as c_int; - while l < ipph { - t1 += idl1; - t2 -= idl1; - ar1h = dcp * ar1 - dsp * ai1; - ai1 = dcp * ai1 + dsp * ar1; - ar1 = ar1h; - t4 = t1; - t5 = t2; - t6 = t3; - t7 = idl1; - ik = 0 as c_int; - while ik < idl1 { - let fresh2 = t7; - t7 = t7 + 1; - let fresh3 = t4; - t4 = t4 + 1; - *ch2.offset(fresh3 as isize) = - *c2.offset(ik as isize) + ar1 * *c2.offset(fresh2 as isize); - let fresh4 = t6; - t6 = t6 + 1; - let fresh5 = t5; - t5 = t5 + 1; - *ch2.offset(fresh5 as isize) = ai1 * *c2.offset(fresh4 as isize); - ik += 1 - } - dc2 = ar1; - ds2 = ai1; - ar2 = ar1; - ai2 = ai1; - t4 = idl1; - t5 = (ipp2 - 1 as c_int) * idl1; - j = 2 as c_int; - while j < ipph { - t4 += idl1; - t5 -= idl1; - ar2h = dc2 * ar2 - ds2 * ai2; - ai2 = dc2 * ai2 + ds2 * ar2; - ar2 = ar2h; - t6 = t1; - t7 = t2; - t8 = t4; - t9 = t5; - ik = 0 as c_int; - while ik < idl1 { - let fresh6 = t8; - t8 = t8 + 1; - let fresh7 = t6; - t6 = t6 + 1; - *ch2.offset(fresh7 as isize) += - ar2 * *c2.offset(fresh6 as isize); - let fresh8 = t9; - t9 = t9 + 1; - let fresh9 = t7; - t7 = t7 + 1; - *ch2.offset(fresh9 as isize) += - ai2 * *c2.offset(fresh8 as isize); - ik += 1 - } - j += 1 - } - l += 1 - } - t1 = 0 as c_int; - j = 1 as c_int; - while j < ipph { - t1 += idl1; - t2 = t1; - ik = 0 as c_int; - while ik < idl1 { - let fresh10 = t2; - t2 = t2 + 1; - *ch2.offset(ik as isize) += *c2.offset(fresh10 as isize); - ik += 1 - } - j += 1 - } - if ido < l1 { - i = 0 as c_int; - while i < ido { - t1 = i; - t2 = i; - k = 0 as c_int; - while k < l1 { - *cc.offset(t2 as isize) = *ch.offset(t1 as isize); - t1 += ido; - t2 += t10; - k += 1 - } - i += 1 - } - } else { - t1 = 0 as c_int; - t2 = 0 as c_int; - k = 0 as c_int; - while k < l1 { - t3 = t1; - t4 = t2; - i = 0 as c_int; - while i < ido { - let fresh11 = t3; - t3 = t3 + 1; - let fresh12 = t4; - t4 = t4 + 1; - *cc.offset(fresh12 as isize) = *ch.offset(fresh11 as isize); - i += 1 - } - t1 += ido; - t2 += t10; - k += 1 - } - } - t1 = 0 as c_int; - t2 = ido << 1 as c_int; - t3 = 0 as c_int; - t4 = ipp2 * t0; - j = 1 as c_int; - while j < ipph { - t1 += t2; - t3 += t0; - t4 -= t0; - t5 = t1; - t6 = t3; - t7 = t4; - k = 0 as c_int; - while k < l1 { - *cc.offset((t5 - 1 as c_int) as isize) = *ch.offset(t6 as isize); - *cc.offset(t5 as isize) = *ch.offset(t7 as isize); - t5 += t10; - t6 += ido; - t7 += ido; - k += 1 + dradf2(ido, l1, c, ch, &mut wa[iw_temp..]); } - j += 1 } - if ido == 1 as c_int { - return; - } - if nbd < l1 { - t1 = -ido; - t3 = 0 as c_int; - t4 = 0 as c_int; - t5 = ipp2 * t0; - j = 1 as c_int; - while j < ipph { - t1 += t2; - t3 += t2; - t4 += t0; - t5 -= t0; - i = 2 as c_int; - while i < ido { - t6 = idp2 + t1 - i; - t7 = i + t3; - t8 = i + t4; - t9 = i + t5; - k = 0 as c_int; - while k < l1 { - *cc.offset((t7 - 1 as c_int) as isize) = *ch - .offset((t8 - 1 as c_int) as isize) - + *ch.offset((t9 - 1 as c_int) as isize); - *cc.offset((t6 - 1 as c_int) as isize) = *ch - .offset((t8 - 1 as c_int) as isize) - - *ch.offset((t9 - 1 as c_int) as isize); - *cc.offset(t7 as isize) = - *ch.offset(t8 as isize) + *ch.offset(t9 as isize); - *cc.offset(t6 as isize) = - *ch.offset(t9 as isize) - *ch.offset(t8 as isize); - t6 += t10; - t7 += t10; - t8 += ido; - t9 += ido; - k += 1 - } - i += 2 as c_int - } - j += 1 - } - return; - } else { - t1 = -ido; - t3 = 0 as c_int; - t4 = 0 as c_int; - t5 = ipp2 * t0; - j = 1 as c_int; - while j < ipph { - t1 += t2; - t3 += t2; - t4 += t0; - t5 -= t0; - t6 = t1; - t7 = t3; - t8 = t4; - t9 = t5; - k = 0 as c_int; - while k < l1 { - i = 2 as c_int; - while i < ido { - ic = idp2 - i; - *cc.offset((i + t7 - 1 as c_int) as isize) = *ch - .offset((i + t8 - 1 as c_int) as isize) - + *ch.offset((i + t9 - 1 as c_int) as isize); - *cc.offset((ic + t6 - 1 as c_int) as isize) = *ch - .offset((i + t8 - 1 as c_int) as isize) - - *ch.offset((i + t9 - 1 as c_int) as isize); - *cc.offset((i + t7) as isize) = *ch - .offset((i + t8) as isize) - + *ch.offset((i + t9) as isize); - *cc.offset((ic + t6) as isize) = *ch - .offset((i + t9) as isize) - - *ch.offset((i + t8) as isize); - i += 2 as c_int - } - t6 += t10; - t7 += t10; - t8 += ido; - t9 += ido; - k += 1 - } - j += 1 - } - return; - }; } -unsafe extern "C" fn drftf1( - mut n: c_int, - mut c: *mut c_float, - mut ch: *mut c_float, - mut wa: *mut c_float, - mut ifac: *mut c_int, + +fn drftf1( + n: i32, + c: &mut [f32], + ch: &mut [f32], + wa: &mut [f32], + ifac: &mut [i32], ) { - let mut i: c_int = 0; - let mut k1: c_int = 0; - let mut l1: c_int = 0; - let mut l2: c_int = 0; - let mut na: c_int = 0; - let mut kh: c_int = 0; - let mut nf: c_int = 0; - let mut ip: c_int = 0; - let mut iw: c_int = 0; - let mut ido: c_int = 0; - let mut idl1: c_int = 0; - let mut ix2: c_int = 0; - let mut ix3: c_int = 0; - nf = *ifac.offset(1 as c_int as isize); - na = 1 as c_int; - l2 = n; - iw = n; - k1 = 0 as c_int; - while k1 < nf { - kh = nf - k1; - ip = *ifac.offset((kh + 1 as c_int) as isize); - l1 = l2 / ip; - ido = n / l2; - idl1 = ido * l1; - iw -= (ip - 1 as c_int) * ido; - na = 1 as c_int - na; - if ip != 4 as c_int { - if ip != 2 as c_int { - if ido == 1 as c_int { - na = 1 as c_int - na - } - if na != 0 as c_int { - dradfg( - ido, - ip, - l1, - idl1, - ch, - ch, - ch, - c, - c, - wa.offset(iw as isize).offset(-(1 as c_int as isize)), - ); - na = 0 as c_int - } else { - dradfg( - ido, - ip, - l1, - idl1, - c, - c, - c, - ch, - ch, - wa.offset(iw as isize).offset(-(1 as c_int as isize)), - ); - na = 1 as c_int - } - } else if na != 0 as c_int { - dradf2( - ido, - l1, - ch, - c, - wa.offset(iw as isize).offset(-(1 as c_int as isize)), - ); - } else { - dradf2( - ido, - l1, - c, - ch, - wa.offset(iw as isize).offset(-(1 as c_int as isize)), - ); - } + let nf = ifac[1]; + let mut na = 1; + let mut l2 = n; + let mut iw = n; + + for k1 in 0..nf { + let kh = nf - k1; + let ip = ifac[kh as usize + 1]; + let l1 = l2 / ip; + let ido = n / l2; + let idl1 = ido * l1; + iw -= (ip - 1) * ido; + na = 1 - na; + + if ip != 4 { + drftf1_l102(ip, ido, l1, idl1, iw, &mut na, c, ch, wa); } else { - ix2 = iw + ido; - ix3 = ix2 + ido; - if na != 0 as c_int { + let ix2 = iw + ido; + let ix3 = ix2 + ido; + if na != 0 { dradf4( ido, l1, ch, c, - wa.offset(iw as isize).offset(-(1 as c_int as isize)), - wa.offset(ix2 as isize).offset(-(1 as c_int as isize)), - wa.offset(ix3 as isize).offset(-(1 as c_int as isize)), + &wa[(iw as usize - 1)..], + &wa[(ix2 as usize - 1)..], + &wa[(ix3 as usize - 1)..], ); } else { dradf4( @@ -1002,1062 +194,173 @@ unsafe extern "C" fn drftf1( l1, c, ch, - wa.offset(iw as isize).offset(-(1 as c_int as isize)), - wa.offset(ix2 as isize).offset(-(1 as c_int as isize)), - wa.offset(ix3 as isize).offset(-(1 as c_int as isize)), + &wa[(iw as usize - 1)..], + &wa[(ix2 as usize - 1)..], + &wa[(ix3 as usize - 1)..], ); } } + l2 = l1; - k1 += 1 } - if na == 1 as c_int { - return; - } - i = 0 as c_int; - while i < n { - *c.offset(i as isize) = *ch.offset(i as isize); - i += 1 - } -} -unsafe extern "C" fn dradb2( - mut ido: c_int, - mut l1: c_int, - mut cc: *mut c_float, - mut ch: *mut c_float, - mut wa1: *mut c_float, -) { - let mut i: c_int = 0; - let mut k: c_int = 0; - let mut t0: c_int = 0; - let mut t1: c_int = 0; - let mut t2: c_int = 0; - let mut t3: c_int = 0; - let mut t4: c_int = 0; - let mut t5: c_int = 0; - let mut t6: c_int = 0; - let mut ti2: c_float = 0.; - let mut tr2: c_float = 0.; - t0 = l1 * ido; - t1 = 0 as c_int; - t2 = 0 as c_int; - t3 = (ido << 1 as c_int) - 1 as c_int; - k = 0 as c_int; - while k < l1 { - *ch.offset(t1 as isize) = - *cc.offset(t2 as isize) + *cc.offset((t3 + t2) as isize); - *ch.offset((t1 + t0) as isize) = - *cc.offset(t2 as isize) - *cc.offset((t3 + t2) as isize); - t1 += ido; - t2 = t1 << 1 as c_int; - k += 1 - } - if ido < 2 as c_int { - return; - } - if !(ido == 2 as c_int) { - t1 = 0 as c_int; - t2 = 0 as c_int; - k = 0 as c_int; - while k < l1 { - t3 = t1; - t4 = t2; - t5 = t4 + (ido << 1 as c_int); - t6 = t0 + t1; - i = 2 as c_int; - while i < ido { - t3 += 2 as c_int; - t4 += 2 as c_int; - t5 -= 2 as c_int; - t6 += 2 as c_int; - *ch.offset((t3 - 1 as c_int) as isize) = *cc - .offset((t4 - 1 as c_int) as isize) - + *cc.offset((t5 - 1 as c_int) as isize); - tr2 = *cc.offset((t4 - 1 as c_int) as isize) - - *cc.offset((t5 - 1 as c_int) as isize); - *ch.offset(t3 as isize) = - *cc.offset(t4 as isize) - *cc.offset(t5 as isize); - ti2 = *cc.offset(t4 as isize) + *cc.offset(t5 as isize); - *ch.offset((t6 - 1 as c_int) as isize) = - *wa1.offset((i - 2 as c_int) as isize) * tr2 - - *wa1.offset((i - 1 as c_int) as isize) * ti2; - *ch.offset(t6 as isize) = - *wa1.offset((i - 2 as c_int) as isize) * ti2 - + *wa1.offset((i - 1 as c_int) as isize) * tr2; - i += 2 as c_int - } - t1 += ido; - t2 = t1 << 1 as c_int; - k += 1 - } - if ido % 2 as c_int == 1 as c_int { - return; - } - } - t1 = ido - 1 as c_int; - t2 = ido - 1 as c_int; - k = 0 as c_int; - while k < l1 { - *ch.offset(t1 as isize) = - *cc.offset(t2 as isize) + *cc.offset(t2 as isize); - *ch.offset((t1 + t0) as isize) = -(*cc - .offset((t2 + 1 as c_int) as isize) - + *cc.offset((t2 + 1 as c_int) as isize)); - t1 += ido; - t2 += ido << 1 as c_int; - k += 1 - } -} -unsafe extern "C" fn dradb3( - mut ido: c_int, - mut l1: c_int, - mut cc: *mut c_float, - mut ch: *mut c_float, - mut wa1: *mut c_float, - mut wa2: *mut c_float, -) { - static mut taur: c_float = -0.5f32; - static mut taui: c_float = 0.8660254037844386f32; - let mut i: c_int = 0; - let mut k: c_int = 0; - let mut t0: c_int = 0; - let mut t1: c_int = 0; - let mut t2: c_int = 0; - let mut t3: c_int = 0; - let mut t4: c_int = 0; - let mut t5: c_int = 0; - let mut t6: c_int = 0; - let mut t7: c_int = 0; - let mut t8: c_int = 0; - let mut t9: c_int = 0; - let mut t10: c_int = 0; - let mut ci2: c_float = 0.; - let mut ci3: c_float = 0.; - let mut di2: c_float = 0.; - let mut di3: c_float = 0.; - let mut cr2: c_float = 0.; - let mut cr3: c_float = 0.; - let mut dr2: c_float = 0.; - let mut dr3: c_float = 0.; - let mut ti2: c_float = 0.; - let mut tr2: c_float = 0.; - t0 = l1 * ido; - t1 = 0 as c_int; - t2 = t0 << 1 as c_int; - t3 = ido << 1 as c_int; - t4 = ido + (ido << 1 as c_int); - t5 = 0 as c_int; - k = 0 as c_int; - while k < l1 { - tr2 = *cc.offset((t3 - 1 as c_int) as isize) - + *cc.offset((t3 - 1 as c_int) as isize); - cr2 = *cc.offset(t5 as isize) + taur * tr2; - *ch.offset(t1 as isize) = *cc.offset(t5 as isize) + tr2; - ci3 = taui * (*cc.offset(t3 as isize) + *cc.offset(t3 as isize)); - *ch.offset((t1 + t0) as isize) = cr2 - ci3; - *ch.offset((t1 + t2) as isize) = cr2 + ci3; - t1 += ido; - t3 += t4; - t5 += t4; - k += 1 - } - if ido == 1 as c_int { - return; - } - t1 = 0 as c_int; - t3 = ido << 1 as c_int; - k = 0 as c_int; - while k < l1 { - t7 = t1 + (t1 << 1 as c_int); - t5 = t7 + t3; - t6 = t5; - t8 = t1; - t9 = t1 + t0; - t10 = t9 + t0; - i = 2 as c_int; - while i < ido { - t5 += 2 as c_int; - t6 -= 2 as c_int; - t7 += 2 as c_int; - t8 += 2 as c_int; - t9 += 2 as c_int; - t10 += 2 as c_int; - tr2 = *cc.offset((t5 - 1 as c_int) as isize) - + *cc.offset((t6 - 1 as c_int) as isize); - cr2 = *cc.offset((t7 - 1 as c_int) as isize) + taur * tr2; - *ch.offset((t8 - 1 as c_int) as isize) = - *cc.offset((t7 - 1 as c_int) as isize) + tr2; - ti2 = *cc.offset(t5 as isize) - *cc.offset(t6 as isize); - ci2 = *cc.offset(t7 as isize) + taur * ti2; - *ch.offset(t8 as isize) = *cc.offset(t7 as isize) + ti2; - cr3 = taui - * (*cc.offset((t5 - 1 as c_int) as isize) - - *cc.offset((t6 - 1 as c_int) as isize)); - ci3 = taui * (*cc.offset(t5 as isize) + *cc.offset(t6 as isize)); - dr2 = cr2 - ci3; - dr3 = cr2 + ci3; - di2 = ci2 + cr3; - di3 = ci2 - cr3; - *ch.offset((t9 - 1 as c_int) as isize) = - *wa1.offset((i - 2 as c_int) as isize) * dr2 - - *wa1.offset((i - 1 as c_int) as isize) * di2; - *ch.offset(t9 as isize) = *wa1.offset((i - 2 as c_int) as isize) - * di2 - + *wa1.offset((i - 1 as c_int) as isize) * dr2; - *ch.offset((t10 - 1 as c_int) as isize) = - *wa2.offset((i - 2 as c_int) as isize) * dr3 - - *wa2.offset((i - 1 as c_int) as isize) * di3; - *ch.offset(t10 as isize) = *wa2.offset((i - 2 as c_int) as isize) - * di3 - + *wa2.offset((i - 1 as c_int) as isize) * dr3; - i += 2 as c_int - } - t1 += ido; - k += 1 - } -} -unsafe extern "C" fn dradb4( - mut ido: c_int, - mut l1: c_int, - mut cc: *mut c_float, - mut ch: *mut c_float, - mut wa1: *mut c_float, - mut wa2: *mut c_float, - mut wa3: *mut c_float, -) { - static mut sqrt2: c_float = std::f32::consts::SQRT_2; - let mut i: c_int = 0; - let mut k: c_int = 0; - let mut t0: c_int = 0; - let mut t1: c_int = 0; - let mut t2: c_int = 0; - let mut t3: c_int = 0; - let mut t4: c_int = 0; - let mut t5: c_int = 0; - let mut t6: c_int = 0; - let mut t7: c_int = 0; - let mut t8: c_int = 0; - let mut ci2: c_float = 0.; - let mut ci3: c_float = 0.; - let mut ci4: c_float = 0.; - let mut cr2: c_float = 0.; - let mut cr3: c_float = 0.; - let mut cr4: c_float = 0.; - let mut ti1: c_float = 0.; - let mut ti2: c_float = 0.; - let mut ti3: c_float = 0.; - let mut ti4: c_float = 0.; - let mut tr1: c_float = 0.; - let mut tr2: c_float = 0.; - let mut tr3: c_float = 0.; - let mut tr4: c_float = 0.; - t0 = l1 * ido; - t1 = 0 as c_int; - t2 = ido << 2 as c_int; - t3 = 0 as c_int; - t6 = ido << 1 as c_int; - k = 0 as c_int; - while k < l1 { - t4 = t3 + t6; - t5 = t1; - tr3 = *cc.offset((t4 - 1 as c_int) as isize) - + *cc.offset((t4 - 1 as c_int) as isize); - tr4 = *cc.offset(t4 as isize) + *cc.offset(t4 as isize); - t4 += t6; - tr1 = *cc.offset(t3 as isize) - *cc.offset((t4 - 1 as c_int) as isize); - tr2 = *cc.offset(t3 as isize) + *cc.offset((t4 - 1 as c_int) as isize); - *ch.offset(t5 as isize) = tr2 + tr3; - t5 += t0; - *ch.offset(t5 as isize) = tr1 - tr4; - t5 += t0; - *ch.offset(t5 as isize) = tr2 - tr3; - t5 += t0; - *ch.offset(t5 as isize) = tr1 + tr4; - t1 += ido; - t3 += t2; - k += 1 - } - if ido < 2 as c_int { + + if na == 1 { return; } - if !(ido == 2 as c_int) { - t1 = 0 as c_int; - k = 0 as c_int; - while k < l1 { - t2 = t1 << 2 as c_int; - t3 = t2 + t6; - t4 = t3; - t5 = t4 + t6; - t7 = t1; - i = 2 as c_int; - while i < ido { - t2 += 2 as c_int; - t3 += 2 as c_int; - t4 -= 2 as c_int; - t5 -= 2 as c_int; - t7 += 2 as c_int; - ti1 = *cc.offset(t2 as isize) + *cc.offset(t5 as isize); - ti2 = *cc.offset(t2 as isize) - *cc.offset(t5 as isize); - ti3 = *cc.offset(t3 as isize) - *cc.offset(t4 as isize); - tr4 = *cc.offset(t3 as isize) + *cc.offset(t4 as isize); - tr1 = *cc.offset((t2 - 1 as c_int) as isize) - - *cc.offset((t5 - 1 as c_int) as isize); - tr2 = *cc.offset((t2 - 1 as c_int) as isize) - + *cc.offset((t5 - 1 as c_int) as isize); - ti4 = *cc.offset((t3 - 1 as c_int) as isize) - - *cc.offset((t4 - 1 as c_int) as isize); - tr3 = *cc.offset((t3 - 1 as c_int) as isize) - + *cc.offset((t4 - 1 as c_int) as isize); - *ch.offset((t7 - 1 as c_int) as isize) = tr2 + tr3; - cr3 = tr2 - tr3; - *ch.offset(t7 as isize) = ti2 + ti3; - ci3 = ti2 - ti3; - cr2 = tr1 - tr4; - cr4 = tr1 + tr4; - ci2 = ti1 + ti4; - ci4 = ti1 - ti4; - t8 = t7 + t0; - *ch.offset((t8 - 1 as c_int) as isize) = - *wa1.offset((i - 2 as c_int) as isize) * cr2 - - *wa1.offset((i - 1 as c_int) as isize) * ci2; - *ch.offset(t8 as isize) = - *wa1.offset((i - 2 as c_int) as isize) * ci2 - + *wa1.offset((i - 1 as c_int) as isize) * cr2; - t8 += t0; - *ch.offset((t8 - 1 as c_int) as isize) = - *wa2.offset((i - 2 as c_int) as isize) * cr3 - - *wa2.offset((i - 1 as c_int) as isize) * ci3; - *ch.offset(t8 as isize) = - *wa2.offset((i - 2 as c_int) as isize) * ci3 - + *wa2.offset((i - 1 as c_int) as isize) * cr3; - t8 += t0; - *ch.offset((t8 - 1 as c_int) as isize) = - *wa3.offset((i - 2 as c_int) as isize) * cr4 - - *wa3.offset((i - 1 as c_int) as isize) * ci4; - *ch.offset(t8 as isize) = - *wa3.offset((i - 2 as c_int) as isize) * ci4 - + *wa3.offset((i - 1 as c_int) as isize) * cr4; - i += 2 as c_int - } - t1 += ido; - k += 1 - } - if ido % 2 as c_int == 1 as c_int { - return; - } - } - t1 = ido; - t2 = ido << 2 as c_int; - t3 = ido - 1 as c_int; - t4 = ido + (ido << 1 as c_int); - k = 0 as c_int; - while k < l1 { - t5 = t3; - ti1 = *cc.offset(t1 as isize) + *cc.offset(t4 as isize); - ti2 = *cc.offset(t4 as isize) - *cc.offset(t1 as isize); - tr1 = *cc.offset((t1 - 1 as c_int) as isize) - - *cc.offset((t4 - 1 as c_int) as isize); - tr2 = *cc.offset((t1 - 1 as c_int) as isize) - + *cc.offset((t4 - 1 as c_int) as isize); - *ch.offset(t5 as isize) = tr2 + tr2; - t5 += t0; - *ch.offset(t5 as isize) = sqrt2 * (tr1 - ti1); - t5 += t0; - *ch.offset(t5 as isize) = ti2 + ti2; - t5 += t0; - *ch.offset(t5 as isize) = -sqrt2 * (tr1 + ti1); - t3 += ido; - t1 += t2; - t4 += t2; - k += 1 - } + + c[..n as usize].copy_from_slice(&ch[..n as usize]); } -unsafe extern "C" fn dradbg( - mut ido: c_int, - mut ip: c_int, - mut l1: c_int, - mut idl1: c_int, - mut cc: *mut c_float, - mut c1: *mut c_float, - mut c2: *mut c_float, - mut ch: *mut c_float, - mut ch2: *mut c_float, - mut wa: *mut c_float, + +#[inline(always)] +fn drftb1_l102( + ip: i32, + ido: i32, + l1: i32, + idl1: i32, + iw: i32, + na: &mut i32, + c: &mut [f32], + ch: &mut [f32], + wa: &mut [f32], ) { - static mut tpi: c_float = 6.283185307179586f32; - let mut idij: c_int = 0; - let mut ipph: c_int = 0; - let mut i: c_int = 0; - let mut j: c_int = 0; - let mut k: c_int = 0; - let mut l: c_int = 0; - let mut ik: c_int = 0; - let mut is: c_int = 0; - let mut t0: c_int = 0; - let mut t1: c_int = 0; - let mut t2: c_int = 0; - let mut t3: c_int = 0; - let mut t4: c_int = 0; - let mut t5: c_int = 0; - let mut t6: c_int = 0; - let mut t7: c_int = 0; - let mut t8: c_int = 0; - let mut t9: c_int = 0; - let mut t10: c_int = 0; - let mut t11: c_int = 0; - let mut t12: c_int = 0; - let mut dc2: c_float = 0.; - let mut ai1: c_float = 0.; - let mut ai2: c_float = 0.; - let mut ar1: c_float = 0.; - let mut ar2: c_float = 0.; - let mut ds2: c_float = 0.; - let mut nbd: c_int = 0; - let mut dcp: c_float = 0.; - let mut arg: c_float = 0.; - let mut dsp: c_float = 0.; - let mut ar1h: c_float = 0.; - let mut ar2h: c_float = 0.; - let mut ipp2: c_int = 0; - t10 = ip * ido; - t0 = l1 * ido; - arg = tpi / ip as c_float; - dcp = f64::cos(arg as c_double) as c_float; - dsp = f64::sin(arg as c_double) as c_float; - nbd = ido - 1 as c_int >> 1 as c_int; - ipp2 = ip; - ipph = ip + 1 as c_int >> 1 as c_int; - if ido < l1 { - t1 = 0 as c_int; - i = 0 as c_int; - while i < ido { - t2 = t1; - t3 = t1; - k = 0 as c_int; - while k < l1 { - *ch.offset(t2 as isize) = *cc.offset(t3 as isize); - t2 += ido; - t3 += t10; - k += 1 - } - t1 += 1; - i += 1 - } - } else { - t1 = 0 as c_int; - t2 = 0 as c_int; - k = 0 as c_int; - while k < l1 { - t3 = t1; - t4 = t2; - i = 0 as c_int; - while i < ido { - *ch.offset(t3 as isize) = *cc.offset(t4 as isize); - t3 += 1; - t4 += 1; - i += 1 + let iw_temp = iw as usize - 1; + if ip != 2 { + if ip != 3 { + if *na != 0 { + dradbg(ido, ip, l1, idl1, ch, c, &wa[iw_temp..]); + } else { + dradbg(ido, ip, l1, idl1, c, ch, &wa[iw_temp..]); } - t1 += ido; - t2 += t10; - k += 1 - } - } - t1 = 0 as c_int; - t2 = ipp2 * t0; - t5 = ido << 1 as c_int; - t7 = t5; - j = 1 as c_int; - while j < ipph { - t1 += t0; - t2 -= t0; - t3 = t1; - t4 = t2; - t6 = t5; - k = 0 as c_int; - while k < l1 { - *ch.offset(t3 as isize) = *cc.offset((t6 - 1 as c_int) as isize) - + *cc.offset((t6 - 1 as c_int) as isize); - *ch.offset(t4 as isize) = - *cc.offset(t6 as isize) + *cc.offset(t6 as isize); - t3 += ido; - t4 += ido; - t6 += t10; - k += 1 - } - t5 += t7; - j += 1 - } - if !(ido == 1 as c_int) { - if nbd < l1 { - t1 = 0 as c_int; - t2 = ipp2 * t0; - t7 = 0 as c_int; - j = 1 as c_int; - while j < ipph { - t1 += t0; - t2 -= t0; - t3 = t1; - t4 = t2; - t7 += ido << 1 as c_int; - t8 = t7; - t9 = t7; - i = 2 as c_int; - while i < ido { - t3 += 2 as c_int; - t4 += 2 as c_int; - t8 += 2 as c_int; - t9 -= 2 as c_int; - t5 = t3; - t6 = t4; - t11 = t8; - t12 = t9; - k = 0 as c_int; - while k < l1 { - *ch.offset((t5 - 1 as c_int) as isize) = *cc - .offset((t11 - 1 as c_int) as isize) - + *cc.offset((t12 - 1 as c_int) as isize); - *ch.offset((t6 - 1 as c_int) as isize) = *cc - .offset((t11 - 1 as c_int) as isize) - - *cc.offset((t12 - 1 as c_int) as isize); - *ch.offset(t5 as isize) = *cc.offset(t11 as isize) - - *cc.offset(t12 as isize); - *ch.offset(t6 as isize) = *cc.offset(t11 as isize) - + *cc.offset(t12 as isize); - t5 += ido; - t6 += ido; - t11 += t10; - t12 += t10; - k += 1 - } - i += 2 as c_int - } - j += 1 + + if ido == 1 { + *na = 1 - *na; } } else { - t1 = 0 as c_int; - t2 = ipp2 * t0; - t7 = 0 as c_int; - j = 1 as c_int; - while j < ipph { - t1 += t0; - t2 -= t0; - t3 = t1; - t4 = t2; - t7 += ido << 1 as c_int; - t8 = t7; - k = 0 as c_int; - while k < l1 { - t5 = t3; - t6 = t4; - t9 = t8; - t11 = t8; - i = 2 as c_int; - while i < ido { - t5 += 2 as c_int; - t6 += 2 as c_int; - t9 += 2 as c_int; - t11 -= 2 as c_int; - *ch.offset((t5 - 1 as c_int) as isize) = *cc - .offset((t9 - 1 as c_int) as isize) - + *cc.offset((t11 - 1 as c_int) as isize); - *ch.offset((t6 - 1 as c_int) as isize) = *cc - .offset((t9 - 1 as c_int) as isize) - - *cc.offset((t11 - 1 as c_int) as isize); - *ch.offset(t5 as isize) = - *cc.offset(t9 as isize) - *cc.offset(t11 as isize); - *ch.offset(t6 as isize) = - *cc.offset(t9 as isize) + *cc.offset(t11 as isize); - i += 2 as c_int - } - t3 += ido; - t4 += ido; - t8 += t10; - k += 1 - } - j += 1 - } - } - } - ar1 = 1.0f32; - ai1 = 0.0f32; - t1 = 0 as c_int; - t2 = ipp2 * idl1; - t9 = t2; - t3 = (ip - 1 as c_int) * idl1; - l = 1 as c_int; - while l < ipph { - t1 += idl1; - t2 -= idl1; - ar1h = dcp * ar1 - dsp * ai1; - ai1 = dcp * ai1 + dsp * ar1; - ar1 = ar1h; - t4 = t1; - t5 = t2; - t6 = 0 as c_int; - t7 = idl1; - t8 = t3; - ik = 0 as c_int; - while ik < idl1 { - let fresh13 = t6; - t6 = t6 + 1; - let fresh14 = t7; - t7 = t7 + 1; - let fresh15 = t4; - t4 = t4 + 1; - *c2.offset(fresh15 as isize) = *ch2.offset(fresh13 as isize) - + ar1 * *ch2.offset(fresh14 as isize); - let fresh16 = t8; - t8 = t8 + 1; - let fresh17 = t5; - t5 = t5 + 1; - *c2.offset(fresh17 as isize) = ai1 * *ch2.offset(fresh16 as isize); - ik += 1 - } - dc2 = ar1; - ds2 = ai1; - ar2 = ar1; - ai2 = ai1; - t6 = idl1; - t7 = t9 - idl1; - j = 2 as c_int; - while j < ipph { - t6 += idl1; - t7 -= idl1; - ar2h = dc2 * ar2 - ds2 * ai2; - ai2 = dc2 * ai2 + ds2 * ar2; - ar2 = ar2h; - t4 = t1; - t5 = t2; - t11 = t6; - t12 = t7; - ik = 0 as c_int; - while ik < idl1 { - let fresh18 = t11; - t11 = t11 + 1; - let fresh19 = t4; - t4 = t4 + 1; - *c2.offset(fresh19 as isize) += - ar2 * *ch2.offset(fresh18 as isize); - let fresh20 = t12; - t12 = t12 + 1; - let fresh21 = t5; - t5 = t5 + 1; - *c2.offset(fresh21 as isize) += - ai2 * *ch2.offset(fresh20 as isize); - ik += 1 + let ix2 = iw + ido - 1; + if *na != 0 { + dradb3(ido, l1, &ch, c, &wa[iw_temp..], &wa[ix2 as usize..]); + } else { + dradb3(ido, l1, &c, ch, &wa[iw_temp..], &wa[ix2 as usize..]); } - j += 1 + *na = 1 - *na; } - l += 1 - } - t1 = 0 as c_int; - j = 1 as c_int; - while j < ipph { - t1 += idl1; - t2 = t1; - ik = 0 as c_int; - while ik < idl1 { - let fresh22 = t2; - t2 = t2 + 1; - *ch2.offset(ik as isize) += *ch2.offset(fresh22 as isize); - ik += 1 - } - j += 1 - } - t1 = 0 as c_int; - t2 = ipp2 * t0; - j = 1 as c_int; - while j < ipph { - t1 += t0; - t2 -= t0; - t3 = t1; - t4 = t2; - k = 0 as c_int; - while k < l1 { - *ch.offset(t3 as isize) = - *c1.offset(t3 as isize) - *c1.offset(t4 as isize); - *ch.offset(t4 as isize) = - *c1.offset(t3 as isize) + *c1.offset(t4 as isize); - t3 += ido; - t4 += ido; - k += 1 - } - j += 1 - } - if !(ido == 1 as c_int) { - if nbd < l1 { - t1 = 0 as c_int; - t2 = ipp2 * t0; - j = 1 as c_int; - while j < ipph { - t1 += t0; - t2 -= t0; - t3 = t1; - t4 = t2; - i = 2 as c_int; - while i < ido { - t3 += 2 as c_int; - t4 += 2 as c_int; - t5 = t3; - t6 = t4; - k = 0 as c_int; - while k < l1 { - *ch.offset((t5 - 1 as c_int) as isize) = *c1 - .offset((t5 - 1 as c_int) as isize) - - *c1.offset(t6 as isize); - *ch.offset((t6 - 1 as c_int) as isize) = *c1 - .offset((t5 - 1 as c_int) as isize) - + *c1.offset(t6 as isize); - *ch.offset(t5 as isize) = *c1.offset(t5 as isize) - + *c1.offset((t6 - 1 as c_int) as isize); - *ch.offset(t6 as isize) = *c1.offset(t5 as isize) - - *c1.offset((t6 - 1 as c_int) as isize); - t5 += ido; - t6 += ido; - k += 1 - } - i += 2 as c_int - } - j += 1 - } + } else { + if *na != 0 { + dradb2(ido, l1, &ch, c, &wa[iw_temp..]); } else { - t1 = 0 as c_int; - t2 = ipp2 * t0; - j = 1 as c_int; - while j < ipph { - t1 += t0; - t2 -= t0; - t3 = t1; - t4 = t2; - k = 0 as c_int; - while k < l1 { - t5 = t3; - t6 = t4; - i = 2 as c_int; - while i < ido { - t5 += 2 as c_int; - t6 += 2 as c_int; - *ch.offset((t5 - 1 as c_int) as isize) = *c1 - .offset((t5 - 1 as c_int) as isize) - - *c1.offset(t6 as isize); - *ch.offset((t6 - 1 as c_int) as isize) = *c1 - .offset((t5 - 1 as c_int) as isize) - + *c1.offset(t6 as isize); - *ch.offset(t5 as isize) = *c1.offset(t5 as isize) - + *c1.offset((t6 - 1 as c_int) as isize); - *ch.offset(t6 as isize) = *c1.offset(t5 as isize) - - *c1.offset((t6 - 1 as c_int) as isize); - i += 2 as c_int - } - t3 += ido; - t4 += ido; - k += 1 - } - j += 1 - } - } - } - if ido == 1 as c_int { - return; - } - ik = 0 as c_int; - while ik < idl1 { - *c2.offset(ik as isize) = *ch2.offset(ik as isize); - ik += 1 - } - t1 = 0 as c_int; - j = 1 as c_int; - while j < ip { - t1 += t0; - t2 = t1; - k = 0 as c_int; - while k < l1 { - *c1.offset(t2 as isize) = *ch.offset(t2 as isize); - t2 += ido; - k += 1 + dradb2(ido, l1, &c, ch, &wa[iw_temp..]); } - j += 1 + *na = 1 - *na; } - if nbd > l1 { - is = -ido - 1 as c_int; - t1 = 0 as c_int; - j = 1 as c_int; - while j < ip { - is += ido; - t1 += t0; - t2 = t1; - k = 0 as c_int; - while k < l1 { - idij = is; - t3 = t2; - i = 2 as c_int; - while i < ido { - idij += 2 as c_int; - t3 += 2 as c_int; - *c1.offset((t3 - 1 as c_int) as isize) = *wa - .offset((idij - 1 as c_int) as isize) - * *ch.offset((t3 - 1 as c_int) as isize) - - *wa.offset(idij as isize) * *ch.offset(t3 as isize); - *c1.offset(t3 as isize) = *wa - .offset((idij - 1 as c_int) as isize) - * *ch.offset(t3 as isize) - + *wa.offset(idij as isize) - * *ch.offset((t3 - 1 as c_int) as isize); - i += 2 as c_int - } - t2 += ido; - k += 1 - } - j += 1 - } - return; - } else { - is = -ido - 1 as c_int; - t1 = 0 as c_int; - j = 1 as c_int; - while j < ip { - is += ido; - t1 += t0; - idij = is; - t2 = t1; - i = 2 as c_int; - while i < ido { - t2 += 2 as c_int; - idij += 2 as c_int; - t3 = t2; - k = 0 as c_int; - while k < l1 { - *c1.offset((t3 - 1 as c_int) as isize) = *wa - .offset((idij - 1 as c_int) as isize) - * *ch.offset((t3 - 1 as c_int) as isize) - - *wa.offset(idij as isize) * *ch.offset(t3 as isize); - *c1.offset(t3 as isize) = *wa - .offset((idij - 1 as c_int) as isize) - * *ch.offset(t3 as isize) - + *wa.offset(idij as isize) - * *ch.offset((t3 - 1 as c_int) as isize); - t3 += ido; - k += 1 - } - i += 2 as c_int - } - j += 1 - } - return; - }; } -unsafe extern "C" fn drftb1( - mut n: c_int, - mut c: *mut c_float, - mut ch: *mut c_float, - mut wa: *mut c_float, - mut ifac: *mut c_int, + +fn drftb1( + n: i32, + c: &mut [f32], + ch: &mut [f32], + wa: &mut [f32], + ifac: &mut [i32], ) { - let mut i: c_int = 0; - let mut k1: c_int = 0; - let mut l1: c_int = 0; - let mut l2: c_int = 0; - let mut na: c_int = 0; - let mut nf: c_int = 0; - let mut ip: c_int = 0; - let mut iw: c_int = 0; - let mut ix2: c_int = 0; - let mut ix3: c_int = 0; - let mut ido: c_int = 0; - let mut idl1: c_int = 0; - nf = *ifac.offset(1 as c_int as isize); - na = 0 as c_int; - l1 = 1 as c_int; - iw = 1 as c_int; - k1 = 0 as c_int; - while k1 < nf { - ip = *ifac.offset((k1 + 2 as c_int) as isize); - l2 = ip * l1; - ido = n / l2; - idl1 = ido * l1; - if ip != 4 as c_int { - if ip != 2 as c_int { - if ip != 3 as c_int { - /* The radix five case can be translated later..... */ - /* if(ip!=5)goto L112; + let nf = ifac[1]; + let mut l1 = 1; + let mut na = 0; + let mut iw = 1; - ix2=iw+ido; - ix3=ix2+ido; - ix4=ix3+ido; - if(na!=0) - dradb5(ido,l1,ch,c,wa+iw-1,wa+ix2-1,wa+ix3-1,wa+ix4-1); - else - dradb5(ido,l1,c,ch,wa+iw-1,wa+ix2-1,wa+ix3-1,wa+ix4-1); - na=1-na; - goto L115; + for k1 in 0..nf { + let ip = ifac[k1 as usize + 2]; + let l2 = ip * l1; + let ido = n / l2; + let idl1 = ido * l1; - L112:*/ - if na != 0 as c_int { - dradbg( - ido, - ip, - l1, - idl1, - ch, - ch, - ch, - c, - c, - wa.offset(iw as isize) - .offset(-(1 as c_int as isize)), - ); - } else { - dradbg( - ido, - ip, - l1, - idl1, - c, - c, - c, - ch, - ch, - wa.offset(iw as isize) - .offset(-(1 as c_int as isize)), - ); - } - if ido == 1 as c_int { - na = 1 as c_int - na - } - } else { - ix2 = iw + ido; - if na != 0 as c_int { - dradb3( - ido, - l1, - ch, - c, - wa.offset(iw as isize) - .offset(-(1 as c_int as isize)), - wa.offset(ix2 as isize) - .offset(-(1 as c_int as isize)), - ); - } else { - dradb3( - ido, - l1, - c, - ch, - wa.offset(iw as isize) - .offset(-(1 as c_int as isize)), - wa.offset(ix2 as isize) - .offset(-(1 as c_int as isize)), - ); - } - na = 1 as c_int - na - } - } else { - if na != 0 as c_int { - dradb2( - ido, - l1, - ch, - c, - wa.offset(iw as isize).offset(-(1 as c_int as isize)), - ); - } else { - dradb2( - ido, - l1, - c, - ch, - wa.offset(iw as isize).offset(-(1 as c_int as isize)), - ); - } - na = 1 as c_int - na - } + if ip != 4 { + drftb1_l102(ip, ido, l1, idl1, iw, &mut na, c, ch, wa); } else { - ix2 = iw + ido; - ix3 = ix2 + ido; - if na != 0 as c_int { + let ix2 = iw + ido; + let ix3 = ix2 + ido; + if na != 0 { dradb4( ido, l1, - ch, + &ch, c, - wa.offset(iw as isize).offset(-(1 as c_int as isize)), - wa.offset(ix2 as isize).offset(-(1 as c_int as isize)), - wa.offset(ix3 as isize).offset(-(1 as c_int as isize)), + &wa[(iw as usize - 1)..], + &wa[(ix2 as usize - 1)..], + &wa[(ix3 as usize - 1)..], ); } else { dradb4( ido, l1, - c, + &c, ch, - wa.offset(iw as isize).offset(-(1 as c_int as isize)), - wa.offset(ix2 as isize).offset(-(1 as c_int as isize)), - wa.offset(ix3 as isize).offset(-(1 as c_int as isize)), + &wa[(iw as usize - 1)..], + &wa[(ix2 as usize - 1)..], + &wa[(ix3 as usize - 1)..], ); } - na = 1 as c_int - na + na = 1 - na } + l1 = l2; - iw += (ip - 1 as c_int) * ido; - k1 += 1 + iw += (ip - 1) * ido; } - if na == 0 as c_int { + + if na == 0 { return; } - i = 0 as c_int; - while i < n { - *c.offset(i as isize) = *ch.offset(i as isize); - i += 1 - } + + c[..n as usize].copy_from_slice(&ch[..n as usize]); } -#[no_mangle] -pub unsafe extern "C" fn spx_drft_forward( - mut l: *mut drft_lookup, - mut data: *mut c_float, -) { - if (*l).n == 1 as c_int { + +pub fn spx_drft_forward(l: &mut DrftLookup, data: &mut [f32]) { + if l.n == 1 { return; } + + let mut trigcache_temp = l.trigcache[l.n as usize..].to_vec(); + drftf1( - (*l).n, + l.n as i32, data, - (*l).trigcache.as_mut_ptr(), - (*l).trigcache.as_mut_ptr().offset((*l).n as isize), - (*l).splitcache.as_mut_ptr(), + &mut l.trigcache, + &mut trigcache_temp, + &mut l.splitcache, ); + + l.trigcache[l.n as usize..].copy_from_slice(&trigcache_temp); } -#[no_mangle] -pub unsafe extern "C" fn spx_drft_backward( - mut l: *mut drft_lookup, - mut data: *mut c_float, -) { - if (*l).n == 1 as c_int { + +pub fn spx_drft_backward(l: &mut DrftLookup, data: &mut [f32]) { + if l.n == 1 { return; } + + let mut trigcache_temp = l.trigcache[l.n as usize..].to_vec(); + drftb1( - (*l).n, + l.n as i32, data, - (*l).trigcache.as_mut_ptr(), - (*l).trigcache.as_mut_ptr().offset((*l).n as isize), - (*l).splitcache.as_mut_ptr(), + &mut l.trigcache, + &mut trigcache_temp, + &mut l.splitcache, ); + + l.trigcache[l.n as usize..].copy_from_slice(&trigcache_temp); } #[cfg(test)] mod tests { use super::*; - use std::os::raw::{c_float, c_int}; - const EPSILON: c_float = 1e-6; + const EPSILON: f32 = 1e-6; #[test] fn fdrffti_simple() { - let mut trigcache = [42. as c_float; 3]; - let mut splitcache = [24 as c_int; 32]; + let mut trigcache = [42.; 3]; + let mut splitcache = [24; 32]; + + fdrffti(1, &mut trigcache, &mut splitcache); - unsafe { - fdrffti(1, &mut trigcache, &mut splitcache); - } assert!(trigcache.iter().all(|&x| (x - 42.).abs() < EPSILON)); assert!(splitcache.iter().all(|&x| x == 24)); } diff --git a/fft/tests/fftwrap.rs b/fft/tests/fftwrap.rs index 2ca73a3..87ee622 100644 --- a/fft/tests/fftwrap.rs +++ b/fft/tests/fftwrap.rs @@ -5,7 +5,7 @@ mod orig; use crate::orig::fftwrap as original; use speexdsp_fft::*; -const INPUT: [std::os::raw::c_float; 64] = [ +const INPUT: [f32; 64] = [ 1.0, 1.999002, 2.832922, @@ -74,62 +74,40 @@ const INPUT: [std::os::raw::c_float; 64] = [ const EPSILON: f32 = 1e-8; -#[test] -fn fft() { - let mut output = [0.; 64]; - unsafe { - let table = spx_fft_init(INPUT.len() as i32); +macro_rules! test_fftwrap { + ($func: ident) => { + let mut output = [0.; 64]; + let mut table = spx_fft_init(INPUT.len()); let mut input = INPUT.clone(); - spx_fft(table, input.as_mut_ptr(), output.as_mut_ptr()); - spx_fft_destroy(table); - }; + $func(&mut table, &mut input, &mut output); - let mut expected_output = [0.; 64]; - unsafe { - let table = original::spx_fft_init(INPUT.len() as i32); - let mut input = INPUT.clone(); + let mut expected_output = [0.; 64]; + unsafe { + let table = original::spx_fft_init(INPUT.len() as i32); + let mut input = INPUT.clone(); - original::spx_fft( - table, - input.as_mut_ptr(), - expected_output.as_mut_ptr(), - ); - original::spx_fft_destroy(table); + original::$func( + table, + input.as_mut_ptr(), + expected_output.as_mut_ptr(), + ); + original::spx_fft_destroy(table); + }; + + assert!(output + .iter() + .zip(expected_output.iter()) + .all(|(a, b)| (a - b).abs() < EPSILON)); }; +} - assert!(output - .iter() - .zip(expected_output.iter()) - .all(|(a, b)| (a - b).abs() < EPSILON)); +#[test] +fn fft() { + test_fftwrap!(spx_fft); } #[test] fn ifft() { - let mut output = [0.; 64]; - unsafe { - let table = spx_fft_init(INPUT.len() as i32); - let mut input = INPUT.clone(); - - spx_ifft(table, input.as_mut_ptr(), output.as_mut_ptr()); - spx_fft_destroy(table); - }; - - let mut expected_output = [0.; 64]; - unsafe { - let table = original::spx_fft_init(INPUT.len() as i32); - let mut input = INPUT.clone(); - - original::spx_ifft( - table, - input.as_mut_ptr(), - expected_output.as_mut_ptr(), - ); - original::spx_fft_destroy(table); - }; - - assert!(output - .iter() - .zip(expected_output.iter()) - .all(|(a, b)| (a - b).abs() < EPSILON)); + test_fftwrap!(spx_ifft); } diff --git a/fft/tests/smallft.rs b/fft/tests/smallft.rs index 3865514..3bacaf8 100644 --- a/fft/tests/smallft.rs +++ b/fft/tests/smallft.rs @@ -14,11 +14,10 @@ const SPLITCACHE_SIZE: usize = 32; macro_rules! drft { ($func: ident) => { - let mut drft_lookup = drft_lookup::new(N); + let mut drft_lookup = DrftLookup::new(N); let mut data = vec![0.; 32]; - unsafe { - $func(&mut drft_lookup, data.as_mut_ptr()); - } + + $func(&mut drft_lookup, &mut data); let mut original_drft_lookup = original::drft_lookup { n: 0, @@ -62,7 +61,7 @@ macro_rules! drft { #[test] fn fdrffti() { - let drft_lookup = drft_lookup::new(N); + let drft_lookup = DrftLookup::new(N); let mut original_drft_lookup = original::drft_lookup { n: 0, From 1377a10659b9ea02b92f6eba25246905b6b29c89 Mon Sep 17 00:00:00 2001 From: Luni-4 Date: Sun, 19 Apr 2020 23:54:03 +0200 Subject: [PATCH 2/2] Simplify fft's API --- fft/src/fftwrap.rs | 114 ++++++++++++++++++++++++++++------------- fft/src/lib.rs | 1 - fft/src/smallft.rs | 65 ++--------------------- fft/tests/fftwrap.rs | 119 ++++++++++++++++++++++++++++++++++++++++--- fft/tests/smallft.rs | 108 --------------------------------------- 5 files changed, 196 insertions(+), 211 deletions(-) delete mode 100644 fft/tests/smallft.rs diff --git a/fft/src/fftwrap.rs b/fft/src/fftwrap.rs index e8773d9..ccd6672 100644 --- a/fft/src/fftwrap.rs +++ b/fft/src/fftwrap.rs @@ -1,5 +1,3 @@ -use crate::smallft::*; - /* Copyright (C) 2005-2006 Jean-Marc Valin File: fftwrap.c @@ -33,46 +31,94 @@ use crate::smallft::*; SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. */ -pub fn spx_fft_init(size: usize) -> DrftLookup { - DrftLookup::new(size) + +use crate::smallft::*; + +#[derive(Clone)] +pub struct DrftLookup { + pub n: usize, + pub trigcache: Vec, + pub splitcache: Vec, } -pub fn spx_fft(table: &mut DrftLookup, in_0: &mut [f32], out: &mut [f32]) { - let scale = (1.0f64 / table.n as f64) as f32; - if in_0 == out { - eprintln!("FFT should not be done in-place"); +impl DrftLookup { + pub fn new(n: usize) -> Self { + let mut drft = Self { + n: n, + trigcache: vec![0.0; 3 * n], + splitcache: vec![0; 32], + }; + + fdrffti(n, &mut drft.trigcache, &mut drft.splitcache); + + drft } - out.iter_mut() - .zip(in_0.iter()) - .take(table.n as usize) - .for_each(|(o, i)| *o = scale * *i); + pub fn spx_fft(&mut self, in_0: &[f32], out: &mut [f32]) { + let scale = (1.0f64 / self.n as f64) as f32; + if in_0 == out { + eprintln!("FFT should not be done in-place"); + } - spx_drft_forward(table, out); -} + out.iter_mut() + .zip(in_0.iter()) + .take(self.n as usize) + .for_each(|(o, i)| *o = scale * *i); -pub fn spx_ifft(table: &mut DrftLookup, in_0: &mut [f32], out: &mut [f32]) { - if in_0 == out { - eprintln!("FFT should not be done in-place"); - } else { - out.copy_from_slice(&in_0[..table.n as usize]); + self.spx_drft_forward(out); } - spx_drft_backward(table, out); -} + pub fn spx_ifft(&mut self, in_0: &[f32], out: &mut [f32]) { + if in_0 == out { + eprintln!("FFT should not be done in-place"); + } else { + out.copy_from_slice(&in_0[..self.n as usize]); + } -pub fn spx_fft_float( - table: &mut DrftLookup, - in_0: &mut [f32], - out: &mut [f32], -) { - spx_fft(table, in_0, out); -} + self.spx_drft_backward(out); + } + + pub fn spx_fft_float(&mut self, in_0: &[f32], out: &mut [f32]) { + self.spx_fft(in_0, out); + } + + pub fn spx_ifft_float(&mut self, in_0: &[f32], out: &mut [f32]) { + self.spx_ifft(in_0, out); + } + + pub fn spx_drft_forward(&mut self, data: &mut [f32]) { + if self.n == 1 { + return; + } + + let mut trigcache_temp = self.trigcache[self.n as usize..].to_vec(); + + drftf1( + self.n as i32, + data, + &mut self.trigcache, + &mut trigcache_temp, + &mut self.splitcache, + ); -pub fn spx_ifft_float( - table: &mut DrftLookup, - in_0: &mut [f32], - out: &mut [f32], -) { - spx_ifft(table, in_0, out); + self.trigcache[self.n as usize..].copy_from_slice(&trigcache_temp); + } + + pub fn spx_drft_backward(&mut self, data: &mut [f32]) { + if self.n == 1 { + return; + } + + let mut trigcache_temp = self.trigcache[self.n as usize..].to_vec(); + + drftb1( + self.n as i32, + data, + &mut self.trigcache, + &mut trigcache_temp, + &mut self.splitcache, + ); + + self.trigcache[self.n as usize..].copy_from_slice(&trigcache_temp); + } } diff --git a/fft/src/lib.rs b/fft/src/lib.rs index 74e6239..ce2b455 100644 --- a/fft/src/lib.rs +++ b/fft/src/lib.rs @@ -4,4 +4,3 @@ mod fftwrap; mod smallft; pub use crate::fftwrap::*; -pub use crate::smallft::*; diff --git a/fft/src/smallft.rs b/fft/src/smallft.rs index 09e6744..e4c850c 100644 --- a/fft/src/smallft.rs +++ b/fft/src/smallft.rs @@ -1,27 +1,6 @@ use crate::dradb::*; use crate::dradf::*; -#[derive(Clone)] -pub struct DrftLookup { - pub n: usize, - pub trigcache: Vec, - pub splitcache: Vec, -} - -impl DrftLookup { - pub fn new(n: usize) -> Self { - let mut drft = Self { - n: n, - trigcache: vec![0.0; 3 * n], - splitcache: vec![0; 32], - }; - - fdrffti(n, &mut drft.trigcache, &mut drft.splitcache); - - drft - } -} - #[inline(always)] fn drfti1_c_10244(ifac: &mut [i32], n: i32, nf: &mut i32) { const NTRYH: [i32; 4] = [4, 2, 3, 5]; @@ -62,7 +41,7 @@ fn drfti1_c_10244(ifac: &mut [i32], n: i32, nf: &mut i32) { } } -fn drfti1(wa: &mut [f32], ifac: &mut [i32]) { +pub(crate) fn drfti1(wa: &mut [f32], ifac: &mut [i32]) { const TPI: f32 = 6.283_185_307_179_586_48; let n = wa.len() as i32; @@ -112,7 +91,7 @@ fn drfti1(wa: &mut [f32], ifac: &mut [i32]) { } } -fn fdrffti(n: usize, wsave: &mut [f32], ifac: &mut [i32]) { +pub(crate) fn fdrffti(n: usize, wsave: &mut [f32], ifac: &mut [i32]) { if n == 1 { return; } @@ -152,7 +131,7 @@ fn drftf1_l102( } } -fn drftf1( +pub(crate) fn drftf1( n: i32, c: &mut [f32], ch: &mut [f32], @@ -254,7 +233,7 @@ fn drftb1_l102( } } -fn drftb1( +pub(crate) fn drftb1( n: i32, c: &mut [f32], ch: &mut [f32], @@ -312,42 +291,6 @@ fn drftb1( c[..n as usize].copy_from_slice(&ch[..n as usize]); } -pub fn spx_drft_forward(l: &mut DrftLookup, data: &mut [f32]) { - if l.n == 1 { - return; - } - - let mut trigcache_temp = l.trigcache[l.n as usize..].to_vec(); - - drftf1( - l.n as i32, - data, - &mut l.trigcache, - &mut trigcache_temp, - &mut l.splitcache, - ); - - l.trigcache[l.n as usize..].copy_from_slice(&trigcache_temp); -} - -pub fn spx_drft_backward(l: &mut DrftLookup, data: &mut [f32]) { - if l.n == 1 { - return; - } - - let mut trigcache_temp = l.trigcache[l.n as usize..].to_vec(); - - drftb1( - l.n as i32, - data, - &mut l.trigcache, - &mut trigcache_temp, - &mut l.splitcache, - ); - - l.trigcache[l.n as usize..].copy_from_slice(&trigcache_temp); -} - #[cfg(test)] mod tests { use super::*; diff --git a/fft/tests/fftwrap.rs b/fft/tests/fftwrap.rs index 87ee622..df7cd31 100644 --- a/fft/tests/fftwrap.rs +++ b/fft/tests/fftwrap.rs @@ -2,9 +2,13 @@ extern crate speexdsp_fft; mod orig; -use crate::orig::fftwrap as original; +use std::ptr::null_mut; + use speexdsp_fft::*; +use crate::orig::fftwrap as original_fftwrap; +use crate::orig::smallft as original_smallft; + const INPUT: [f32; 64] = [ 1.0, 1.999002, @@ -73,26 +77,29 @@ const INPUT: [f32; 64] = [ ]; const EPSILON: f32 = 1e-8; +const N: usize = 2; +const SPLITCACHE_SIZE: usize = 32; macro_rules! test_fftwrap { ($func: ident) => { let mut output = [0.; 64]; - let mut table = spx_fft_init(INPUT.len()); + let mut drft_lookup = DrftLookup::new(INPUT.len()); let mut input = INPUT.clone(); - $func(&mut table, &mut input, &mut output); + drft_lookup.$func(&mut input, &mut output); let mut expected_output = [0.; 64]; unsafe { - let table = original::spx_fft_init(INPUT.len() as i32); + let drft_lookup = + original_fftwrap::spx_fft_init(INPUT.len() as i32); let mut input = INPUT.clone(); - original::$func( - table, + original_fftwrap::$func( + drft_lookup, input.as_mut_ptr(), expected_output.as_mut_ptr(), ); - original::spx_fft_destroy(table); + original_fftwrap::spx_fft_destroy(drft_lookup); }; assert!(output @@ -102,6 +109,56 @@ macro_rules! test_fftwrap { }; } +macro_rules! drft { + ($func: ident) => { + let mut drft_lookup = DrftLookup::new(N); + let mut data = vec![0.; 32]; + + drft_lookup.$func(&mut data); + + let mut original_drft_lookup = original_smallft::drft_lookup { + n: 0, + trigcache: null_mut(), + splitcache: null_mut(), + }; + let mut data_original = vec![0.; 32]; + unsafe { + original_smallft::spx_drft_init( + &mut original_drft_lookup, + N as i32, + ); + original_smallft::$func( + &mut original_drft_lookup, + data_original.as_mut_ptr(), + ); + } + + let expected_trigcache = unsafe { + Vec::from_raw_parts( + original_drft_lookup.trigcache as *mut f32, + 3 * N, + 3 * N, + ) + }; + + let expected_splitcache = unsafe { + Vec::from_raw_parts( + original_drft_lookup.splitcache as *mut i32, + SPLITCACHE_SIZE, + SPLITCACHE_SIZE, + ) + }; + + assert!(drft_lookup + .trigcache + .iter() + .zip(expected_trigcache.iter()) + .all(|(&a, &b)| (a - b).abs() < EPSILON)); + + assert_eq!(&drft_lookup.splitcache, &expected_splitcache); + }; +} + #[test] fn fft() { test_fftwrap!(spx_fft); @@ -111,3 +168,51 @@ fn fft() { fn ifft() { test_fftwrap!(spx_ifft); } + +#[test] +fn fdrffti() { + let drft_lookup = DrftLookup::new(N); + + let mut original_drft_lookup = original_smallft::drft_lookup { + n: 0, + trigcache: null_mut(), + splitcache: null_mut(), + }; + unsafe { + original_smallft::spx_drft_init(&mut original_drft_lookup, N as i32); + } + + let expected_trigcache = unsafe { + Vec::from_raw_parts( + original_drft_lookup.trigcache as *mut f32, + 3 * N, + 3 * N, + ) + }; + + let expected_splitcache = unsafe { + Vec::from_raw_parts( + original_drft_lookup.splitcache as *mut i32, + SPLITCACHE_SIZE, + SPLITCACHE_SIZE, + ) + }; + + assert!(drft_lookup + .trigcache + .iter() + .zip(expected_trigcache.iter()) + .all(|(&a, &b)| (a - b).abs() < EPSILON)); + + assert_eq!(&drft_lookup.splitcache, &expected_splitcache); +} + +#[test] +fn drftf1() { + drft!(spx_drft_forward); +} + +#[test] +fn drftb1() { + drft!(spx_drft_backward); +} diff --git a/fft/tests/smallft.rs b/fft/tests/smallft.rs deleted file mode 100644 index 3bacaf8..0000000 --- a/fft/tests/smallft.rs +++ /dev/null @@ -1,108 +0,0 @@ -extern crate speexdsp_fft; - -mod orig; - -use crate::orig::smallft as original; -use speexdsp_fft::*; - -use std::os::raw::{c_float, c_int}; -use std::ptr::null_mut; - -const EPSILON: c_float = 1e-6; -const N: usize = 2; -const SPLITCACHE_SIZE: usize = 32; - -macro_rules! drft { - ($func: ident) => { - let mut drft_lookup = DrftLookup::new(N); - let mut data = vec![0.; 32]; - - $func(&mut drft_lookup, &mut data); - - let mut original_drft_lookup = original::drft_lookup { - n: 0, - trigcache: null_mut(), - splitcache: null_mut(), - }; - let mut data_original = vec![0.; 32]; - unsafe { - original::spx_drft_init(&mut original_drft_lookup, N as c_int); - original::$func( - &mut original_drft_lookup, - data_original.as_mut_ptr(), - ); - } - - let expected_trigcache = unsafe { - Vec::from_raw_parts( - original_drft_lookup.trigcache as *mut c_float, - 3 * N, - 3 * N, - ) - }; - - let expected_splitcache = unsafe { - Vec::from_raw_parts( - original_drft_lookup.splitcache as *mut c_int, - SPLITCACHE_SIZE, - SPLITCACHE_SIZE, - ) - }; - - assert!(drft_lookup - .trigcache - .iter() - .zip(expected_trigcache.iter()) - .all(|(&a, &b)| (a - b).abs() < EPSILON)); - - assert_eq!(&drft_lookup.splitcache, &expected_splitcache); - }; -} - -#[test] -fn fdrffti() { - let drft_lookup = DrftLookup::new(N); - - let mut original_drft_lookup = original::drft_lookup { - n: 0, - trigcache: null_mut(), - splitcache: null_mut(), - }; - unsafe { - original::spx_drft_init(&mut original_drft_lookup, N as c_int); - } - - let expected_trigcache = unsafe { - Vec::from_raw_parts( - original_drft_lookup.trigcache as *mut c_float, - 3 * N, - 3 * N, - ) - }; - - let expected_splitcache = unsafe { - Vec::from_raw_parts( - original_drft_lookup.splitcache as *mut c_int, - SPLITCACHE_SIZE, - SPLITCACHE_SIZE, - ) - }; - - assert!(drft_lookup - .trigcache - .iter() - .zip(expected_trigcache.iter()) - .all(|(&a, &b)| (a - b).abs() < EPSILON)); - - assert_eq!(&drft_lookup.splitcache, &expected_splitcache); -} - -#[test] -fn drftf1() { - drft!(spx_drft_forward); -} - -#[test] -fn drftb1() { - drft!(spx_drft_backward); -}