CMSIS DSP library for ARM and MATLAB gives two different results of FFT

[langwadt]

ARM code uses 32-bit floats, Matlab is likely to use 80-bit long doubles or so. So what you see may quite well be rounding errors.

JW

I used "single()" casting to input signal of MATLAB to make the MATLAB data format to match 32 bit floating point format of STM32. Result is same.

But the intermediate calculations are still done with float64. 

are you sure?

x=single(rand(1024,1));
plot(abs(fft(double(x)))-abs(fft(x)))

is not zeros

[Nominal Animal]

DataINFFT.txt file contains input samples for both MATLAB fft function and STM32 PreFFTProcess function shown in my original post. STM32fftvsMATLABfft.PNG shows MATLAB fft result for these samples in orange. STM32 fft result for these samples is blue.

Now we're cooking!

The rounding error between single- and double-precision DFT for these samples is less than ±0.000050, and magnitude differences less than ±0.000053.
See the attached dft.txt file for the float (single-precision IEEE 754 Binary32) and double (Binary64) magnitudes and real and imaginary components, listed at sufficient precision to be unambiguous/exact when converted to their respective types.  (I used a dedicated C program and FFTW3 library for this.)
As expected, the Matlab graph agrees with these results.

Now, if I only use 128 of your samples, inserting a zero after each sample, I get your "bad" results pretty exactly; see the attached bad.txt.

Thus, your error is interleaving zeros ("imaginary components") between the consecutive samples.
The DFT your CMSIS DSP library is doing a real-to-complex Discrete Fourier Transform (r2c in FFTW3 planning system); you simply supply it with the real components only, as it does the computations assuming all imaginary components are zeroes (which also speeds things up considerably, which is what you want). This is a computational optimization, and not an approximation or "shortcut", by the way: its results are the same (to within rounding error, i.e. somewhere in the fifth decimal digit, at most, for this sample sequence, using Binary32/float/f32) as doing the computations the "hard way".

[Nominal Animal]

