From c98b10b83a361d7fa4c9537888390fe0694b071b Mon Sep 17 00:00:00 2001
From: Casper <c.jansen@ascee.nl>
Date: Fri, 13 Aug 2021 11:40:25 +0200
Subject: [PATCH] Fixed scaling issue in FFT

---
 lasp/c/lasp_ps.c | 29 ++++++++++++++++-------------
 1 file changed, 16 insertions(+), 13 deletions(-)

diff --git a/lasp/c/lasp_ps.c b/lasp/c/lasp_ps.c
index bc0a244..1d6018b 100644
--- a/lasp/c/lasp_ps.c
+++ b/lasp/c/lasp_ps.c
@@ -1,7 +1,7 @@
 // lasp_ps.c
 //
 // Author: J.A. de Jong -ASCEE
-// 
+//
 // Description:
 //
 //////////////////////////////////////////////////////////////////////
@@ -21,7 +21,7 @@ typedef struct PowerSpectra_s {
 
 PowerSpectra* PowerSpectra_alloc(const us nfft,
                                  const WindowType wt) {
-    
+
     fsTRACE(15);
     int rv;
 
@@ -48,7 +48,7 @@ PowerSpectra* PowerSpectra_alloc(const us nfft,
 
     /* Allocate vectors and matrices */
     ps->window = vd_alloc(nfft);
-    
+
     rv = window_create(wt,&ps->window,&ps->win_pow);
     check_overflow_vx(ps->window);
     if(rv!=0) {
@@ -70,11 +70,11 @@ void PowerSpectra_free(PowerSpectra* ps) {
 void PowerSpectra_compute(const PowerSpectra* ps,
                          const dmat * timedata,
                          cmat * result) {
-    
+
     fsTRACE(15);
 
     dbgassert(ps && timedata && result,NULLPTRDEREF);
-    
+
     const us nchannels = timedata->n_cols;
     const us nfft = Fft_nfft(ps->fft);
     uVARTRACE(15,nchannels);
@@ -100,13 +100,13 @@ void PowerSpectra_compute(const PowerSpectra* ps,
         vd column_work = dmat_column(&timedata_work,i);
 
         vd_elem_prod(&column_work,&column,&ps->window);
-        
+
         vd_free(&column);
         vd_free(&column_work);
     }
     check_overflow_xmat(timedata_work);
     /* print_dmat(&timedata_work); */
-    
+
     /* Compute fft of the time data */
     cmat fft_work = cmat_alloc(nfft/2+1,nchannels);
     Fft_fft(ps->fft,
@@ -114,19 +114,22 @@ void PowerSpectra_compute(const PowerSpectra* ps,
             &fft_work);
 
     dmat_free(&timedata_work);
-    
+
     TRACE(15,"fft done");
-    
-    /* Scale fft such that power is easily computed */
-    const c scale_fac = d_sqrt(2/win_pow)/nfft;
+
+    /* Scale fft such that power is easily computed
+     * - Multiply power spectral density by 2 except at f=0 and f=fNq
+     * - Divide by energy of window function = nfft * window_power
+     * - .. sqrt(factors) because it is applied to output fft instead of psd */
+    const c scale_fac = d_sqrt(2/(nfft*win_pow));
     cmat_scale(&fft_work,scale_fac);
     TRACE(15,"scale done");
 
     for(us i=0;i< nchannels;i++) {
-        /* Multiply DC term with 1/sqrt(2) */
+        /* Multiply DC term by 1/sqrt(2) */
         *getcmatval(&fft_work,0,i) *= 1/d_sqrt(2.)+0*I;
 
-        /* Multiply Nyquist term with 1/sqrt(2) */
+        /* Multiply Nyquist term by 1/sqrt(2) */
         *getcmatval(&fft_work,nfft/2,i) *= 1/d_sqrt(2.)+0*I;
     }
     check_overflow_xmat(fft_work);