Different outputs than FFTW

[TheSirC]

Thank you for your library !

I am comparing the output of this library from the one from FFTW to see if I can port some code using the latter and using the following code :

use rustfft::{num_complex::Complex64, num_traits::Zero, FftPlanner};

fn fft_rustfft(input: Vec<Complex64>) -> Vec<Complex64> {
    let mut planner = FftPlanner::<f64>::new();
    let fft = planner.plan_fft_forward(input.len());
    let mut binding = input;
    let mut buffer = binding.as_mut_slice();
    fft.process(buffer);
    buffer.to_vec()
}

fn fft_fftw(input: Vec<Complex64>) -> Vec<Complex64> {
    use concrete_fftw::array::AlignedVec;
    use concrete_fftw::plan::*;
    use concrete_fftw::types::*;
    let n = input.len();
    let plan: C2CPlan64 = C2CPlan::aligned(&[n], Sign::Forward, Flag::MEASURE).unwrap();
    let mut binding = input.clone();
    let mut a = binding.as_mut_slice();
    let mut b = AlignedVec::new(n);
    plan.c2c(a, &mut b).unwrap();
    b.as_slice().to_owned()
}

fn main() {
    const size: usize = 750;

    let mut test_array = [Complex64::default(); size];

    // Zero vectors
    assert_eq!(
        fft_fftw(test_array.to_vec()),
        fft_rustfft(test_array.to_vec())
    );

    for (i, value) in test_array.iter_mut().enumerate() {
        *value = Complex64 {
            re: (-((i as f64 - size as f64 / 2.0) / (size / 10) as f64).powi(2)).exp(),
            im: 0.0,
        };
    }

    // Gaussian vectors
    assert_ne!(
        fft_fftw(test_array.to_vec()),
        fft_rustfft(test_array.to_vec())
    );
}

I am getting different outputs on the second assertion (no panics if you run it).

Is this behavior expected or am I doing something unexpected ?

[WalterSmuts]

https://docs.rs/rustfft/latest/rustfft/#normalization

[TheSirC]

docs.rs/rustfft/latest/rustfft/#normalization

By normalizing the difference is even worse , by sending back

    buffer
        .iter_mut()
        .map(|value| {
            *value /= (length as f64).sqrt();
            *value
        })
        .collect::<Vec<_>>()

[HEnquist]

I would guess that the input vector is modified by fftw, so that you give different input to the two ffts.

[ejmahler]

Is this just a floating point precision issue? I see it’s doing raw equality checks on floating point values, which is a red flag to me.

It would help to have more information about which FFT sizes this happens at, exactly how they are different, etc. at the moment we’re just throwing out guesses.

On Tue, Jan 3, 2023 at 9:36 AM Henrik Enquist @.***> wrote:

I would guess that the input vector is modified by fftw, so that you give different input to the two ffts.

— Reply to this email directly, view it on GitHub https://github.com/ejmahler/RustFFT/issues/106#issuecomment-1370046299, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAI2M6SKTPVI2SI7NOQ2JN3WQRPR5ANCNFSM6AAAAAATOV6W2Q . You are receiving this because you are subscribed to this thread.Message ID: @.***>

[HEnquist]

@TheSirC Can this be closed?

[macthecadillac]

I would guess that the input vector is modified by fftw, so that you give different input to the two ffts.

The input vector was cloned through to_vec so whatever fft_fftw modified wouldn't propagate to fft_rustfft.

I ran the code by @TheSirC but with proper numerical equality tests. Both functions give identical outputs up to 10^-13

[cpmech]

Hi, I've run this example with my thin wrapper to FFTW and the results are OK

use num_complex::Complex64;
use russell_lab::{complex_vec_approx_eq, complex_vec_copy, ComplexVector, FFTw};
use rustfft::FftPlanner;

fn fft_rfft(u: &mut ComplexVector) {
    let mut planner = FftPlanner::<f64>::new();
    let fft = planner.plan_fft_forward(u.dim());
    fft.process(u.as_mut_data());
}

fn fft_fftw(u: &mut ComplexVector) {
    let mut v = ComplexVector::new(u.dim());
    let mut fft = FFTw::new();
    fft.dft_1d(&mut v, u, false).unwrap();
    complex_vec_copy(u, &v).unwrap();
}

fn main() {
    const SIZE: usize = 750;

    let mut test_rfft = ComplexVector::new(SIZE);
    let mut test_fftw = ComplexVector::new(SIZE);

    // Zero vectors
    fft_rfft(&mut test_rfft);
    fft_fftw(&mut test_fftw);
    complex_vec_approx_eq(&test_rfft, &test_fftw, 1e-15);

    for i in 0..SIZE {
        let ui = Complex64 {
            re: (-((i as f64 - SIZE as f64 / 2.0) / (SIZE / 10) as f64).powi(2)).exp(),
            im: 0.0,
        };
        test_rfft[i] = ui;
        test_fftw[i] = ui;
    }

    // Gaussian vectors
    fft_rfft(&mut test_rfft);
    fft_fftw(&mut test_fftw);
    complex_vec_approx_eq(&test_rfft, &test_fftw, 1e-13);
}

[ejmahler]

Closing for now. Feel free to reopen this if it's still a problem.