Merge pull request #53 from wantysal/freq_correction

Frequency domain corrections for nth octave spectrum and screening for tones functions

Author: Martin Glesser

mosqito/sound_level_meter/__init__.py

@@ -1,3 +1,3 @@
-from mosqito.sound_level_meter.noct_spectrum.noct_spectrum import (
-    noct_spectrum,
-)
+from mosqito.sound_level_meter.noct_spectrum.noct_spectrum import noct_spectrum
+from mosqito.sound_level_meter.noct_spectrum.noct_synthesis import noct_synthesis
+

mosqito/sound_level_meter/noct_spectrum/_filter_bandwidth.py

@@ -30,6 +30,11 @@ def _filter_bandwidth(fc, n=3, N=3):
     alpha : numpy.ndarray
         Ratio of the upper and lower band-edge frequencies to the mid-band
         frequency
+    f1 : float
+        Lower frequency of the band in Hz
+    f2 : float
+        Upper frequency of the band in Hz
+
     """
 
     """
@@ -54,7 +59,7 @@ def _filter_bandwidth(fc, n=3, N=3):
     # Ratio of the upper and lower band-edge frequencies to the mid-band frequency
     alpha = (1 + np.sqrt(1 + 4 * qd ** 2)) / 2 / qd
 
-    return alpha
+    return alpha, f1, f2
 
 
 if __name__ == "__main__":

mosqito/sound_level_meter/noct_spectrum/_n_oct_freq_filter.py

@@ -1,15 +1,11 @@
 # -*- coding: utf-8 -*-
-"""
-Created on Thu Mar 10 10:00:03 2022
 
-@author: wantysal
-"""
+
 # Standard library imports
 import numpy as np
-from scipy.signal import butter, freqz
-
+from scipy.signal import butter, sosfreqz
 
-def _n_oct_freq_filter(spectrum, fs, fc, alpha, n=3):
+def _n_oct_freq_filter(spectrum, fs, fc, alpha, n=3):  
     """ n-th octave filtering in frequency domain
 
     Designs a digital 1/3-octave filter with center frequency fc for
@@ -26,7 +22,7 @@ def _n_oct_freq_filter(spectrum, fs, fc, alpha, n=3):
     alpha : float
         Ratio of the upper and lower band-edge frequencies to the mid-band
         frequency
-    N : int, optional
+    n : int, optional
         Filter order. Default to 3
 
     Outputs
@@ -35,40 +31,17 @@ def _n_oct_freq_filter(spectrum, fs, fc, alpha, n=3):
         Rms level of sig in the third octave band centered on fc
     """
 
-    # Check for Nyquist-Shannon criteria
-    if fc > 0.88 * (fs / 2):
-        raise ValueError(
-            """ERROR: Design not possible. Filter center frequency shall
-            verify: fc <= 0.88 * (fs / 2)"""
-        )
-
     # Normalized cutoff frequencies
-    w1 = fc / (fs / 2) / alpha
-    w2 = fc / (fs / 2) * alpha
-
-    # define filter coefficient
-    b, a = butter(n, [w1, w2], "bandpass", analog=False)
-    
+    w1 = fc / (fs  / 2) / alpha
+    w2 = fc / (fs  / 2) * alpha
+
+    # Define filter coefficient
+    sos = butter(n, [w1, w2], "bandpass", analog=False, output ='sos')  
+    # Get FRF and apply it
+    w, h = sosfreqz(sos, worN=len(spectrum))
+    spec_filt = np.multiply(h, spectrum.T).T
 
-    # filter signal
-    if len(spectrum.shape)>1:
-        level = []
-        # go into frequency domain
-        w, h = freqz(b, a, worN=spectrum.shape[1])
-        for i in range(spectrum.shape[0]):
-            spec_filt = spectrum[i,:] * h
-            # Compute overall rms level
-            level.append(np.sqrt(np.sum(np.abs(spec_filt) ** 2)))
-        level = np.array(level)
-    else:
-        # go into frequency domain
-        w, h = freqz(b, a, worN=len(spectrum))
-        spec_filt = spectrum * h
-        # Compute overall rms level
-        level = np.sqrt(np.sum(np.abs(spec_filt) ** 2)) 
+    # Compute overall rms level
+    level = np.sqrt(np.sum(np.abs(spec_filt) ** 2, axis=0)) 
 
     return level
-
-
-
-

mosqito/sound_level_meter/noct_spectrum/_n_oct_time_filter.py

