From bccaee48218c850e6d305e619a2fd38052b79981 Mon Sep 17 00:00:00 2001
From: dwarning <dwarning>
Date: Wed, 15 Feb 2023 15:24:55 +0100
Subject: [PATCH] correct fft dc scaling bug #620

---
 src/frontend/com_fft.c    |  6 ++++--
 src/maths/cmaths/cmath4.c | 14 ++++++++++----
 2 files changed, 14 insertions(+), 6 deletions(-)

diff --git a/src/frontend/com_fft.c b/src/frontend/com_fft.c
index 0dc42ddc6..2d592525b 100644
--- a/src/frontend/com_fft.c
+++ b/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;
diff --git a/src/maths/cmaths/cmath4.c b/src/maths/cmaths/cmath4.c
index b3dd42f9a..11d7bf8f7 100644
--- a/src/maths/cmaths/cmath4.c
+++ b/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;