diff --git a/src/lib.rs b/src/lib.rs index 6d67a14..1db875d 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -27,6 +27,63 @@ pub fn forward(data: &mut [c64]) { perform(data, n, false); } +/// Perform the forward transform of real data. +/// +/// The number of points should be a power of two. The data are replaced by the +/// positive frequency half of their complex Fourier transform. The real-valued +/// first and last components of the complex transform are returned as elements +/// data[0] and data[1], respectively. +/// +/// ## References +/// +/// 1. William H. Press, Saul A. Teukolsky, William T. Vetterling, Brian P. +/// Flannery, “Numerical Recipes 3rd Edition: The Art of Scientific +/// Computing,” Cambridge University Press, 2007. +pub fn forward_real(data: &mut [f64]) { + use std::f64::consts::PI; + use std::slice::from_raw_parts_mut; + + let n = data.len(); + power_of_two!(n); + + forward(unsafe { + from_raw_parts_mut(data.as_mut_ptr() as *mut _, n / 2) + }); + + let (c1, c2) = (0.5, -0.5); + let (wpr, wpi) = { + let theta = -PI / (n / 2) as f64; + let temp = (0.5 * theta).sin(); + (-2.0 * temp * temp, theta.sin()) + }; + let mut wr = 1.0 + wpr; + let mut wi = wpi; + for i in 1..(n / 2 / 2) { + let i1 = i + i; + let i2 = 1 + i1; + let i3 = n - i1; + let i4 = 1 + i3; + + let h1r = c1 * (data[i1] + data[i3]); + let h1i = c1 * (data[i2] - data[i4]); + let h2r = -c2 * (data[i2] + data[i4]); + let h2i = c2 * (data[i1] - data[i3]); + + data[i1] = h1r + wr * h2r - wi * h2i; + data[i2] = h1i + wr * h2i + wi * h2r; + data[i3] = h1r - wr * h2r + wi * h2i; + data[i4] = -h1i + wr * h2i + wi * h2r; + + let temp = wr; + wr = wr * wpr - wi * wpi + wr; + wi = wi * wpr + temp * wpi + wi; + } + data[n / 2 + 1] *= -1.0; + let temp = data[0]; + data[0] = temp + data[1]; + data[1] = temp - data[1]; +} + /// Perform the inverse transform. /// /// The number of points should be a power of two. @@ -40,7 +97,26 @@ pub fn inverse(data: &mut [c64], scaling: bool) { } } -#[inline(always)] +/// Expand a compressed representation produced by `forward_real`. +pub fn unpack_real(data: &[f64]) -> Vec { + let n = data.len(); + power_of_two!(n); + + let mut cdata = Vec::with_capacity(n); + unsafe { cdata.set_len(n) }; + + cdata[0] = c64(data[0], 0.0); + for i in 1..(n / 2) { + cdata[i] = c64(data[2 * i], data[2 * i + 1]); + } + cdata[n / 2] = c64(data[1], 0.0); + for i in (n / 2 + 1)..n { + cdata[i] = c64(data[2 * (n - i)], -data[2 * (n - i) + 1]); + } + + cdata +} + fn rearrange(data: &mut [c64], n: usize) { let mut target = 0; for position in 0..n { @@ -56,7 +132,6 @@ fn rearrange(data: &mut [c64], n: usize) { } } -#[inline(always)] fn perform(data: &mut [c64], n: usize, inverse: bool) { use std::f64::consts::PI; @@ -83,10 +158,28 @@ fn perform(data: &mut [c64], n: usize, inverse: bool) { } } -#[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; + + #[test] + fn unpack_real() { + let data = (0..4).map(|i| (i + 1) as f64).collect::>(); + assert_eq!(::unpack_real(&data), vec![ + c64(1.0, 0.0), c64(3.0, 4.0), c64(2.0, 0.0), c64(3.0, -4.0), + ]); + + let data = (0..8).map(|i| (i + 1) as f64).collect::>(); + assert_eq!(::unpack_real(&data), vec![ + c64(1.0, 0.0), c64(3.0, 4.0), c64(5.0, 6.0), c64(7.0, 8.0), + c64(2.0, 0.0), c64(7.0, -8.0), c64(5.0, -6.0), c64(3.0, -4.0), + ]); + } +} diff --git a/tests/fixtures.rs b/tests/fixtures.rs index cc8793a..12c7147 100644 --- a/tests/fixtures.rs +++ b/tests/fixtures.rs @@ -1,4 +1,4 @@ -pub const TIME_DATA: [f64; 256] = [ +pub const TIME_DATA: [f64; 2 * 128] = [ 0.242578298897752, 0.013469574513598, 0.383138850044058, 0.414652690484492, 0.067768972864267, 0.993126929734427, @@ -129,7 +129,7 @@ pub const TIME_DATA: [f64; 256] = [ 0.452619465744411, 0.175360766321123, ]; -pub const FREQUENCY_DATA: [f64; 256] = [ +pub const FREQUENCY_DATA_FROM_COMPLEX: [f64; 2 * 128] = [ 6.750656266720340e+01, 6.279874768750683e+01, 7.714838837855442e-02, 5.819389938527752e+00, -5.303978155417409e+00, -2.887549849606836e-01, @@ -259,3 +259,262 @@ pub const FREQUENCY_DATA: [f64; 256] = [ -1.213411748796469e+00, -2.824305441697994e+00, -2.916435892389174e+00, 2.258137395096740e+00, ]; + +pub const FREQUENCY_DATA_FROM_REAL: [f64; 2 * 256] = [ + 1.303053103547102e+02, 0.000000000000000e+00, + 2.581170396731584e+00, 1.851687150427199e-01, + -4.712992966498973e+00, 3.386970116076343e+00, + -3.526081640989057e+00, -6.639837287112855e-01, + 2.470203860537831e+00, 1.880268489878626e+00, + -6.963879517225685e+00, -1.359220918694192e+00, + -3.729301277306694e+00, 1.188929523990150e+00, + 9.470909676459618e-01, -1.133537735465806e+00, + -5.515738367157939e-01, -7.532580478755269e-01, + -3.549075039135673e+00, -7.723961671240310e+00, + -4.554423293554938e+00, 3.928695392899054e+00, + -3.376141283756011e-01, 3.873091478779083e-01, + -2.239363597649238e+00, -2.622333161216238e+00, + -1.209463421383755e+00, 2.081549368638949e+00, + 2.787306134089560e+00, -4.406059514292826e-01, + -5.867907573418165e-01, -1.188914452808486e+00, + -1.143465327856217e+00, 3.381662456994260e+00, + 3.502762058115909e-01, 1.939344791605027e+00, + -3.966269981614962e+00, 2.578889617469152e+00, + 3.152848007718109e+00, -1.882857879463846e+00, + -3.804426226346088e+00, -2.360628066063453e+00, + -1.645310465960079e+00, 1.955745233096744e-01, + 8.727435711779785e-02, 5.200789512413666e+00, + 1.120142690086254e+00, -1.554454960578805e+00, + -3.568747367377906e-01, -1.708135182867247e+00, + -5.343117610163509e+00, 5.248479799162289e+00, + 5.372292947253837e-01, -4.845214246729507e+00, + 6.450571433636654e-02, 1.463660986831282e+00, + 1.455587976273307e+00, -5.712855301747824e+00, + -1.042998871212450e-01, -1.101626615691207e+00, + -4.844351629884981e+00, 2.524159012105676e+00, + -6.773746016629306e+00, -1.640287552626490e+00, + -3.831851193928844e+00, 3.944330086139181e+00, + -1.492165597423655e+00, -5.490578426839169e+00, + 1.689911156425059e+00, 8.268603075350915e-01, + 1.497677302367308e+00, 4.248256795204092e+00, + -3.579498902285684e+00, -1.349669158977704e+00, + 3.044651859669117e+00, 1.746167472193509e+00, + 7.601446066545735e-01, -3.263458881514898e+00, + 5.150135192916220e+00, 2.544519592918017e+00, + -7.037646240419698e-01, 2.824304911746514e+00, + 4.772947230001412e+00, 5.943735485336072e-01, + 3.480135296267813e+00, -1.372346334547339e+00, + 4.028770258877906e+00, 7.235631972573928e+00, + 2.312768769350037e+00, 3.796429682741911e+00, + -1.782620750385735e-01, 6.191786592630609e+00, + -3.719950405898309e+00, 3.654794911963335e+00, + -2.543476747185558e+00, 1.221476001211498e+00, + 1.118615654278817e+00, -1.238455571070932e+00, + -9.959356075165555e-01, -2.054128218014691e+00, + 2.639682132748301e+00, -1.701633154204415e+00, + 2.173264828603092e+00, 8.248948121629960e-01, + -2.798464518993120e+00, 2.355688055046297e+00, + -8.937461326272578e-01, 7.588584122960702e-01, + 2.021391810490292e-02, -1.805385996360099e+00, + -4.681753970646017e+00, -6.445270984734710e+00, + -3.901614280067265e+00, -6.215718415805858e+00, + -4.468302161079752e+00, 6.560811803588420e-01, + -8.057225542680848e-01, 1.502186010252237e+00, + -3.574988339289840e+00, -4.774591638105163e-01, + -2.126408608816460e+00, 6.496133889843878e-01, + 3.088495481715207e+00, -4.835922591012088e+00, + 3.119768676762977e+00, -2.647386306731720e+00, + -7.608791527035966e+00, 1.199605557244359e+00, + -3.991693675979821e+00, 2.395612192990072e+00, + -3.174227772880593e+00, 3.445449943922659e+00, + -1.555563215691436e+00, -4.167697436792233e+00, + -2.644960216780778e+00, 1.689639671556237e+00, + 2.017147084475464e+00, -3.909911542141400e+00, + -1.572925906246435e+00, 2.734328259634699e+00, + 2.238378162968156e+00, -3.333393703264438e+00, + -1.135715370257396e+00, 5.432062981893455e+00, + 3.141032040974913e+00, 4.233118739006381e+00, + -2.244646753286113e+00, -8.133506606215240e+00, + 6.183938698585155e-01, 9.413936219073080e-01, + 5.074411884642281e+00, 2.371137223330093e+00, + 1.508419785650108e+00, 7.575939593955644e-01, + -3.307128656712724e-01, 1.491213952741355e+00, + -1.664463833645840e+00, -7.738272374453387e+00, + 1.691355121676642e+00, -5.149258565119261e+00, + 3.234961768517784e+00, -2.208252455058217e+00, + -3.050344297876335e-01, 5.928476092041176e+00, + -2.585342599501457e+00, 1.416415611654855e-01, + -1.623242793328371e+00, -7.524531068476575e-01, + 3.370243371985926e-02, 3.533526922763765e+00, + 5.114525363780107e-01, 3.372272082430229e+00, + 3.270734032557803e+00, -4.681342238696377e+00, + -2.751931334150639e+00, -2.651729205325561e+00, + 3.840150506416120e+00, -1.699490542466287e+00, + -3.775415231602213e+00, 2.718975366537348e-01, + 2.794228129172385e-01, -5.617586502190708e-01, + 4.044899973907196e+00, 3.871341157699510e+00, + -5.667678788988395e+00, 2.665378472115218e+00, + 2.032171253558434e+00, 5.788737425987882e+00, + 3.842811024005547e+00, 2.182079534549411e+00, + 5.732968806672320e+00, 6.229199225361914e+00, + 2.016971203832223e+00, -3.163328419816367e+00, + 1.125523162622212e+00, -1.354631807011866e+00, + -2.586009495883930e+00, 2.331792485469409e+00, + -2.717460886180259e+00, 1.425577896601245e+00, + 6.404452401830746e+00, -1.278706155866979e+00, + 1.963612499177015e+00, -7.149444742400535e-01, + 2.854512196711251e+00, -1.617625082200842e-01, + -6.170863635190793e-01, 4.618556574885894e+00, + 9.087219524177454e-01, -2.422226101938949e+00, + -1.786745383326708e+00, -1.276937137201355e+00, + 5.230935547165801e+00, -1.293868043086650e+00, + 3.687342499608545e+00, 3.666842103111971e-01, + 4.690446349542567e-01, 3.354686164448831e+00, + 1.459522623361330e+00, -7.836282087117322e-01, + -1.556476042488200e+00, 9.977746254405300e+00, + 3.615874982837306e+00, 7.587231541051001e-01, + -1.663885294732292e+00, -4.893547445581650e+00, + 6.195623059940027e+00, 6.233312200748062e-01, + -7.389630549527272e+00, -1.774569316216395e+00, + -9.065368226930496e-01, -5.980324595351290e+00, + 3.272997662902595e+00, -8.599788876820805e-01, + 2.931287676659427e+00, -4.844159203841096e+00, + -6.026900117940132e+00, -5.010921049858612e+00, + 7.092481670250197e-01, -2.703433464882523e+00, + -3.310419708767940e+00, 3.606335895631941e+00, + -1.115479477600923e+01, -4.306474671574700e+00, + -7.448658248653059e-01, 5.402306741679653e+00, + 7.270283940301631e+00, 1.866320798541251e+00, + 2.470666768058556e+00, 3.302494992111579e+00, + 3.949947040413504e+00, -6.999949435933990e+00, + -1.804396937714905e+00, 8.514196593390326e-01, + -5.420457900742202e+00, -3.376083828388293e+00, + 4.707814979696565e+00, 0.000000000000000e+00, + -5.420457900742202e+00, 3.376083828388293e+00, + -1.804396937714905e+00, -8.514196593390326e-01, + 3.949947040413504e+00, 6.999949435933990e+00, + 2.470666768058556e+00, -3.302494992111579e+00, + 7.270283940301631e+00, -1.866320798541251e+00, + -7.448658248653059e-01, -5.402306741679653e+00, + -1.115479477600923e+01, 4.306474671574700e+00, + -3.310419708767940e+00, -3.606335895631941e+00, + 7.092481670250197e-01, 2.703433464882523e+00, + -6.026900117940132e+00, 5.010921049858612e+00, + 2.931287676659427e+00, 4.844159203841096e+00, + 3.272997662902595e+00, 8.599788876820805e-01, + -9.065368226930496e-01, 5.980324595351290e+00, + -7.389630549527272e+00, 1.774569316216395e+00, + 6.195623059940027e+00, -6.233312200748062e-01, + -1.663885294732292e+00, 4.893547445581650e+00, + 3.615874982837306e+00, -7.587231541051001e-01, + -1.556476042488200e+00, -9.977746254405300e+00, + 1.459522623361330e+00, 7.836282087117322e-01, + 4.690446349542567e-01, -3.354686164448831e+00, + 3.687342499608545e+00, -3.666842103111971e-01, + 5.230935547165801e+00, 1.293868043086650e+00, + -1.786745383326708e+00, 1.276937137201355e+00, + 9.087219524177454e-01, 2.422226101938949e+00, + -6.170863635190793e-01, -4.618556574885894e+00, + 2.854512196711251e+00, 1.617625082200842e-01, + 1.963612499177015e+00, 7.149444742400535e-01, + 6.404452401830746e+00, 1.278706155866979e+00, + -2.717460886180259e+00, -1.425577896601245e+00, + -2.586009495883930e+00, -2.331792485469409e+00, + 1.125523162622212e+00, 1.354631807011866e+00, + 2.016971203832223e+00, 3.163328419816367e+00, + 5.732968806672320e+00, -6.229199225361914e+00, + 3.842811024005547e+00, -2.182079534549411e+00, + 2.032171253558434e+00, -5.788737425987882e+00, + -5.667678788988395e+00, -2.665378472115218e+00, + 4.044899973907196e+00, -3.871341157699510e+00, + 2.794228129172385e-01, 5.617586502190708e-01, + -3.775415231602213e+00, -2.718975366537348e-01, + 3.840150506416120e+00, 1.699490542466287e+00, + -2.751931334150639e+00, 2.651729205325561e+00, + 3.270734032557803e+00, 4.681342238696377e+00, + 5.114525363780107e-01, -3.372272082430229e+00, + 3.370243371985926e-02, -3.533526922763765e+00, + -1.623242793328371e+00, 7.524531068476575e-01, + -2.585342599501457e+00, -1.416415611654855e-01, + -3.050344297876335e-01, -5.928476092041176e+00, + 3.234961768517784e+00, 2.208252455058217e+00, + 1.691355121676642e+00, 5.149258565119261e+00, + -1.664463833645840e+00, 7.738272374453387e+00, + -3.307128656712724e-01, -1.491213952741355e+00, + 1.508419785650108e+00, -7.575939593955644e-01, + 5.074411884642281e+00, -2.371137223330093e+00, + 6.183938698585155e-01, -9.413936219073080e-01, + -2.244646753286113e+00, 8.133506606215240e+00, + 3.141032040974913e+00, -4.233118739006381e+00, + -1.135715370257396e+00, -5.432062981893455e+00, + 2.238378162968156e+00, 3.333393703264438e+00, + -1.572925906246435e+00, -2.734328259634699e+00, + 2.017147084475464e+00, 3.909911542141400e+00, + -2.644960216780778e+00, -1.689639671556237e+00, + -1.555563215691436e+00, 4.167697436792233e+00, + -3.174227772880593e+00, -3.445449943922659e+00, + -3.991693675979821e+00, -2.395612192990072e+00, + -7.608791527035966e+00, -1.199605557244359e+00, + 3.119768676762977e+00, 2.647386306731720e+00, + 3.088495481715207e+00, 4.835922591012088e+00, + -2.126408608816460e+00, -6.496133889843878e-01, + -3.574988339289840e+00, 4.774591638105163e-01, + -8.057225542680848e-01, -1.502186010252237e+00, + -4.468302161079752e+00, -6.560811803588420e-01, + -3.901614280067265e+00, 6.215718415805858e+00, + -4.681753970646017e+00, 6.445270984734710e+00, + 2.021391810490292e-02, 1.805385996360099e+00, + -8.937461326272578e-01, -7.588584122960702e-01, + -2.798464518993120e+00, -2.355688055046297e+00, + 2.173264828603092e+00, -8.248948121629960e-01, + 2.639682132748301e+00, 1.701633154204415e+00, + -9.959356075165555e-01, 2.054128218014691e+00, + 1.118615654278817e+00, 1.238455571070932e+00, + -2.543476747185558e+00, -1.221476001211498e+00, + -3.719950405898309e+00, -3.654794911963335e+00, + -1.782620750385735e-01, -6.191786592630609e+00, + 2.312768769350037e+00, -3.796429682741911e+00, + 4.028770258877906e+00, -7.235631972573928e+00, + 3.480135296267813e+00, 1.372346334547339e+00, + 4.772947230001412e+00, -5.943735485336072e-01, + -7.037646240419698e-01, -2.824304911746514e+00, + 5.150135192916220e+00, -2.544519592918017e+00, + 7.601446066545735e-01, 3.263458881514898e+00, + 3.044651859669117e+00, -1.746167472193509e+00, + -3.579498902285684e+00, 1.349669158977704e+00, + 1.497677302367308e+00, -4.248256795204092e+00, + 1.689911156425059e+00, -8.268603075350915e-01, + -1.492165597423655e+00, 5.490578426839169e+00, + -3.831851193928844e+00, -3.944330086139181e+00, + -6.773746016629306e+00, 1.640287552626490e+00, + -4.844351629884981e+00, -2.524159012105676e+00, + -1.042998871212450e-01, 1.101626615691207e+00, + 1.455587976273307e+00, 5.712855301747824e+00, + 6.450571433636654e-02, -1.463660986831282e+00, + 5.372292947253837e-01, 4.845214246729507e+00, + -5.343117610163509e+00, -5.248479799162289e+00, + -3.568747367377906e-01, 1.708135182867247e+00, + 1.120142690086254e+00, 1.554454960578805e+00, + 8.727435711779785e-02, -5.200789512413666e+00, + -1.645310465960079e+00, -1.955745233096744e-01, + -3.804426226346088e+00, 2.360628066063453e+00, + 3.152848007718109e+00, 1.882857879463846e+00, + -3.966269981614962e+00, -2.578889617469152e+00, + 3.502762058115909e-01, -1.939344791605027e+00, + -1.143465327856217e+00, -3.381662456994260e+00, + -5.867907573418165e-01, 1.188914452808486e+00, + 2.787306134089560e+00, 4.406059514292826e-01, + -1.209463421383755e+00, -2.081549368638949e+00, + -2.239363597649238e+00, 2.622333161216238e+00, + -3.376141283756011e-01, -3.873091478779083e-01, + -4.554423293554938e+00, -3.928695392899054e+00, + -3.549075039135673e+00, 7.723961671240310e+00, + -5.515738367157939e-01, 7.532580478755269e-01, + 9.470909676459618e-01, 1.133537735465806e+00, + -3.729301277306694e+00, -1.188929523990150e+00, + -6.963879517225685e+00, 1.359220918694192e+00, + 2.470203860537831e+00, -1.880268489878626e+00, + -3.526081640989057e+00, 6.639837287112855e-01, + -4.712992966498973e+00, -3.386970116076343e+00, + 2.581170396731584e+00, -1.851687150427199e-01, +]; diff --git a/tests/lib.rs b/tests/lib.rs index 5920722..8754697 100644 --- a/tests/lib.rs +++ b/tests/lib.rs @@ -6,21 +6,44 @@ mod fixtures; #[test] fn forward() { let mut data = fixtures::TIME_DATA.to_vec(); - fft::forward(reinterpret(&mut data)); - assert::close(&data, &fixtures::FREQUENCY_DATA[..], 1e-14); + fft::forward(as_c64_mut(&mut data)); + assert::close(&data, &fixtures::FREQUENCY_DATA_FROM_COMPLEX[..], 1e-14); +} + +#[test] +fn forward_real() { + let mut data = fixtures::TIME_DATA.to_vec(); + { + let mut data = to_c64(&data); + fft::forward(&mut data); + assert::close(as_f64(&data), &fixtures::FREQUENCY_DATA_FROM_REAL[..], 1e-13); + } + { + fft::forward_real(&mut data); + let data = fft::unpack_real(&data); + assert::close(as_f64(&data), &fixtures::FREQUENCY_DATA_FROM_REAL[..], 1e-13); + } } #[test] fn inverse() { - let mut data = fixtures::FREQUENCY_DATA.to_vec(); - fft::inverse(reinterpret(&mut data), true); + let mut data = fixtures::FREQUENCY_DATA_FROM_COMPLEX.to_vec(); + fft::inverse(as_c64_mut(&mut data), true); assert::close(&data, &fixtures::TIME_DATA[..], 1e-14); } -fn reinterpret<'l>(slice: &'l mut [f64]) -> &'l mut [fft::c64] { +fn as_f64<'l>(slice: &'l [fft::c64]) -> &'l [f64] { unsafe { - let length = slice.len(); - assert!(length % 2 == 0, "the number of elements should be even"); - std::slice::from_raw_parts_mut(slice.as_mut_ptr() as *mut _, length / 2) + std::slice::from_raw_parts(slice.as_ptr() as *const _, 2 * slice.len()) } } + +fn as_c64_mut<'l>(slice: &'l mut [f64]) -> &'l mut [fft::c64] { + unsafe { + std::slice::from_raw_parts_mut(slice.as_mut_ptr() as *mut _, slice.len() / 2) + } +} + +fn to_c64(slice: &[f64]) -> Vec { + slice.iter().map(|&re| fft::c64(re, 0.0)).collect() +}