correct fft dc scaling bug #620

src/frontend/com_fft.c

@@ -179,7 +179,9 @@ com_fft(wordlist *wl)
         fftw_execute(plan_forward);
 
         scale = (double) fpts - 1.0;
-        for (j = 0; j < fpts; j++) {
+        fdvec[i][0].cx_real = out[0][0]/scale/2.0;
+        fdvec[i][0].cx_imag = 0.0;
+        for (j = 1; j < fpts; j++) {
             fdvec[i][j].cx_real = out[j][0]/scale;
             fdvec[i][j].cx_imag = out[j][1]/scale;
         }
@@ -212,7 +214,7 @@ com_fft(wordlist *wl)
 
         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]). */
-        fdvec[i][0].cx_real = in[0]/scale;
+        fdvec[i][0].cx_real = in[0]/scale/2.0;
         fdvec[i][0].cx_imag = 0.0;
         for (j = 1; j < fpts-1; j++) {
             fdvec[i][j].cx_real = in[2*j]/scale;

src/maths/cmaths/cmath4.c

@@ -713,7 +713,9 @@ cx_fft(void *data, short int type, int length, int *newlength, short int *newtyp
         outdata = alloc_c(fpts);
 
         scale = (double) fpts;
-        for (i = 0; i < fpts; i++) {
+        outdata[0].cx_real = out[0][0]/scale/2.0;
+        outdata[0].cx_imag = out[0][1]/scale/2.0;
+        for (i = 1; i < fpts; i++) {
             outdata[i].cx_real = out[i][0]/scale;
             outdata[i].cx_imag = out[i][1]/scale;
         }
@@ -744,7 +746,9 @@ cx_fft(void *data, short int type, int length, int *newlength, short int *newtyp
         outdata = alloc_c(N);
 
         scale = (double) N;
-        for (i = 0; i < N; i++) {
+        outdata[0].cx_real = datax[0]/scale/2.0;
+        outdata[0].cx_imag = datax[1]/scale/2.0;
+        for (i = 1; i < N; i++) {
             outdata[i].cx_real = datax[2*i]/scale;
             outdata[i].cx_imag = datax[2*i+1]/scale;
         }
@@ -773,7 +777,9 @@ cx_fft(void *data, short int type, int length, int *newlength, short int *newtyp
         fftw_execute(plan_forward);
 
         scale = (double) fpts - 1.0;
-        for (i = 0; i < fpts; i++) {
+        outdata[0].cx_real = out[0][0]/scale/2.0;
+        outdata[0].cx_imag = 0.0;
+        for (i = 1; i < fpts; i++) {
             outdata[i].cx_real = out[i][0]/scale;
             outdata[i].cx_imag = out[i][1]/scale;
         }
@@ -798,7 +804,7 @@ cx_fft(void *data, short int type, int length, int *newlength, short int *newtyp
 
         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]). */
-        outdata[0].cx_real = datax[0]/scale;
+        outdata[0].cx_real = datax[0]/scale/2.0;
         outdata[0].cx_imag = 0.0;
         for (i = 1; i < fpts-1; i++) {
             outdata[i].cx_real = datax[2*i]/scale;