Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Correctness tests? #65

Open
mcourteaux opened this issue Sep 29, 2017 · 9 comments
Open

Correctness tests? #65

mcourteaux opened this issue Sep 29, 2017 · 9 comments

Comments

@mcourteaux
Copy link

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.

@g40
Copy link

g40 commented Sep 29, 2017

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.

@mcourteaux
Copy link
Author

mcourteaux commented Sep 29, 2017

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.

@mcourteaux
Copy link
Author

Also figured the last problem out: 2D FFTS only works on sizes of multiples of 8 complex numbers. This definitely should be documented somewhere.

@sevagh
Copy link

sevagh commented Aug 12, 2020

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.

@g40
Copy link

g40 commented Aug 12, 2020

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.

@sevagh
Copy link

sevagh commented Aug 12, 2020

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))

@g40
Copy link

g40 commented Aug 12, 2020

@sevagh can you share a minimal test example? I will try it here.

@sevagh
Copy link

sevagh commented Aug 12, 2020

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

@sevagh
Copy link

sevagh commented Aug 12, 2020

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)

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

3 participants