(In more simple terms, syntax333's bug is in using PreFFTProcess() at all; just like Kalvin suspected, earlier in the thread.  The proof is in the two data sets, bad.txt duplicating the error, with dft.txt showing that when omitted/done correctly, the f32/IEEE 754 Binary32 DFT implementation yields results that only minimally differ from the f64/IEEE Binary64 DFT, using a similar DFT-on-real-data-yielding-complex-results from FFTW3 library.)

[turny]

The DFT your CMSIS DSP library is doing a real-to-complex Discrete Fourier Transform (r2c in FFTW3 planning system);

Except it does not. arm_cfft_f32 is complex to complex, arm_rfft_f32 is real to complex.
But I agree that most likely PreFFTProcess is the culprit. One thing syntax333 could do is to take the CMSIS FFT output and do inverse FFT of it in MATLAB and compare the result with the original input. Or just post the code here.

[Nominal Animal]

The DFT your CMSIS DSP library is doing a real-to-complex Discrete Fourier Transform (r2c in FFTW3 planning system);

Except it does not. arm_cfft_f32 is complex to complex, arm_rfft_f32 is real to complex.

Ah, I didn't check that; thanks.

In any case, the 512-sample DFT of

Code: [Select]

     12, 0,  10, 0,  13, 0,   7, 0,  11, 0,  11, 0,  18, 0,  20, 0,  16, 0,  12, 0,   8, 0,   6, 0,   2, 0,   4, 0,   0, 0,   2, 0,
      1, 0,  -1, 0,  -4, 0,  -3, 0,  -9, 0, -12, 0, -11, 0, -14, 0, -10, 0,  -8, 0,  -4, 0,  -2, 0,   4, 0,   8, 0,  11, 0,  14, 0,
     14, 0,  16, 0,  20, 0,  17, 0,  18, 0,  16, 0,  13, 0,  12, 0,   8, 0,   4, 0,  -3, 0,  -4, 0,  -5, 0, -12, 0, -15, 0, -20, 0,
    -23, 0, -25, 0, -24, 0, -21, 0, -23, 0, -18, 0, -16, 0, -11, 0,  -3, 0,   0, 0,   4, 0,  10, 0,  11, 0,  15, 0,  13, 0,  14, 0,
     14, 0,  16, 0,  19, 0,  15, 0,  20, 0,  21, 0,  22, 0,  20, 0,  17, 0,  17, 0,  15, 0,  13, 0,  14, 0,   8, 0,  11, 0,  11, 0,
     11, 0,   9, 0,   1, 0,  -2, 0,  -7, 0,  -8, 0,  -9, 0, -14, 0, -14, 0, -17, 0, -13, 0, -11, 0, -11, 0,  -7, 0,  -9, 0,  -8, 0,
    -11, 0, -10, 0, -10, 0,  -9, 0,  -6, 0, -10, 0,  -8, 0,  -3, 0,  -4, 0,  -3, 0,  -4, 0,  -6, 0,  -8, 0,  -4, 0,  -2, 0,   0, 0,
      3, 0,   9, 0,  14, 0,  22, 0,  25, 0,  24, 0,  25, 0,  24, 0,  27, 0,  23, 0,  19, 0,  10, 0,   2, 0,   2, 0,   1, 0,   0, 0,
      0, 0,   0, 0,   0, 0,   0, 0,   0, 0,   0, 0,   0, 0,   0, 0,   0, 0,   0, 0,   0, 0,   0, 0,   0, 0,   0, 0,   0, 0,   0, 0,
      0, 0,   0, 0,   0, 0,   0, 0,   0, 0,   0, 0,   0, 0,   0, 0,   0, 0,   0, 0,   0, 0,   0, 0,   0, 0,   0, 0,   0, 0,   0, 0,
      0, 0,   0, 0,   0, 0,   0, 0,   0, 0,   0, 0,   0, 0,   0, 0,   0, 0,   0, 0,   0, 0,   0, 0,   0, 0,   0, 0,   0, 0,   0, 0,
      0, 0,   0, 0,   0, 0,   0, 0,   0, 0,   0, 0,   0, 0,   0, 0,   0, 0,   0, 0,   0, 0,   0, 0,   0, 0,   0, 0,   0, 0,   0, 0,
      0, 0,   0, 0,   0, 0,   0, 0,   0, 0,   0, 0,   0, 0,   0, 0,   0, 0,   0, 0,   0, 0,   0, 0,   0, 0,   0, 0,   0, 0,   0, 0,
      0, 0,   0, 0,   0, 0,   0, 0,   0, 0,   0, 0,   0, 0,   0, 0,   0, 0,   0, 0,   0, 0,   0, 0,   0, 0,   0, 0,   0, 0,   0, 0,
      0, 0,   0, 0,   0, 0,   0, 0,   0, 0,   0, 0,   0, 0,   0, 0,   0, 0,   0, 0,   0, 0,   0, 0,   0, 0,   0, 0,   0, 0,   0, 0,
      0, 0,   0, 0,   0, 0,   0, 0,   0, 0,   0, 0,   0, 0,   0, 0,   0, 0,   0, 0,   0, 0,   0, 0,   0, 0,   0, 0,   0, 0,   0, 0where each value is a complex number if one does a complex DFT, and real if one does a real DFT –– i.e., first 128 samples from DataINFFT.txt with an extra zero after each sample, plus an additional 256 zeroes –– reproduces OP's blue curve.

The expected 256-sample DFT of (complex samples if one does a complex DFT, and real if one does a real DFT)

Code: [Select]

    12, 10, 13,  7, 11, 11, 18, 20, 16, 12,  8,  6,  2,  4,  0,  2,
     1, -1, -4, -3, -9,-12,-11,-14,-10, -8, -4, -2,  4,  8, 11, 14,
    14, 16, 20, 17, 18, 16, 13, 12,  8,  4, -3, -4, -5,-12,-15,-20,
   -23,-25,-24,-21,-23,-18,-16,-11, -3,  0,  4, 10, 11, 15, 13, 14,
    14, 16, 19, 15, 20, 21, 22, 20, 17, 17, 15, 13, 14,  8, 11, 11,
    11,  9,  1, -2, -7, -8, -9,-14,-14,-17,-13,-11,-11, -7, -9, -8,
   -11,-10,-10, -9, -6,-10, -8, -3, -4, -3, -4, -6, -8, -4, -2,  0,
     3,  9, 14, 22, 25, 24, 25, 24, 27, 23, 19, 10,  2,  2,  1,  0,
    -2, -2,  0,  5, 11, 11,  7,  6,  5,  5,  1, -2, -2, -6, -5, -9,
   -16,-21,-22,-27,-36,-37,-46,  0,  0,  0,  0,  0,  0,  0,  0,  0,
     0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,
     0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,
     0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,
     0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,
     0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,
     0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0is shown in red:


This matches the STM32fftvsMATLABfft.PNG OP attached to reply #17 exactly:


I did overlap the two curves in Inkscape (my Gnuplot original is in SVG), and they do match precisely.

[turny]

That to me indicates that the values that syntax333 uses to call arm_cfft_f32 are 12, 0, 0, 0, 10, 0, 0, 0, 13, 0, 0, 0 etc (12+0i, 0+0i, 10+0i, 0+0i, 13+0i, 0+0i, ...) instead of 12, 0, 10, 0, 13, 0 etc (12+0i, 10+0i, 13+0i, ...).

[Nominal Animal]

That to me indicates that the values that syntax333 uses to call arm_cfft_f32 are 12, 0, 0, 0, 10, 0, 0, 0, 13, 0, 0, 0 etc (12+0i, 0+0i, 10+0i, 0+0i, 13+0i, 0+0i, ...) instead of 12, 0, 10, 0, 13, 0 etc (12+0i, 10+0i, 13+0i, ...).

I modified my wording to hopefully make it clear: by "samples", I mean logical samples as in complex numbers if doing a complex DFT, and real numbers if doing a real DFT; not elements of the array supplied to some specific arm*fft*() call.  (Remember, we want the discussion to encompass the Matlab or other-tool derived DFTs for comparison.)

(In the frequency domain, I prefer to use the term 'bin', as in 'frequency bin'.  These are always complex numbers, and only the first half are useful, since the latter half corresponds to negative frequencies and thus contain no additional information, just a mirrored copy.  The very first bin is the DC component, i.e. average of the input samples.)

Assuming OP is using github.com/ARM-software/CMSIS/blob/master/CMSIS/DSP_Lib/Source/TransformFunctions/arm_cfft_f32.c:arm_cfft_f32() (or an earlier version), then there are several issues in OP's approach.  One is the extra zero insertion, the other is using complex FFT when the samples are purely real.

As I showed above, OP's current ARM code uses a DFT over 512 complex samples, where every odd complex number is zero, 0.0+0.0i, and the last 256 complex numbers are also zero, 0.0+0.0i.  The first 128 even complex samples are the real values copied from the data OP supplied.
For the arm_cfft_f32(), the input array contents are indeed { 12.0f, 0.0f, 0.0f, 0.0f, 10.0f, 0.0f, 0.0f, 0.0f, 13.0f, 0.0f, 0.0f, 0.0f, ... }, because for arm_cfft_f32(), each consecutive pair of input f32's forms one complex sample.

The proper functions to use with real samples is arm_rfft_fast_f32() (which supersedes the older arm_rfft_f32() and related functions).  The input array to this is simply the float samples; the output array is an array of floats, with each bin described by two consecutive elements.  The first pair is the DC bin, with imaginary part always zero (because the DC component for real samples is always real).

Exactly where OP's error is, I cannot tell.  It could be that tempBuf_1[] already contains complex samples, i.e. all odd indexes (imaginary parts) are already zero, and the call to PreFFTProcess() expands that so that every other complex sample, i.e. every other f32 pair, supplied to arm_cfft_f32(), ends up being zero.  Or, it could be a bug in PreFFTProcess() itself.  Also, we haven't seen the contents of the arm_cfft_instance_f32 fftInstance; its fftLen is the number of samples, not the number of f32's, supplied to arm_cfft_f32(); here, it should be 256, not 512.