@@ -2,7 +2,7 @@
 
 # Standard library imports
 import numpy as np
-from scipy.signal import decimate, butter, lfilter
+from scipy.signal import decimate, butter, sosfilt
 
 
 def _n_oct_time_filter(sig, fs, fc, alpha, N=3):
@@ -29,7 +29,7 @@ def _n_oct_time_filter(sig, fs, fc, alpha, N=3):
     alpha : float
         Ratio of the upper and lower band-edge frequencies to the mid-band
         frequency
-    N : int, optional
+    n : int, optional
         Filter order. Default to 3
 
     Outputs
@@ -59,10 +59,10 @@ def _n_oct_time_filter(sig, fs, fc, alpha, N=3):
     w2 = fc / (fs / 2) * alpha
 
     # define filter coefficient
-    b, a = butter(N, [w1, w2], "bandpass", analog=False)
+    sos = butter(int(N), (w1, w2), "bandpass", analog=False, output='sos')
 
     # filter signal
-    sig_filt = lfilter(b, a, sig, axis=0)
+    sig_filt = sosfilt(sos, sig, axis=0)
 
     # Compute overall rms level
     level = np.sqrt(np.mean(sig_filt ** 2, axis=0))

mosqito/sound_level_meter/noct_spectrum/noct_spectrum.py

@@ -53,7 +53,7 @@ def noct_spectrum(sig, fs, fmin, fmax, n=3, G=10, fr=1000):
     fc_vec, fpref = _center_freq(fmin=fmin, fmax=fmax, n=n, G=G, fr=fr)
 
     # Compute the filters bandwidth
-    alpha_vec = _filter_bandwidth(fc_vec, n=n)
+    alpha_vec, _, _ = _filter_bandwidth(fc_vec, n=n)
 
     # Calculation of the rms level of the signal in each band
     spec = []

mosqito/sound_level_meter/noct_spectrum/noct_synthesis.py

@@ -5,14 +5,14 @@
 
 # Local application imports
 from mosqito.sound_level_meter.noct_spectrum._center_freq import _center_freq
+from mosqito.sound_level_meter.noct_spectrum._filter_bandwidth import _filter_bandwidth
+from mosqito.sound_level_meter.noct_spectrum._n_oct_freq_filter import _n_oct_freq_filter
 
 
 def noct_synthesis(spectrum, freqs, fmin, fmax, n=3, G=10, fr=1000):
     """Adapt input spectrum to nth-octave band spectrum
-
     Convert the input spectrum to third-octave band spectrum
     between "fc_min" and "fc_max".
-
     Parameters
     ----------
     spectrum : numpy.ndarray
@@ -32,48 +32,47 @@ def noct_synthesis(spectrum, freqs, fmin, fmax, n=3, G=10, fr=1000):
         Reference frequency. Shall be set to 1 kHz for audible frequency
         range, to 1 Hz for infrasonic range (f < 20 Hz) and to 1 MHz for
         ultrasonic range (f > 31.5 kHz).
-
     Outputs
     -------
     spec : numpy.ndarray
         Third octave band spectrum of signal sig [dB re.2e-5 Pa], size (nbands, nseg).
     fpref : numpy.ndarray
         Corresponding preferred third octave band center frequencies, size (nbands).
     """
+    
+    # Deduce sampling frequency
+    fs = np.mean(freqs[1:] - freqs[:-1]) * 2 * len(spectrum)
 
     # Get filters center frequencies
     fc_vec, fpref = _center_freq(fmin=fmin, fmax=fmax, n=n, G=G, fr=fr)
+    
+    # Compute the filters bandwidth
+    alpha_vec, f_low, f_high = _filter_bandwidth(fc_vec, n=n)
+    
+    # Delete ends frequencies to prevent aliasing
+    idx = np.asarray(np.where(f_high > fs / 2))
+    if any(idx[0]):
+        fc_vec = np.delete(fc_vec, idx)
+        f_low = np.delete(f_low, idx)
+        f_high = np.delete(f_high, idx)
 
-    nband = len(fpref)
-
+    # Number of nth bands
+    nband = len(fc_vec)
+    
+    # Results array initialization
     if len(spectrum.shape) > 1:
         nseg = spectrum.shape[1]
         spec = np.zeros((nband, nseg))
+        # If only one axis is given, it is used for all the spectra
         if len(freqs.shape) == 1:
             freqs = np.tile(freqs, (nseg, 1)).T
