From f36132e84be6d10e1c0aeb5071a431a9e5f5c863 Mon Sep 17 00:00:00 2001 From: daniellga Date: Sun, 14 Jul 2024 15:18:45 -0300 Subject: [PATCH] fixed fft and starting build stft --- Cargo.toml | 2 + harmonium-core/src/array.rs | 9 +- harmonium-core/src/audioop.rs | 15 +- harmonium-fft/Cargo.toml | 1 - harmonium-fft/src/fft.rs | 646 +++++++++++------- harmonium-stft/Cargo.toml | 9 + harmonium-stft/src/lib.rs | 1 + harmonium-stft/src/stft.rs | 101 +++ harmonium-window/src/windows.rs | 8 +- r-harmonium/.fft_bench.R | 10 +- r-harmonium/DESCRIPTION | 2 +- r-harmonium/NAMESPACE | 6 +- r-harmonium/R/000-wrappers.R | 167 +++-- r-harmonium/R/structs.R | 5 +- r-harmonium/R/zzz.R | 4 +- r-harmonium/src/init.c | 90 ++- r-harmonium/src/rust/Cargo.toml | 21 +- r-harmonium/src/rust/api.h | 29 +- r-harmonium/src/rust/src/conversions.rs | 6 +- r-harmonium/src/rust/src/harray.rs | 2 +- r-harmonium/src/rust/src/haudiosink.rs | 2 +- r-harmonium/src/rust/src/hfft.rs | 632 +++++++++++------ .../{test_hfftplanner.R => test_hfft.R} | 24 +- 23 files changed, 1153 insertions(+), 639 deletions(-) create mode 100644 harmonium-stft/Cargo.toml create mode 100644 harmonium-stft/src/lib.rs create mode 100644 harmonium-stft/src/stft.rs rename r-harmonium/tests/testthat/{test_hfftplanner.R => test_hfft.R} (82%) diff --git a/Cargo.toml b/Cargo.toml index 3e7eb57..7f345be 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -6,6 +6,7 @@ members = [ "harmonium-fft", "harmonium-resample", "harmonium-window", + "harmonium-stft", ] [workspace.dependencies] @@ -23,6 +24,7 @@ harmonium-io = { path = "harmonium-io", default-features = false } harmonium-fft = { path = "harmonium-fft", default-features = false } harmonium-resample = { path = "harmonium-resample", default-features = false } harmonium-window = { path = "harmonium-window", default-features = false } +harmonium-stft = { path = "harmonium-stft", default-features = false } [profile.release] opt-level = 3 diff --git a/harmonium-core/src/array.rs b/harmonium-core/src/array.rs index c1a7749..007f900 100644 --- a/harmonium-core/src/array.rs +++ b/harmonium-core/src/array.rs @@ -14,8 +14,13 @@ where where Sh: Into>, { - let ndarray = - ArcArray::from_shape_vec(shape, v).map_err(|_| HError::OutOfSpecError("shape does not correspond to the number of elements in v or if the shape/strides would result in overflowing isize".to_string()))?; + let ndarray = ArcArray::from_shape_vec(shape, v).map_err(|_| { + HError::OutOfSpecError( + "shape does not correspond to the number of elements + in v or if the shape/strides would result in overflowing isize" + .to_string(), + ) + })?; Ok(HArray(ndarray)) } diff --git a/harmonium-core/src/audioop.rs b/harmonium-core/src/audioop.rs index 2d3eb91..6ba492f 100644 --- a/harmonium-core/src/audioop.rs +++ b/harmonium-core/src/audioop.rs @@ -2,19 +2,18 @@ use crate::{ array::HArray, errors::{HError, HResult}, }; -use ndarray::{Axis, Dimension, Ix1, Ix2, IxDyn}; +use ndarray::{Axis, Dimension, Ix0, Ix1, Ix2, IxDyn}; use num_traits::{Float, FloatConst, FromPrimitive}; -pub trait AudioOp +pub trait AudioOp where T: Float + FloatConst + FromPrimitive, D: Dimension, - U: Dimension, { fn nchannels(&self) -> usize; fn nframes(&self) -> usize; fn db_to_amplitude(&mut self, reference: T, power: T); - fn to_mono(&self) -> HResult>; + fn to_mono(&self) -> HResult>; } pub enum Audio<'a, T> @@ -26,7 +25,7 @@ where Dyn(&'a HArray), } -impl AudioOp for HArray +impl AudioOp for HArray where T: Float + FloatConst + FromPrimitive, { @@ -50,7 +49,7 @@ where .mapv_inplace(|x| reference * a.powf(b * x).powf(power)); } - fn to_mono(&self) -> HResult> { + fn to_mono(&self) -> HResult> { // To return an error is a design choice. This wasn't supposed to error. Err(HError::OutOfSpecError( "The length of the axis is zero.".into(), @@ -58,7 +57,7 @@ where } } -impl AudioOp for HArray +impl AudioOp for HArray where T: Float + FloatConst + FromPrimitive, { @@ -89,7 +88,7 @@ where } } -impl AudioOp for HArray +impl AudioOp for HArray where T: Float + FloatConst + FromPrimitive, { diff --git a/harmonium-fft/Cargo.toml b/harmonium-fft/Cargo.toml index 2e61e12..38c1b9f 100644 --- a/harmonium-fft/Cargo.toml +++ b/harmonium-fft/Cargo.toml @@ -7,7 +7,6 @@ edition = "2021" [dependencies] harmonium-core = { workspace = true } -harmonium-io = { workspace = true } rustfft = { workspace = true } realfft = { workspace = true } ndarray = { workspace = true } diff --git a/harmonium-fft/src/fft.rs b/harmonium-fft/src/fft.rs index 7c74489..d6302a8 100644 --- a/harmonium-fft/src/fft.rs +++ b/harmonium-fft/src/fft.rs @@ -1,262 +1,367 @@ +use std::sync::Arc; + use harmonium_core::{array::HArray, errors::HError, errors::HResult}; use ndarray::{ArcArray1, ArcArray2, Axis, Dimension, Ix1, Ix2, IxDyn, Zip}; -use realfft::RealFftPlanner; +use realfft::{ComplexToReal, RealFftPlanner, RealToComplex}; use rustfft::{ num_complex::Complex, num_traits::{Float, FloatConst}, - FftPlanner, + FftNum, FftPlanner, }; +#[derive(Clone)] +pub struct Fft { + pub fft: Arc>, + pub scratch_buffer: Arc<[Complex]>, +} + +#[derive(Clone)] +pub struct RealFftForward { + fft: Arc>, + scratch_buffer: Arc<[Complex]>, +} + +#[derive(Clone)] +pub struct RealFftInverse { + fft: Arc>, + scratch_buffer: Arc<[Complex]>, +} + +impl Fft { + pub fn new_fft_forward(length: usize) -> Self { + let mut planner = FftPlanner::new(); + let fft = planner.plan_fft_forward(length); + let scratch_len = fft.get_inplace_scratch_len(); + let zero = T::zero(); + let scratch_buffer = vec![Complex::new(zero, zero); scratch_len]; + let scratch_buffer: Arc<[Complex]> = Arc::from(scratch_buffer); + + Self { + fft, + scratch_buffer, + } + } + + pub fn new_fft_inverse(length: usize) -> Self { + let mut planner = FftPlanner::new(); + let fft = planner.plan_fft_inverse(length); + let scratch_len = fft.get_inplace_scratch_len(); + let zero = T::zero(); + let scratch_buffer = vec![Complex::new(zero, zero); scratch_len]; + let scratch_buffer: Arc<[Complex]> = Arc::from(scratch_buffer); + + Self { + fft, + scratch_buffer, + } + } +} + +impl RealFftForward { + pub fn new_real_fft_forward(length: usize) -> Self { + let mut planner = RealFftPlanner::new(); + let fft = planner.plan_fft_forward(length); + let zero = T::zero(); + let scratch_len = fft.get_scratch_len(); + let scratch_buffer = vec![Complex::new(zero, zero); scratch_len]; + let scratch_buffer: Arc<[Complex]> = Arc::from(scratch_buffer); + + Self { + fft, + scratch_buffer, + } + } +} + +impl RealFftInverse { + pub fn new_real_fft_inverse(length: usize) -> Self { + let mut planner = RealFftPlanner::new(); + let fft = planner.plan_fft_inverse(length); + let zero = T::zero(); + let scratch_len = fft.get_scratch_len(); + let scratch_buffer = vec![Complex::new(zero, zero); scratch_len]; + let scratch_buffer: Arc<[Complex]> = Arc::from(scratch_buffer); + + Self { + fft, + scratch_buffer, + } + } +} + pub trait ProcessFft where - T: Float + FloatConst, + T: FftNum + Float + FloatConst, D: Dimension, { - fn fft(&mut self, harray: &mut HArray, D>) -> HResult<()>; - fn ifft(&mut self, harray: &mut HArray, D>) -> HResult<()>; + fn process(&mut self, harray: &mut HArray, D>) -> HResult<()>; } -pub trait ProcessRealFft +pub trait ProcessRealFftForward where - T: Float + FloatConst, + T: FftNum + Float + FloatConst, D: Dimension, { - fn rfft(&mut self, harray: &mut HArray) -> HResult, D>>; - fn irfft(&mut self, harray: &mut HArray, D>, length: usize) - -> HResult>; + fn process(&mut self, harray: &mut HArray) -> HResult, D>>; } -macro_rules! impl_fft { - ($($t: ty),+) => { - $( - impl ProcessFft<$t, Ix1> for FftPlanner<$t> { - fn fft(&mut self, harray: &mut HArray, Ix1>) -> HResult<()> { - let length = harray.len(); - let fft = self.plan_fft_forward(length); - let scratch_len = fft.get_inplace_scratch_len(); - let mut scratch_buffer = vec![Complex::new(0., 0.); scratch_len]; - fft.process_with_scratch(harray.as_slice_mut().unwrap(), &mut scratch_buffer); - Ok(()) - } - - fn ifft(&mut self, harray: &mut HArray, Ix1>) -> HResult<()> { - let length = harray.len(); - let fft = self.plan_fft_inverse(length); - let scratch_len = fft.get_inplace_scratch_len(); - let mut scratch_buffer = vec![Complex::new(0., 0.); scratch_len]; - fft.process_with_scratch(harray.as_slice_mut().unwrap(), &mut scratch_buffer); - Ok(()) - } - } +pub trait ProcessRealFftInverse +where + T: FftNum + Float + FloatConst, + D: Dimension, +{ + fn process(&mut self, harray: &mut HArray, D>) -> HResult>; +} - impl ProcessRealFft<$t, Ix1> for RealFftPlanner<$t> { - fn rfft(&mut self, harray: &mut HArray<$t, Ix1>) -> HResult, Ix1>> { - let length = harray.len(); - let r2c = self.plan_fft_forward(length); - let scratch_len = r2c.get_scratch_len(); - let mut scratch_buffer = vec![Complex::new(0., 0.); scratch_len]; - let mut ndarray = ArcArray1::from_elem(length / 2 + 1, Complex::new(0., 0.)); - r2c.process_with_scratch(harray.as_slice_mut().unwrap(), ndarray.as_slice_mut().unwrap(), &mut scratch_buffer).unwrap(); - Ok(HArray(ndarray)) - } - - fn irfft(&mut self, harray: &mut HArray, Ix1>, length: usize) -> HResult> { - let c2r = self.plan_fft_inverse(length); - let scratch_len = c2r.get_scratch_len(); - let mut scratch_buffer = vec![Complex::new(0., 0.); scratch_len]; - let mut ndarray = ArcArray1::from_elem(length, 0.); - c2r.process_with_scratch(harray.as_slice_mut().unwrap(), ndarray.as_slice_mut().unwrap(), &mut scratch_buffer).unwrap(); - Ok(HArray(ndarray)) - } - } +impl ProcessFft for Fft +where + T: FftNum + Float + FloatConst, +{ + fn process(&mut self, harray: &mut HArray, Ix1>) -> HResult<()> { + let scratch_buffer = make_mut_slice(&mut self.scratch_buffer); + self.fft + .process_with_scratch(harray.as_slice_mut().unwrap(), scratch_buffer); + Ok(()) + } +} - impl ProcessFft<$t, Ix2> for FftPlanner<$t> { - fn fft(&mut self, harray: &mut HArray, Ix2>) -> HResult<()> { - let ncols = harray.0.ncols(); - let fft = self.plan_fft_forward(ncols); - let scratch_len = fft.get_inplace_scratch_len(); - let mut scratch_buffer = vec![Complex::new(0., 0.); scratch_len]; +impl ProcessRealFftForward for RealFftForward +where + T: FftNum + Float + FloatConst, +{ + fn process(&mut self, harray: &mut HArray) -> HResult, Ix1>> { + let zero = T::zero(); + let length = harray.len(); + let mut ndarray = ArcArray1::from_elem(length / 2 + 1, Complex::new(zero, zero)); + let scratch_buffer = make_mut_slice(&mut self.scratch_buffer); + self.fft + .process_with_scratch( + harray.as_slice_mut().unwrap(), + ndarray.as_slice_mut().unwrap(), + scratch_buffer, + ) + .unwrap(); + Ok(HArray(ndarray)) + } +} - Zip::from(harray.0.lanes_mut(Axis(1))).for_each(|mut row| { - fft.process_with_scratch(row.as_slice_mut().unwrap(), &mut scratch_buffer); - }); - Ok(()) - } +impl ProcessRealFftInverse for RealFftInverse +where + T: FftNum + Float + FloatConst, +{ + fn process(&mut self, harray: &mut HArray, Ix1>) -> HResult> { + let zero = T::zero(); + let length = self.fft.len(); + let scratch_buffer = make_mut_slice(&mut self.scratch_buffer); + let mut ndarray = ArcArray1::from_elem(length, zero); + self.fft + .process_with_scratch( + harray.as_slice_mut().unwrap(), + ndarray.as_slice_mut().unwrap(), + scratch_buffer, + ) + .unwrap(); + Ok(HArray(ndarray)) + } +} + +impl ProcessFft for Fft +where + T: FftNum + Float + FloatConst, +{ + fn process(&mut self, harray: &mut HArray, Ix2>) -> HResult<()> { + let scratch_buffer = make_mut_slice(&mut self.scratch_buffer); + + Zip::from(harray.0.lanes_mut(Axis(1))).for_each(|mut row| { + self.fft + .process_with_scratch(row.as_slice_mut().unwrap(), scratch_buffer); + }); + Ok(()) + } +} - fn ifft(&mut self, harray: &mut HArray, Ix2>) -> HResult<()> { - let ncols = harray.0.ncols(); - let fft = self.plan_fft_inverse(ncols); - let scratch_len = fft.get_inplace_scratch_len(); - let mut scratch_buffer = vec![Complex::new(0., 0.); scratch_len]; +impl ProcessRealFftForward for RealFftForward +where + T: FftNum + Float + FloatConst, +{ + fn process(&mut self, harray: &mut HArray) -> HResult, Ix2>> { + let zero = T::zero(); + let nrows = harray.0.nrows(); + let ncols = harray.0.ncols(); + let scratch_buffer = make_mut_slice(&mut self.scratch_buffer); + let mut ndarray = ArcArray2::from_elem((nrows, ncols / 2 + 1), Complex::new(zero, zero)); + + Zip::from(ndarray.lanes_mut(Axis(1))) + .and(harray.0.lanes_mut(Axis(1))) + .for_each(|mut row_output, mut row_input| { + self.fft + .process_with_scratch( + row_input.as_slice_mut().unwrap(), + row_output.as_slice_mut().unwrap(), + scratch_buffer, + ) + .unwrap(); + }); + + Ok(HArray(ndarray)) + } +} - Zip::from(harray.0.lanes_mut(Axis(1))).for_each(|mut row| { - fft.process_with_scratch(row.as_slice_mut().unwrap(), &mut scratch_buffer); - }); - Ok(()) - } +impl ProcessRealFftInverse for RealFftInverse +where + T: FftNum + Float + FloatConst, +{ + fn process(&mut self, harray: &mut HArray, Ix2>) -> HResult> { + let zero = T::zero(); + let length = self.fft.len(); + let nrows = harray.0.nrows(); + let scratch_buffer = make_mut_slice(&mut self.scratch_buffer); + let mut ndarray = ArcArray2::from_elem((nrows, length), zero); + + Zip::from(ndarray.lanes_mut(Axis(1))) + .and(harray.0.lanes_mut(Axis(1))) + .for_each(|mut row_output, mut row_input| { + self.fft + .process_with_scratch( + row_input.as_slice_mut().unwrap(), + row_output.as_slice_mut().unwrap(), + scratch_buffer, + ) + .unwrap(); + }); + + Ok(HArray(ndarray)) + } +} + +impl ProcessFft for Fft +where + T: FftNum + Float + FloatConst, +{ + fn process(&mut self, harray: &mut HArray, IxDyn>) -> HResult<()> { + let scratch_buffer = make_mut_slice(&mut self.scratch_buffer); + match harray.ndim() { + 1 => { + self.fft + .process_with_scratch(harray.as_slice_mut().unwrap(), scratch_buffer); + Ok(()) + } + 2 => { + Zip::from(harray.0.lanes_mut(Axis(1))).for_each(|mut row| { + self.fft + .process_with_scratch(row.as_slice_mut().unwrap(), scratch_buffer); + }); + Ok(()) } + _ => Err(HError::OutOfSpecError( + "The HArray's ndim should be 1 or 2.".into(), + )), + } + } +} - impl ProcessRealFft<$t, Ix2> for RealFftPlanner<$t> { - fn rfft(&mut self, harray: &mut HArray<$t, Ix2>) -> HResult, Ix2>> { - let nrows = harray.0.nrows(); - let ncols = harray.0.ncols(); - let r2c = self.plan_fft_forward(ncols); - let scratch_len = r2c.get_scratch_len(); - let mut scratch_buffer = vec![Complex::new(0., 0.); scratch_len]; - let mut ndarray = ArcArray2::from_elem((nrows, ncols / 2 + 1), Complex::new(0., 0.)); - - Zip::from(ndarray.lanes_mut(Axis(1))) - .and(harray.0.lanes_mut(Axis(1))) - .for_each(|mut row_output, mut row_input| { - r2c - .process_with_scratch(row_input.as_slice_mut().unwrap(), row_output.as_slice_mut().unwrap(), &mut scratch_buffer) - .unwrap(); - }); - - Ok(HArray(ndarray)) - } - - fn irfft(&mut self, harray: &mut HArray, Ix2>, length: usize) -> HResult> { - let nrows = harray.0.nrows(); - let c2r = self.plan_fft_inverse(length); - let scratch_len = c2r.get_scratch_len(); - let mut scratch_buffer = vec![Complex::new(0., 0.); scratch_len]; - let mut ndarray = ArcArray2::from_elem((nrows, length), 0.); - - Zip::from(ndarray.lanes_mut(Axis(1))) - .and(harray.0.lanes_mut(Axis(1))) - .for_each(|mut row_output, mut row_input| { - c2r - .process_with_scratch(row_input.as_slice_mut().unwrap(), row_output.as_slice_mut().unwrap(), &mut scratch_buffer) - .unwrap(); - }); - - Ok(HArray(ndarray)) - } +impl ProcessRealFftForward for RealFftForward +where + T: FftNum + Float + FloatConst, +{ + fn process(&mut self, harray: &mut HArray) -> HResult, IxDyn>> { + let zero = T::zero(); + let scratch_buffer = make_mut_slice(&mut self.scratch_buffer); + match harray.ndim() { + 1 => { + let length = harray.len(); + let mut ndarray = ArcArray1::from_elem(length / 2 + 1, Complex::new(zero, zero)); + self.fft + .process_with_scratch( + harray.as_slice_mut().unwrap(), + ndarray.as_slice_mut().unwrap(), + scratch_buffer, + ) + .unwrap(); + Ok(HArray(ndarray.into_dyn())) } + 2 => { + let nrows = harray.0.len_of(Axis(0)); + let ncols = harray.0.len_of(Axis(1)); + let mut ndarray = + ArcArray2::from_elem((nrows, ncols / 2 + 1), Complex::new(zero, zero)) + .into_dyn(); + + Zip::from(ndarray.lanes_mut(Axis(1))) + .and(harray.0.lanes_mut(Axis(1))) + .for_each(|mut row_output, mut row_input| { + self.fft + .process_with_scratch( + row_input.as_slice_mut().unwrap(), + row_output.as_slice_mut().unwrap(), + scratch_buffer, + ) + .unwrap(); + }); + + Ok(HArray(ndarray)) + } + _ => Err(HError::OutOfSpecError( + "The HArray's ndim should be 1 or 2.".into(), + )), + } + } +} - impl ProcessFft<$t, IxDyn> for FftPlanner<$t> { - fn fft(&mut self, harray: &mut HArray, IxDyn>) -> HResult<()> { - match harray.ndim() { - 1 => { - let length = harray.len(); - let fft = self.plan_fft_forward(length); - let scratch_len = fft.get_inplace_scratch_len(); - let mut scratch_buffer = vec![Complex::new(0., 0.); scratch_len]; - fft.process_with_scratch(harray.as_slice_mut().unwrap(), &mut scratch_buffer); - Ok(()) - }, - 2 => { - let ncols = harray.0.len_of(Axis(1)); - let fft = self.plan_fft_forward(ncols); - let scratch_len = fft.get_inplace_scratch_len(); - let mut scratch_buffer = vec![Complex::new(0., 0.); scratch_len]; - - Zip::from(harray.0.lanes_mut(Axis(1))).for_each(|mut row| { - fft.process_with_scratch(row.as_slice_mut().unwrap(), &mut scratch_buffer); - }); - Ok(()) - }, - _ => Err(HError::OutOfSpecError("The HArray's ndim should be 1 or 2.".into())), - } - } - - fn ifft(&mut self, harray: &mut HArray, IxDyn>) -> HResult<()> { - match harray.ndim() { - 1 => { - let length = harray.len(); - let fft = self.plan_fft_inverse(length); - let scratch_len = fft.get_inplace_scratch_len(); - let mut scratch_buffer = vec![Complex::new(0., 0.); scratch_len]; - fft.process_with_scratch(harray.as_slice_mut().unwrap(), &mut scratch_buffer); - Ok(()) - }, - 2 => { - let ncols = harray.0.len_of(Axis(1)); - let fft = self.plan_fft_inverse(ncols); - let scratch_len = fft.get_inplace_scratch_len(); - let mut scratch_buffer = vec![Complex::new(0., 0.); scratch_len]; - - Zip::from(harray.0.lanes_mut(Axis(1))).for_each(|mut row| { - fft.process_with_scratch(row.as_slice_mut().unwrap(), &mut scratch_buffer); - }); - Ok(()) - }, - _ => Err(HError::OutOfSpecError("The HArray's ndim should be 1 or 2.".into())), - } - } +impl ProcessRealFftInverse for RealFftInverse +where + T: FftNum + Float + FloatConst, +{ + fn process(&mut self, harray: &mut HArray, IxDyn>) -> HResult> { + let zero = T::zero(); + let length = self.fft.len(); + let scratch_buffer = make_mut_slice(&mut self.scratch_buffer); + match harray.ndim() { + 1 => { + let mut ndarray = ArcArray1::from_elem(length, zero); + self.fft + .process_with_scratch( + harray.as_slice_mut().unwrap(), + ndarray.as_slice_mut().unwrap(), + scratch_buffer, + ) + .unwrap(); + Ok(HArray(ndarray.into_dyn())) } + 2 => { + let nrows = harray.0.len_of(Axis(0)); + let mut ndarray = ArcArray2::from_elem((nrows, length), zero).into_dyn(); + + Zip::from(ndarray.lanes_mut(Axis(1))) + .and(harray.0.lanes_mut(Axis(1))) + .for_each(|mut row_output, mut row_input| { + self.fft + .process_with_scratch( + row_input.as_slice_mut().unwrap(), + row_output.as_slice_mut().unwrap(), + scratch_buffer, + ) + .unwrap(); + }); - impl ProcessRealFft<$t, IxDyn> for RealFftPlanner<$t> { - fn rfft(&mut self, harray: &mut HArray<$t, IxDyn>) -> HResult, IxDyn>> { - match harray.ndim() { - 1 => { - let length = harray.len(); - let r2c = self.plan_fft_forward(length); - let scratch_len = r2c.get_scratch_len(); - let mut scratch_buffer = vec![Complex::new(0., 0.); scratch_len]; - let mut ndarray = ArcArray1::from_elem(length / 2 + 1, Complex::new(0., 0.)); - r2c.process_with_scratch(harray.as_slice_mut().unwrap(), ndarray.as_slice_mut().unwrap(), &mut scratch_buffer).unwrap(); - Ok(HArray(ndarray.into_dyn())) - }, - 2 => { - let nrows = harray.0.len_of(Axis(0)); - let ncols = harray.0.len_of(Axis(1)); - let r2c = self.plan_fft_forward(ncols); - let scratch_len = r2c.get_scratch_len(); - let mut scratch_buffer = vec![Complex::new(0., 0.); scratch_len]; - let mut ndarray = ArcArray2::from_elem((nrows, ncols / 2 + 1), Complex::new(0., 0.)).into_dyn(); - - Zip::from(ndarray.lanes_mut(Axis(1))) - .and(harray.0.lanes_mut(Axis(1))) - .for_each(|mut row_output, mut row_input| { - r2c - .process_with_scratch(row_input.as_slice_mut().unwrap(), row_output.as_slice_mut().unwrap(), &mut scratch_buffer) - .unwrap(); - }); - - Ok(HArray(ndarray)) - } - _ => Err(HError::OutOfSpecError("The HArray's ndim should be 1 or 2.".into())), - } - } - - fn irfft(&mut self, harray: &mut HArray, IxDyn>, length: usize) -> HResult> { - match harray.ndim() { - 1 => { - let c2r = self.plan_fft_inverse(length); - let scratch_len = c2r.get_scratch_len(); - let mut scratch_buffer = vec![Complex::new(0., 0.); scratch_len]; - let mut ndarray = ArcArray1::from_elem(length, 0.); - c2r.process_with_scratch(harray.as_slice_mut().unwrap(), ndarray.as_slice_mut().unwrap(), &mut scratch_buffer).unwrap(); - Ok(HArray(ndarray.into_dyn())) - }, - 2 => { - let nrows = harray.0.len_of(Axis(0)); - let c2r = self.plan_fft_inverse(length); - let scratch_len = c2r.get_scratch_len(); - let mut scratch_buffer = vec![Complex::new(0., 0.); scratch_len]; - let mut ndarray = ArcArray2::from_elem((nrows, length), 0.).into_dyn(); - - Zip::from(ndarray.lanes_mut(Axis(1))) - .and(harray.0.lanes_mut(Axis(1))) - .for_each(|mut row_output, mut row_input| { - c2r - .process_with_scratch(row_input.as_slice_mut().unwrap(), row_output.as_slice_mut().unwrap(), &mut scratch_buffer) - .unwrap(); - }); - - Ok(HArray(ndarray)) - }, - _ => Err(HError::OutOfSpecError("The HArray's ndim should be 1 or 2.".into())), - } - } + Ok(HArray(ndarray)) } - )+ - }; + _ => Err(HError::OutOfSpecError( + "The HArray's ndim should be 1 or 2.".into(), + )), + } + } } -impl_fft!(f32, f64); +// replace this function by make_mut when in stable (it is currently, but doesn't work for slices.) +fn make_mut_slice(arc: &mut Arc<[T]>) -> &mut [T] { + if Arc::get_mut(arc).is_none() { + *arc = Arc::from(&arc[..]); + } + // Replace by get_mut_unchecked when available in stable. This can't failed since get_mut was + // checked above. + unsafe { Arc::get_mut(arc).unwrap_unchecked() } +} #[cfg(test)] mod tests { @@ -278,8 +383,9 @@ mod tests { Complex::new(11_f32, 12_f32), ]; let mut lhs = HArray::new_from_shape_vec(6, v).unwrap(); - let mut planner = FftPlanner::new(); - planner.fft(&mut lhs).unwrap(); + let length = lhs.len(); + let mut fft = Fft::new_fft_forward(length); + fft.process(&mut lhs).unwrap(); let result = vec![ Complex::new(36.0, 42.0), Complex::new(-16.392305, 4.392305), @@ -302,8 +408,9 @@ mod tests { Complex::new(11_f32, 12_f32), ]; let mut lhs = HArray::new_from_shape_vec((3, 2), v).unwrap(); - let mut planner = FftPlanner::new(); - planner.fft(&mut lhs).unwrap(); + let ncols = lhs.0.ncols(); + let mut fft = Fft::new_fft_forward(ncols); + fft.process(&mut lhs).unwrap(); let result = vec![ Complex::new(4_f32, 6_f32), Complex::new(-2_f32, -2_f32), @@ -328,8 +435,9 @@ mod tests { let mut lhs = HArray::new_from_shape_vec((3, 2), v) .unwrap() .into_dynamic(); - let mut planner = FftPlanner::new(); - planner.fft(&mut lhs).unwrap(); + let ncols = lhs.0.len_of(Axis(1)); + let mut fft = Fft::new_fft_forward(ncols); + fft.process(&mut lhs).unwrap(); let result = vec![ Complex::new(4_f32, 6_f32), Complex::new(-2_f32, -2_f32), @@ -355,8 +463,9 @@ mod tests { Complex::new(11_f32, 12_f32), ]; let mut lhs = HArray::new_from_shape_vec(6, v.clone()).unwrap(); - let mut planner = FftPlanner::new(); - planner.ifft(&mut lhs).unwrap(); + let length = lhs.len(); + let mut fft = Fft::new_fft_inverse(length); + fft.process(&mut lhs).unwrap(); let result = vec![ Complex::new(36.0, 42.0), Complex::new(4.392305, -16.392305), @@ -380,8 +489,9 @@ mod tests { Complex::new(11_f32, 12_f32), ]; let mut lhs = HArray::new_from_shape_vec((3, 2), v).unwrap(); - let mut planner = FftPlanner::new(); - planner.ifft(&mut lhs).unwrap(); + let ncols = lhs.0.ncols(); + let mut fft = Fft::new_fft_forward(ncols); + fft.process(&mut lhs).unwrap(); let result = vec![ Complex::new(4.0, 6.0), Complex::new(-2.0, -2.0), @@ -407,8 +517,9 @@ mod tests { let mut lhs = HArray::new_from_shape_vec((3, 2), v) .unwrap() .into_dynamic(); - let mut planner = FftPlanner::new(); - planner.ifft(&mut lhs).unwrap(); + let ncols = lhs.0.len_of(Axis(1)); + let mut fft = Fft::new_fft_forward(ncols); + fft.process(&mut lhs).unwrap(); let result = vec![ Complex::new(4.0, 6.0), Complex::new(-2.0, -2.0), @@ -427,8 +538,9 @@ mod tests { fn rfft_1d_test() { let v = vec![1_f32, 2_f32, 3_f32, 4_f32, 5_f32, 6_f32]; let mut harray = HArray::new_from_shape_vec(6, v).unwrap(); - let mut planner = RealFftPlanner::new(); - let lhs = planner.rfft(&mut harray).unwrap(); + let length = harray.len(); + let mut rfft = RealFftForward::new_real_fft_forward(length); + let lhs = rfft.process(&mut harray).unwrap(); let result = vec![ Complex::new(21.0, 0.0), Complex::new(-3.0, 5.196152), @@ -445,8 +557,9 @@ mod tests { 1_f32, 2_f32, 3_f32, 4_f32, 5_f32, 6_f32, 7_f32, 8_f32, 9_f32, 10_f32, 11_f32, 12_f32, ]; let mut harray = HArray::new_from_shape_vec((3, 4), v).unwrap(); - let mut planner = RealFftPlanner::new(); - let lhs = planner.rfft(&mut harray).unwrap(); + let ncols = harray.0.ncols(); + let mut rfft = RealFftForward::new_real_fft_forward(ncols); + let lhs = rfft.process(&mut harray).unwrap(); let result = vec![ Complex::new(10_f32, 0_f32), Complex::new(-2_f32, 2_f32), @@ -470,8 +583,9 @@ mod tests { let mut harray = HArray::new_from_shape_vec((3, 4), v) .unwrap() .into_dynamic(); - let mut planner = RealFftPlanner::new(); - let lhs = planner.rfft(&mut harray).unwrap(); + let ncols = harray.0.len_of(Axis(1)); + let mut rfft = RealFftForward::new_real_fft_forward(ncols); + let lhs = rfft.process(&mut harray).unwrap(); let result = vec![ Complex::new(10_f32, 0_f32), Complex::new(-2_f32, 2_f32), @@ -495,9 +609,10 @@ mod tests { let v = vec![1_f32, 2., 3., 4., 5., 6.]; let mut harray = HArray::new_from_shape_vec(6, v).unwrap(); let length = harray.len(); - let mut planner = RealFftPlanner::new(); - let mut spectrum = planner.rfft(&mut harray).unwrap(); - let lhs = planner.irfft(&mut spectrum, length).unwrap(); + let mut rfft = RealFftForward::new_real_fft_forward(length); + let mut spectrum = rfft.process(&mut harray).unwrap(); + let mut irfft = RealFftInverse::new_real_fft_inverse(length); + let lhs = irfft.process(&mut spectrum).unwrap(); let result = vec![6_f32, 12., 18., 24., 30., 36.]; let rhs = HArray::new_from_shape_vec(length, result).unwrap(); assert!(compare_harray(&lhs, &rhs)); @@ -505,9 +620,10 @@ mod tests { let v = vec![1_f32, 2., 3., 4., 5., 6.]; let mut harray = HArray::new_from_shape_vec(6, v).unwrap(); let length = harray.len() + 1; - let mut planner = RealFftPlanner::new(); - let mut spectrum = planner.rfft(&mut harray).unwrap(); - let lhs = planner.irfft(&mut spectrum, length).unwrap(); + let mut rfft = RealFftForward::new_real_fft_forward(length - 1); + let mut spectrum = rfft.process(&mut harray).unwrap(); + let mut irfft = RealFftInverse::new_real_fft_inverse(length); + let lhs = irfft.process(&mut spectrum).unwrap(); let result = vec![ 3.000000000000007_f32, 12.497_72, @@ -527,9 +643,10 @@ mod tests { let v = vec![1_f32, 2., 3., 4., 5., 6.]; let mut harray = HArray::new_from_shape_vec((3, 2), v).unwrap(); let length = harray.0.ncols(); - let mut planner = RealFftPlanner::new(); - let mut spectrum = planner.rfft(&mut harray).unwrap(); - let lhs = planner.irfft(&mut spectrum, length).unwrap(); + let mut rfft = RealFftForward::new_real_fft_forward(length); + let mut spectrum = rfft.process(&mut harray).unwrap(); + let mut irfft = RealFftInverse::new_real_fft_inverse(length); + let lhs = irfft.process(&mut spectrum).unwrap(); let result = vec![2., 4., 6., 8., 10., 12.]; let rhs = HArray::new_from_shape_vec((3, length), result).unwrap(); assert!(compare_harray(&lhs, &rhs)); @@ -537,9 +654,10 @@ mod tests { let v = vec![1_f32, 2., 3., 4., 5., 6.]; let mut harray = HArray::new_from_shape_vec((3, 2), v).unwrap(); let length = harray.0.ncols() + 1; - let mut planner = RealFftPlanner::new(); - let mut spectrum = planner.rfft(&mut harray).unwrap(); - let lhs = planner.irfft(&mut spectrum, length).unwrap(); + let mut rfft = RealFftForward::new_real_fft_forward(length - 1); + let mut spectrum = rfft.process(&mut harray).unwrap(); + let mut irfft = RealFftInverse::new_real_fft_inverse(length); + let lhs = irfft.process(&mut spectrum).unwrap(); let result = vec![1., 4., 4., 5., 8., 8., 9., 12., 12.]; let rhs = HArray::new_from_shape_vec((3, length), result).unwrap(); assert!(compare_harray(&lhs, &rhs)); @@ -553,9 +671,10 @@ mod tests { .unwrap() .into_dynamic(); let length = harray.0.len_of(Axis(1)); - let mut planner = RealFftPlanner::new(); - let mut spectrum = planner.rfft(&mut harray).unwrap(); - let lhs = planner.irfft(&mut spectrum, length).unwrap(); + let mut rfft = RealFftForward::new_real_fft_forward(length); + let mut spectrum = rfft.process(&mut harray).unwrap(); + let mut irfft = RealFftInverse::new_real_fft_inverse(length); + let lhs = irfft.process(&mut spectrum).unwrap(); let result = vec![2., 4., 6., 8., 10., 12.]; let rhs = HArray::new_from_shape_vec((3, length), result) .unwrap() @@ -567,9 +686,10 @@ mod tests { .unwrap() .into_dynamic(); let length = harray.0.len_of(Axis(1)) + 1; - let mut planner = RealFftPlanner::new(); - let mut spectrum = planner.rfft(&mut harray).unwrap(); - let lhs = planner.irfft(&mut spectrum, length).unwrap(); + let mut rfft = RealFftForward::new_real_fft_forward(length - 1); + let mut spectrum = rfft.process(&mut harray).unwrap(); + let mut irfft = RealFftInverse::new_real_fft_inverse(length); + let lhs = irfft.process(&mut spectrum).unwrap(); let result = vec![1., 4., 4., 5., 8., 8., 9., 12., 12.]; let rhs = HArray::new_from_shape_vec((3, length), result) .unwrap() diff --git a/harmonium-stft/Cargo.toml b/harmonium-stft/Cargo.toml new file mode 100644 index 0000000..9f1e3bd --- /dev/null +++ b/harmonium-stft/Cargo.toml @@ -0,0 +1,9 @@ +[package] +name = "harmonium-stft" +version = "0.1.0" +edition = "2021" + +[dependencies] +harmonium-core = { workspace = true } +rustfft = { workspace = true } +ndarray = { workspace = true } diff --git a/harmonium-stft/src/lib.rs b/harmonium-stft/src/lib.rs new file mode 100644 index 0000000..0a6ba6c --- /dev/null +++ b/harmonium-stft/src/lib.rs @@ -0,0 +1 @@ +mod stft; diff --git a/harmonium-stft/src/stft.rs b/harmonium-stft/src/stft.rs new file mode 100644 index 0000000..2f81467 --- /dev/null +++ b/harmonium-stft/src/stft.rs @@ -0,0 +1,101 @@ +use harmonium_core::array::HArray; +use ndarray::{ArcArray2, Dimension, Ix1, Ix2, IxDyn}; +use rustfft::{ + num_complex::{Complex, ComplexFloat}, + num_traits::Float, + Fft, FftNum, FftPlanner, +}; +use std::sync::Arc; + +struct StftPlanner +where + T: FftNum, +{ + fft: Arc>, + buffer: Vec, +} + +impl StftPlanner +where + T: FftNum + ComplexFloat, +{ + fn new(fft_length: usize) -> Self { + let mut fft_planner = FftPlanner::new(); + let fft = fft_planner.plan_fft_forward(fft_length); + let buffer = vec![T::zero(); fft_length]; + StftPlanner { fft, buffer } + } + + fn len(&self) -> usize { + self.fft.len() + } +} + +// The input buffer is used as scratch space, so the contents of input should be considered garbage after calling. +pub trait Stft +where + T: FftNum + ComplexFloat, + D: Dimension, +{ + fn process( + &mut self, + harray: HArray, + hop_length: usize, + window_length: usize, + window: Option<&[T]>, + ) -> HArray; +} + +impl Stft for StftPlanner +where + T: FftNum + ComplexFloat, +{ + fn process( + &mut self, + harray: HArray, + hop_length: usize, + window_length: usize, + window: Option<&[T]>, + ) -> HArray { + assert!(hop_length > 0); + + let fft_length = self.len(); + let n_fft = (harray.len() - window_length) / hop_length + 1; + let stft_ndarray = ArcArray2::zeros((fft_length, n_fft)); + todo!() + } +} + +trait ApplyWindow { + fn apply_window(&mut self, window: &[T]); +} + +impl ApplyWindow for [Complex] +where + T: Float, +{ + fn apply_window(&mut self, window: &[T]) { + for (i, w) in self.iter_mut().zip(window.iter()) { + *i = (*i).scale(*w); + } + } +} + +impl ApplyWindow for [T] +where + T: Float, +{ + fn apply_window(&mut self, window: &[T]) { + for (i, w) in self.iter_mut().zip(window.iter()) { + *i = *i * *w; + } + } +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn it_works() {} +} diff --git a/harmonium-window/src/windows.rs b/harmonium-window/src/windows.rs index ff244e7..48e7608 100644 --- a/harmonium-window/src/windows.rs +++ b/harmonium-window/src/windows.rs @@ -225,7 +225,7 @@ where let pi = T::PI(); let zero = T::zero(); - let one = T::from(1.0).unwrap(); + let one = T::one(); let two = T::from(2.0).unwrap(); let step = two / (np_f); let mut fac = -one; @@ -255,7 +255,7 @@ pub fn boxcar(npoints: usize) -> HArray where T: Float + FloatConst, { - let one = T::from(1.0).unwrap(); + let one = T::one(); let window: Vec = (0..npoints).map(|_| one).collect(); HArray::new_from_shape_vec(npoints, window).unwrap() @@ -336,7 +336,7 @@ where let pi = T::PI(); let zero = T::zero(); - let one = T::from(1.0).unwrap(); + let one = T::one(); let two = T::from(2.0).unwrap(); let ten = T::from(10.0).unwrap(); let twenty = ten + ten; @@ -512,7 +512,7 @@ where T: Float + FloatConst, { let np_f = T::from(npoints).unwrap(); - let one = T::from(1.0).unwrap(); + let one = T::one(); let two = T::from(2.0).unwrap(); let mut window: Vec = Vec::with_capacity(npoints); diff --git a/r-harmonium/.fft_bench.R b/r-harmonium/.fft_bench.R index e432e84..80378bb 100644 --- a/r-harmonium/.fft_bench.R +++ b/r-harmonium/.fft_bench.R @@ -5,18 +5,18 @@ devtools::load_all(".", export_all = FALSE) # fft_matrix with complexes results <- bench::press( - n = seq(30, 400000, 30000), + n = seq.int(30, 400000, 30000L), { r = as.double(sample(100, n, replace = TRUE)) i = as.double(sample(100, n, replace = TRUE)) x = complex(real=r, imaginary=i) x = matrix(x, ncol = 10) - fft_planner = HFftPlanner$new(HDataType$Complex64) + hfft = HFft$new_fft_forward(as.integer(n/10), HDataType$Complex64) mark( torch = as_array(torch::torch_fft_fft(torch_tensor(x, dtype = torch_cfloat64()), dim = 1)), harmonium_fft = { harray = HArray$new_from_values(x, HDataType$Complex64) - fft_planner$fft(harray) + hfft$process(harray) harray$collect() }, base_r = stats::mvfft(x), @@ -33,12 +33,12 @@ results <- bench::press( { x = as.double(sample(100, n, replace = TRUE)) x = matrix(x, ncol = 10) - real_fft_planner = HRealFftPlanner$new(HDataType$Float64) + hfft = HRealFft$new_real_fft(as.integer(n/10), HDataType$Float64) mark( torch = as_array(torch::torch_fft_rfft(torch_tensor(x, dtype = torch_float64()), dim = 1)), harmonium_fft = { harray = HArray$new_from_values(x, HDataType$Float64) - real_fft_planner$rfft(harray) + hfft$process(harray) harray$collect() }, iterations = 50, diff --git a/r-harmonium/DESCRIPTION b/r-harmonium/DESCRIPTION index 962f51f..2474cc1 100644 --- a/r-harmonium/DESCRIPTION +++ b/r-harmonium/DESCRIPTION @@ -1,6 +1,6 @@ Package: harmonium Title: Audio analysis and I/O -Version: 0.1.1 +Version: 0.2.0 Authors@R: person("Daniel", "Gurgel", email = "daniellga@gmail.com", role = c("aut", "cre")) Description: Audio analysis and I/O. diff --git a/r-harmonium/NAMESPACE b/r-harmonium/NAMESPACE index e7389b0..44a7a58 100644 --- a/r-harmonium/NAMESPACE +++ b/r-harmonium/NAMESPACE @@ -11,8 +11,8 @@ export(HPolynomialDegree) export(HWindow) export(HMetadataType) export(HFile) -export(HFftPlanner) -export(HRealFftPlanner) +export(HFft) +export(HRealFft) export(HAudioOp) export(hdocs) export(HDecoderStream) @@ -43,5 +43,7 @@ S3method(print,HPolynomialDegree) S3method(print,HWindowType) S3method(print,HInterpolationType) S3method(print,HSincInterpolationParameters) +S3method(print,HFft) +S3method(print,HRealFft) useDynLib(harmonium, .registration = TRUE) diff --git a/r-harmonium/R/000-wrappers.R b/r-harmonium/R/000-wrappers.R index c17025b..cf1279d 100644 --- a/r-harmonium/R/000-wrappers.R +++ b/r-harmonium/R/000-wrappers.R @@ -663,86 +663,104 @@ class(`HDecoderStream`) <- "HDecoderStream__bundle" #' @export `[[<-.HDecoderStream__bundle` <- function(x, i, value) stop("HDecoderStream cannot be modified", call. = FALSE) -### wrapper functions for HFftPlanner +### wrapper functions for HFft -`HFftPlanner_fft` <- function(self) { +`HFft_process` <- function(self) { function(`harray`) { `harray` <- .savvy_extract_ptr(`harray`, "HArray") - invisible(.Call(savvy_HFftPlanner_fft__impl, `self`, `harray`)) + invisible(.Call(savvy_HFft_process__impl, `self`, `harray`)) } } -`HFftPlanner_ifft` <- function(self) { - function(`harray`) { - `harray` <- .savvy_extract_ptr(`harray`, "HArray") - invisible(.Call(savvy_HFftPlanner_ifft__impl, `self`, `harray`)) +`HFft_dtype` <- function(self) { + function() { + .savvy_wrap_HDataType(.Call(savvy_HFft_dtype__impl, `self`)) + } +} + +`HFft_print` <- function(self) { + function() { + invisible(.Call(savvy_HFft_print__impl, `self`)) + } +} + +`HFft_clone` <- function(self) { + function() { + .savvy_wrap_HFft(.Call(savvy_HFft_clone__impl, `self`)) } } -`HFftPlanner_dtype` <- function(self) { +`HFft_is_shared` <- function(self) { function() { - .savvy_wrap_HDataType(.Call(savvy_HFftPlanner_dtype__impl, `self`)) + .Call(savvy_HFft_is_shared__impl, `self`) } } -`HFftPlanner_print` <- function(self) { +`HFft_invalidate` <- function(self) { function() { - invisible(.Call(savvy_HFftPlanner_print__impl, `self`)) + invisible(.Call(savvy_HFft_invalidate__impl, `self`)) } } -`.savvy_wrap_HFftPlanner` <- function(ptr) { +`.savvy_wrap_HFft` <- function(ptr) { e <- new.env(parent = emptyenv()) e$.ptr <- ptr - e$`fft` <- `HFftPlanner_fft`(ptr) - e$`ifft` <- `HFftPlanner_ifft`(ptr) - e$`dtype` <- `HFftPlanner_dtype`(ptr) - e$`print` <- `HFftPlanner_print`(ptr) - - class(e) <- "HFftPlanner" + e$`process` <- `HFft_process`(ptr) + e$`dtype` <- `HFft_dtype`(ptr) + e$`print` <- `HFft_print`(ptr) + e$`clone` <- `HFft_clone`(ptr) + e$`is_shared` <- `HFft_is_shared`(ptr) + e$`invalidate` <- `HFft_invalidate`(ptr) + + class(e) <- "HFft" e } #' @export -`$<-.HFftPlanner` <- function(x, name, value) stop("HFftPlanner cannot be modified", call. = FALSE) +`$<-.HFft` <- function(x, name, value) stop("HFft cannot be modified", call. = FALSE) #' @export -`[[<-.HFftPlanner` <- function(x, i, value) stop("HFftPlanner cannot be modified", call. = FALSE) +`[[<-.HFft` <- function(x, i, value) stop("HFft cannot be modified", call. = FALSE) -#' HFftPlanner -#' A planner is used to create FFTs. It caches results internally, so when making more than one FFT it is advisable to reuse the same planner. +#' HFft +#' An `HFft` is used to create FFTs. It caches results internally, so when making more than one FFT it is advisable to reuse the same `HFft` instance. #' #' # Methods #' -`HFftPlanner` <- new.env(parent = emptyenv()) +`HFft` <- new.env(parent = emptyenv()) #' @export -`$<-.HFftPlanner` <- function(x, name, value) stop("HFftPlanner cannot be modified", call. = FALSE) +`$<-.HFft` <- function(x, name, value) stop("HFft cannot be modified", call. = FALSE) #' @export -`[[<-.HFftPlanner` <- function(x, i, value) stop("HFftPlanner cannot be modified", call. = FALSE) +`[[<-.HFft` <- function(x, i, value) stop("HFft cannot be modified", call. = FALSE) -### associated functions for HFftPlanner +### associated functions for HFft -`HFftPlanner`$`new` <- function(`dtype`) { +`HFft`$`new_fft_forward` <- function(`length`, `dtype`) { `dtype` <- .savvy_extract_ptr(`dtype`, "HDataType") - .savvy_wrap_HFftPlanner(.Call(savvy_HFftPlanner_new__impl, `dtype`)) + .savvy_wrap_HFft(.Call(savvy_HFft_new_fft_forward__impl, `length`, `dtype`)) } +`HFft`$`new_fft_inverse` <- function(`length`, `dtype`) { + `dtype` <- .savvy_extract_ptr(`dtype`, "HDataType") + .savvy_wrap_HFft(.Call(savvy_HFft_new_fft_inverse__impl, `length`, `dtype`)) +} -class(`HFftPlanner`) <- "HFftPlanner__bundle" + +class(`HFft`) <- "HFft__bundle" #' @export -`print.HFftPlanner__bundle` <- function(x, ...) { - cat('HFftPlanner') +`print.HFft__bundle` <- function(x, ...) { + cat('HFft') } #' @export -`$<-.HFftPlanner__bundle` <- function(x, name, value) stop("HFftPlanner cannot be modified", call. = FALSE) +`$<-.HFft__bundle` <- function(x, name, value) stop("HFft cannot be modified", call. = FALSE) #' @export -`[[<-.HFftPlanner__bundle` <- function(x, i, value) stop("HFftPlanner cannot be modified", call. = FALSE) +`[[<-.HFft__bundle` <- function(x, i, value) stop("HFft cannot be modified", call. = FALSE) ### wrapper functions for HFile @@ -1160,88 +1178,95 @@ class(`HPolynomialDegree`) <- "HPolynomialDegree__bundle" #' @export `[[<-.HPolynomialDegree__bundle` <- function(x, i, value) stop("HPolynomialDegree cannot be modified", call. = FALSE) -### wrapper functions for HRealFftPlanner +### wrapper functions for HRealFft -`HRealFftPlanner_rfft` <- function(self) { +`HRealFft_process` <- function(self) { function(`harray`) { `harray` <- .savvy_extract_ptr(`harray`, "HArray") - invisible(.Call(savvy_HRealFftPlanner_rfft__impl, `self`, `harray`)) + invisible(.Call(savvy_HRealFft_process__impl, `self`, `harray`)) } } -`HRealFftPlanner_irfft` <- function(self) { - function(`harray`, `length`) { - `harray` <- .savvy_extract_ptr(`harray`, "HArray") - invisible(.Call(savvy_HRealFftPlanner_irfft__impl, `self`, `harray`, `length`)) +`HRealFft_dtype` <- function(self) { + function() { + .savvy_wrap_HDataType(.Call(savvy_HRealFft_dtype__impl, `self`)) + } +} + +`HRealFft_print` <- function(self) { + function() { + invisible(.Call(savvy_HRealFft_print__impl, `self`)) } } -`HRealFftPlanner_dtype` <- function(self) { +`HRealFft_clone` <- function(self) { function() { - .savvy_wrap_HDataType(.Call(savvy_HRealFftPlanner_dtype__impl, `self`)) + .savvy_wrap_HRealFft(.Call(savvy_HRealFft_clone__impl, `self`)) } } -`HRealFftPlanner_print` <- function(self) { +`HRealFft_is_shared` <- function(self) { function() { - invisible(.Call(savvy_HRealFftPlanner_print__impl, `self`)) + .Call(savvy_HRealFft_is_shared__impl, `self`) } } -`.savvy_wrap_HRealFftPlanner` <- function(ptr) { +`HRealFft_invalidate` <- function(self) { + function() { + invisible(.Call(savvy_HRealFft_invalidate__impl, `self`)) + } +} + +`.savvy_wrap_HRealFft` <- function(ptr) { e <- new.env(parent = emptyenv()) e$.ptr <- ptr - e$`rfft` <- `HRealFftPlanner_rfft`(ptr) - e$`irfft` <- `HRealFftPlanner_irfft`(ptr) - e$`dtype` <- `HRealFftPlanner_dtype`(ptr) - e$`print` <- `HRealFftPlanner_print`(ptr) - - class(e) <- "HRealFftPlanner" + e$`process` <- `HRealFft_process`(ptr) + e$`dtype` <- `HRealFft_dtype`(ptr) + e$`print` <- `HRealFft_print`(ptr) + e$`clone` <- `HRealFft_clone`(ptr) + e$`is_shared` <- `HRealFft_is_shared`(ptr) + e$`invalidate` <- `HRealFft_invalidate`(ptr) + + class(e) <- "HRealFft" e } #' @export -`$<-.HRealFftPlanner` <- function(x, name, value) stop("HRealFftPlanner cannot be modified", call. = FALSE) +`$<-.HRealFft` <- function(x, name, value) stop("HRealFft cannot be modified", call. = FALSE) #' @export -`[[<-.HRealFftPlanner` <- function(x, i, value) stop("HRealFftPlanner cannot be modified", call. = FALSE) +`[[<-.HRealFft` <- function(x, i, value) stop("HRealFft cannot be modified", call. = FALSE) -#' HRealFftPlanner -#' A planner is used to create FFTs. It caches results internally, so when making more than one FFT it is advisable to reuse the same planner. -#' -#' This planner is used to calculate FFTs of real valued inputs and its inverse operation. -#' -#' # Methods -#' -`HRealFftPlanner` <- new.env(parent = emptyenv()) + +`HRealFft` <- new.env(parent = emptyenv()) #' @export -`$<-.HRealFftPlanner` <- function(x, name, value) stop("HRealFftPlanner cannot be modified", call. = FALSE) +`$<-.HRealFft` <- function(x, name, value) stop("HRealFft cannot be modified", call. = FALSE) #' @export -`[[<-.HRealFftPlanner` <- function(x, i, value) stop("HRealFftPlanner cannot be modified", call. = FALSE) +`[[<-.HRealFft` <- function(x, i, value) stop("HRealFft cannot be modified", call. = FALSE) -### associated functions for HRealFftPlanner +### associated functions for HRealFft -`HRealFftPlanner`$`new` <- function(`dtype`) { +`HRealFft`$`new_real_fft` <- function(`length`, `dtype`) { `dtype` <- .savvy_extract_ptr(`dtype`, "HDataType") - .savvy_wrap_HRealFftPlanner(.Call(savvy_HRealFftPlanner_new__impl, `dtype`)) + .savvy_wrap_HRealFft(.Call(savvy_HRealFft_new_real_fft__impl, `length`, `dtype`)) } -class(`HRealFftPlanner`) <- "HRealFftPlanner__bundle" +class(`HRealFft`) <- "HRealFft__bundle" #' @export -`print.HRealFftPlanner__bundle` <- function(x, ...) { - cat('HRealFftPlanner') +`print.HRealFft__bundle` <- function(x, ...) { + cat('HRealFft') } #' @export -`$<-.HRealFftPlanner__bundle` <- function(x, name, value) stop("HRealFftPlanner cannot be modified", call. = FALSE) +`$<-.HRealFft__bundle` <- function(x, name, value) stop("HRealFft cannot be modified", call. = FALSE) #' @export -`[[<-.HRealFftPlanner__bundle` <- function(x, i, value) stop("HRealFftPlanner cannot be modified", call. = FALSE) +`[[<-.HRealFft__bundle` <- function(x, i, value) stop("HRealFft cannot be modified", call. = FALSE) ### wrapper functions for HResampler diff --git a/r-harmonium/R/structs.R b/r-harmonium/R/structs.R index 2538561..4ed2ae1 100644 --- a/r-harmonium/R/structs.R +++ b/r-harmonium/R/structs.R @@ -50,10 +50,11 @@ print.HInterpolationType = function(x, ...) { "==.HInterpolationType" <- function(e1,e2) e1$eq(e2) "!=.HInterpolationType" <- function(e1,e2) e1$ne(e2) -print.HFftPlanner = function(x, ...) { +print.HFft = function(x, ...) { x$print() } -print.HRealFftPlanner = function(x, ...) { +print.HRealFft = function(x, ...) { x$print() } + diff --git a/r-harmonium/R/zzz.R b/r-harmonium/R/zzz.R index 43f22c0..22eb4c1 100644 --- a/r-harmonium/R/zzz.R +++ b/r-harmonium/R/zzz.R @@ -10,8 +10,8 @@ lockEnvironment(HWindow, bindings = TRUE) lockEnvironment(HFile, bindings = TRUE) lockEnvironment(HResampler, bindings = TRUE) - lockEnvironment(HFftPlanner, bindings = TRUE) - lockEnvironment(HRealFftPlanner, bindings = TRUE) + lockEnvironment(HFft, bindings = TRUE) + lockEnvironment(HRealFft, bindings = TRUE) lockEnvironment(HDecoderStream, bindings = TRUE) } diff --git a/r-harmonium/src/init.c b/r-harmonium/src/init.c index 87d17ee..35d5749 100644 --- a/r-harmonium/src/init.c +++ b/r-harmonium/src/init.c @@ -275,28 +275,43 @@ SEXP savvy_HDecoderStream_stream__impl(SEXP self__) { return handle_result(res); } -SEXP savvy_HFftPlanner_new__impl(SEXP dtype) { - SEXP res = savvy_HFftPlanner_new__ffi(dtype); +SEXP savvy_HFft_new_fft_forward__impl(SEXP length, SEXP dtype) { + SEXP res = savvy_HFft_new_fft_forward__ffi(length, dtype); return handle_result(res); } -SEXP savvy_HFftPlanner_fft__impl(SEXP self__, SEXP harray) { - SEXP res = savvy_HFftPlanner_fft__ffi(self__, harray); +SEXP savvy_HFft_new_fft_inverse__impl(SEXP length, SEXP dtype) { + SEXP res = savvy_HFft_new_fft_inverse__ffi(length, dtype); return handle_result(res); } -SEXP savvy_HFftPlanner_ifft__impl(SEXP self__, SEXP harray) { - SEXP res = savvy_HFftPlanner_ifft__ffi(self__, harray); +SEXP savvy_HFft_process__impl(SEXP self__, SEXP harray) { + SEXP res = savvy_HFft_process__ffi(self__, harray); return handle_result(res); } -SEXP savvy_HFftPlanner_dtype__impl(SEXP self__) { - SEXP res = savvy_HFftPlanner_dtype__ffi(self__); +SEXP savvy_HFft_dtype__impl(SEXP self__) { + SEXP res = savvy_HFft_dtype__ffi(self__); return handle_result(res); } -SEXP savvy_HFftPlanner_print__impl(SEXP self__) { - SEXP res = savvy_HFftPlanner_print__ffi(self__); +SEXP savvy_HFft_print__impl(SEXP self__) { + SEXP res = savvy_HFft_print__ffi(self__); + return handle_result(res); +} + +SEXP savvy_HFft_clone__impl(SEXP self__) { + SEXP res = savvy_HFft_clone__ffi(self__); + return handle_result(res); +} + +SEXP savvy_HFft_is_shared__impl(SEXP self__) { + SEXP res = savvy_HFft_is_shared__ffi(self__); + return handle_result(res); +} + +SEXP savvy_HFft_invalidate__impl(SEXP self__) { + SEXP res = savvy_HFft_invalidate__ffi(self__); return handle_result(res); } @@ -370,28 +385,38 @@ SEXP savvy_HPolynomialDegree_ne__impl(SEXP self__, SEXP other) { return handle_result(res); } -SEXP savvy_HRealFftPlanner_new__impl(SEXP dtype) { - SEXP res = savvy_HRealFftPlanner_new__ffi(dtype); +SEXP savvy_HRealFft_new_real_fft__impl(SEXP length, SEXP dtype) { + SEXP res = savvy_HRealFft_new_real_fft__ffi(length, dtype); + return handle_result(res); +} + +SEXP savvy_HRealFft_process__impl(SEXP self__, SEXP harray) { + SEXP res = savvy_HRealFft_process__ffi(self__, harray); + return handle_result(res); +} + +SEXP savvy_HRealFft_dtype__impl(SEXP self__) { + SEXP res = savvy_HRealFft_dtype__ffi(self__); return handle_result(res); } -SEXP savvy_HRealFftPlanner_rfft__impl(SEXP self__, SEXP harray) { - SEXP res = savvy_HRealFftPlanner_rfft__ffi(self__, harray); +SEXP savvy_HRealFft_print__impl(SEXP self__) { + SEXP res = savvy_HRealFft_print__ffi(self__); return handle_result(res); } -SEXP savvy_HRealFftPlanner_irfft__impl(SEXP self__, SEXP harray, SEXP length) { - SEXP res = savvy_HRealFftPlanner_irfft__ffi(self__, harray, length); +SEXP savvy_HRealFft_clone__impl(SEXP self__) { + SEXP res = savvy_HRealFft_clone__ffi(self__); return handle_result(res); } -SEXP savvy_HRealFftPlanner_dtype__impl(SEXP self__) { - SEXP res = savvy_HRealFftPlanner_dtype__ffi(self__); +SEXP savvy_HRealFft_is_shared__impl(SEXP self__) { + SEXP res = savvy_HRealFft_is_shared__ffi(self__); return handle_result(res); } -SEXP savvy_HRealFftPlanner_print__impl(SEXP self__) { - SEXP res = savvy_HRealFftPlanner_print__ffi(self__); +SEXP savvy_HRealFft_invalidate__impl(SEXP self__) { + SEXP res = savvy_HRealFft_invalidate__ffi(self__); return handle_result(res); } @@ -576,11 +601,14 @@ static const R_CallMethodDef CallEntries[] = { {"savvy_HDecodedAudio_sr__impl", (DL_FUNC) &savvy_HDecodedAudio_sr__impl, 1}, {"savvy_HDecodedAudio_invalidate__impl", (DL_FUNC) &savvy_HDecodedAudio_invalidate__impl, 1}, {"savvy_HDecoderStream_stream__impl", (DL_FUNC) &savvy_HDecoderStream_stream__impl, 1}, - {"savvy_HFftPlanner_new__impl", (DL_FUNC) &savvy_HFftPlanner_new__impl, 1}, - {"savvy_HFftPlanner_fft__impl", (DL_FUNC) &savvy_HFftPlanner_fft__impl, 2}, - {"savvy_HFftPlanner_ifft__impl", (DL_FUNC) &savvy_HFftPlanner_ifft__impl, 2}, - {"savvy_HFftPlanner_dtype__impl", (DL_FUNC) &savvy_HFftPlanner_dtype__impl, 1}, - {"savvy_HFftPlanner_print__impl", (DL_FUNC) &savvy_HFftPlanner_print__impl, 1}, + {"savvy_HFft_new_fft_forward__impl", (DL_FUNC) &savvy_HFft_new_fft_forward__impl, 2}, + {"savvy_HFft_new_fft_inverse__impl", (DL_FUNC) &savvy_HFft_new_fft_inverse__impl, 2}, + {"savvy_HFft_process__impl", (DL_FUNC) &savvy_HFft_process__impl, 2}, + {"savvy_HFft_dtype__impl", (DL_FUNC) &savvy_HFft_dtype__impl, 1}, + {"savvy_HFft_print__impl", (DL_FUNC) &savvy_HFft_print__impl, 1}, + {"savvy_HFft_clone__impl", (DL_FUNC) &savvy_HFft_clone__impl, 1}, + {"savvy_HFft_is_shared__impl", (DL_FUNC) &savvy_HFft_is_shared__impl, 1}, + {"savvy_HFft_invalidate__impl", (DL_FUNC) &savvy_HFft_invalidate__impl, 1}, {"savvy_HFile_decode__impl", (DL_FUNC) &savvy_HFile_decode__impl, 2}, {"savvy_HFile_decode_stream__impl", (DL_FUNC) &savvy_HFile_decode_stream__impl, 3}, {"savvy_HFile_metadata__impl", (DL_FUNC) &savvy_HFile_metadata__impl, 2}, @@ -595,11 +623,13 @@ static const R_CallMethodDef CallEntries[] = { {"savvy_HPolynomialDegree_print__impl", (DL_FUNC) &savvy_HPolynomialDegree_print__impl, 1}, {"savvy_HPolynomialDegree_eq__impl", (DL_FUNC) &savvy_HPolynomialDegree_eq__impl, 2}, {"savvy_HPolynomialDegree_ne__impl", (DL_FUNC) &savvy_HPolynomialDegree_ne__impl, 2}, - {"savvy_HRealFftPlanner_new__impl", (DL_FUNC) &savvy_HRealFftPlanner_new__impl, 1}, - {"savvy_HRealFftPlanner_rfft__impl", (DL_FUNC) &savvy_HRealFftPlanner_rfft__impl, 2}, - {"savvy_HRealFftPlanner_irfft__impl", (DL_FUNC) &savvy_HRealFftPlanner_irfft__impl, 3}, - {"savvy_HRealFftPlanner_dtype__impl", (DL_FUNC) &savvy_HRealFftPlanner_dtype__impl, 1}, - {"savvy_HRealFftPlanner_print__impl", (DL_FUNC) &savvy_HRealFftPlanner_print__impl, 1}, + {"savvy_HRealFft_new_real_fft__impl", (DL_FUNC) &savvy_HRealFft_new_real_fft__impl, 2}, + {"savvy_HRealFft_process__impl", (DL_FUNC) &savvy_HRealFft_process__impl, 2}, + {"savvy_HRealFft_dtype__impl", (DL_FUNC) &savvy_HRealFft_dtype__impl, 1}, + {"savvy_HRealFft_print__impl", (DL_FUNC) &savvy_HRealFft_print__impl, 1}, + {"savvy_HRealFft_clone__impl", (DL_FUNC) &savvy_HRealFft_clone__impl, 1}, + {"savvy_HRealFft_is_shared__impl", (DL_FUNC) &savvy_HRealFft_is_shared__impl, 1}, + {"savvy_HRealFft_invalidate__impl", (DL_FUNC) &savvy_HRealFft_invalidate__impl, 1}, {"savvy_HResampler_new_fft__impl", (DL_FUNC) &savvy_HResampler_new_fft__impl, 7}, {"savvy_HResampler_new_sinc__impl", (DL_FUNC) &savvy_HResampler_new_sinc__impl, 7}, {"savvy_HResampler_new_fast__impl", (DL_FUNC) &savvy_HResampler_new_fast__impl, 7}, diff --git a/r-harmonium/src/rust/Cargo.toml b/r-harmonium/src/rust/Cargo.toml index 35fe3a5..b75502d 100644 --- a/r-harmonium/src/rust/Cargo.toml +++ b/r-harmonium/src/rust/Cargo.toml @@ -1,6 +1,7 @@ [package] name = 'r-harmonium' edition = '2021' +version = '0.2.0' [lib] crate-type = [ 'staticlib' ] @@ -16,16 +17,16 @@ panic = "unwind" # otherwise savvy crashes R when panic. # prevents package from thinking it's in the workspace. [dependencies] -# harmonium-core = { path = "../../../harmonium-core" } -# harmonium-io = { path = "../../../harmonium-io" } -# harmonium-resample = { path = "../../../harmonium-resample" } -# harmonium-fft = { path = "../../../harmonium-fft" } -# harmonium-window = { path = "../../../harmonium-window" } -harmonium-core = { git = "https://github.com/daniellga/harmonium.git" } -harmonium-io = { git = "https://github.com/daniellga/harmonium.git", features = ["all", "opt-simd"] } -harmonium-resample = { git = "https://github.com/daniellga/harmonium.git" } -harmonium-fft = { git = "https://github.com/daniellga/harmonium.git" } -harmonium-window = { git = "https://github.com/daniellga/harmonium.git" } +harmonium-core = { path = "../../../harmonium-core" } +harmonium-io = { path = "../../../harmonium-io" } +harmonium-resample = { path = "../../../harmonium-resample" } +harmonium-fft = { path = "../../../harmonium-fft" } +harmonium-window = { path = "../../../harmonium-window" } +#harmonium-core = { git = "https://github.com/daniellga/harmonium.git" } +#harmonium-io = { git = "https://github.com/daniellga/harmonium.git", features = ["all", "opt-simd"] } +#harmonium-resample = { git = "https://github.com/daniellga/harmonium.git" } +#harmonium-fft = { git = "https://github.com/daniellga/harmonium.git" } +#harmonium-window = { git = "https://github.com/daniellga/harmonium.git" } rubato = "0.15.0" rustfft = "6.2" realfft = "3.3" diff --git a/r-harmonium/src/rust/api.h b/r-harmonium/src/rust/api.h index 064c20a..e0fad9f 100644 --- a/r-harmonium/src/rust/api.h +++ b/r-harmonium/src/rust/api.h @@ -60,12 +60,15 @@ SEXP savvy_HDecodedAudio_invalidate__ffi(SEXP self__); // methods and associated functions for HDecoderStream SEXP savvy_HDecoderStream_stream__ffi(SEXP self__); -// methods and associated functions for HFftPlanner -SEXP savvy_HFftPlanner_new__ffi(SEXP dtype); -SEXP savvy_HFftPlanner_fft__ffi(SEXP self__, SEXP harray); -SEXP savvy_HFftPlanner_ifft__ffi(SEXP self__, SEXP harray); -SEXP savvy_HFftPlanner_dtype__ffi(SEXP self__); -SEXP savvy_HFftPlanner_print__ffi(SEXP self__); +// methods and associated functions for HFft +SEXP savvy_HFft_new_fft_forward__ffi(SEXP length, SEXP dtype); +SEXP savvy_HFft_new_fft_inverse__ffi(SEXP length, SEXP dtype); +SEXP savvy_HFft_process__ffi(SEXP self__, SEXP harray); +SEXP savvy_HFft_dtype__ffi(SEXP self__); +SEXP savvy_HFft_print__ffi(SEXP self__); +SEXP savvy_HFft_clone__ffi(SEXP self__); +SEXP savvy_HFft_is_shared__ffi(SEXP self__); +SEXP savvy_HFft_invalidate__ffi(SEXP self__); // methods and associated functions for HFile SEXP savvy_HFile_decode__ffi(SEXP fpath, SEXP dtype); @@ -89,12 +92,14 @@ SEXP savvy_HPolynomialDegree_print__ffi(SEXP self__); SEXP savvy_HPolynomialDegree_eq__ffi(SEXP self__, SEXP other); SEXP savvy_HPolynomialDegree_ne__ffi(SEXP self__, SEXP other); -// methods and associated functions for HRealFftPlanner -SEXP savvy_HRealFftPlanner_new__ffi(SEXP dtype); -SEXP savvy_HRealFftPlanner_rfft__ffi(SEXP self__, SEXP harray); -SEXP savvy_HRealFftPlanner_irfft__ffi(SEXP self__, SEXP harray, SEXP length); -SEXP savvy_HRealFftPlanner_dtype__ffi(SEXP self__); -SEXP savvy_HRealFftPlanner_print__ffi(SEXP self__); +// methods and associated functions for HRealFft +SEXP savvy_HRealFft_new_real_fft__ffi(SEXP length, SEXP dtype); +SEXP savvy_HRealFft_process__ffi(SEXP self__, SEXP harray); +SEXP savvy_HRealFft_dtype__ffi(SEXP self__); +SEXP savvy_HRealFft_print__ffi(SEXP self__); +SEXP savvy_HRealFft_clone__ffi(SEXP self__); +SEXP savvy_HRealFft_is_shared__ffi(SEXP self__); +SEXP savvy_HRealFft_invalidate__ffi(SEXP self__); // methods and associated functions for HResampler SEXP savvy_HResampler_new_fft__ffi(SEXP sr_in, SEXP sr_out, SEXP chunk_size, SEXP sub_chunks, SEXP nchannels, SEXP res_type, SEXP dtype); diff --git a/r-harmonium/src/rust/src/conversions.rs b/r-harmonium/src/rust/src/conversions.rs index 59af9ac..e701a1a 100644 --- a/r-harmonium/src/rust/src/conversions.rs +++ b/r-harmonium/src/rust/src/conversions.rs @@ -29,7 +29,7 @@ impl ToScalar for Sexp { Ok(integer_sexp.as_slice()[0]) } _ => { - let err = "Argument must be a string of length 1.".to_string(); + let err = "Argument must be an integer of length 1.".to_string(); Err(err.into()) } } @@ -41,7 +41,7 @@ impl ToScalar for Sexp { match self.into_typed() { TypedSexp::Real(real_sexp) if real_sexp.len() == 1 => Ok(real_sexp.as_slice()[0]), _ => { - let err = "Argument must be a string of length 1.".to_string(); + let err = "Argument must be a double of length 1.".to_string(); Err(err.into()) } } @@ -55,7 +55,7 @@ impl ToScalar for Sexp { Ok(logical_sexp.as_slice_raw()[0] == 1) } _ => { - let err = "Argument must be a string of length 1.".to_string(); + let err = "Argument must be a logical of length 1.".to_string(); Err(err.into()) } } diff --git a/r-harmonium/src/rust/src/harray.rs b/r-harmonium/src/rust/src/harray.rs index 950aa30..5d649f6 100644 --- a/r-harmonium/src/rust/src/harray.rs +++ b/r-harmonium/src/rust/src/harray.rs @@ -487,7 +487,7 @@ impl HArray { /// harray1$is_shared() # FALSE. /// /// harray2 = harray1$clone() - /// harray$is_shared() # TRUE, HArray object shared with harray2. + /// harray1$is_shared() # TRUE, HArray object shared with harray2. /// ``` /// /// _________ diff --git a/r-harmonium/src/rust/src/haudiosink.rs b/r-harmonium/src/rust/src/haudiosink.rs index 5bbd98f..570d9f4 100644 --- a/r-harmonium/src/rust/src/haudiosink.rs +++ b/r-harmonium/src/rust/src/haudiosink.rs @@ -264,7 +264,7 @@ impl HAudioSink { /// /// Returns the position of the sound that’s being played. /// This takes into account any speedup or delay applied. - /// Example: if you apply a speedup of 2 to an mp3 decoder source and get_pos() returns 5s then the position in the mp3 recording is 10s from its start. + /// Example: if you apply a speedup of 2 to an mp3 decoder source and `get_pos()` returns 5s then the position in the mp3 recording is 10s from its start. /// /// #### Returns /// diff --git a/r-harmonium/src/rust/src/hfft.rs b/r-harmonium/src/rust/src/hfft.rs index 23d1f2f..5ee77f4 100644 --- a/r-harmonium/src/rust/src/hfft.rs +++ b/r-harmonium/src/rust/src/hfft.rs @@ -1,51 +1,50 @@ use crate::{conversions::ToScalar, errors::HErrorR, harray::HArray, hdatatype::HDataType}; +use harmonium_fft::fft::{Fft, RealFftForward, RealFftInverse}; use ndarray::IxDyn; use num_complex::Complex; -use realfft::RealFftPlanner; -use rustfft::FftPlanner; -use savvy::{r_println, savvy, Sexp}; +use savvy::{r_println, savvy, OwnedLogicalSexp, Sexp}; use std::sync::Arc; -pub trait HFftPlannerR: Send { - fn fft(&mut self, harray: &mut HArray) -> savvy::Result<()>; - fn ifft(&mut self, harray: &mut HArray) -> savvy::Result<()>; +pub trait HFftR: Send + Sync { + fn process(&mut self, harray: &mut HArray) -> savvy::Result<()>; fn dtype(&self) -> savvy::Result; - fn print(&self) -> savvy::Result<()>; + fn clone_inner(&self) -> Arc; } -pub trait HRealFftPlannerR: Send { - fn rfft(&mut self, harray: &mut HArray) -> savvy::Result<()>; - fn irfft(&mut self, harray: &mut HArray, length: usize) -> savvy::Result<()>; +pub trait HRealFftR: Send + Sync { + fn process(&mut self, harray: &mut HArray) -> savvy::Result<()>; fn dtype(&self) -> savvy::Result; - fn print(&self) -> savvy::Result<()>; + fn clone_inner(&self) -> Arc; } -/// HFftPlanner -/// A planner is used to create FFTs. It caches results internally, so when making more than one FFT it is advisable to reuse the same planner. +/// HFft +/// An `HFft` is used to create FFTs. It caches results internally, so when making more than one FFT it is advisable to reuse the same `HFft` instance. /// /// # Methods /// #[savvy] -pub struct HFftPlanner(pub Box); +#[derive(Clone)] +pub struct HFft(pub Arc); -/// HRealFftPlanner -/// A planner is used to create FFTs. It caches results internally, so when making more than one FFT it is advisable to reuse the same planner. -/// -/// This planner is used to calculate FFTs of real valued inputs and its inverse operation. -/// -/// # Methods -/// +///// HRealFft +///// An HRealFft is used to create real FFTs. It caches results internally, so when making more than one FFT it is advisable to reuse the same planner. +///// +///// This planner is used to calculate FFTs of real valued inputs and its inverse operation. +///// +///// # Methods +///// #[savvy] -pub struct HRealFftPlanner(pub Box); +#[derive(Clone)] +pub struct HRealFft(pub Arc); #[savvy] -impl HFftPlanner { - /// HFftPlanner - /// ## new +impl HFft { + /// HFft + /// ## new_fft_forward /// - /// `new(dtype: HDataType) -> HFftPlanner` + /// `new_fft_forward(length: integer, dtype: HDataType) -> HFft` /// - /// Creates a new `HFftPlanner` instance. + /// Creates a new `HFft` instance which will be used to calculate forward FFTs. /// /// If you plan on creating multiple FFT instances, it is recommended to reuse the same planner for all of them. This is because the planner re-uses internal data /// across FFT instances wherever possible, saving memory and reducing setup time (FFT instances created with one planner will never re-use data and buffers with @@ -56,13 +55,18 @@ impl HFftPlanner { /// /// #### Arguments /// + /// - `length` + /// + /// An integer denoting the length of the input. For 2D `HArray`'s, nrows must + /// be provided. + /// /// - `dtype` /// - /// A complex `HDataType` to indicate the dtype that the `HFftPlanner` will be working with. + /// A complex `HDataType` to indicate the dtype that the `HFft` will be working with. /// /// #### Returns /// - /// An `HFftPlanner`. + /// An `HFft`. /// /// Will return an error if dtype is of a float type. /// @@ -70,79 +74,90 @@ impl HFftPlanner { /// /// ```r /// library(harmonium) - /// arr = array(c(1,2,3,4,5,6,7,8,9,10,11,12), c(3,4)) - /// dtype = HDataType$Float32 + /// arr = array(c(1+1i,2+2i,3+3i,4+4i,5+5i,6+6i), c(3,2)) + /// dtype = HDataType$Complex32 /// harray = HArray$new_from_values(arr, dtype) - /// fft_planner = HFftPlanner$new(harray$dtype()) + /// hfft = HFft$new_fft_forward(3L, harray$dtype()) /// ``` /// /// _________ /// - fn new(dtype: &HDataType) -> savvy::Result { + fn new_fft_forward(length: Sexp, dtype: &HDataType) -> savvy::Result { + let length: i32 = length.to_scalar()?; + let length: usize = length + .try_into() + .map_err(|_| savvy::Error::new("Cannot convert i32 to usize."))?; match dtype { - HDataType::Float32 => Err("The HFftPlanner is for Complex dtypes.".into()), - HDataType::Float64 => Err("The HFftPlanner is for Complex dtypes.".into()), - HDataType::Complex32 => Ok(HFftPlanner(Box::new(FftPlanner::::new()))), - HDataType::Complex64 => Ok(HFftPlanner(Box::new(FftPlanner::::new()))), + HDataType::Float32 => Err("The HFft is for Complex dtypes.".into()), + HDataType::Float64 => Err("The HFft is for Complex dtypes.".into()), + HDataType::Complex32 => Ok(HFft(Arc::new(Fft::::new_fft_forward(length)))), + HDataType::Complex64 => Ok(HFft(Arc::new(Fft::::new_fft_forward(length)))), } } - /// HFftPlanner - /// ## fft + /// HFft + /// ## new_fft_inverse /// - /// `fft(harray: HArray)` + /// `new_fft_inverse(length: integer, dtype: HDataType) -> HFft` /// - /// Computes the fast fourier transform of a complex `HArray`. + /// Creates a new `HFft` instance which will be used to calculate inverse FFTs. /// - /// The operation is done in-place. + /// If you plan on creating multiple FFT instances, it is recommended to reuse the same planner for all of them. This is because the planner re-uses internal data + /// across FFT instances wherever possible, saving memory and reducing setup time (FFT instances created with one planner will never re-use data and buffers with + /// FFT instances created by a different planner). /// - /// FFT (Fast Fourier Transform) refers to a way the discrete Fourier Transform (DFT) can be calculated efficiently, - /// by using symmetries in the calculated terms. The symmetry is highest when n is a power of 2, and the transform - /// is therefore most efficient for these sizes. + /// In the constructor, the FftPlanner will detect available CPU features. If AVX, SSE, Neon, or WASM SIMD are available, it will set itself up to plan FFTs with + /// the fastest available instruction set. If no SIMD instruction sets are available, the planner will seamlessly fall back to planning non-SIMD FFTs. /// - /// The function does not normalize outputs. Callers must manually normalize the results by scaling each element by - /// `1/sqrt(n)`. Multiple normalization steps can be merged into one via pairwise multiplication, so when doing - /// a forward FFT followed by an inverse callers can normalize once by scaling each element by `1/n`. + /// #### Arguments /// - /// Elements in the output are ordered by ascending frequency, with the first element corresponding to frequency 0. + /// - `length` /// - /// #### Arguments + /// An integer denoting the length of the input. For 2D `HArray`'s, nrows must + /// be provided. /// - /// - `harray` + /// - `dtype` /// - /// A complex `HArray`. + /// A complex `HDataType` to indicate the dtype that the `HFft` will be working with. /// /// #### Returns /// - /// Will return an error if: - /// - /// - The `HArray`'s dtype is incompatible with the `HFftPlanner`'s dtype. + /// An `HFft`. /// - /// - The `HArray`'s `ndim` is greater than 2. + /// Will return an error if dtype is of a float type. /// /// #### Examples /// /// ```r /// library(harmonium) - /// arr = array(c(1,2,3,4,5,6,7,8,9,10,11,12), c(3,4)) - /// dtype = HDataType$Float32 + /// arr = array(c(1+1i,2+2i,3+3i,4+4i,5+5i,6+6i), c(3,2)) + /// dtype = HDataType$Complex32 /// harray = HArray$new_from_values(arr, dtype) - /// fft_planner = HFftPlanner$new(harray$dtype()) - /// fft_planner$fft(harray) + /// hfft = HFft$new_fft_inverse(3L, harray$dtype()) /// ``` /// /// _________ /// - fn fft(&mut self, harray: &mut HArray) -> savvy::Result<()> { - self.0.fft(harray) + fn new_fft_inverse(length: Sexp, dtype: &HDataType) -> savvy::Result { + let length: i32 = length.to_scalar()?; + let length: usize = length + .try_into() + .map_err(|_| savvy::Error::new("Cannot convert i32 to usize."))?; + match dtype { + HDataType::Float32 => Err("The HFft is for Complex dtypes.".into()), + HDataType::Float64 => Err("The HFft is for Complex dtypes.".into()), + HDataType::Complex32 => Ok(HFft(Arc::new(Fft::::new_fft_inverse(length)))), + HDataType::Complex64 => Ok(HFft(Arc::new(Fft::::new_fft_inverse(length)))), + } } - /// HFftPlanner - /// ## ifft + /// HFft + /// ## process /// - /// `ifft(harray: HArray)` + /// `process(harray: HArray)` /// - /// Computes the inverse fast fourier transform of a complex `HArray`. + /// Computes the fast fourier transform of a complex `HArray`. + /// The FFT computed may be forward or inverse, depending on the `HFFT` created. /// /// The operation is done in-place. /// @@ -154,6 +169,8 @@ impl HFftPlanner { /// `1/sqrt(n)`. Multiple normalization steps can be merged into one via pairwise multiplication, so when doing /// a forward FFT followed by an inverse callers can normalize once by scaling each element by `1/n`. /// + /// Elements in the output are ordered by ascending frequency, with the first element corresponding to frequency 0. + /// /// #### Arguments /// /// - `harray` @@ -164,33 +181,42 @@ impl HFftPlanner { /// /// Will return an error if: /// - /// - The `HArray`'s dtype is incompatible with the `HFftPlanner`'s dtype. + /// - The `HArray`'s dtype is incompatible with the `HFft`'s dtype. /// /// - The `HArray`'s `ndim` is greater than 2. /// /// #### Examples /// /// ```r + /// # Forward FFT. /// library(harmonium) /// arr = array(c(1+1i,2+2i,3+3i,4+4i,5+5i,6+6i), c(3,2)) /// dtype = HDataType$Complex32 /// harray = HArray$new_from_values(arr, dtype) - /// fft_planner = HFftPlanner$new(harray$dtype()) - /// fft_planner$ifft(harray) + /// hfft = HFft$new_fft_forward(3L, harray$dtype()) + /// hfft$process(harray) + /// + /// # Inverse FFT. + /// arr = array(c(1+1i,2+2i,3+3i,4+4i,5+5i,6+6i), c(3,2)) + /// dtype = HDataType$Complex32 + /// harray = HArray$new_from_values(arr, dtype) + /// hfft = HFft$new_fft_inverse(3L, harray$dtype()) + /// hfft$process(harray) /// ``` /// /// _________ /// - fn ifft(&mut self, harray: &mut HArray) -> savvy::Result<()> { - self.0.ifft(harray) + fn process(&mut self, harray: &mut HArray) -> savvy::Result<()> { + let inner_mut = self.get_inner_mut(); + inner_mut.process(harray) } - /// HFftPlanner + /// HFft /// ## dtype /// /// `dtype() -> HDataType` /// - /// Gets the `HFftPlanner`'s dtype. + /// Gets the `HFft`'s dtype. /// /// #### Returns /// @@ -203,8 +229,8 @@ impl HFftPlanner { /// arr = array(c(1+1i,2+2i,3+3i,4+4i,5+5i,6+6i), c(3,2)) /// dtype = HDataType$Complex32 /// harray = HArray$new_from_values(arr, dtype) - /// fft_planner = HFftPlanner$new(harray$dtype()) - /// fft_planner$dtype() + /// hfft = HFft$new_fft_forward(3L, harray$dtype()) + /// hfft$dtype() /// ``` /// /// _________ @@ -213,12 +239,12 @@ impl HFftPlanner { self.0.dtype() } - /// HFftPlanner + /// HFft /// ## print /// /// `print()` /// - /// Print the `HFftPlanner`. + /// Prints the `HFft`. /// /// Differently from R's normal behaviour, `print` doesn't return the value invisibly. /// @@ -229,28 +255,122 @@ impl HFftPlanner { /// arr = array(c(1+1i,2+2i,3+3i,4+4i,5+5i,6+6i), c(3,2)) /// dtype = HDataType$Complex32 /// harray = HArray$new_from_values(arr, dtype) - /// fft_planner = HFftPlanner$new(harray$dtype()) - /// fft_planner$print() + /// hfft = HFft$new_fft_forward(3L, harray$dtype()) + /// hfft$print() /// /// # or similarly: - /// print(fft_planner) + /// print(hfft) /// ``` /// /// _________ /// fn print(&self) -> savvy::Result<()> { - self.0.print() + r_println!("HFft"); + Ok(()) + } + + /// HFft + /// ## clone + /// + /// `clone() -> HFft` + /// + /// Clones the `HFft`. + /// + /// Creates a new `HFft`, with the underlying data pointing to the same place in memory. + /// When `HFFT` is cloned, thus having more than one reference to the same internal struct, and `process` is run, + /// it uses the same cached `Fft` instance, but a new scratch buffer will have to be allocated. + /// + /// #### Returns + /// + /// An `HFft`. + /// + /// #### Examples + /// + /// ```r + /// library(harmonium) + /// arr = array(c(1+1i,2+2i,3+3i,4+4i,5+5i,6+6i), c(3,2)) + /// dtype = HDataType$Complex32 + /// harray = HArray$new_from_values(arr, dtype) + /// hfft = HFft$new_fft_forward(3L, harray$dtype()) + /// hfft$clone() + /// ``` + /// + /// _________ + /// + fn clone(&self) -> savvy::Result { + Ok(std::clone::Clone::clone(self)) + } + + /// HFft + /// ## is_shared + /// + /// `is_shared() -> bool` + /// + /// Checks if the object is shared. + /// + /// Since `HFft` has a COW ([clone-on-write](https://doc.rust-lang.org/std/borrow/enum.Cow.html)) behaviour, this function is useful to check if a new + /// object will be created or if the change will be done in-place. + /// + /// #### Returns + /// + /// A bool. + /// + /// #### Examples + /// + /// ```r + /// library(harmonium) + /// arr = array(c(1+1i,2+2i,3+3i,4+4i,5+5i,6+6i), c(3,2)) + /// dtype = HDataType$Complex32 + /// harray = HArray$new_from_values(arr, dtype) + /// hfft = HFft$new_fft_forward(3L, harray$dtype()) + /// hfft$is_shared() # FALSE. + /// + /// hfft2 = hfft$clone() + /// hfft$is_shared() # TRUE, hfft shares the same inner object with hfft2. + /// ``` + /// + /// _________ + /// + fn is_shared(&self) -> savvy::Result { + let bool = Arc::weak_count(&self.0) + Arc::strong_count(&self.0) != 1; + let logical_sexp: OwnedLogicalSexp = bool.try_into()?; + logical_sexp.into() + } + + /// HFft + /// ## invalidate + /// + /// `invalidate()` + /// + /// Replaces the inner value of the external pointer, invalidating it. + /// This function is useful to remove one of the shared references of the inner pointer in rust. + /// + /// #### Examples + /// + /// ```r + /// library(harmonium) + /// arr = array(c(1+1i,2+2i,3+3i,4+4i,5+5i,6+6i), c(3,2)) + /// dtype = HDataType$Complex32 + /// harray = HArray$new_from_values(arr, dtype) + /// hfft = HFft$new_fft_forward(3L, harray$dtype()) + /// hfft$invalidate() + /// ``` + /// + /// _________ + /// + pub fn invalidate(self) -> savvy::Result<()> { + Ok(()) } } #[savvy] -impl HRealFftPlanner { - /// HRealFftPlanner - /// ## new +impl HRealFft { + /// HRealFft + /// ## new_real_fft /// - /// `new(dtype: HDataType) -> HRealFftPlanner` + /// `new_real_fft(length: integer, dtype: HDataType) -> HRealFft` /// - /// Creates a new `HRealFftPlanner` instance. + /// Creates a new `HRealFft` instance which will be used to calculate forward FFTs. /// /// If you plan on creating multiple FFT instances, it is recommended to reuse the same planner for all of them. This is because the planner re-uses internal data /// across FFT instances wherever possible, saving memory and reducing setup time (FFT instances created with one planner will never re-use data and buffers with @@ -261,15 +381,21 @@ impl HRealFftPlanner { /// /// #### Arguments /// + /// - `length` + /// + /// An integer denoting the length of the input for forward FFTs and the length of the output for inverse FFTs. For 2D `HArray`'s, nrows must + /// be provided. + /// /// - `dtype` /// - /// A float `HDataType` to indicate the dtype that the `HFftPlanner` will be working with. + /// An `HDataType` to indicate the dtype that the `HRealFft` will be working with. If float, + /// will calculate the forward FFT. If complex, will calculate the inverse FFT. /// /// #### Returns /// - /// An `HRealFftPlanner`. + /// An `HRealFft`. /// - /// Will return an error if dtype is of a complex type. + /// Will return an error if dtype is of a float type. /// /// #### Examples /// @@ -278,26 +404,41 @@ impl HRealFftPlanner { /// arr = array(c(1,2,3,4,5,6,7,8,9,10,11,12), c(3,4)) /// dtype = HDataType$Float32 /// harray = HArray$new_from_values(arr, dtype) - /// real_fft_planner = HRealFftPlanner$new(harray$dtype()) + /// hfft = HRealFft$new_real_fft(3L, harray$dtype()) /// ``` /// /// _________ /// - fn new(dtype: &HDataType) -> savvy::Result { + fn new_real_fft(length: Sexp, dtype: &HDataType) -> savvy::Result { + let length: i32 = length.to_scalar()?; + let length: usize = length + .try_into() + .map_err(|_| savvy::Error::new("Cannot convert i32 to usize."))?; match dtype { - HDataType::Float32 => Ok(HRealFftPlanner(Box::new(RealFftPlanner::::new()))), - HDataType::Float64 => Ok(HRealFftPlanner(Box::new(RealFftPlanner::::new()))), - HDataType::Complex32 => panic!(), - HDataType::Complex64 => panic!(), + HDataType::Float32 => Ok(HRealFft(Arc::new( + RealFftForward::::new_real_fft_forward(length), + ))), + HDataType::Float64 => Ok(HRealFft(Arc::new( + RealFftForward::::new_real_fft_forward(length), + ))), + HDataType::Complex32 => Ok(HRealFft(Arc::new( + RealFftInverse::::new_real_fft_inverse(length), + ))), + HDataType::Complex64 => Ok(HRealFft(Arc::new( + RealFftInverse::::new_real_fft_inverse(length), + ))), } } - /// HRealFftPlanner - /// ## rfft + /// HRealFft + /// ## process /// - /// `rfft(harray: HArray)` + /// `process(harray: HArray)` /// - /// Computes the fast fourier transform of a float `HArray`. Transforms a real signal of length `N` to a complex-valued spectrum of length `N/2+1` (with `N/2` rounded down). + /// Computes the fast fourier transform of a float `HArray` or the inverse fast fourier transform of a complex `HArray`. + /// For a forward FFT, transforms a real signal of length `N` to a complex-valued spectrum of length `N/2+1` (with `N/2` rounded down). + /// For an inverse FFT, transforms a complex spectrum of length `N/2+1` (with `N/2` rounded down) to a real-valued + /// signal of length `N`. /// /// The operation is not done in-place, although the same external pointer is used to store the new HArray. /// @@ -318,13 +459,13 @@ impl HRealFftPlanner { /// /// - `harray` /// - /// A float `HArray`. + /// An `HArray`. /// /// #### Returns /// /// Will return an error if: /// - /// - The `HArray`'s dtype is incompatible with the `HFftPlanner`'s dtype. + /// - The `HArray`'s dtype is incompatible with the `HRealFft`'s dtype. /// /// - The `HArray`'s `ndim` is greater than 2. /// @@ -335,194 +476,250 @@ impl HRealFftPlanner { /// arr = array(c(1,2,3,4,5,6,7,8,9,10,11,12), c(3,4)) /// dtype = HDataType$Float32 /// harray = HArray$new_from_values(arr, dtype) - /// real_fft_planner = HRealFftPlanner$new(harray$dtype()) - /// real_fft_planner$rfft(harray) + /// # Forward fft + /// fft = HRealFft$new_real_fft(3L, harray$dtype()) + /// fft$process(harray) + /// # Inverse fft + /// ifft = HRealFft$new_real_fft(3L, HDataType$Complex32) + /// ifft$process(harray) /// ``` /// /// _________ /// - fn rfft(&mut self, harray: &mut HArray) -> savvy::Result<()> { - self.0.rfft(harray) + fn process(&mut self, harray: &mut HArray) -> savvy::Result<()> { + let inner_mut = self.get_inner_mut(); + inner_mut.process(harray) } - /// HRealFftPlanner - /// ## irfft + /// HRealFft + /// ## dtype /// - /// `irfft(harray: HArray, length: integer)` + /// `dtype() -> HDataType` /// - /// Computes the inverse fast fourier transform of a complex `HArray`. Transforms a complex spectrum of length `N/2+1` (with `N/2` rounded down) to a real-valued - /// signal of length `N`. + /// Gets the `HRealFft`'s dtype. /// - /// The operation is not done in-place, although the same external pointer is used to store the new HArray. - /// The FFT of a real signal is Hermitian-symmetric, X[i] = conj(X[-i]) so the output contains only the positive frequencies - /// below the Nyquist frequency. + /// #### Returns /// - /// FFT (Fast Fourier Transform) refers to a way the discrete Fourier Transform (DFT) can be calculated efficiently, - /// by using symmetries in the calculated terms. The symmetry is highest when n is a power of 2, and the transform - /// is therefore most efficient for these sizes. + /// An `HDataType`. /// - /// The function does not normalize outputs. Callers must manually normalize the results by scaling each element by - /// `1/sqrt(n)`. Multiple normalization steps can be merged into one via pairwise multiplication, so when doing - /// a forward FFT followed by an inverse callers can normalize once by scaling each element by `1/n`. + /// #### Examples /// - /// #### Arguments + /// ```r + /// library(harmonium) + /// arr = array(c(1+1i,2+2i,3+3i,4+4i,5+5i,6+6i), c(3,2)) + /// dtype = HDataType$Complex32 + /// harray = HArray$new_from_values(arr, dtype) + /// hfft = HRealFft$new_real_fft(3L, HDataType$Complex32) + /// hfft$dtype() + /// ``` /// - /// - `harray` + /// _________ /// - /// A complex `HArray`. The `HArray`'s dtype must be the complex equivalent of the `HRealFftPlanner`'s dtype. For example if `HRealFftPlanner` dtype is `Float64`, - /// the `HArray`'s dtype must be `Complex64`. + fn dtype(&self) -> savvy::Result { + self.0.dtype() + } + + /// HRealFft + /// ## print /// - /// - `length` + /// `print()` /// - /// An integer. The output length of the signal. Since the spectrum is `N/2+1`, the length can be `N` and `N+1`, if `N` is even, or can be `N` and `N-1` if `N` is odd. + /// Prints the `HRealFft`. /// - /// #### Returns + /// Differently from R's normal behaviour, `print` doesn't return the value invisibly. /// - /// Will return an error if: + /// #### Examples /// - /// - The `HArray`'s dtype is incompatible with the `HFftPlanner`'s dtype. + /// ```r + /// library(harmonium) + /// arr = array(c(1+1i,2+2i,3+3i,4+4i,5+5i,6+6i), c(3,2)) + /// dtype = HDataType$Complex32 + /// harray = HArray$new_from_values(arr, dtype) + /// hfft = HRealFft$new_real_fft(3L, HDataType$Complex32) + /// hfft$print() /// - /// - The `HArray`'s `ndim` is greater than 2. + /// # or similarly: + /// print(hfft) + /// ``` + /// + /// _________ + /// + fn print(&self) -> savvy::Result<()> { + r_println!("HRealFft"); + Ok(()) + } + + /// HRealFft + /// ## clone + /// + /// `clone() -> HRealFft` /// - /// - The `length` argument is not compatible with the spectrum length. + /// Clones the `HRealFft`. + /// + /// Creates a new `HRealFft`, with the underlying data pointing to the same place in memory. + /// When `HFFT` is cloned, having more than one reference to the same internal struct, and `process` is run, it uses the same cached `Fft` instance, but a new + /// scratch buffer will have to be allocated. + /// + /// #### Returns + /// + /// An `HRealFft`. /// /// #### Examples /// /// ```r /// library(harmonium) - /// r = as.double(sample(100, 4, replace = TRUE)) - /// i = as.double(sample(100, 3, replace = TRUE)) - /// arr = array(complex(real=r, imaginary=c(0,i)), c(4,1)) + /// arr = array(c(1+1i,2+2i,3+3i,4+4i,5+5i,6+6i), c(3,2)) /// dtype = HDataType$Complex32 /// harray = HArray$new_from_values(arr, dtype) - /// real_fft_planner = HRealFftPlanner$new(HDataType$Float32) - /// real_fft_planner$irfft(harray, 7L) + /// hfft = HRealFft$new_real_fft(3L, HDataType$Complex32) + /// hfft$clone() /// ``` /// /// _________ /// - fn irfft(&mut self, harray: &mut HArray, length: Sexp) -> savvy::Result<()> { - let length: i32 = length.to_scalar()?; - let length: usize = length - .try_into() - .map_err(|_| savvy::Error::new("Cannot convert i32 to usize."))?; - self.0.irfft(harray, length) + fn clone(&self) -> savvy::Result { + Ok(std::clone::Clone::clone(self)) } - /// HRealFftPlanner - /// ## dtype + /// HRealFft + /// ## is_shared /// - /// `dtype() -> HDataType` + /// `is_shared() -> bool` /// - /// Gets the `HRealFftPlanner`'s dtype. + /// Checks if the object is shared. + /// + /// Since `HRealFft` has a COW ([clone-on-write](https://doc.rust-lang.org/std/borrow/enum.Cow.html)) behaviour, this function is useful to check if a new + /// object will be created or if the change will be done in-place. /// /// #### Returns /// - /// An `HDataType`. + /// A bool. /// /// #### Examples /// /// ```r /// library(harmonium) - /// arr = array(c(1,2,3,4,5,6,7,8,9,10,11,12), c(3,4)) - /// dtype = HDataType$Float32 + /// arr = array(c(1+1i,2+2i,3+3i,4+4i,5+5i,6+6i), c(3,2)) + /// dtype = HDataType$Complex32 /// harray = HArray$new_from_values(arr, dtype) - /// real_fft_planner = HRealFftPlanner$new(harray$dtype()) - /// real_fft_planner$dtype() + /// hfft = HRealFft$new_real_fft(3L, HDataType$Complex32) + /// hfft$is_shared() # FALSE. + /// + /// hfft2 = hfft$clone() + /// hfft$is_shared() # TRUE, HRealFft object shared with hfft2. /// ``` /// /// _________ /// - fn dtype(&self) -> savvy::Result { - self.0.dtype() + fn is_shared(&self) -> savvy::Result { + let bool = Arc::weak_count(&self.0) + Arc::strong_count(&self.0) != 1; + let logical_sexp: OwnedLogicalSexp = bool.try_into()?; + logical_sexp.into() } - /// HRealFftPlanner - /// ## print - /// - /// `print()` + /// HRealFft + /// ## invalidate /// - /// Print the `HRealFftPlanner`. + /// `invalidate()` /// - /// Differently from R's normal behaviour, `print` doesn't return the value invisibly. + /// Replaces the inner value of the external pointer, invalidating it. + /// This function is useful to remove one of the shared references of the inner pointer in rust. /// /// #### Examples /// /// ```r /// library(harmonium) - /// arr = array(c(1,2,3,4,5,6,7,8,9,10,11,12), c(3,4)) - /// dtype = HDataType$Float32 + /// arr = array(c(1+1i,2+2i,3+3i,4+4i,5+5i,6+6i), c(3,2)) + /// dtype = HDataType$Complex32 /// harray = HArray$new_from_values(arr, dtype) - /// real_fft_planner = HRealFftPlanner$new(harray$dtype()) - /// real_fft_planner$print() - /// - /// # or similarly: - /// print(real_fft_planner) + /// hfft = HRealFft$new_real_fft(3L, HDataType$Complex32) + /// hfft$invalidate() /// ``` /// /// _________ /// - fn print(&self) -> savvy::Result<()> { - self.0.print() + pub fn invalidate(self) -> savvy::Result<()> { + Ok(()) } } -macro_rules! impl_hfftplannerr { - ($(($t1:ty, $t2:ty, $e1:expr, $e2: expr)),+) => { +macro_rules! impl_hfft { + ($(($t1:ty, $t2:ty, $e1:expr)),+) => { $( - impl HFftPlannerR for $t1 { - fn fft(&mut self, harray: &mut HArray) -> savvy::Result<()> { - let harray = harray.get_inner_mut().as_any_mut().downcast_mut::<$t2>().ok_or_else(|| savvy::Error::new("HArray and HFftPlanner must have the same HDataType."))?; - harmonium_fft::fft::ProcessFft::fft(self, harray).map_err(|err| savvy::Error::from(HErrorR::from(err))) - } - - fn ifft(&mut self, harray: &mut HArray) -> savvy::Result<()> { - let harray = harray.get_inner_mut().as_any_mut().downcast_mut::<$t2>().ok_or_else(|| savvy::Error::new("HArray and HFftPlanner must have the same HDataType."))?; - harmonium_fft::fft::ProcessFft::ifft(self, harray).map_err(|err| savvy::Error::from(HErrorR::from(err))) + impl HFftR for $t1 { + fn process(&mut self, harray: &mut HArray) -> savvy::Result<()> { + let harray_inner = harray.get_inner_mut().as_any_mut().downcast_mut::<$t2>().ok_or_else(|| savvy::Error::new("HArray and HFft must have the same HDataType."))?; + harmonium_fft::fft::ProcessFft::process(self, harray_inner).map_err(|err| savvy::Error::from(HErrorR::from(err))) } fn dtype(&self) -> savvy::Result { Ok($e1) } - fn print(&self) -> savvy::Result<()> { - r_println!($e2); - Ok(()) + fn clone_inner(&self) -> Arc { + Arc::new(self.clone()) } } )+ }; } -impl_hfftplannerr!( +impl_hfft!( ( - FftPlanner, + Fft, harmonium_core::array::HArray, IxDyn>, - HDataType::Complex32, - "FftPlanner" + HDataType::Complex32 ), ( - FftPlanner, + Fft, harmonium_core::array::HArray, IxDyn>, - HDataType::Complex64, - "FftPlanner" + HDataType::Complex64 ) ); -macro_rules! impl_hrealfftplannerr { - ($(($t1:ty, $t2:ty, $t3:ty, $e1:expr, $e2: expr)),+) => { +macro_rules! impl_hrealfftforward { + ($(($t1:ty, $t2:ty, $e1:expr)),+) => { $( - impl HRealFftPlannerR for $t1 { - fn rfft(&mut self, harray: &mut HArray) -> savvy::Result<()> { - let harray_downcasted = harray.get_inner_mut().as_any_mut().downcast_mut::<$t2>().ok_or_else(|| savvy::Error::new("HArray and HFftPlanner must have the same HDataType."))?; - let result = harmonium_fft::fft::ProcessRealFft::rfft(self, harray_downcasted).map_err(|err| savvy::Error::from(HErrorR::from(err)))?; + impl HRealFftR for $t1 { + fn process(&mut self, harray: &mut HArray) -> savvy::Result<()> { + let harray_inner = harray.get_inner_mut().as_any_mut().downcast_mut::<$t2>().ok_or_else(|| savvy::Error::new("HArray and HFft must have the same HDataType."))?; + let result = harmonium_fft::fft::ProcessRealFftForward::process(self, harray_inner).map_err(|err| savvy::Error::from(HErrorR::from(err)))?; *harray = HArray(Arc::new(result)); Ok(()) } - fn irfft(&mut self, harray: &mut HArray, length: usize) -> savvy::Result<()> { - let harray_downcasted = harray.get_inner_mut().as_any_mut().downcast_mut::<$t3>().ok_or_else(|| savvy::Error::new("HArray must be the complex equivalent of HFftPlanner's HDataType."))?; - let result = harmonium_fft::fft::ProcessRealFft::irfft(self, harray_downcasted, length).map_err(|err| savvy::Error::from(HErrorR::from(err)))?; + fn dtype(&self) -> savvy::Result { + Ok($e1) + } + + fn clone_inner(&self) -> Arc { + Arc::new(self.clone()) + } + } + )+ + }; +} + +impl_hrealfftforward!( + ( + RealFftForward, + harmonium_core::array::HArray, + HDataType::Float32 + ), + ( + RealFftForward, + harmonium_core::array::HArray, + HDataType::Float64 + ) +); + +macro_rules! impl_hrealfftinverse { + ($(($t1:ty, $t2:ty, $e1:expr)),+) => { + $( + impl HRealFftR for $t1 { + fn process(&mut self, harray: &mut HArray) -> savvy::Result<()> { + let harray_inner = harray.get_inner_mut().as_any_mut().downcast_mut::<$t2>().ok_or_else(|| savvy::Error::new("HArray and HFft must have the same HDataType."))?; + let result = harmonium_fft::fft::ProcessRealFftInverse::process(self, harray_inner).map_err(|err| savvy::Error::from(HErrorR::from(err)))?; *harray = HArray(Arc::new(result)); Ok(()) } @@ -531,28 +728,45 @@ macro_rules! impl_hrealfftplannerr { Ok($e1) } - fn print(&self) -> savvy::Result<()> { - r_println!($e2); - Ok(()) + fn clone_inner(&self) -> Arc { + Arc::new(self.clone()) } } )+ }; } -impl_hrealfftplannerr!( +impl_hrealfftinverse!( ( - RealFftPlanner, - harmonium_core::array::HArray, + RealFftInverse, harmonium_core::array::HArray, IxDyn>, - HDataType::Float32, - "RealFftPlanner" + HDataType::Complex32 ), ( - RealFftPlanner, - harmonium_core::array::HArray, + RealFftInverse, harmonium_core::array::HArray, IxDyn>, - HDataType::Float64, - "RealFftPlanner" + HDataType::Complex64 ) ); + +impl HFft { + #[doc(hidden)] + pub fn get_inner_mut(&mut self) -> &mut dyn HFftR { + if Arc::get_mut(&mut self.0).is_none() { + self.0 = self.0.clone_inner(); + } + // Safe to unwrap_unchecked since get_mut was checked above. + unsafe { Arc::get_mut(&mut self.0).unwrap_unchecked() } + } +} + +impl HRealFft { + #[doc(hidden)] + pub fn get_inner_mut(&mut self) -> &mut dyn HRealFftR { + if Arc::get_mut(&mut self.0).is_none() { + self.0 = self.0.clone_inner(); + } + // Safe to unwrap_unchecked since get_mut was checked above. + unsafe { Arc::get_mut(&mut self.0).unwrap_unchecked() } + } +} diff --git a/r-harmonium/tests/testthat/test_hfftplanner.R b/r-harmonium/tests/testthat/test_hfft.R similarity index 82% rename from r-harmonium/tests/testthat/test_hfftplanner.R rename to r-harmonium/tests/testthat/test_hfft.R index 3f92ff0..43828cb 100644 --- a/r-harmonium/tests/testthat/test_hfftplanner.R +++ b/r-harmonium/tests/testthat/test_hfft.R @@ -8,8 +8,8 @@ test_that( dtype = HDataType$Complex64 harray = HArray$new_from_values(x, dtype) - fft_planner = HFftPlanner$new(dtype) - fft_planner$fft(harray) + hfft = HFft$new_fft_forward(30L, dtype) + hfft$process(harray) expect_equal(stats::fft(x), harray$collect(), tolerance = 1.0e-06) } @@ -20,8 +20,8 @@ test_that( dtype = HDataType$Complex64 harray = HArray$new_from_values(x, dtype) - fft_planner = HFftPlanner$new(dtype) - fft_planner$ifft(harray) + hfft = HFft$new_fft_inverse(30L, dtype) + hfft$process(harray) expect_equal(stats::fft(x, inverse = TRUE), harray$collect(), tolerance = 1.0e-06) } @@ -33,8 +33,8 @@ test_that( dtype = HDataType$Complex64 harray = HArray$new_from_values(x, dtype) - fft_planner = HFftPlanner$new(dtype) - fft_planner$fft(harray) + hfft = HFft$new_fft_forward(as.integer(30/10), dtype) + hfft$process(harray) expect_equal(stats::mvfft(x), harray$collect(), tolerance = 1.0e-06) } @@ -46,8 +46,8 @@ test_that( dtype = HDataType$Complex64 harray = HArray$new_from_values(x, dtype) - fft_planner = HFftPlanner$new(dtype) - fft_planner$ifft(harray) + hfft = HFft$new_fft_inverse(as.integer(30/10), dtype) + hfft$process(harray) expect_equal(stats::mvfft(x, inverse = TRUE), harray$collect(), tolerance = 1.0e-06) } @@ -55,8 +55,8 @@ test_that( x = as.array(as.double(1:6)) dtype = HDataType$Float64 harray = HArray$new_from_values(x, dtype) - real_fft_planner = HRealFftPlanner$new(dtype) - real_fft_planner$rfft(harray) + hfft = HRealFft$new_real_fft(as.integer(6), dtype) + hfft$process(harray) result = as.array(c(21+0i, -3+5.196152i, -3+1.732051i, -3+0i)) expect_equal(harray$collect(), result, tolerance = 1.0e-06) } @@ -66,8 +66,8 @@ test_that( x = matrix(x, ncol = 3) dtype = HDataType$Float64 harray = HArray$new_from_values(x, dtype) - real_fft_planner = HRealFftPlanner$new(dtype) - real_fft_planner$rfft(harray) + hfft = HRealFft$new_real_fft(as.integer(4), dtype) + hfft$process(harray) result = matrix(c(10+0i, -2+2i, -2+0i, 26+0i, -2+2i, -2+0i, 42+0i, -2+2i, -2+0i), ncol = 3) expect_equal(harray$collect(), result, tolerance = 1.0e-06) }