From 9c8c546566112ebb73d42c3165361125a734a64b Mon Sep 17 00:00:00 2001
From: Casper <c.jansen@ascee.nl>
Date: Fri, 22 Oct 2021 13:36:23 +0200
Subject: [PATCH] Restored FFT scaling to power/bin instead of PSD

---
 lasp/c/lasp_ps.c    | 8 ++++++--
 lasp/tools/tools.py | 5 +++--
 2 files changed, 9 insertions(+), 4 deletions(-)

diff --git a/lasp/c/lasp_ps.c b/lasp/c/lasp_ps.c
index 1d6018b..c8e664d 100644
--- a/lasp/c/lasp_ps.c
+++ b/lasp/c/lasp_ps.c
@@ -117,11 +117,15 @@ void PowerSpectra_compute(const PowerSpectra* ps,
 
     TRACE(15,"fft done");
 
-    /* Scale fft such that power is easily computed
+    /* Scale fft such that power is easily comxputed */
+    const c scale_fac = d_sqrt(2/win_pow)/nfft;
+
+    /* POWER SPECTRAL DENSITY? 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));
+    // const c scale_fac = d_sqrt(2/(nfft*win_pow));
+
     cmat_scale(&fft_work,scale_fac);
     TRACE(15,"scale done");
 
diff --git a/lasp/tools/tools.py b/lasp/tools/tools.py
index e9bbf15..befb0fb 100644
--- a/lasp/tools/tools.py
+++ b/lasp/tools/tools.py
@@ -48,6 +48,7 @@ class SmoothingWidth(Enum):
     twelve = (12, '1/12th octave smoothing')
     twfo = (24, '1/24th octave smoothing')
     ftei = (48, '1/48th octave smoothing')
+#    hundred = (100, '1/100th octave smoothing')  # useful for removing 'grass'
 
     @staticmethod
     def fillComboBox(cb):
@@ -136,7 +137,6 @@ def smoothSpectralData(freq, M, sw: SmoothingWidth,
     x0 = 1 if freq[0] == 0 else 0  # Skip first data point if zero frequency
     Psm[0] = P[0]                  # Reuse old value in case first data..
                                    # ..point is skipped. Not plotted any way.
-    P[-1] = P[-2]                  # Last data point (Nyquist f) can be invalid
 
     # Loop through data points
     for x in range(x0, L):
@@ -162,6 +162,7 @@ def smoothSpectralData(freq, M, sw: SmoothingWidth,
             fu = fNq  # no data beyond fNq
             fl = fc**2 / fu  # keep window symmetric
 
+        # Find indices corresponding to frequencies
         xl = bisect.bisect_left(freq, fl)  # index corresponding to fl
         xu = bisect.bisect_left(freq, fu)
 
@@ -178,7 +179,7 @@ def smoothSpectralData(freq, M, sw: SmoothingWidth,
 
         # Apply smoothing
         Psm[x] = np.dot(g, P[xl:xu])
-    
+
     if st == SmoothingType.levels:
         Psm = 10*np.log10(Psm)