FFTS output difference with fft in Matlab

[dzungpv]

Hello, I am converting code from Matlab to java/c to use in Android device. I can do this with pure java fft lib like Jtransform and Jwave but it slow. I want to use FFTS but the out put not the same, i don't know why.

[linkotec]

Hi! It cannot be a big problem, but it would help if you could provide code snippet how you have used FFTS? And have you looked at the examples and test cases?

[dzungpv]

@linkotec I use java jni call:

                       FFTS fft = FFTS.real(FFTS.FORWARD, FFT_size);
			float[] Y1 = new float[FFT_size*2];
			float[] y1 = System.arraycopy(displaySamples, 0, y1, 0, FFT_size);
                        fft.execute(y1, Y1);

And this code use JWave https://github.com/cscheiblich/JWave (output same with matlab fft)

			Transform t = new Transform( new DiscreteFourierTransform( ) );
			double[] Y1 = new double[FFT_size];
			double[] y1 = System.arraycopy(displaySamples, 0, y1, 0, FFT_size);
			Y1 = t.forward(y1); // 1-D DFT forward

But it give difference output. I guess it may has the same spectrum but i am not sure.

[linkotec]

I assume that "displaySamples" contains complex numbers, in-order re1 im1 re2 im2... And the actual FFT length is FFT_size/2, as FFT_size is the size of float array.

Right away I don't see anything wrong, but what is magnitude of difference between JWave and FFTS? FFTS uses internally 32 bit float, and JWave uses Java native 64 bit double all the way, so that will result in loss of accuracy when using FFTS.

[dzungpv]

@linkotec I have edited the code above to use real, not complex, it is 1d array with float or double. And then i get magnitude by those code:

			// Calculate the Real and imaginary and Magnitude.
			for (int i = 0; i < FFT_size/2 -1; i++) {
				// real is stored in first part of array
				re[i] = Y1[i*2];
				// imaginary is stored in the sequential part
				im[i] = Y1[i*2+1];
				// magnitude is calculated by the square root of (imaginary^2 + real^2)
				YY1[i] = Math.sqrt(Math.pow(re[i], 2.0) + Math.pow(im[i], 2.0));
			}

But it far difference, i don't know if i miss understand about ffts. I tested with simple array and they are not the same.

[dzungpv]

This is the code i test:
Matlab:

x = [0,1,2,3,4,5,6,7]
X=fft(x,8);

28.0000000000000 + 0.00000000000000i    
-4.00000000000000 + 9.65685424949238i
-4.00000000000000 + 4.00000000000000i
-4.00000000000000 + 1.65685424949238i
-4.00000000000000 + 0.00000000000000i
-4.00000000000000 - 1.65685424949238i
-4.00000000000000 - 4.00000000000000i
-4.00000000000000 - 9.65685424949238i

FFTS:

FFTS fft = FFTS.real(FFTS.FORWARD, 8);
fft.execute(x,output);

28.000000 + 0.000000i
-4.000000 + 9.656855i
-4.000000 + 4.000000i
-4.000000 + 1.656854i
0.000000 + 0.000000i
0.000000 + 0.000000i
0.000000 + 0.000000i
0.000000 + 0.000000i

[linkotec]

I am sure that the problem is lack of FFTS documentation #12, but there might also be a bug in JNI. I am not using Java or JNI interface so I can only provide you with C versions and explanation.

Here is MSVC C code version of your complex problem:

float __declspec(align(32)) x[16] = {0,0,1,0,2,0,3,0,4,0,5,0,6,0,7,0};
float __declspec(align(32)) X[16];
ffts_plan_t *plan;
plan = ffts_init_1d(8, FFTS_FORWARD);
ffts_execute(plan, x, X);
ffts_free(plan);

Input x has 8 complex numbers, where imaginary parts are set to zero, thus real numbers. After execute output X contains 8 complex numbers. The result is "same" as with Matlab:

0. 28.000000 + 0.00000000i
1. -4.0000000 + 9.6568546i
2. -4.0000000 + 4.0000000i
3. -4.0000000 + 1.6568542i
4. -4.0000000 + 0.00000000i
5. -4.0000000 - 1.6568542i
6. -4.0000000 - 4.0000000i
7. -4.0000000 - 9.6568546i

But your input is pure real numbers, so you should use real valued transform. Now if you look at https://github.com/anthonix/ffts/blob/master/include/ffts.h , comment says "The output of a real-to-complex transform is N/2+1 complex numbers, where the redundant outputs have been omitted."
This is reason why you see different results compared to Matlab.

When input data is purely real numbers, DFT output satisfies the Hermitian redundancy; output[n-i] is complex conjugate of output[i]. Output[0] is the DC, output[N/2] is the Nyquist.

Here is real version in C:

float __declspec(align(32)) x[8] = {0,1,2,3,4,5,6,7};
float __declspec(align(32)) X[8+2];
ffts_plan_t *plan;
plan = ffts_init_1d_real(8, FFTS_FORWARD);
ffts_execute(plan, x, X);
ffts_free(plan);

Input x has 8 real numbers, and output is (8/2)+1=5 complex numbers,

0. 28.000000 + 0.00000000i
1. -4.0000000 + 9.6568546i
2. -4.0000000 + 4.0000000i
3. -4.0000000 + 1.6568542i
4. -4.0000000 + 0.00000000i

I didn't find any open issues or commits regarding JNI, but for some reason your test is missing the Nyquist frequency, last component -4.0000000 + 0.00000000i

This needs more testing.

[dzungpv]

Yes you are right, JNI code has some bug (maybe JNI author not understand real DFT in FFTS ) https://github.com/anthonix/ffts/blob/master/java/jni/ffts_jni.c he use the same size, it must be the length of input and output array. After fix this, the result the same with you.
And the off topic, i don't know why can not build your https://github.com/linkotec/ffts for android? I must use this repo which merge some class: https://github.com/billthefarmer/ffts-android. And your repo don't include issue tab, it will more easy to report bug in your repo?

[linkotec]

If you can push your fix, I would like to include that to my repo. I have totally discarded GNU autoconf configure scripts, so there may be a lot of things to fix. I cannot make any promises for Android support, but I have now enabled the issue tab.