ejmahler / RustFFT

RustFFT is a high-performance FFT library written in pure Rust.
Apache License 2.0
678 stars 47 forks source link

Wrong results with plan_fft_inverse #94

Closed nilgoyette closed 1 year ago

nilgoyette commented 1 year ago

Hi. I'm porting an application from C++ to Rust and I obtain strange results from rustfft.

I have a vector v of 112 (nb_rows) elements, already filled with complex data.

let fft_row = planner.plan_fft_forward(nb_rows);
let mut scratch = vec![Complex::default(); fft_row.get_inplace_scratch_len()];
fft_row.process_with_scratch(&mut v, &mut scratch);

This is working perfectly. I obtain the same results with NumPy and Eigen.

When I try to use plan_fft_inverse however...

let ifft_row = planner.plan_fft_inverse(nb_rows);
let mut v = vec![
Complex::new(5832073.1052856445, 0.0),
Complex::new(-3049261.8609998133, -2288111.7097871057),
Complex::new(-129389.3950178826, 2258022.797078504),
Complex::new(438990.11088426923, -631833.654686513),
Complex::new(144056.06915030538, -4551.628559942569),
Complex::new(-242084.7372503671, -443724.7675956719),
Complex::new(-503436.3796742589, 802147.7260735594),
Complex::new(1038730.916229399, -597.1541990541982),
Complex::new(-489654.58492202143, -499897.2565499133),
Complex::new(66369.17936628888, 56452.586691041215),
Complex::new(-258348.94420511462, 169146.02619753848),
Complex::new(358939.57714444696, 152139.69157675526),
Complex::new(-123420.64002263799, -199606.9070699016),
Complex::new(-173873.36311204438, 15967.91987310239),
Complex::new(118928.29534449398, 148326.80456079103),
Complex::new(26263.691789717763, -37610.4786907175),
Complex::new(108502.46167855956, -68682.65915318011),
Complex::new(-246467.22838259608, -50882.52874742227),
Complex::new(124749.60138293434, 238415.1794321786),
Complex::new(-16181.09092776988, -210212.7023419696),
Complex::new(-13153.314898120372, 92133.82420056402),
Complex::new(16290.17449797971, -69207.44392908699),
Complex::new(-64927.019606064154, 115692.51751086662),
Complex::new(108763.65461218973, -68479.59069933197),
Complex::new(-109675.06337793106, -21128.88014280519),
Complex::new(79694.24825696828, 75696.52528360071),
Complex::new(-48355.84542114351, -65826.23394951625),
Complex::new(36583.79929987297, 79843.17523646136),
Complex::new(-25408.950291951496, -98219.97228495279),
Complex::new(-40628.24404739229, 91102.80096031739),
Complex::new(64330.71353648114, -38000.65567398356),
Complex::new(-43610.30287604069, 997.3307829497949),
Complex::new(17555.081180546793, 23197.711745610766),
Complex::new(-12090.060161055571, -27428.11386058035),
Complex::new(6787.734442903624, 22435.151098983228),
Complex::new(-6391.645054227046, -24579.30598717564),
Complex::new(-17916.607126975552, 27779.677629061836),
Complex::new(25341.348840095903, -4852.624238974698),
Complex::new(-14858.975303855528, -14800.93986824459),
Complex::new(3216.3683765778937, 15122.883218650342),
Complex::new(-3881.78222042828, -11838.377440844368),
Complex::new(1639.7748166714841, 11958.31102811436),
Complex::new(-2377.785347294428, -10778.76950716625),
Complex::new(-2380.8275263176415, 10309.355428373097),
Complex::new(4490.237013621637, -5217.682973700906),
Complex::new(-4584.819612130072, 1105.4652113463499),
Complex::new(3912.482593593805, 3.5625575710960287),
Complex::new(-3179.4295863959014, -1960.6974920150055),
Complex::new(1362.5368858986083, 2433.6735357510292),
Complex::new(-1971.5755762742274, -1098.5451515243583),
Complex::new(962.5664958224329, 1469.0557363163402),
Complex::new(-333.85384962638176, -858.594021963471),
Complex::new(-21.523455285972137, 92.72770171267979),
Complex::new(-123.43172661820809, 171.5767408920077),
Complex::new(81.85566733637236, -97.71373324257463),
Complex::new(-39.92986915222728, 11.38712786275599),
Complex::new(0.0, 0.0),
Complex::new(-39.9298691522274, -11.387127862755955),
Complex::new(81.85566733637381, 97.71373324257651),
Complex::new(-123.43172661820837, -171.57674089200773),
Complex::new(-21.523455285972474, -92.72770171268014),
Complex::new(-333.85384962638074, 858.5940219634657),
Complex::new(962.5664958224345, -1469.0557363163398),
Complex::new(-1971.5755762742235, 1098.545151524355),
Complex::new(1362.536885898608, -2433.6735357510242),
Complex::new(-3179.4295863959032, 1960.6974920150026),
Complex::new(3912.4825935937997, -3.5625575710997817),
Complex::new(-4584.819612130066, -1105.4652113463503),
Complex::new(4490.237013621644, 5217.682973700908),
Complex::new(-2380.8275263176342, -10309.355428373094),
Complex::new(-2377.7853472944303, 10778.769507166238),
Complex::new(1639.7748166714691, -11958.311028114347),
Complex::new(-3881.782220428261, 11838.377440844351),
Complex::new(3216.3683765779033, -15122.883218650355),
Complex::new(-14858.975303855526, 14800.939868244594),
Complex::new(25341.348840095907, 4852.624238974722),
Complex::new(-17916.607126975545, -27779.677629061785),
Complex::new(-6391.645054227022, 24579.305987175616),
Complex::new(6787.734442903583, -22435.1510989832),
Complex::new(-12090.060161055544, 27428.11386058033),
Complex::new(17555.081180546804, -23197.711745610766),
Complex::new(-43610.30287604063, -997.3307829498154),
Complex::new(64330.71353648113, 38000.655673983565),
Complex::new(-40628.24404739228, -91102.80096031731),
Complex::new(-25408.950291951493, 98219.97228495278),
Complex::new(36583.79929987292, -79843.1752364613),
Complex::new(-48355.84542114344, 65826.23394951622),
Complex::new(79694.24825696833, -75696.52528360074),
Complex::new(-109675.06337793102, 21128.880142805196),
Complex::new(108763.6546121897, 68479.59069933198),
Complex::new(-64927.01960606412, -115692.51751086648),
Complex::new(16290.174497979733, 69207.44392908698),
Complex::new(-13153.314898120445, -92133.82420056402),
Complex::new(-16181.090927769885, 210212.70234196956),
Complex::new(124749.60138293437, -238415.17943217856),
Complex::new(-246467.22838259605, 50882.5287474222),
Complex::new(108502.46167855953, 68682.65915318011),
Complex::new(26263.691789717745, 37610.47869071766),
Complex::new(118928.29534449399, -148326.80456079103),
Complex::new(-173873.3631120445, -15967.919873102443),
Complex::new(-123420.64002263799, 199606.9070699016),
Complex::new(358939.57714444696, -152139.69157675517),
Complex::new(-258348.9442051145, -169146.02619753854),
Complex::new(66369.17936628876, -56452.586691041244),
Complex::new(-489654.58492202166, 499897.25654991344),
Complex::new(1038730.9162293993, 597.1541990543087),
Complex::new(-503436.3796742589, -802147.7260735597),
Complex::new(-242084.7372503671, 443724.76759567193),
Complex::new(144056.06915030538, 4551.628559942616),
Complex::new(438990.1108842691, 631833.6546865131),
Complex::new(-129389.39501788229, -2258022.797078504),
Complex::new(-3049261.860999814, 2288111.7097871047)
];
let mut scratch = vec![Complex::default(); ifft_row.get_inplace_scratch_len()];
ifft_row.process_with_scratch(&mut v, &mut scratch);

