diff --git a/src/complex.rs b/src/complex.rs new file mode 100644 index 0000000..1e07b21 --- /dev/null +++ b/src/complex.rs @@ -0,0 +1,42 @@ +use std::ops::{Add, Mul, Sub}; + +/// A complex number. +#[allow(non_camel_case_types)] +#[derive(Clone, Copy, Debug)] +pub struct c64(pub f64, pub f64); + +impl Add for c64 { + type Output = Self; + + #[inline(always)] + fn add(self, rhs: c64) -> c64 { + c64(self.0 + rhs.0, self.1 + rhs.1) + } +} + +impl Mul for c64 { + type Output = Self; + + #[inline(always)] + fn mul(self, rhs: c64) -> c64 { + c64(self.0 * rhs.0 - self.1 * rhs.1, self.0 * rhs.1 + self.1 * rhs.0) + } +} + +impl Mul for c64 { + type Output = Self; + + #[inline(always)] + fn mul(self, rhs: f64) -> c64 { + c64(self.0 * rhs, self.1 * rhs) + } +} + +impl Sub for c64 { + type Output = Self; + + #[inline(always)] + fn sub(self, rhs: c64) -> c64 { + c64(self.0 - rhs.0, self.1 - rhs.1) + } +} diff --git a/src/lib.rs b/src/lib.rs index 6334796..7f57964 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -1,96 +1,139 @@ -#![feature(step_by)] +use std::slice; + +mod complex; +pub use complex::c64; + +/// A means of obtaining a slice of mutable complex numbers. +pub trait AsMutComplex<'l> { + fn as_mut_complex(self) -> &'l mut [c64]; +} + +impl<'l> AsMutComplex<'l> for &'l mut [c64] { + #[inline(always)] + fn as_mut_complex(self) -> &'l mut [c64] { + self + } +} + +impl<'l> AsMutComplex<'l> for &'l mut [f64] { + /// Treat the slice as a collection of pairs of real and imaginary parts and + /// reinterpret it as a slice of complex numbers. + /// + /// ## Panics + /// + /// The function panics if the number of elements is not even. + #[inline] + fn as_mut_complex(self) -> &'l mut [c64] { + unsafe { + let length = self.len(); + assert!(length % 2 == 0, "the number of elements should be even"); + slice::from_raw_parts_mut(self.as_mut_ptr() as *mut _, length / 2) + } + } +} + +impl<'l> AsMutComplex<'l> for &'l mut Vec { + #[inline] + fn as_mut_complex(self) -> &'l mut [c64] { + <&mut [f64]>::as_mut_complex(&mut *self) + } +} /// Perform the Fourier transform. /// -/// `data` should contain `n` complex numbers where `n` is a power of two. Each -/// complex number is stored as a pair of `f64`s so that the first is the real -/// part and the second is the corresponding imaginary part. Hence, the total -/// length of `data` should be `2 × n`. -/// -/// # Panics -/// -/// The function panics if `data.len()` is not even or `data.len() / 2` is not a -/// power of two. +/// The number of points should be a power of two. #[inline(always)] -pub fn forward(data: &mut [f64]) { - transform(data, 1.0); +pub fn forward<'l, T: AsMutComplex<'l>>(data: T) { + let data = data.as_mut_complex(); + + let n = data.len(); + if n < 1 || n & (n - 1) != 0 { + panic!("expected the number of points to be a power of two"); + } + + rearrange(data, n); + perform(data, n, false); } /// Perform the inverse Fourier transform. /// -/// `data` should contain `n` complex numbers where `n` is a power of two. Each -/// complex number is stored as a pair of `f64`s so that the first is the real -/// part and the second is the corresponding imaginary part. Hence, the total -/// length of `data` should be `2 × n`. -/// -/// # Panics -/// -/// The function panics if `data.len()` is not even or `data.len() / 2` is not a -/// power of two. +/// The number of points should be a power of two. #[inline(always)] -pub fn inverse(data: &mut [f64]) { - transform(data, -1.0); -} - -fn transform(data: &mut [f64], isign: f64) { - use std::f64::consts::PI; - - let l = data.len(); - if l % 2 != 0 { - panic!("expected the length of the data to be even"); - } +pub fn inverse<'l, T: AsMutComplex<'l>>(data: T, scaling: bool) { + let data = data.as_mut_complex(); - let n = l / 2; + let n = data.len(); if n < 1 || n & (n - 1) != 0 { panic!("expected the number of points to be a power of two"); } - let nn = n << 1; + rearrange(data, n); + perform(data, n, true); + if scaling { + scale(data, n); + } +} - let mut j = 1; - for i in (1..nn).step_by(2) { - if j > i { - data.swap(j - 1, i - 1); - data.swap(j, i); +#[inline(always)] +fn rearrange(data: &mut [c64], n: usize) { + let mut target = 0; + for position in 0..n { + if target > position { + data.swap(position, target); } - let mut m = n; - while m >= 2 && j > m { - j -= m; - m >>= 1; + let mut mask = n >> 1; + while target & mask != 0 { + target &= !mask; + mask >>= 1; } - j += m; + target |= mask; } +} - let mut mmax = 2; - while nn > mmax { - let istep = mmax << 1; - let theta = isign * (2.0 * PI / mmax as f64); - let wtemp = (0.5 * theta).sin(); - let wpr = -2.0 * wtemp * wtemp; - let wpi = theta.sin(); - let mut wr = 1.0; - let mut wi = 0.0; - for m in (1..mmax).step_by(2) { - for i in (m..(nn + 1)).step_by(istep) { - let j = i + mmax; - let tempr = wr * data[j - 1] - wi * data[j]; - let tempi = wr * data[j] + wi * data[j - 1]; - data[j - 1] = data[i - 1] - tempr; - data[j] = data[i] - tempi; - data[i - 1] += tempr; - data[i] += tempi; +#[inline(always)] +fn perform(data: &mut [c64], n: usize, inverse: bool) { + use std::f64::consts::PI; + + let pi = if inverse { PI } else { -PI }; + + let mut step = 1; + while step < n { + let jump = step << 1; + let delta = pi / step as f64; + let sine = (0.5 * delta).sin(); + let multiplier = c64(-2.0 * sine * sine, delta.sin()); + let mut factor = c64(1.0, 0.0); + for group in 0..step { + let mut pair = group; + while pair < n { + let match_pair = pair + step; + let product = factor * data[match_pair]; + data[match_pair] = data[pair] - product; + data[pair] = data[pair] + product; + pair += jump; } - let wtemp = wr; - wr = wr * wpr - wi * wpi + wr; - wi = wi * wpr + wtemp * wpi + wi; + factor = multiplier * factor + factor; } - mmax = istep; + step <<= 1; } +} - if isign == -1.0 { - let scale = 1.0 / n as f64; - for i in 0..l { - data[i] *= scale; - } +#[inline(always)] +fn scale(data: &mut [c64], n: usize) { + let factor = 1.0 / n as f64; + for position in 0..n { + data[position] = data[position] * factor; + } +} + +#[cfg(test)] +mod tests { + use c64; + use std::mem; + + #[test] + fn size_of() { + assert_eq!(mem::size_of::(), 2 * mem::size_of::()); + assert_eq!(mem::size_of::<[c64; 42]>(), 2 * mem::size_of::<[f64; 42]>()); } } diff --git a/tests/fixtures.rs b/tests/fixtures.rs index 43a464a..cc8793a 100644 --- a/tests/fixtures.rs +++ b/tests/fixtures.rs @@ -130,132 +130,132 @@ pub const TIME_DATA: [f64; 256] = [ ]; pub const FREQUENCY_DATA: [f64; 256] = [ - 67.506562667203411, 62.798747687506840, - -2.916435892389177, 2.258137395096745, - -1.213411748796466, -2.824305441697990, - -3.884636396852776, -6.613971933163763, - 5.049316157478026, 0.456883084114653, - -0.466361441955149, -5.481824322589378, - 0.803910676769211, 0.147054274059631, - -6.749332581521054, 4.840407656923564, - -0.262756329494097, 3.254410744830774, - -6.973431206318843, 1.575164194822261, - -5.636666719031166, -3.624154259337717, - -1.286794839414061, -3.596629452384591, - -1.949442402686498, -1.250892366300883, - -2.956509867164284, -3.563268353779529, - -1.629751417352084, 4.497170521285449, - 1.320095075374985, -2.156073489574812, - -2.002497059929997, -3.607915605047517, - 2.554697843376545, -2.629716248632329, - 2.398999344465010, -0.074112453018754, - 1.495969367959889, 1.905315704750475, - -2.236600417683027, 0.738926830782050, - -0.048604481554162, -2.372949299902079, - 3.012458946272938, -6.457540983540477, - -0.751768100557872, 2.124100304630884, - -1.792775263985138, 0.264154836959835, - -0.307176194080528, -5.087687483603743, - -1.005151880556302, 2.902422614798144, - 0.724987856799877, -2.068196780758698, - -0.342018034885117, 2.522012423589103, - -0.434808099354017, 2.147158479497558, - -2.674474141084141, -2.563372692380105, - -6.632075418392929, -1.685124294109070, - -2.699185197567189, -5.897826029871508, - -0.241335159541328, 3.101405290332420, - 2.979104548159058, -1.160028660033103, - 4.840508528923102, -3.204368559739159, - -3.399157829092539, 2.161358937350348, - 4.878539287722978, -1.459376706664861, - -0.426496319168839, 3.030254700260016, - 5.146915245537944, 0.282061559224328, - -0.008414044119338, -3.991754204800235, - 3.638877999600389, 1.258926814764483, - 1.909130437842078, 0.995539872732337, - 6.414605960628886, -5.679739832207080, - 3.905872946117078, -2.826501461749101, - 0.967386492473095, -5.576552424605286, - -2.853890281653033, -3.715099310600197, - -0.998730499887771, -1.368332585830762, - 0.539665253171003, 0.702327760413163, - -2.202133267022090, 1.329217526808806, - 0.923776499928182, 2.150735454872208, - 2.473335187775022, -0.373719659978438, - -2.253868990511847, -2.913770876062729, - -0.368246167526135, -1.498059975229275, - -0.075789241425793, 1.719766675858378, - -6.249259838595410, 6.001167756824215, - -4.027346047616200, 5.509694798065765, - -3.923350081315725, -0.896137928694981, - -0.923595644189471, -1.735427210626094, - -3.429327683732347, 0.363408609675781, - -2.276214983007463, -0.860532742279556, - 2.965000493451818, 5.042549679265357, - 2.949752700932401, 2.757985604650260, - -7.551126038911468, -1.253320876553430, - -3.991693675979835, -2.395612192990079, - -3.231893261005089, -3.499165263231712, - -1.385547239860866, 4.278296734710755, - -2.521465228517418, -1.483012583302978, - 2.166953458666462, 3.698992188846230, - -1.718586561803937, -2.848378813769423, - 2.356251252889540, 3.100152502890587, - -1.680667450021412, -5.672119730229592, - 3.266763808523838, -4.939142356746459, - -0.677140885336755, 7.689403378304775, - 0.714397029389201, -1.027012942409025, - 4.548911919541158, -3.110338786263284, - 0.963824257168846, -1.315676780411989, - -0.630783224843193, -1.040038800556799, - 0.051441799174254, 8.187374675121179, - 2.897552781182147, 4.424347873913396, - 3.813912169625594, 1.672124644400444, - -1.849780677085411, -6.075332676660434, - -3.451402723746721, -0.201945959802337, - -2.768891360840025, 1.367687274872978, - -1.559401743047184, -2.563598701770972, - -1.874383165372955, -1.816379942063397, - 4.841738890983525, 4.304535776881384, - -1.617862103749635, 4.505029568623651, - 3.144799926493496, 0.532041249412581, - -3.772195284223951, 2.554683615488583, - 1.466063738740647, 0.328554468964208, - 2.211012545853337, -3.584550392170861, - -5.848019862181538, -1.853688693742572, - -1.310659972997355, -4.744849190522960, - 2.553617632271563, -2.515247887047404, - 4.482138368790012, -8.618372361868635, - 0.884305207470568, 1.209832476084043, - 0.983852564385845, -1.970780039723666, - -4.755886984684762, -2.371006165743844, - -2.386952673947499, -0.380046032794893, - 8.202058412989162, -1.912136722291743, - 1.303130356713501, 0.110408680312637, - 4.396893371992936, -1.781029123711262, - -5.653027779602041, -4.457764259327377, - 2.344622479665089, 0.978245756031540, - 0.085165407317421, 1.846582481253430, - 2.305750958010686, 0.037116571959843, - 2.090636515202623, -2.544058986903603, - -1.098781173708802, -4.976387399730235, - 3.116401263119525, 0.806086033998370, - -7.921745368568139, -7.472969089954924, - 1.411453345272358, -1.449094611132407, - -0.804853562658514, 4.667294297528407, - 4.288737227223235, -3.968319162458097, - -2.972572998085656, 5.831133886072563, - 0.840509623087476, 4.498605610210712, - 2.983076467939862, -3.013246639835036, - 3.880468387697881, 1.634838899334426, - -4.944656692463916, 5.315462183419951, - 4.133604334208203, -3.445364011535503, - -3.599237215989622, -1.105183198676706, - -3.458371226842255, 8.013344593032475, - -5.278077778941197, -4.066322943629892, - 0.772765865031132, -8.707366039824819, - -0.108445528881629, -0.965343418118287, - 4.308501796277244, -0.278006225941031, - -5.303978155417413, -0.288754984960703, - 0.077148388378536, 5.819389938527753, + 6.750656266720340e+01, 6.279874768750683e+01, + 7.714838837855442e-02, 5.819389938527752e+00, + -5.303978155417409e+00, -2.887549849606836e-01, + 4.308501796277235e+00, -2.780062259410463e-01, + -1.084455288816293e-01, -9.653434181182980e-01, + 7.727658650311002e-01, -8.707366039824819e+00, + -5.278077778941212e+00, -4.066322943629877e+00, + -3.458371226842232e+00, 8.013344593032468e+00, + -3.599237215989637e+00, -1.105183198676696e+00, + 4.133604334208186e+00, -3.445364011535524e+00, + -4.944656692463901e+00, 5.315462183419959e+00, + 3.880468387697888e+00, 1.634838899334416e+00, + 2.983076467939860e+00, -3.013246639835038e+00, + 8.405096230874929e-01, 4.498605610210712e+00, + -2.972572998085648e+00, 5.831133886072569e+00, + 4.288737227223233e+00, -3.968319162458111e+00, + -8.048535626585105e-01, 4.667294297528397e+00, + 1.411453345272351e+00, -1.449094611132408e+00, + -7.921745368568175e+00, -7.472969089954909e+00, + 3.116401263119543e+00, 8.060860339983653e-01, + -1.098781173708808e+00, -4.976387399730230e+00, + 2.090636515202629e+00, -2.544058986903604e+00, + 2.305750958010683e+00, 3.711657195983209e-02, + 8.516540731741751e-02, 1.846582481253431e+00, + 2.344622479665092e+00, 9.782457560315421e-01, + -5.653027779602048e+00, -4.457764259327358e+00, + 4.396893371992927e+00, -1.781029123711273e+00, + 1.303130356713504e+00, 1.104086803126383e-01, + 8.202058412989162e+00, -1.912136722291750e+00, + -2.386952673947501e+00, -3.800460327948856e-01, + -4.755886984684770e+00, -2.371006165743835e+00, + 9.838525643858385e-01, -1.970780039723673e+00, + 8.843052074705682e-01, 1.209832476084040e+00, + 4.482138368789980e+00, -8.618372361868662e+00, + 2.553617632271545e+00, -2.515247887047421e+00, + -1.310659972997361e+00, -4.744849190522960e+00, + -5.848019862181545e+00, -1.853688693742562e+00, + 2.211012545853336e+00, -3.584550392170881e+00, + 1.466063738740644e+00, 3.285544689641922e-01, + -3.772195284223943e+00, 2.554683615488591e+00, + 3.144799926493497e+00, 5.320412494125712e-01, + -1.617862103749616e+00, 4.505029568623644e+00, + 4.841738890983537e+00, 4.304535776881372e+00, + -1.874383165372953e+00, -1.816379942063393e+00, + -1.559401743047192e+00, -2.563598701770973e+00, + -2.768891360840021e+00, 1.367687274872984e+00, + -3.451402723746730e+00, -2.019459598023322e-01, + -1.849780677085415e+00, -6.075332676660432e+00, + 3.813912169625591e+00, 1.672124644400447e+00, + 2.897552781182169e+00, 4.424347873913383e+00, + 5.144179917427660e-02, 8.187374675121173e+00, + -6.307832248431945e-01, -1.040038800556793e+00, + 9.638242571688409e-01, -1.315676780411990e+00, + 4.548911919541162e+00, -3.110338786263292e+00, + 7.143970293892079e-01, -1.027012942409028e+00, + -6.771408853367402e-01, 7.689403378304775e+00, + 3.266763808523832e+00, -4.939142356746464e+00, + -1.680667450021413e+00, -5.672119730229586e+00, + 2.356251252889546e+00, 3.100152502890580e+00, + -1.718586561803939e+00, -2.848378813769422e+00, + 2.166953458666474e+00, 3.698992188846235e+00, + -2.521465228517406e+00, -1.483012583302976e+00, + -1.385547239860872e+00, 4.278296734710759e+00, + -3.231893261005089e+00, -3.499165263231714e+00, + -3.991693675979821e+00, -2.395612192990072e+00, + -7.551126038911470e+00, -1.253320876553413e+00, + 2.949752700932412e+00, 2.757985604650246e+00, + 2.965000493451835e+00, 5.042549679265349e+00, + -2.276214983007470e+00, -8.605327422795508e-01, + -3.429327683732336e+00, 3.634086096757940e-01, + -9.235956441894744e-01, -1.735427210626094e+00, + -3.923350081315737e+00, -8.961379286949720e-01, + -4.027346047616183e+00, 5.509694798065775e+00, + -6.249259838595390e+00, 6.001167756824247e+00, + -7.578924142578991e-02, 1.719766675858380e+00, + -3.682461675261361e-01, -1.498059975229270e+00, + -2.253868990511852e+00, -2.913770876062722e+00, + 2.473335187775014e+00, -3.737196599784345e-01, + 9.237764999281854e-01, 2.150735454872202e+00, + -2.202133267022083e+00, 1.329217526808811e+00, + 5.396652531710098e-01, 7.023277604131626e-01, + -9.987304998877746e-01, -1.368332585830754e+00, + -2.853890281653036e+00, -3.715099310600182e+00, + 9.673864924730763e-01, -5.576552424605284e+00, + 3.905872946117088e+00, -2.826501461749119e+00, + 6.414605960628869e+00, -5.679739832207091e+00, + 1.909130437842080e+00, 9.955398727323348e-01, + 3.638877999600390e+00, 1.258926814764474e+00, + -8.414044119346986e-03, -3.991754204800230e+00, + 5.146915245537949e+00, 2.820615592243083e-01, + -4.264963191688312e-01, 3.030254700260019e+00, + 4.878539287722978e+00, -1.459376706664878e+00, + -3.399157829092533e+00, 2.161358937350359e+00, + 4.840508528923103e+00, -3.204368559739170e+00, + 2.979104548159062e+00, -1.160028660033101e+00, + -2.413351595413149e-01, 3.101405290332420e+00, + -2.699185197567189e+00, -5.897826029871508e+00, + -6.632075418392930e+00, -1.685124294109047e+00, + -2.674474141084143e+00, -2.563372692380101e+00, + -4.348080993540033e-01, 2.147158479497566e+00, + -3.420180348851098e-01, 2.522012423589093e+00, + 7.249878567998773e-01, -2.068196780758695e+00, + -1.005151880556292e+00, 2.902422614798151e+00, + -3.071761940805398e-01, -5.087687483603752e+00, + -1.792775263985136e+00, 2.641548369598412e-01, + -7.517681005578722e-01, 2.124100304630884e+00, + 3.012458946272915e+00, -6.457540983540484e+00, + -4.860448155416430e-02, -2.372949299902082e+00, + -2.236600417683022e+00, 7.389268307820533e-01, + 1.495969367959897e+00, 1.905315704750479e+00, + 2.398999344465012e+00, -7.411245301876057e-02, + 2.554697843376545e+00, -2.629716248632333e+00, + -2.002497059929998e+00, -3.607915605047513e+00, + 1.320095075374978e+00, -2.156073489574819e+00, + -1.629751417352065e+00, 4.497170521285455e+00, + -2.956509867164297e+00, -3.563268353779526e+00, + -1.949442402686502e+00, -1.250892366300881e+00, + -1.286794839414063e+00, -3.596629452384589e+00, + -5.636666719031170e+00, -3.624154259337708e+00, + -6.973431206318839e+00, 1.575164194822264e+00, + -2.627563294940978e-01, 3.254410744830774e+00, + -6.749332581521038e+00, 4.840407656923577e+00, + 8.039106767692121e-01, 1.470542740596267e-01, + -4.663614419551561e-01, -5.481824322589377e+00, + 5.049316157478017e+00, 4.568830841146547e-01, + -3.884636396852787e+00, -6.613971933163752e+00, + -1.213411748796469e+00, -2.824305441697994e+00, + -2.916435892389174e+00, 2.258137395096740e+00, ]; diff --git a/tests/lib.rs b/tests/lib.rs index 453e8c1..1aa55e3 100644 --- a/tests/lib.rs +++ b/tests/lib.rs @@ -13,6 +13,6 @@ fn forward() { #[test] fn inverse() { let mut data = fixtures::FREQUENCY_DATA.to_vec(); - fft::inverse(&mut data); + fft::inverse(&mut data, true); assert::close(&data, &fixtures::TIME_DATA[..], 1e-13); }