-
     else:
         nseg = 1
         spec = np.zeros((nband))
 
-    # Frequency resolution
-    # df = freqs[1:] - freqs[:-1]
-    # df = np.concatenate((df, [df[-1]]))
-
-    # Get upper and lower frequencies
-    fu = fc_vec * 2**(1/(2*n))
-    fl = fc_vec / 2**(1/(2*n))
-
-    for s in range(nseg):
-        for i in range(nband):
-            if len(spectrum.shape) > 1:
-                # index of the frequencies within the band
-                idx = np.where((freqs[:, s] >= fl[i]) & (freqs[:, s] < fu[i]))
-                spec[i, s] = np.sqrt(
-                    np.sum(np.power(np.abs(spectrum[idx,s]), 2)))
-            else:
-                # index of the frequencies within the band
-                idx = np.where((freqs >= fl[i]) & (freqs < fu[i]))
-                spec[i] = np.sqrt(np.sum(np.abs(spectrum[idx])**2))
+    # Calculation of the rms level of the signal in each band
+    spec = []
+    for fc, alpha in zip(fc_vec, alpha_vec):
+        spec.append(_n_oct_freq_filter(spectrum, fs, fc, alpha))
 
-    return spec, fpref
+    return np.array(spec), fpref

mosqito/sq_metrics/tonality/prominence_ratio_ecma/_pr_main_calc.py

@@ -205,9 +205,7 @@ def _pr_main_calc(spectrum_db, freq_axis):
 
             # suppression from the list of tones of all the candidates belonging to the
             # same critical band
-            sup = np.where((peaks >= low_limit_idx) & (peaks <= high_limit_idx))[
-                0
-            ]
+            sup = np.where((peaks >= low_limit_idx) & (peaks <= high_limit_idx))[0]
             peaks = np.delete(peaks, sup)
             nb_tones -= len(sup)
 

mosqito/sq_metrics/tonality/tone_to_noise_ecma/_screening_for_tones.py

@@ -8,9 +8,7 @@
 import numpy as np
 
 # Mosqito functions import
-from mosqito.sq_metrics.tonality.tone_to_noise_ecma._spectrum_smoothing import (
-    _spectrum_smoothing,
-)
+from mosqito.sq_metrics.tonality.tone_to_noise_ecma._spectrum_smoothing import _spectrum_smoothing
 from mosqito.sq_metrics.tonality.tone_to_noise_ecma._LTH import _LTH
 from mosqito.sq_metrics.tonality.tone_to_noise_ecma._critical_band import _critical_band
 
@@ -54,7 +52,7 @@ def _screening_for_tones(freqs, spec_db, method, low_freq, high_freq):
     # Detection of the tonal candidates according to their level
 
     # Creation of the smoothed spectrum
-    smooth_spec = _spectrum_smoothing(freqs, spec_db, 24, low_freq, high_freq, freqs)
+    smooth_spec = _spectrum_smoothing(freqs, spec_db.T, 24, low_freq, high_freq, freqs).T
 
 
     n = spec_db.shape[0]
@@ -149,7 +147,7 @@ def _screening_for_tones(freqs, spec_db, method, low_freq, high_freq):
         # Screen the left points of the peak
         temp = peak_index - 1
         # As long as the level decreases,
-        while (spec_db[temp] > smooth_spec[temp] + 6) and (temp + 1 < (block+1)*m):
+        while (spec_db[temp] > smooth_spec[temp] + 6) and (temp +1 > (block)*m):
             # if a highest spectral line is found, it becomes the candidate
             if spec_db[temp] > spec_db[peak_index]:
                 peak_index = temp

mosqito/sq_metrics/tonality/tone_to_noise_ecma/_spectrum_smoothing.py

@@ -18,9 +18,9 @@ def _spectrum_smoothing(freqs_in, spec, noct, low_freq, high_freq, freqs_out):
     Parameters
     ----------
     freqs : numpy.array
-        frequency axis (frequency axis)
+        frequency axis dim (nperseg)
     spec : numpy.array
-        spectrum in dB (n block x spectrum)
+        spectrum in dB (nperseg, nseg)
     noct : integer
         n-th octave-band according to which smooth the spectrum
     low_freq : float
