espressif / esp-dsp

DSP library for ESP-IDF
Apache License 2.0
474 stars 88 forks source link

FFT is wrong by factor 2. Violation of the Parseval–Plancherel identity (DSP-130) #81

Open Bobobel opened 8 months ago

Bobobel commented 8 months ago

Answers checklist.

General issue report

I know that the prefactor of FFT/DFT is under heavy discussion. But for me the Parseval–Plancherel identity has to be obeyed! Please check your example dsps_fft_main.c without windowing and you will see, that the absolute magnitudes of x1 and x2 are not represented in the frequency domain but with a factor of 2. Their relative value is OK, so you did not notice a problem. This was hidden because you used the windowing facility (thanks for showing that anyway!). Thanks for that DSP component!

dmitry1945 commented 8 months ago

Dear @Bobobel,

thank you for the issue.

We will change the documentation (or the code). The factor of 2 is related to calculation performance.

Regards, Dmitry

Bobobel commented 8 months ago

Dear dmitry1945,

I am not shure if my last mail arrived this way. I have a more serious problem with fft_main_sc16 (see my version in last mail from Feb. 3rd: When I scale sin signals up to an amplitude of 32767 everything worked as expected. When I scale a rectangluar signal (puls 0.5 and amplitude 1) upto 25000 : results OK, with little deviations But when I scale higher, e.g. 28000 or 32767 (max), then I have got lower results. Here is an ouput of two rectangular signals: x1 scaled to 16383 and x2 to 32766 (frequencies 1/16 resp. 1/32,amplitudes 1 and 2 with a scaling factor of 16383):

I (1577) main: FFT for 256 complex points took 23239 cycles I (1587) main: End Example. Start example with N=256, ifactor=16383.      i,    x1,      x2,      R(y1),   I(y1),   R(y2),   I(y2), y1p,   y2p       0,       0,  -32766,      57,       0,    -260, 0,-49.17037,-35.98840       1,   16383,  -32766,     120,     -51,    -503, -100,-41.98312,-30.08816       2,   16383,  -32766,     107,     -98,    -474, -195,-41.05473,-30.09325       3,   16383,  -32766,      83,    -142,    -426, -283,-39.96571,-30.11207       4,   16383,  -32766,      51,    -181,    -363, -361,-38.80250,-30.10336       5,   16383,  -32766,      13,    -213,    -285, -425,-37.70413,-30.10734       6,   16383,  -32766,     -32,    -236,    -198, -472,-36.75051,-30.10505       7,   16383,  -32766,     -80,    -250,    -102, -502,-35.90568,-30.09810       8,   16383,  -32766,    -131,    -257,   23959, 3577,-35.08611, 3.39724       9,  -16383,   32766,    -180,    -250,      98, -502,-34.51521,-30.11136      10,  -16383,   32766,    -228,    -236,     194, -472,-33.96652,-30.13113      11,  -16383,   32766,    -271,    -213,     281, -425,-33.53934,-30.14506      12,  -16383,   32766,    -311,    -181,     359, -361,-33.16578,-30.15148      13,  -16383,   32766,    -342,    -143,     423, -284,-32.90759,-30.14507      14,  -16383,   32766,    -366,     -98,     470, -196,-32.71753,-30.14958      15,  -16383,   32766,    -380,     -50,     500, -100,-32.61765,-30.13813      16,   16383,   32766,    3707,  -20586,     510,      -1, 2.12216,-30.13645      17,   16383,   32766,    -380,      50,     500, 100,-32.61765,-30.13813

As you see for i=8 in signal 2 R(y2)=23959 I(y2)=3577 while for i=16: R(y1)=3707, I(y1)=-20585 (phase was -90° with x2, so we have switched R and I) Signal never went beyond +-32767 but the spectral components display an overflow:

With a scaling factor of 8191 I have got R(y2)=-20783, I(y2)=1785 and therefor a doubling gives R=-41566  and overflow +2^16 = 23970 (almost the number above, deviation due to 8191 and 16383 not doubling).

As you see: the problem willnot occur with simpletests using sin, but with real complex signals there will be! Is there a chance to handle overflow and "underrun" in your assembler codings? Otherwise you should tell people in the docs that there is a severe problem because aou have missed the "factor 2 because of efficiency" as you explained before.

With other tests I had an irregular phase change in spectrum in complex summed signals. With the above foundings I have an explanation now as well.

Sorry for my insisting on that subject, but I plan a frequency analysis for real signals in sc16, 16 bit, and with that I have investigated these problems.

Thank you for your work and kind regards from

Jürgen/Bobobel