make fft scaling independent from rounding behaviour for odd data

This commit is contained in:
dwarning 2026-03-05 09:30:55 +01:00
parent b522bcc6c0
commit a31ac54200
2 changed files with 4 additions and 4 deletions

View File

@ -63,6 +63,7 @@ com_fft(wordlist *wl)
#ifdef HAVE_LIBFFTW3 #ifdef HAVE_LIBFFTW3
fpts = length/2 + 1; fpts = length/2 + 1;
scale = ((double)length)/2.0;
#else #else
/* size of fft input vector is power of two and larger or equal than spice vector */ /* size of fft input vector is power of two and larger or equal than spice vector */
N = 1; N = 1;
@ -72,6 +73,7 @@ com_fft(wordlist *wl)
M++; M++;
} }
fpts = N/2 + 1; fpts = N/2 + 1;
scale = ((double)N)/2;
#endif #endif
win = TMALLOC(double, length); win = TMALLOC(double, length);
@ -178,7 +180,6 @@ com_fft(wordlist *wl)
fftw_execute(plan_forward); fftw_execute(plan_forward);
scale = (double) fpts - 1.0;
fdvec[i][0].cx_real = out[0][0]/scale/2.0; fdvec[i][0].cx_real = out[0][0]/scale/2.0;
fdvec[i][0].cx_imag = 0.0; fdvec[i][0].cx_imag = 0.0;
for (j = 1; j < fpts; j++) { for (j = 1; j < fpts; j++) {
@ -215,7 +216,6 @@ com_fft(wordlist *wl)
rffts(in, M, 1); rffts(in, M, 1);
fftFree(); fftFree();
scale = (double) fpts - 1.0;
/* Re(x[0]), Re(x[N/2]), Re(x[1]), Im(x[1]), Re(x[2]), Im(x[2]), ... Re(x[N/2-1]), Im(x[N/2-1]). */ /* Re(x[0]), Re(x[N/2]), Re(x[1]), Im(x[1]), Re(x[2]), Im(x[2]), ... Re(x[N/2-1]), Im(x[N/2-1]). */
fdvec[i][0].cx_real = in[0]/scale/2.0; fdvec[i][0].cx_real = in[0]/scale/2.0;
fdvec[i][0].cx_imag = 0.0; fdvec[i][0].cx_imag = 0.0;

View File

@ -774,7 +774,7 @@ cx_fft(void *data, short int type, int length, int *newlength, short int *newtyp
fftw_execute(plan_forward); fftw_execute(plan_forward);
scale = (double) fpts - 1.0; scale = ((double)length)/2.0;
outdata[0].cx_real = out[0][0]/scale/2.0; outdata[0].cx_real = out[0][0]/scale/2.0;
outdata[0].cx_imag = 0.0; outdata[0].cx_imag = 0.0;
for (i = 1; i < fpts; i++) { for (i = 1; i < fpts; i++) {
@ -803,7 +803,7 @@ cx_fft(void *data, short int type, int length, int *newlength, short int *newtyp
rffts(datax, M, 1); rffts(datax, M, 1);
fftFree(); fftFree();
scale = (double) fpts - 1.0; scale = ((double)N)/2;
/* Re(x[0]), Re(x[N/2]), Re(x[1]), Im(x[1]), Re(x[2]), Im(x[2]), ... Re(x[N/2-1]), Im(x[N/2-1]). */ /* Re(x[0]), Re(x[N/2]), Re(x[1]), Im(x[1]), Re(x[2]), Im(x[2]), ... Re(x[N/2-1]), Im(x[N/2-1]). */
outdata[0].cx_real = datax[0]/scale/2.0; outdata[0].cx_real = datax[0]/scale/2.0;
outdata[0].cx_imag = 0.0; outdata[0].cx_imag = 0.0;