Complete tool for detection of prominent tones according to ISO1996-2-K + wiki

.gitignore

@@ -3,3 +3,4 @@ __pycache__/
 .vscode/
 .ipynb_checkpoints/
 tests/output
+prueba

docs/readme.md

@@ -18,3 +18,4 @@ In this folder, you will find the documentation and validation of the different
 - [Sharpness according to DIN standard](./sharpness_din.md)
 - [Roughness (Daniel and Weber method)](./roughness_dw.md)
 - [Tone-to-noise ratio and Prominence ratio of stationary and time-varying signals](./tone-to-noise-prominence-ratio.md)
+- [Third octave band tonality according to ISO 1996 annex K](./tonality_ISO1996K.md)

docs/tonality_ISO1996K.md

@@ -0,0 +1,16 @@
+# MOSQITO Documentation
+## Tone detection according to ISO 1996-2 Annex K
+
+### Introduction
+ISO 1996-2 introduces different methods for determining the presence of a prominent discrete frequency spectral component (tone).  
+
+In particular, Annex K describes the calculations for applying the "Objective method for the assessment of the audibility of tones in noise", also known as the "Inspection method".  
+
+These calculations are based on the comparison of time-averaged sound pressure levels in one-third octave bands taking into account the frequency range in which these levels are found.
+
+### Validation
+The validation of the function has been performed by testing the function to perform tone detection using different signals: pure tones at low, medium and high frequencies, multi-tone signals and atonal signals.
+By means of manual calculations, it is verified that the calculations obtained with the function are correct.
+
+### References
+    ISO 1996. Acoustics. Description, measurement and assessment of environmental noise. Part 2: Determination of sound pressure levels. (October 2020 version)

mosqito/sound_level_meter/noct_spectrum/noct_spectrum.py

@@ -40,7 +40,15 @@ def noct_spectrum(sig, fs, fmin, fmax, n=3, G=10, fr=1000):
     Outputs
     -------
     spec : numpy.ndarray
+<<<<<<< HEAD
+<<<<<<< HEAD
+        Third octave band spectrum of signal sig (units: Pascals [Pa]).
+=======
         The third octave band spectrum of signal sig with size (nfreq, nseg)
