xfft 0.2.1

A library of Fast Fourier Transforms
Documentation
// +--------------------------+
// |                          |
// |    Copyright (c) 2023    |
// |       Keith Cullen       |
// |                          |
// +--------------------------+

use super::*;
use std::f64::consts::PI;

fn compute_snr(s1: &[f64], s2: &[f64]) -> f64
{
    let mut sss: f64 = 0.0;
    let mut ses: f64 = 0.0;
    let mut err: f64;

    for i in 0..s1.len() {
        println!("s1: {}. s2: {}", s1[i], s2[i]);
        sss += s1[i] * s1[i];
        err  = s1[i] - s2[i];
        ses += err * err;
    }
    10.0 * (sss / ses).log10()
}

#[test]
fn forward() {
    const THRESH_SNR: f64 = 300.0;
    const SIZE: usize = 32;
    const PHI: f64 = PI / 3.5;
    const W: f64 = PI / 177.0;
    let mut x1: [f64; SIZE] = [0.0; SIZE];
    let x2: [f64; SIZE] = [
        12.239078880915997161,
        0.258284950888739928,
        0.160212836566969219,
        -2.653745085949461124,
        0.234387948037571653,
        -1.305931990926573638,
        0.248021462284764932,
        -0.855361295555976286,
        0.252782547758519494,
        -0.626169715592798770,
        0.254981195766465940,
        -0.485154368539490011,
        0.256171008871032502,
        -0.388061892290399391,
        0.256883731967520179,
        -0.315932673182518720,
        0.257341246814549507,
        -0.259269344179668071,
        0.257649369632107328,
        -0.212771380102750429,
        0.257863640596742494,
        -0.173230684205333330,
        0.258015347770489067,
        -0.138574554934469352,
        0.258123047432126773,
        -0.107385871351043577,
        0.258198118195194226,
        -0.078642871187757701,
        0.258247614995273378,
        -0.051567977843477927,
        0.258275797148043718,
        -0.025533809816936967];
    let ctx = Ctx::new(SIZE);
    for i in 0..SIZE {
        x1[i] = (W * i as f64 + PHI).cos();
    }
    ctx.fwd(&mut x1[..]);
    let snr = compute_snr(&x1[..], &x2[..]);
    if snr < THRESH_SNR {
        panic!("SNR below threshold");
    }
}

#[test]
fn forward_and_backward() {
    const THRESH_SNR: f64 = 300.0;
    const SIZE: usize = 32;
    const PHI: f64 = PI / 3.5;
    const W: f64 = PI / 177.0;
    let mut x1: [f64; SIZE] = [0.0; SIZE];
    let mut x2: [f64; SIZE] = [0.0; SIZE];
    let ctx = Ctx::new(SIZE);
    for i in 0..SIZE {
        x1[i] = (W * i as f64 + PHI).cos();
        x2[i] = x1[i];
    }
    ctx.fwd(&mut x1[..]);
    ctx.bwd(&mut x1[..]);
    for i in 0..SIZE {
        x1[i] *= 2.0 / SIZE as f64;
    }
    let snr = compute_snr(&x1[..], &x2[..]);
    if snr < THRESH_SNR {
        panic!("SNR below threshold");
    }
}

#[test]
fn unpack() {
    const THRESH_SNR: f64 = 300.0;
    const SIZE: usize = 32;
    const PHI: f64 = PI / 3.5;
    const W: f64 = PI / 177.0;
    let mut x1: [f64; SIZE] = [0.0; SIZE];
    let mut x2: [f64; 2 * SIZE] = [0.0; 2 * SIZE];
    let x3: [f64; 2 * SIZE] = [
        12.239078880915997161,
        0.000000000000000000,
        0.160212836566969219,
        -2.653745085949461124,
        0.234387948037571653,
        -1.305931990926573638,
        0.248021462284764932,
        -0.855361295555976286,
        0.252782547758519494,
        -0.626169715592798770,
        0.254981195766465940,
        -0.485154368539490011,
        0.256171008871032502,
        -0.388061892290399391,
        0.256883731967520179,
        -0.315932673182518720,
        0.257341246814549507,
        -0.259269344179668071,
        0.257649369632107328,
        -0.212771380102750429,
        0.257863640596742494,
        -0.173230684205333330,
        0.258015347770489067,
        -0.138574554934469352,
        0.258123047432126773,
        -0.107385871351043577,
        0.258198118195194226,
        -0.078642871187757701,
        0.258247614995273378,
        -0.051567977843477927,
        0.258275797148043718,
        -0.025533809816936967,
        0.258284950888739928,
        0.000000000000000000,
        0.258275797148043718,
        0.025533809816936967,
        0.258247614995273378,
        0.051567977843477927,
        0.258198118195194226,
        0.078642871187757701,
        0.258123047432126773,
        0.107385871351043577,
        0.258015347770489067,
        0.138574554934469352,
        0.257863640596742494,
        0.173230684205333330,
        0.257649369632107328,
        0.212771380102750429,
        0.257341246814549507,
        0.259269344179668071,
        0.256883731967520179,
        0.315932673182518720,
        0.256171008871032502,
        0.388061892290399391,
        0.254981195766465940,
        0.485154368539490011,
        0.252782547758519494,
        0.626169715592798770,
        0.248021462284764932,
        0.855361295555976286,
        0.234387948037571653,
        1.305931990926573638,
        0.160212836566969219,
        2.653745085949461124];
    let ctx = Ctx::new(SIZE);
    for i in 0..SIZE {
        x1[i] = (W * i as f64 + PHI).cos();
    }
    ctx.fwd(&mut x1[..]);
    ctx.unpack(&mut x2[..], &x1[..]);
    let snr = compute_snr(&x2[..], &x3[..]);
    if snr < THRESH_SNR {
        panic!("SNR below threshold");
    }
}

#[test]
fn pack() {
    const THRESH_SNR: f64 = 300.0;
    const SIZE: usize = 32;
    const PHI: f64 = PI / 3.5;
    const W: f64 = PI / 177.0;
    let mut x1: [f64; SIZE] = [0.0; SIZE];
    let mut x2: [f64; 2 * SIZE] = [0.0; 2 * SIZE];
    let mut x3: [f64; SIZE] = [0.0; SIZE];
    let mut x4: [f64; SIZE] = [0.0; SIZE];
    let ctx = Ctx::new(SIZE);
    for i in 0..SIZE {
        x1[i] = (W * i as f64 + PHI).cos();
        x4[i] = x1[i];
    }
    ctx.fwd(&mut x1[..]);
    ctx.unpack(&mut x2[..], &x1[..]);
    ctx.pack(&mut x3[..], &x2[..]);
    ctx.bwd(&mut x3[..]);
    for i in 0..SIZE {
        x3[i] *= 2.0 / SIZE as f64;
    }
    let snr = compute_snr(&x3[..], &x4[..]);
    if snr < THRESH_SNR {
        panic!("SNR below threshold");
    }
}