From 59f82ae14cb9ee54acf7d7e506ce36171d487de4 Mon Sep 17 00:00:00 2001
From: "J.A. de Jong - ASCEE" <j.a.dejong@ascee.nl>
Date: Sun, 27 Oct 2019 14:54:28 +0100
Subject: [PATCH] Added tools to bin narrow-band data into frequency bands

---
 lasp/tools/bin_narrow.py | 86 ++++++++++++++++++++++++++++++++++++++++
 1 file changed, 86 insertions(+)
 create mode 100644 lasp/tools/bin_narrow.py

diff --git a/lasp/tools/bin_narrow.py b/lasp/tools/bin_narrow.py
new file mode 100644
index 0000000..e98d6ff
--- /dev/null
+++ b/lasp/tools/bin_narrow.py
@@ -0,0 +1,86 @@
+#!/usr/bin/python
+"""
+Bin narrow band power in octave/third octave band data
+"""
+from lasp.filter.bandpass_fir import (OctaveBankDesigner, 
+                                      ThirdOctaveBankDesigner)
+import numpy as np
+import warnings
+
+
+def binPower(freq, narrow_power, band=3, start_band=-16, stop_band=13):
+    """
+    Apply binning to narrow band frequency domain power results
+    
+    Args:
+        freq: Array of frequency indices
+        narrow_power: narrow-band power values in units of [W] or [Pa^2]
+        band: 1, or 3
+    
+    Returns:
+        ( ['25', '31.5', '40', '50', ... ],
+          [float(power_25), float(power_31p5), ...]) putting NAN values where
+           inaccurate.
+    """
+
+    if band == 3:
+        designer = ThirdOctaveBankDesigner()
+    elif band == 1:
+        designer = OctaveBankDesigner()
+    else:
+        raise ValueError("Parameter 'Band' should be either '1', or '3'")
+
+    freq = np.copy(freq)
+    narrow_power = np.copy(narrow_power)
+
+    
+    # Exact midband, lower and upper frequency limit of each band
+    fm = [designer.fm(x) for x in range(start_band, stop_band+1)]
+    fl = [designer.fl(x) for x in range(start_band, stop_band+1)]
+    fu = [designer.fu(x) for x in range(start_band, stop_band+1)]
+    fex = [designer.nominal_txt(x) for x in range(start_band, stop_band+1)]
+#    print(fl)
+    binned_power = np.zeros(len(fm), dtype=float)
+    ## Start: linear interpolation between bins while Parseval is conserved
+
+    # current frequency resolution
+    df_old = freq[1]-freq[0]
+    
+    # preferred new frequency resolution
+    df_new = .1
+
+    # ratio of resolutions
+
+    ratio = int(df_old/df_new)
+    # calculate new frequency bins
+    freq_new = np.linspace(freq[0],freq[-1],(len(freq)-1)*ratio+1)         
+    
+    # calculate the new bin data
+    interp_power = np.interp(freq_new, freq, narrow_power)/ratio
+    
+    # adapt first and last bin values so that Parseval still holds
+    interp_power[0] = binned_power[0]*(1.+1./ratio)/2.
+    interp_power[-1] = binned_power[-1]*(1.+1./ratio)/2.    
+
+    # check if Parseval still holds        
+    # print(np.sum(y, axis=0))
+    # print(np.sum(y_new, axis=0))
+    
+    ## Stop: linear interpolation between bins while Parseval is conserved        
+    binned_power = np.zeros(len(fm), dtype=float)
+    for k in range(len(fm)):
+#        print(k)
+        # find the bins which are in the corresponding band        
+        bins = (fl[k] <= freq_new) & (freq_new < fu[k])
+#        print(bins)
+        # sum the output values of these bins to obtain the band value
+        binned_power[k] = np.sum(interp_power[bins], axis=0)
+        # if no frequency bin falls in a certain band, skip previous bands 
+#        if not any(bins):
+#            binned_power[0:k+1] = np.nan
+
+    # check if data is valid
+    if(np.isnan(binned_power).all()):
+        warnings.warn('Invalid frequency array, we cannot bin these values')
+
+    return fm, fex, binned_power
\ No newline at end of file