+>>>>>>> ad054c17c80ccea5874344f38a135c8fe1315a12
+=======
+        The third octave band spectrum of signal sig with size (nfreq, nseg)
+>>>>>>> ad054c17c80ccea5874344f38a135c8fe1315a12
     fpref : numpy.ndarray
         Corresponding preferred third octave band center frequencies
     """

mosqito/sq_metrics/tonality/tonality_iso1996K/comp_tonality.py

@@ -0,0 +1,119 @@
+# -*- coding: utf-8 -*-
+"""
+    Author: Cristina Taboada (TinaTabo)
+    Start date: 17/11/2021 
+"""
+# Standard library imports
+import math
+
+# Local imports
+from mosqito.sound_level_meter.noct_spectrum.noct_spectrum import noct_spectrum
+
+
+def comp_tonality(sig, fs):
+    """
+    <Function to obtain the prominent tones of a given signal
+    according to the procedure described in ISO 1996-2 Annex K.>
+
+    Parameters
+    ----------
+    sig : numpy.array
+        time signal values
+    fs : integer
+        sampling frequency
+
+    Outputs
+    -------
+    prominent_tones: dictionary
+        dictionary with {fc:Lp_dB} pairs where prominent tones are detected.
+    """
+
+
+    #-- As the tonality is studied for the audible frequency range, 
+    #-- we set the minimum and maximum frequencies at 25 Hz and 20 kHz.
+    fmin = 25
+    fmax = 10000
+
+    # -- we obtain the data of the Lp in Pa in thirds of octave of the signal of which
+    # -- we want to know the prominent tones
+    third_spec = noct_spectrum(sig=sig, fs=fs, fmin=fmin, fmax=fmax)
+
+    #-- Returns a tuple with two arrays, one with the Lp_Pa of each third octave band 
+    #-- and the other with the center frequencies, fc, of each band.
+    #-- Convert tuple to list for further processing
+    third_spec = list(third_spec)
+
+    # -- Obtain the lists of the central frequencies and the average Lp
+    fc = third_spec[1].tolist()
+    Lp_Pa = third_spec[0].tolist()
+
+    #-- Create a list with the Lp conversion in dB.
+    Lp_dB = []
+    P_ref = 2e-05
+    for i in range(0, len(Lp_Pa)):
+        P = Lp_Pa[i][0]
+        level = 20*math.log10(P/P_ref)
+        Lp_dB.append(level)
+
+    # -- List where the indexes corresponding to the positions where there is
+    # -- a prominent tone will be stored.
+    index_tone_list = []
+
+    # -- length of the lists
+    idx = len(fc) - 1
+
+    # -- level differences depending on the frequency range to determine if there
+    # -- is a prominent pitch.
+    diff_low_freqs = 15.0
+    diff_medium_freqs = 8.0
+    diff_high_freqs = 5.0
+
+    for x in range(0, idx):
+        # -- Variables to compare
+        if x == 0:
+            Lp_prev = 0
+            Lp_central = abs(Lp_dB[x])
+            Lp_post = abs(Lp_dB[x + 1])
+        else:
+            Lp_prev = abs(Lp_dB[x - 1])
+            Lp_central = abs(Lp_dB[x])
+            Lp_post = abs(Lp_dB[x + 1])
+
+        # -- calculate the difference
+        Lp_diff_prev = Lp_central - Lp_prev
+        Lp_diff_post = Lp_central - Lp_post
+
+        # -- if the value of the difference is constant with respect to the bands below and above
+        # -- the one studied, we obtain the value of the difference and proceed to check if we have
+        # -- found a prominent tone.
+
+        if x > fc.index(25.0) and x < fc.index(160.0):
+            # -- "LOW FREQUENCY --> difference 15 dB".
+            if Lp_diff_prev >= diff_low_freqs and Lp_diff_post >= diff_low_freqs:
+                # -- there is a tone in x, we store its value.
+                index_tone_list.append(x)
+
+        elif x > fc.index(125.0) and x < fc.index(500.0):
+            # -- "HALF FREQUENCY --> difference 8 dB"
+            if Lp_diff_prev >= diff_medium_freqs and Lp_diff_post >= diff_medium_freqs:
+                # -- there is a tone in x, we store its value.
+                index_tone_list.append(x)
+        elif x > fc.index(400.0) and x < fc.index(10000.0):
+            # -- "HIGH FREQUENCY --> 5 dB difference".
+            if Lp_diff_prev >= diff_high_freqs and Lp_diff_post >= diff_high_freqs:
+                # -- there is a tone in x, we store its value.
+                index_tone_list.append(x)
+
+    # -- Dictionary in which the data corresponding to the prominent tones
+    # -- will be stored, i.e. {fc : Lp_dB}
+    prominent_tones = {}
+
+    for i in range(0, len(index_tone_list)):
+        index = index_tone_list[i]
+        key = fc[index]
+        value = Lp_dB[index]
+        prominent_tones[key] = value
+
+    # -- Return of the function.
+    return prominent_tones
+

tests/sq_metrics/tonality/test_ISO1996K.py

@@ -0,0 +1,102 @@
+import pytest
+
+from mosqito.utils.load import load
+from mosqito.sq_metrics.tonality.tonality_iso1996K.comp_tonality import comp_tonality
+
+
+@pytest.mark.tonality_1996k  # to skip or run PR test
+def test_ISO1996K():
+    """Test function for the detection of a prominent tone 
+    using the inspection method described in ISO 1996-2 Annex K.
+
+    Validation of the function "comp_tonality" with signal array
+    as input. The input signals have been obtained from the website https://freesound.org/.
+
+    Parameters
+    ----------
+    None
+
+    Outputs
+    -------
+    None
+    """
+
+    """-----TONE--100-Hz----"""
+    sig, fs = load("tests\input\TONE100HZ.wav")
+    tones = comp_tonality(sig, fs)
+    print("----RESULT-----")
+    print(tones)
+    print("---------------")
+    assert tones == {100.0: 84.85146128088218}
+
+
+    """-----TONE--200-Hz----"""
+    sig, fs = load("tests\input\TONE200HZ.wav")
+    tones = comp_tonality(sig, fs)
+    print("----RESULT-----")
+    print(tones)
+    print("---------------")
+    assert tones == {200.0: 84.85835085762665}
+
+
+    """-----TONE--1-KHz----"""
+    sig, fs = load("tests/input/TONE1000HZ.WAV")
+    tones = comp_tonality(sig, fs)
+    print("----RESULT-----")
+    print(tones)
+    print("---------------")
+    assert tones == {1000.0: 63.962211067656554}
+
+
+    """-----TONE--2-KHz----"""
+    sig, fs = load("tests\input\TONE2000HZ.wav")
+    tones = comp_tonality(sig, fs)
+    print("----RESULT-----")
+    print(tones)
+    print("---------------")
+    assert tones == {2000.0: 84.95769475723816}
+
+
+    """-----TONE--4-KHz----"""
+    sig, fs = load("tests\input\TONE4000HZ.wav")
+    tones = comp_tonality(sig, fs)
+    print("----RESULT-----")
+    print(tones)
+    print("---------------")
+    assert tones == {4000.0: 84.95769513303429}
+
+
+    """-----TONE--5-kHz----"""
+    sig, fs = load("tests\input\TONE5000HZ.wav")
+    tones = comp_tonality(sig, fs)
+    print("----RESULT-----")
+    print(tones)
+    print("---------------")
+    assert tones == {5000.0: 84.95769534548218}
+
+
+    """-----MULTITONE-ALARM----"""
+    sig, fs = load("tests\input\MULTITONE_ALARM.wav")
+    tones = comp_tonality(sig, fs)
+    print("----RESULT-----")
+    print(tones)
+    print("---------------")
+    assert tones == {800.0: 53.685430439067396, 2000.0: 55.97882939521264, 3150.0: 64.05469163499629, 5000.0: 74.90550800435562}
+
+
+    """-----MUTITONE-SIREN----"""
+    sig, fs = load("tests\input\MULTITONE_SIREN.wav")
+    tones = comp_tonality(sig, fs)
+    print("----RESULT-----")
+    print(tones)
+    print("---------------")
+    assert tones == {1600.0: 61.618174324534216, 2500.0: 56.86145723559116, 8000.0: 38.87188014184665}
+
+
+    """-----ATONAL-SIGNAL----"""
+    sig, fs = load("tests\input\WHITE_NOISE.wav")
+    tones = comp_tonality(sig, fs)
+    print("----RESULT-----")
+    print(tones)
+    print("---------------")
+    assert tones == {}