Correctness tests?
Opened this issue · 9 comments
Is ffts ever checked for correctness, because I'm getting highly suspicious about the correctness. Comparing results from FFTS and matlab, leads me to believe that FFTS is broken for 2D non-square transforms.
Given the (ahem) sparse nature of the example, the complete lack of any test cases and the fact that it won't compile out of the box on non SSE architectures (https://github.com/anthonix/ffts/blob/master/src/ffts_trig.c#L884) I'd say you have cause for concern.
That said non-square may simply be an undocumented limitation. Albeit an unexpected one in what is very useful for image processing etc.
I figured it out (wasted literally days to it)... ffts_init_2d() takes first the height and as second parameter the width...
So, this fixes the non-square issue I had. However, when using non power of two dimensions in the 2D case, it is still broken. In case one dimension is odd (2N+1), transpose segfaults, because it assumes that it can transpose nice blocks.
Also figured the last problem out: 2D FFTS only works on sizes of multiples of 8 complex numbers. This definitely should be documented somewhere.
Also figured the last problem out: 2D FFTS only works on sizes of multiples of 8 complex numbers. This definitely should be documented somewhere.
I might be observing something similar for 1D FFTS. I'm trying to determine when the output of FFTS matches that of MATLAB with a basic real to complex FFT.
If I pass an nfft which is a multiple of 8 (the first parameter of ffts_plan_1d), the result of FFTS comes out strangely different from MATLAB. For other values, FFTS matches MATLAB.
I've been trying to debug a piece of my code for a few days now, which uses the following default sizes:
- input signal size 512
- nfft = 512*2 = 1024
The output ends up wrong and totally different from similar examples in MATLAB and Python (using numpy + scipy).
If I use an input signal size of 513 and nfft = 1026 (i.e. non multiples of 8), FFTS == MATLAB == numpy/scipy.
edit: next thing I'll be verifying if this is related to memory alignment for sizes that are multiples of 8.
also worth checking against https://github.com/linkotec/ffts as that fork seems to have had more input over the last few years than this repo.
I compiled the above fork and still, using 1022 and 1026 produce correct results, and 1024 produces bad results. It also seems unrelated to memory alignment (I used the alignment code here: #55 (comment))
Sure. Here's an example showing several similar tests with n = 1022, 1023, 1024, 1025 - all except 1024 have similar values:
#include <iostream>
#include <complex>
#include <vector>
#include <numeric>
#include "ffts/ffts.h"
void print_fft(std::vector<float> &arr) {
int n = arr.size();
int nfft = 2*n;
ffts_plan_t *p = ffts_init_1d(nfft, FFTS_FORWARD);
std::vector<std::complex<float>> result(nfft);
// fill real input into first half of result complex array
for (int i = 0; i < n; ++i) {
result[i] = std::complex<float>(arr[i], 0.0);
}
// explicitly fill second half with zeros
std::fill(result.begin()+n, result.end(), std::complex<float>(0.0, 0.0));
// fft in-place
ffts_execute(p, result.data(), result.data());
std::cout << "fft result for n = " << n << ", nfft = " << nfft << ", first 8 entries" << std::endl;
for (std::size_t i = 0; i < 8; ++i) {
std::cout << result[i] << std::endl;
}
std::cout << std::endl;
std::cout << std::endl;
}
int
main(int argc, char **argv)
{
std::vector<float> test0(1022);
std::iota(test0.begin(), test0.end(), 1.0);
std::vector<float> test1(1023);
std::iota(test1.begin(), test1.end(), 1.0);
std::vector<float> test2(1024);
std::iota(test2.begin(), test2.end(), 1.0);
std::vector<float> test3(1025);
std::iota(test3.begin(), test3.end(), 1.0);
std::vector<float> test4(1025);
std::iota(test4.begin(), test4.end(), 1.0);
print_fft(test0);
print_fft(test1);
print_fft(test2);
print_fft(test3);
print_fft(test4);
return 0;
}Output:
fft result for n = 1022, nfft = 2044, first 8 entries
(522753,-0.00390625)
(-211145,-333120)
(-510.998,166234)
(-23005.6,-111039)
(-510.996,83116.3)
(-7954.44,-66622.7)
(-511,55410)
(-3807.7,-47586.8)
fft result for n = 1023, nfft = 2046, first 8 entries
(523776,0)
(-211559,-333771)
(-511.499,166560)
(-23051.1,-111256)
(-511.497,83279.1)
(-7970.51,-66753)
(-511.5,55518.5)
(-3815.65,-47679.8)
fft result for n = 1024, nfft = 2048, first 8 entries
(222958,314511)
(71223.6,-126610)
(9896.8,43777.8)
(87232.7,34883.6)
(59082.2,-94674.4)
(-84855.2,182826)
(-185408,9416.42)
(-61190.4,234029)
fft result for n = 1025, nfft = 2050, first 8 entries
(525825,-0.00156403)
(-212388,-335077)
(-512.474,167212)
(-23142.4,-111691)
(-512.498,83605)
(-8002.71,-67014.1)
(-512.498,55735.8)
(-3831.59,-47866.3)
fft result for n = 1025, nfft = 2050, first 8 entries
(525825,-0.00156403)
(-212388,-335077)
(-512.474,167212)
(-23142.4,-111691)
(-512.498,83605)
(-8002.71,-67014.1)
(-512.498,55735.8)
(-3831.59,-47866.3)
MATLAB has values similar to 1022,1023,1025,1026. Seems like 1024 on FFTS is the outlier:
>> x = [1:1024];
>> X = fft(x, 2048);
>> display(X(1:9));
1.0e+05 *
Columns 1 through 3
5.2480 + 0.0000i -2.1197 - 3.3442i -0.0051 + 1.6689i
Columns 4 through 6
-0.2310 - 1.1147i -0.0051 + 0.8344i -0.0799 - 0.6688i
Columns 7 through 9
-0.0051 + 0.5563i -0.0382 - 0.4777i -0.0051 + 0.4172i
>>
>> x = [1:1022];
>> X = fft(x, 2042);
>> display(X(1:9));
1.0e+05 *
Columns 1 through 3
5.2275 + 0.0000i -2.1175 - 3.3247i 0.0051 + 1.6591i
Columns 4 through 6
-0.2398 - 1.1082i 0.0051 + 0.8295i -0.0896 - 0.6649i
Columns 7 through 9
0.0051 + 0.5530i -0.0482 - 0.4749i 0.0051 + 0.4148i
If I modify nfft = 2*n + 1 or nfft = 2*n - 1, the result seems to approach correctness once more:
fft result for n = 1024, nfft = 2047, first 8 entries
(524800,-0.015625)
(-212278,-334098)
(0.375061,166723)
(-23586,-111366)
(0.380371,83362.1)
(-8490.68,-66820.2)
(0.37793,55575.2)
(-4331.76,-47729.2)
fft result for n = 1024, nfft = 2049, first 8 entries
(524800,-0.00390625)
(-211668,-334749)
(-1024.62,167046)
(-22607.2,-111579)
(-1024.61,83518.8)
(-7482.32,-66942.9)
(-1024.58,55674.4)
(-3315.28,-47811.4)