The first 5 elements look like that.

rustfft                     NumPy && Eigen
129099.64136202168, -0.0    1152.6753-0i
-22150.385800415184, -0.0   -197.7713+0i
3802.6734404931776, 0.0     33.9524-0i
-665.6548425187357, -0.0    -5.9433-0i
191.2556146387942, 0.0      1.7076-0i
...

TBH, I don't know what the exact result should be, but I would have believed rustfft, NumPy and Eigen should be mostly identical ± epsilon. So I'm led to believe that rustfft has a bug.

HEnquist commented 1 year ago

The vector you posted here is symmetric around the middle element, not conjugate symmetric. This won't give you a real result. For example, since the second element is Complex::new(-3049261.861, 2288111.7098) I would expect the last element to be Complex::new(-3049261.861, -2288111.7098). Are you comparing with the right values?

For the data above, Numpy gives me this for the first five elements: 1152.6753875 +166290.37557857j 2160.10844305+138306.04022324j 3010.74789284 +98521.72931408j 2787.28305285 +84694.76667457j 2783.03673224 +79827.58606261j That's not what you get from rustfft, so there may be somethng strange here. Do you get the same result if you use the out-of-place transform?

nilgoyette commented 1 year ago

@HEnquist Thank you for checking. You're right, I'm sorry, I used find & replace to format the list of complex numbers for github and I destroyed all negatives :| It's fixed now! I copied my code snippet from the first post and now it gives what I said.

HEnquist commented 1 year ago

Alright, then it's all clear. It's the normalization. Rustfft does not normalize the result, while the others do. See https://docs.rs/rustfft/latest/rustfft/#normalization Divide your result by 112 (the length) and you should have matching numbers.

nilgoyette commented 1 year ago

Oh! This is good to know. I lost several hours debugging this hahaha. It was in the documentation. No excuse :)

Thank you!