@@ -36,14 +36,13 @@ def _spectrum_smoothing(freqs_in, spec, noct, low_freq, high_freq, freqs_out):
         smoothed spectrum along the given frequency axis
 
     """
-    n = spec.shape[0]
+    nperseg = spec.shape[0]
     if len(spec.shape) > 1:
-        m = spec.shape[1]
+        nseg = spec.shape[1]
     else:
-        m = spec.shape[0]
-        n = 1
+        nseg = 1
     spec = spec.ravel()
-    stop = np.arange(1, n + 1) * m
+    stop = np.arange(1, nseg + 1) * nperseg
     freqs_in = freqs_in.ravel()
 
     # n-th octave bands filter
@@ -53,14 +52,12 @@ def _spectrum_smoothing(freqs_in, spec, noct, low_freq, high_freq, freqs_out):
 
     # Smoothed spectrum creation
     nb_bands = filter_freqs.shape[0]
-    smoothed_spectrum = np.zeros((n, nb_bands))
+    smoothed_spectrum = np.zeros((nb_bands, nseg))
     i = 0
     # Each band is considered individually until all of them have been treated
     while nb_bands > 0:
         # Find the index of the spectral components within the frequency bin
-        bin_index = np.where(
-            (freqs_in >= filter_freqs[i, 0]) & (freqs_in <= filter_freqs[i, 2])
-        )[0]
+        bin_index = np.where((freqs_in >= filter_freqs[i, 0]) & (freqs_in <= filter_freqs[i, 2]))[0]
 
         # If the frequency bin is empty, it is deleted from the list
         if len(bin_index) == 0:
@@ -70,45 +67,32 @@ def _spectrum_smoothing(freqs_in, spec, noct, low_freq, high_freq, freqs_out):
 
         else:
             # The spectral components within the frequency bin are averaged on an energy basis
-            spec_sum = np.zeros((n))
-            for j in range(n):
-                if (
-                    len(bin_index[(bin_index < stop[j]) & (bin_index > (stop[j] - m))])
-                    != 0
-                ):
-                    spec_sum[j] = np.mean(
-                        10
-                        ** (
-                            spec[
-                                bin_index[
-                                    (bin_index < stop[j]) & (bin_index > (stop[j] - m))
-                                ]
-                            ]
-                            / 10
-                        )
-                    )
+            spec_sum = np.zeros((nseg))
+            for j in range(nseg):
+                if len(bin_index[(bin_index < stop[j]) & (bin_index > (stop[j] - nperseg))])!= 0:
+                    spec_sum[j] = np.mean(10** (spec[bin_index[(bin_index < stop[j]) & (bin_index > (stop[j] - nperseg))]]/ 10))
                 else:
                     spec_sum[j] = 1e-12
-            smoothed_spectrum[:, i] = 10 * np.log10(
-                spec_sum
-                # / len(bin_index[(bin_index < stop[j]) & (bin_index > (stop[j] - m))])
-            )
+            smoothed_spectrum[i, :] = 10 * np.log10(spec_sum)# / len(bin_index[(bin_index < stop[j]) & (bin_index > (stop[j] - m))]))
         nb_bands -= 1
         i += 1
     # Pose of the smoothed spectrum on the frequency-axis
-    low = np.zeros((n, filter_freqs.shape[0]))
-    high = np.zeros((n, filter_freqs.shape[0]))
+    low = np.zeros((filter_freqs.shape[0], nseg))
+    high = np.zeros((filter_freqs.shape[0], nseg))
 
     # Index of the lower and higher limit of each frequency bin into the original spectrum
     for i in range(len(filter_freqs)):
-        low[:, i] = np.argmin(np.abs(freqs_out - filter_freqs[i, 0]))
-        high[:, i] = np.argmin(np.abs(freqs_out - filter_freqs[i, 2]))
+        low[i,:] = np.argmin(np.abs(freqs_out - filter_freqs[i, 0]))
+        high[i,:] = np.argmin(np.abs(freqs_out - filter_freqs[i, 2]))
     low = low.astype(int)
     high = high.astype(int)
 
-    smooth_spec = np.zeros((n, m))
-    for i in range(n):
+    smooth_spec = np.zeros((nperseg, nseg))
+    for i in range(nseg):
         for j in range(filter_freqs.shape[0]):
-            smooth_spec[i, low[i, j] : high[i, j]] = smoothed_spectrum[i, j]
+            smooth_spec[low[j,i] : high[j,i], i] = smoothed_spectrum[j,i]
+    
+    if smooth_spec.shape[1] == 1:
+        smooth_spec = np.squeeze(smooth_spec)
 
     return smooth_spec