ARM-software / CMSIS_5

CMSIS Version 5 Development Repository
http://arm-software.github.io/CMSIS_5/index.html
Apache License 2.0
1.33k stars 1.08k forks source link

issue with arm_rfft_q31 #219

Closed escherstair closed 4 years ago

escherstair commented 7 years ago

I have a buffer of 1024 samples in q31_t format and I need to do a FFT calculation, so I tried arm_rfft_q31() but the results are far from waht I expected to have. Probably either I missed something or I misunderstood the documentation, but I can't fin an explanation of what I see. This is my snippet of code:

#include "arm_math.h"

q31_t input[1024];
q31_t output[1024];
arm_rfft_instance_q31   S;
arm_status stat;

/* fill the input[] buffer with two cycles of a sine wave */
for(int i=0; i<512; i++) {
    input[i] = arm_sin_q31 (0x00400000*i);
    input[i+512] = input[i];
}

/* init the instance structure */
stat = arm_rfft_init_q31 (&S, 1024, 0, 0);

/* calculate FFT of the input[] buffer
 * the expected ideal output is the whole buffer set to 0, except for real[2] and/or imag[2]
 * because the input[] buffer has been filled in with exactly two cycles of sine wave */
arm_rfft_q31 (&S, input, output);

I found two strange behaviors:

Could you clarify this, please?

escherstair commented 7 years ago

I suspect the problem in the output buffer has to do with bitReverseFlag input parameter to arm_rfft_init_q31. Setting it to 1 seems to produce a better output[] buffer. input[] buffer is modified, anyway. I'm going to further investigate.

CCyros commented 7 years ago

Your right, bitReverseFlag = 1 means that normal bit order is used.

In this video the strange ordering of the result data is explained. The guy from the video also mentioned that the input data buffer is altered.

https://www.youtube.com/watch?v=gIh0pPsyGpI

escherstair commented 7 years ago

Hi @CCyros thank you for your link. It explains the meaning of bitReverseFlag input parameters and it points out that input buffer is altered.

I suggest to add to the documentation a proper advice to the user. As an example:

arm_fft_q31 performs some in-place calculations on the input buffer that is altered. Make sure you don't need anymore the data in this buffer after the execution of arm_fft_q31

One more thig: probably arm_fft_q31 can operate in-place. If this is true, I suggest to add this information to the documentation too.

JonatanAntoni commented 7 years ago

Hi @escherstair,

Thank you for pointing out possible improvements in detail. Please feel free to contribute using pull requests as well.

Cheers, Jonatan

escherstair commented 7 years ago

I've just found that documentation for arm_rfft_q31 needs improvement on another side too: the length for the several buffers involved is not clear enough.

The FFT of a real N-point sequence has even symmetry in the frequency domain. The second half of the data equals the conjugate of the first half flipped in frequency. Looking at the data, we see that we can uniquely represent the FFT using only N/2 complex numbers. These are packed into the output array in alternating real and imaginary components:

X = { real[0], imag[0], real[1], imag[1], real[2], imag[2] ... real[(N/2)-1], imag[(N/2)-1 }

and so it seems that N/2 complex numbers are enough. If this is true, the destination buffer can have a length of N q31_t numbers. But it seems that arm_rfft_q31 writes data to a region of memory with a length of 2*N q31_t elements.

Could you clarify these requirements for the buffers in the documentation, please?

xizhizhang commented 7 years ago

@escherstair you are right, arm_rfft_q31 need 2*N q31_t length destination buffer, the code is here:

https://github.com/ARM-software/CMSIS_5/blob/develop/CMSIS/DSP/Source/TransformFunctions/arm_rfft_q31.c#L137

it access 2*N-1 q31_t length address.

By the way, accurately, the length of a rfft result is N+2, for example, if you have 1024 samples inputs, you will get a 1026 length output.

xizhizhang commented 7 years ago

I think we may delete line 191 and 192 in arm_rfft_q31.c

https://github.com/ARM-software/CMSIS_5/blob/develop/CMSIS/DSP/Source/TransformFunctions/arm_rfft_q31.c#L191

because, if someone is using rfft, they may never interesting those results. It will help to save destination buffer size, and runs a little faster.

escherstair commented 7 years ago

@xizhizhang

I think we may delete line 191 and 192 in arm_rfft_q31.c

I'm not 100% sure this a safe approach. rfft functions work with real data sequences as input, but produces complex sequences as output. If I need to do some kind of operation on the output, then perform the inverse-FFT to get a waveform back, all the complex components are necessary.

xizhizhang commented 7 years ago

@escherstair in the arm_rfft_q31, it direct calls arm_cfft_q31, and direct puts all real input in cfft. the output of cfft is also N length. Then arm_split_rfft_q31 transforms N length of cfft output to 2*N length rfft output.

we can also find, in arm_rfft_q31 inverse-FFT operation, calls arm_split_rifft_q31, but arm_split_rifft_q31 only use N length input, and then puts N length input in inverse cfft function.

so, I think if we only want get back the waveform, it is OK to only use N length input.

christophe0606 commented 4 years ago

It is now covered in issue https://github.com/ARM-software/CMSIS_5/issues/810