DIVSCALAR off by 1

[fbarchard]

The current fixed point DIVSCALAR multiplies by SAMP_MAX which is off by 1

#   define DIVSCALAR(x,k) \
    (x) = sround( smul(  x, SAMP_MAX/k ) )

which produces 0.999938965 / k instead of 1 / k. The value is rounded down, so the larger the k value, the higher the error

it should be

#   define DIVSCALAR(x,k) \
    (x) = sround( smul(  x, (1<<(FRACBITS))/k ) )

It is used in bfly functions like this:

C_FIXDIV(*Fout,2); C_FIXDIV(*Fout2,2);

These are all the instances

FIXDIV(c,div) \
FIXDIV( fk , 2 );
FIXDIV( fnkc , 2 );
FIXDIV( *Fout0,5); C_FIXDIV( *Fout1,5); C_FIXDIV( *Fout2,5); C_FIXDIV( *Fout3,5); C_FIXDIV( *Fout4,5);
FIXDIV(*Fout,2); C_FIXDIV(*Fout2,2);
FIXDIV(*Fout,3); C_FIXDIV(Fout[m],3); C_FIXDIV(Fout[m2],3);
FIXDIV(*Fout,4); C_FIXDIV(Fout[m],4); C_FIXDIV(Fout[m2],4); C_FIXDIV(Fout[m3],4);
FIXDIV(fpk,2);
FIXDIV(fpnk,2);
FIXDIV(scratchbuf[q1],p);
FIXDIV(scratch[q1],p);
FIXDIV(st->tmpbuf[0],2);
FIXDIV(tdc,2);

So 4 constants - 2,3,4,5 and p which is a factor, typically a prime number larger than 5. FIXDIV(Fout,2); in bfly2 becomes (v 16387 + 16384) / 32768 which is accurate for values up to 16384, but off by 1 on larger values e.g. FIXDIV(20000,2) = (20000 × 16383 + 16384) ÷ 32768 = 9999.889648438 which truncates to 9999. With k = 4 32767/4 = 8191. e.g. FIXDIV(20000,4) = (20000 × 8191 + 16384) ÷ 32768 = 4999.889648437 which truncates to 4999. it should be (20000 × 8192 + 16384) ÷ 32768 = 5000.5 which truncates to 5000.

[fbarchard]

To work with FIXED_POINT==32 as well as FIXED_POINT==16 this is the change

==== kissfft/_kiss_fft_guts.h#1 - kissfft/_kiss_fft_guts.h ====
74c74
<   (x) = sround( smul(  x, SAMP_MAX/k ) )
---
>   (x) = sround( smul(  x, (SAMPPROD)(1u<<FRACBITS)/k ) )

[JulienMaille]

What about submitting a PR?