Repository: Eomys/MoSQITo (language: python)

Issue:
Loudness according to Annex F of ECMA-74
Thanks to @djcaminero work, we have a first draft of [ECMA-74](https://www.ecma-international.org/publications-and-standards/standards/ecma-74/) annex F hearing model to compute loudness in MOSQITO. However  the model doesn't pass the validation stage yet (see figure below with the loudness contour from ISO 226 in dotted line and the loudness contour from MOSQITO hearing model in solid lines).

![valid_loudness_ecma](https://user-images.githubusercontent.com/60893296/120209654-38600380-c22f-11eb-98e7-cc239bda9cf9.png)

Any help from the community to review @djcaminero [code ](https://github.com/djcaminero/MoSQITo)will be welcome. You can reproduce the figure above by running mosqito/validation/hearing_model/hearing_model_validation.py script. That's a good starting point for a review. 

If you would like to contribute but are not familiar with github, please contact me. I will be happy to help.




Discussion / hints:
I'm begining a step by step validation of the implementation. Step 1 is the outer and middle/inner ear filtering (section F 3.2 of the ECMA). By running mosqito/validations/hearing_model/valid_ear_filter_design.py, i get a response compliant to the standard, see figure below (ECMA-74 on the left,  MOSQITO on the right)
![ecma-74_fig_F3_ear_filter_cmp](https://user-images.githubusercontent.com/60893296/120796856-9e080480-c53b-11eb-9ad0-6f48edec8acc.png)

Second step is the gammatone filters design. I have created a new 'gammatone' function from 'afb_coefficients' and part of 'afb_ppal_parameters' the objective is to be able to easily compare MOSQITO implementation of the gammatone filter with the one of scipy. The function is not used in the main code yet, but it should be eventually. 

The filters from mosqito match the one from scipy enough to consider that it is, at least, not totally wrong considering that it is not necessarily the same auditory model that is implemented (run mosqito/validations/hearing_model/valid_gammatone.py to get the plots):
![gammatone_400Hz](https://user-images.githubusercontent.com/60893296/120829776-141d6300-c55e-11eb-8c06-c48f1bddbc5e.png)
![gammatone_4000Hz](https://user-images.githubusercontent.com/60893296/120829815-20a1bb80-c55e-11eb-962d-f2ed4ab7b7bd.png)

~~I checked the validation script and realized that it was not giving a correct evaluation of the implementation. I've programmed an [alternative validation script](https://github.com/djcaminero/MoSQITo/blob/master/mosqito/validations/hearing_model/hearing_model_validation_alt.py). The principle is the following : for each frequency, the Loudness is calculated for 4 tones with an SPL amplitude determined by the isophone curves 20, 40, 60 and 80. A satisfying implementation should return values in phon close to 20, 40, 60 and 80 respectively. The result of the current implementation is the following (theoretical isophon in dotted lines, results from MOSQITO in solid lines):~~

Update 13/08/2021
I've programmed an [alternative validation script](https://github.com/djcaminero/MoSQITo/blob/master/mosqito/validations/hearing_model/hearing_model_validation_alt.py). The principle is the following : for each frequency, the Loudness is calculated for 4 tones with an SPL amplitude determined by the isophone curves 20, 40, 60 and 80. A satisfying implementation should return values in phon close to 20, 40, 60 and 80 for each frequency (dotted lines in the figure below). The current implementation is in solid lines.

![image](https://user-images.githubusercontent.com/60893296/129394204-7a0c1773-6b02-4b0c-841f-1e6107c20fc0.png)

I changed variable _centre_freq_array_ in the call of function `annex_f_validation()` in `comp_loudness()` to variable _centre_freq_ for working diagnostics: (`comp_loudness(signal, validation=True)`).

More detailed explanation:
If the validation option of function `comp_loudness(signal, validation=False)` is set to _True_, function `annex_f_validation()` is called at the end with one of the arguments being _centre_freq_array_. In the current version, variable _centre_freq_array_ is an empty list, thus, interrupting the validation process. The appropriate variable to finish the validation is called _centre_freq_:

[...]
    if validation:
        annex_f_validation(
            fs,
            signal,
            signal_filtered,
            sos_ear,
            **centre_freq,**
            f_bandwidth_array,
            t_delay_array,
            d_coefficients_array,
            sb_array,
            sh_array,
            am_mod_coefficient_array,
            bm_mod_coefficient_array,
            band_pass_signal_array,
            n_array,
            t_array,
        )

I've been reviewing the whole code. I took this opportunity to implement an alternative version of the [comp_loudness function](https://github.com/djcaminero/MoSQITo/blob/master/mosqito/functions/hearing_model/comp_loudness_alt.py). Some of the steps have been simplified to allow for easyer debugging. However, the results remains the same.
I noticed that if I do not apply the Outer and middle/inner ear filtering (section 5.1.2 of the ECMA-418-2), I get a overall trend of the hearing model that is closer to the theoretical equal loudness contours. It suggests that the filter is not applied correctly to the signal.
 
![image](https://user-images.githubusercontent.com/60893296/127649012-658c918e-6744-42e8-b342-2d03dd6bf66f.png)

I've finally spotted ~~two~~ one errors in the "outer and midle ear filtering" section of the code.

1. The sosfiltfilt function was used to apply the filter. 
```Python
signal_filtered = sp.signal.sosfiltfilt(sos_ear, signal, axis=0)
```
As stated in the doc : "This function applies a linear digital filter twice, once forward and once backwards. The combined filter has zero phase and a filter order twice that of the original." The filter order was then twice too high.

~~2. The filter predict the response of the ear. In order to "compensate" this response and go from physical amplitude in dB to physiological amplitude in sone, one should apply the inverse filter:~~

~~We are now getting dangerously close to the theoretical equal loudness contours. Next mistakes are to be search in a portion of the code that affect differently the low and high amplitudes. I would guess the "Non-linerarity" section.~~

![image](https://user-images.githubusercontent.com/60893296/129394761-73583543-07ae-4fce-b78e-4c94b0f9a64c.png)




Fix patch:
diff --git a/README.md b/README.md
index 68dcafa3..b31a7027 100644
--- a/README.md
+++ b/README.md
@@ -61,7 +61,9 @@ If you need to cite the current release of MOSQITO, please use the "Cite this re
 
 ## Publications citing MOSQITO
 
-Glesser, M., Ni, S., Degrendele, K., Wanty, S., & Le Besnerais, J. (2021, October 13-14). *Sound quality analysis of Electric Drive Units under different switching control strategies* [Paper presentation]. Automotive NVH Comfort Congress 2021, Le Mans , France.
+M. Glesser, S. Wanty, K. Degrendele, J. L. Besnerais, and S. Ni, *“Perceived sound quality analysis of Electric Drive Units under different switching control strategies,”* in 12th Aachen Acoustic Colloquium, Aachen, Nov. 2021.
+
+M. Glesser, S. Ni, K. Degrendele, S. Wanty, and J. L. Besnerais, *“Sound quality analysis of Electric Drive Units under differentswitching control strategies,”* in SIA Automotive NVH comfort, Le Mans, 2021.
 
 San Millán-Castillo, R., Latorre-Iglesias, E., Glesser, M., Wanty, S., Jiménez-Caminero, D., & Álvarez-Jimeno, J.M. (2021). MOSQITO: an open-source and free toolbox for sound quality metrics in the industry and education. *INTER-NOISE and NOISE-CON Congress and Conference Proceedings*, 12, 1164-1175. https://doi.org/10.3397/IN-2021-1767
 
diff --git a/documentation/loudness-ecma.md b/documentation/loudness-ecma.md
new file mode 100644
index 00000000..f443831e
--- /dev/null
+++ b/documentation/loudness-ecma.md
@@ -0,0 +1,23 @@
+# MOSQITO Documentation
+## Loudness according to ECMA 418-2
+
+### Introduction
+The acoustic loudness calculation according to the Hearing Model of Sottek is described on the ECMA 418-2 satndard. The method is in accordance with ISO 226:2003 equal loudness contours (instead of older ISO 226:1987 definition used for the Zwicker Loudness).
+
+The result of the hearing model can be used as the basis for further psychoacoustic analyses (Psychoacoustic tonality according to ECMA-418-2 section 6 or roughness according to ECMA-418-2 section 7).
+
+### Validation of the implementation
+The ECMA loudness calculation is evaluated as proposed in the annex A of the standard. The loudness is calculated for sinusoidal signals with a frequency of 1000 Hz and a sound pressure level of 20 to 80 dB with a step size of 20 dB. For other frequencies, the level is varied to match the loudness calculated for the 1000 Hz tone. The results are shown in the figure below (to be compared to the Figure A.1 of the standard. The target equal-loudness contours are emulated well by the results of the hearing model.
+
+![](../validations/loudness_ecma/output/ecma_hearing_model_validation.png)
+
+*Results for the equal-loudness contours. The dotted lines show the target equal-loudness contours, the solid lines are the equal-loudness contours obtained with the MOSQITO hearing model implementation*
+
+The validation plots and scripts can be found in [this folder](../validations/loudness_ecma).
+
+
+### References
+
+ECMA-418-2, Psychoacoustic metrics for ITT equipment — Part 2 (models based on human perception), 1st edition (December 2020) [[available online](https://www.ecma-international.org/wp-content/uploads/ECMA-418-2_1st_edition_december_2020.pdf)]
+
+ISO 226:2003, Acoustics — Normal equal-loudness-level contours 
diff --git a/documentation/readme.md b/documentation/readme.md
index f3eb4492..7acb6bfe 100644
--- a/documentation/readme.md
+++ b/documentation/readme.md
@@ -13,6 +13,7 @@ In this folder, you will find the documentation and validation of the different
 
 - [Loudness of stationary signals (Zwicker method)](./loudness-stationary.md)
 - [Loudness of time varying signals (Zwicker method)](./loudness-time-varying.md)
+- [Loudness of time varying signals (ECMA method)](./loudness-ecma.md)
 - [Sharpness of stationary and time-varying signals](./sharpness.md)
 - [Roughness of stationary and time-varying signals (Daniel and Weber method)](./roughness.md)
 - [Tone-to-noise ratio and Prominence ratio of stationary and time-varying signals](./tone-to-noise-prominence-ratio.md)
diff --git a/documentation/scope.md b/documentation/scope.md
index 25854a64..5a59d042 100644
--- a/documentation/scope.md
+++ b/documentation/scope.md
@@ -5,15 +5,15 @@ The scope of the project is to implement the following first set of
 metrics:
 
 |                                                    | Reference                                            | Validated                                          | Available                                     | Under dev. | To do |
-|:-------------------------------------------------- |:---------------------------------------------------- |:--------------------------------------------------:|:---------------------------------------------:|:----------:|:-----:|
-| Loudness for<br>steady signals<br>(Zwicker method) | ISO 532B:1975<br>DIN 45631:1991<br>ISO 532-1:2017 §5 | [x](../validations/loudness_zwicker/output) | [x](./loudness-stationary.md)   |            |       |
-| Loudness for non-stationary<br>(Zwicker method)    | DIN 45631/A1:2010<br>ISO 532-1:2017 §6               | [x](../validations/loudness_zwicker/output) | [x](./loudness-time-varying.md) |            |       |
-| Loudness for non-stationary<br>(ECMA-74 method)    | ECMA-74:2019, annex F<br>Sottek, 2016               |  |  |     x       |       |
-| Roughness                                          | Daniel and Weber, 1997                               | [x](../validations/roughness_danielweber)   | [x](./roughness.md)             |            |       |
-| Roughness | ECMA-418-2:2020 | |                                               |            | x     |
-| Fluctuation Strength                               | To be defined                                        |                                                    |                                               |            | x     |
-| Sharpness                                          | DIN 45692:2009                                       | [x](../validations/sharpness/output)        | [x](./sharpness.md)             |            |       |
-| Tonality (Hearing model)                           | ECMA-74:2019 annex G                                 |                                                    |                                               | x          |       |
+|:-------------------------------------------------- |:---------------------------------------------------- |:-------------------------------------------:|:---------------------------------------------:|:----------:|:-----:|
+| Loudness for<br>steady signals<br>(Zwicker method) | ISO 532B:1975<br>DIN 45631:1991<br>ISO 532-1:2017 §5 | [x](./loudness-stationary.md)               | x |  |  |
+| Loudness for non-stationary<br>(Zwicker method)    | DIN 45631/A1:2010<br>ISO 532-1:2017 §6               | [x](./loudness-time-varying.md)             | x |  |  |
+| Loudness for non-stationary<br>(ECMA method)       | ECMA 418-2:2019 section 5<br>Sottek, 2016            | [x](./loudness-ecma.md)                     | x |  |  |
+| Roughness                                          | Daniel and Weber, 1997                               | [x](../validations/roughness_danielweber)   | x |  |  |
+| Roughness                                          | ECMA-418-2:2020                                      |                                             |  |  | x |
+| Fluctuation Strength                               | To be defined                                        |                                             |  |  | x |
+| Sharpness                                          | DIN 45692:2009                                       | [x](./sharpness.md)                         | x |  |  |
+| Tonality (Hearing model)                           | ECMA-74:2019 annex G                                 |                                             |  | x |  |
 
 As a second priority, the project could address the following metrics:
 
diff --git a/mosqito/functions/loudness_ecma_spain/__init__.py b/mosqito/functions/loudness_ecma/__init__.py
similarity index 100%
rename from mosqito/functions/loudness_ecma_spain/__init__.py
rename to mosqito/functions/loudness_ecma/__init__.py
diff --git a/mosqito/functions/loudness_ecma/comp_loudness.py b/mosqito/functions/loudness_ecma/comp_loudness.py
new file mode 100644
index 00000000..71fc836d
--- /dev/null
+++ b/mosqito/functions/loudness_ecma/comp_loudness.py
@@ -0,0 +1,84 @@
+# -*- coding: utf-8 -*-
+"""
+@author: Daniel Jiménez-Caminero Costa
+"""
+import numpy as np
+import matplotlib.pyplot as plt
+
+# Project Imports
+from mosqito.functions.loudness_ecma.rectified_band_pass_signals import (
+    rectified_band_pass_signals,
+)
+from mosqito.functions.loudness_ecma.nonlinearity import nonlinearity
+
+# Data import
+# Threshold in quiet
+from mosqito.functions.loudness_ecma.loudness_ecma_data import ltq_z
+
+
+def comp_loudness(signal, sb=2048, sh=1024):
+    """Calculation of the specific and total loudness according to ECMA-418-2 section 5
+
+    Parameters
+    ----------
+    signal: numpy.array
+        time signal values in 'Pa'. The sampling frequency of the signal must be 48000 Hz.
+    sb: int or list of int
+        block size.
+    sh: int or list of int
+        Hop size.
+
+    Returns
+    -------
+    n_specific: list of numpy.array
+        Specific Loudness [sone_HMS per Bark]. Each of the 53 element of the list corresponds to the time-dependant
+        specific loudness for a given bark band. Can be a ragged array if a different sb/sh are used for each band.
+
+    bark_axis: numpy.array
+        Bark axis
+
+    """
+
+    # TODO: return time axis (list or list of list)
+
+    # Computaton of rectified band-pass signals
+    # (section 5.1.2 to 5.1.5 of the standard)
+    block_array_rect = rectified_band_pass_signals(signal, sb, sh)
+
+    # # sb and sh for Tonality
+    # z = np.linspace(0.5, 26.5, num=53, endpoint=True)
+    # sb = np.ones(53, dtype="int")
+    # sh = np.ones(53, dtype="int")
+    # sb[z <= 1.5] = 8192
+    # sh[z <= 1.5] = 2048
+    # sb[np.all([z >= 2, z <= 8], axis=0)] = 4096
+    # sh[np.all([z >= 2, z <= 8], axis=0)] = 1024
+    # sb[np.all([z >= 8.5, z <= 12.5], axis=0)] = 2048
+    # sh[np.all([z >= 8.5, z <= 12.5], axis=0)] = 512
+    # sb[z >= 13] = 1024
+    # sh[z >= 13] = 256
+
+    n_specific = []
+    for band_number in range(53):
+        # ROOT-MEAN-SQUARE (section 5.1.6)
+        # After the segmentation of the signal into blocks, root-mean square values of each block are calculated
+        # according to Formula 17.
+        rms_block_value = np.sqrt(
+            2 * np.mean(np.array(block_array_rect[band_number]) ** 2, axis=1)
+        )
+
+        # NON-LINEARITY (section 5.1.7)
+        # This section covers the other part of the calculations needed to consider the non-linear transformation
+        # of sound pressure to specific loudness that does the the auditory system. After this point, the
+        # computation is done equally to every block in which we have divided our signal.
+        a_prime = nonlinearity(rms_block_value)
+
+        # SPECIFIC LOUDNESS CONSIDERING THE THRESHOLD IN QUIET (section 5.1.8)
+        # The next calculation helps us obtain the result for the specific loudness - specific loudness with
+        # consideration of the lower threshold of hearing.
+        a_prime[a_prime < ltq_z[band_number]] = ltq_z[band_number]
+        N_prime = a_prime - ltq_z[band_number]
+        n_specific.append(N_prime)
+
+    bark_axis = np.linspace(0.5, 26.5, num=53, endpoint=True)
+    return n_specific, bark_axis
diff --git a/mosqito/functions/loudness_ecma_spain/ear_filter_design.py b/mosqito/functions/loudness_ecma/ear_filter_design.py
similarity index 100%
rename from mosqito/functions/loudness_ecma_spain/ear_filter_design.py
rename to mosqito/functions/loudness_ecma/ear_filter_design.py
diff --git a/mosqito/functions/loudness_ecma_spain/equal_loudness_contours.py b/mosqito/functions/loudness_ecma/equal_loudness_contours.py
similarity index 100%
rename from mosqito/functions/loudness_ecma_spain/equal_loudness_contours.py
rename to mosqito/functions/loudness_ecma/equal_loudness_contours.py
diff --git a/mosqito/functions/loudness_ecma_spain/gammatone.py b/mosqito/functions/loudness_ecma/gammatone.py
similarity index 100%
rename from mosqito/functions/loudness_ecma_spain/gammatone.py
rename to mosqito/functions/loudness_ecma/gammatone.py
diff --git a/mosqito/functions/loudness_ecma_spain/gen_auditory_filters_centre_freq.py b/mosqito/functions/loudness_ecma/gen_auditory_filters_centre_freq.py
similarity index 100%
rename from mosqito/functions/loudness_ecma_spain/gen_auditory_filters_centre_freq.py
rename to mosqito/functions/loudness_ecma/gen_auditory_filters_centre_freq.py
diff --git a/mosqito/functions/loudness_ecma_spain/loudness_ecma_data.py b/mosqito/functions/loudness_ecma/loudness_ecma_data.py
similarity index 100%
rename from mosqito/functions/loudness_ecma_spain/loudness_ecma_data.py
rename to mosqito/functions/loudness_ecma/loudness_ecma_data.py
diff --git a/mosqito/functions/loudness_ecma/nonlinearity.py b/mosqito/functions/loudness_ecma/nonlinearity.py
new file mode 100644
index 00000000..3301d723
--- /dev/null
+++ b/mosqito/functions/loudness_ecma/nonlinearity.py
@@ -0,0 +1,34 @@
+import numpy as np
+
+
+def nonlinearity(p):
+    """Apply the compressive nonlinearity of the auditory system to the
+    rectified band pass signal rms values according to ECMA 418-2 section 5.1.7
+
+    Parameters
+    ----------
+    p: numpy.array
+        RMS value of the rectified band pass signal blocks
+
+    Returns
+    -------
+    a_prime: numpy.array
+        Specific Loudness [sone_HMS per Bark].
+
+    """
+    p_0 = 2e-5
+    # c_N: In sones/bark
+    c_N = 0.0217406
+    alpha = 1.50
+    v_i = np.array(
+        [1.0, 0.6602, 0.0864, 0.6384, 0.0328, 0.4068, 0.2082, 0.3994, 0.6434]
+    )
+    thresh = np.array([15.0, 25.0, 35.0, 45.0, 55.0, 65.0, 75.0, 85.0])
+    p_ti = p_0 * 10 ** (thresh / 20)
+
+    a_prime = np.ones(p.shape)
+    for i in range(8):
+        a_prime *= (1 + (p / p_ti[i]) ** alpha) ** ((v_i[i + 1] - v_i[i]) / alpha)
+    a_prime *= c_N * p / p_0
+
+    return a_prime
\ No newline at end of file
diff --git a/mosqito/functions/loudness_ecma_spain/phone2spl.py b/mosqito/functions/loudness_ecma/phone2spl.py
similarity index 100%
rename from mosqito/functions/loudness_ecma_spain/phone2spl.py
rename to mosqito/functions/loudness_ecma/phone2spl.py
diff --git a/mosqito/functions/loudness_ecma_spain/rectified_band_pass_signals.py b/mosqito/functions/loudness_ecma/rectified_band_pass_signals.py
similarity index 50%
rename from mosqito/functions/loudness_ecma_spain/rectified_band_pass_signals.py
rename to mosqito/functions/loudness_ecma/rectified_band_pass_signals.py
index 60443832..83d6bc44 100644
--- a/mosqito/functions/loudness_ecma_spain/rectified_band_pass_signals.py
+++ b/mosqito/functions/loudness_ecma/rectified_band_pass_signals.py
@@ -3,23 +3,23 @@
 import numpy as np
 import scipy.signal as sp_signal
 
-from mosqito.functions.loudness_ecma_spain.ear_filter_design import ear_filter_design
-from mosqito.functions.loudness_ecma_spain.gen_auditory_filters_centre_freq import (
+from mosqito.functions.loudness_ecma.ear_filter_design import ear_filter_design
+from mosqito.functions.loudness_ecma.gen_auditory_filters_centre_freq import (
     gen_auditory_filters_centre_freq,
 )
-from mosqito.functions.loudness_ecma_spain.gammatone import gammatone
-from mosqito.functions.loudness_ecma_spain.segmentation_blocks import (
+from mosqito.functions.loudness_ecma.gammatone import gammatone
+from mosqito.functions.loudness_ecma.segmentation_blocks import (
     segmentation_blocks,
 )
 
 
 def rectified_band_pass_signals(sig, sb=2048, sh=1024):
-    """Compute the rectified band-pass signals as per Clause 5.1.5 of ECMA-418-2:2020
+    """Compute the rectified band-pass signals as per Clause 5.1.2 to 5.1.5 of
+    ECMA-418-2:2020
 
     Calculation of the rectified band-pass signals along the 53 critical band rates
     scale. Each band pass signal is segmented into time blocks according to sb and sh
 
-
     Parameters
     ----------
     signal: numpy.array
@@ -60,14 +60,7 @@ def rectified_band_pass_signals(sig, sb=2048, sh=1024):
     block_array_rect = []
     for band_number in range(53):
         bm_mod, am_mod = gammatone(centre_freq[band_number], k=filter_order_k, fs=fs)
-        # bm_mod, am_mod = sp_signal.gammatone(centre_freq[band_number], "fir", fs=fs)
-
-        """ 
-        "scipy.signal.lfilter" instead of "scipy.signal.filtfilt" in order to maintain consistency. That process 
-        makes possible to obtain a signal "band_pass_signal" that does not line up in time with the original signal 
-        because of the non zero-phase filtering of "lfilter", but it has a more appropriate slope than filtfilt. 
-        By using filtfilt the slope is that high that filters too much the signal. 
-        """
+
         band_pass_signal = (
             2.0
             * (
@@ -80,26 +73,26 @@ def rectified_band_pass_signals(sig, sb=2048, sh=1024):
             ).real
         )
 
-        """SEGMENTATION OF THE SIGNAL INTO BLOCKS (5.1.4)
+        # SEGMENTATION OF THE SIGNAL INTO BLOCKS (5.1.4)
+
+        # The segmentation of the signal is done in order to obtain results for intervals of time, not for the whole
+        # duration of the signal. The reason behind this decision resides in the fact that processing the signal in its
+        # full length at one time could end up in imprecise results. By using a "for loop", we are able to decompose the
+        # signal array "band_pass_signal_hr" into blocks. "sb_array" is the block size which changes depending on the
+        # "band_number" in which we are processing the signal. "sh_array" is the step size, the time shift to the next
+        # block.
 
-        The segmentation of the signal is done in order to obtain results for intervals of time, not for the whole
-        duration of the signal. The reason behind this decision resides in the fact that processing the signal in its
-        full length at one time could end up in imprecise results. By using a "for loop", we are able to decompose the
-        signal array "band_pass_signal_hr" into blocks. "sb_array" is the block size which changes depending on the
-        "band_number" in which we are processing the signal. "sh_array" is the step size, the time shift to the next
-        block.
-        """
         block_array = segmentation_blocks(
-            band_pass_signal, sb[band_number], sh[band_number], dim=1
+            band_pass_signal, sb[band_number], sh[band_number]
         )
 
-        """RECTIFICATION (5.1.5)
+        # RECTIFICATION (5.1.5)
+
+        # This part acts as the activation of the auditory nerves when the basilar membrane vibrates in a certain
+        # direction. In order to rectify the signal we are using "np.clip" which establish a minimum and a maximum value
+        # for the signal. "a_min" is set to 0 float, while "a_max" is set to "None" in order to consider the positive
+        # value of the signal.
 
-        This part acts as the activation of the auditory nerves when the basilar membrane vibrates in a certain
-        direction. In order to rectify the signal we are using "np.clip" which establish a minimum and a maximum value
-        for the signal. "a_min" is set to 0 float, while "a_max" is set to "None" in order to consider the positive
-        value of the signal.
-        """
         block_array_rect.append(np.clip(block_array, a_min=0.00, a_max=None))
 
     return block_array_rect
\ No newline at end of file
diff --git a/mosqito/functions/loudness_ecma/segmentation_blocks.py b/mosqito/functions/loudness_ecma/segmentation_blocks.py
new file mode 100644
index 00000000..f32afe85
--- /dev/null
+++ b/mosqito/functions/loudness_ecma/segmentation_blocks.py
@@ -0,0 +1,45 @@
+# -*- coding: utf-8 -*-
+"""
+@author: Daniel Jiménez-Caminero Costa
+"""
+
+import numpy as np
+
+
+def segmentation_blocks(band_pass_signal_hr, sb, sh):
+    """Function used for the segmentation of the signal into smaller parts of audio (blocks).
+    This has been implemented as described in Formula 16 (section 5.1.4) of ECMA-418-2.
+
+    Parameters
+    ----------
+    band_pass_signal_hr: numpy.array
+        time signal values of the input signal
+    sb: int
+        block size
+    sh: int
+        Hop size.
+
+    Returns
+    -------
+    block_array: List
+        Array list of blocks in which the signal has been segmented.
+    """
+
+    # TODO: check first blocks
+
+    l = 0
+    block_array = []
+    while l * sh <= len(band_pass_signal_hr):
+        i_begin = l * sh - sb
+        i_end = l * sh
+        if i_begin < 0 and i_end <= 0:
+            block_array.append(np.zeros(sb))
+        elif i_begin < 0:
+            block_array.append(
+                np.concatenate((np.zeros(-i_begin), band_pass_signal_hr[0:i_end]))
+            )
+        else:
+            block_array.append(band_pass_signal_hr[i_begin:i_end])
+        l += 1
+
+    return block_array
diff --git a/mosqito/functions/loudness_ecma_spain/sone2phone.py b/mosqito/functions/loudness_ecma/sone2phone.py
similarity index 100%
rename from mosqito/functions/loudness_ecma_spain/sone2phone.py
rename to mosqito/functions/loudness_ecma/sone2phone.py
diff --git a/mosqito/functions/loudness_ecma_spain/afb_coefficients.py b/mosqito/functions/loudness_ecma_spain/afb_coefficients.py
deleted file mode 100644
index 70cb1398..00000000
--- a/mosqito/functions/loudness_ecma_spain/afb_coefficients.py
+++ /dev/null
@@ -1,94 +0,0 @@
-# -*- coding: utf-8 -*-
-"""
-@author: Daniel Jiménez-Caminero Costa
-"""
-import cmath
-import numpy as np
-import math
-import scipy as sp
-from scipy.special import comb
-
-
-def afb_coefficients(fs, filter_order_k, centre_freq, d_coefficients):
-    """Function for the calculation of the filter coefficients in the Auditory Filtering Bank section (5.1.3). Here, as
-    it has been mentioned in the principal function for the specific loudness, it is calculated the band-pass equivalent
-    of the low-pass filter. The typo/error found in the formulas (13 and 14) presented in the ECMA-418-2 has been
-    already corrected.
-
-    Parameters
-    ----------
-    fs: float
-        'Hz', sampling frequency.
-
-    filter_order_k: int
-        Order of the auditory filtering bank.
-
-    centre_freq: float
-        'Hz', central frequency of the filter.
-
-    d_coefficients: float
-        "d" coefficient, related with sampling rate and delay.
-
-    Returns
-    -------
-    am_mod_coefficient_band_array: numpy.array
-        "am" coefficient to filter a signal.
-
-    bm_mod_coefficient_band_array: numpy.array
-        "bm" coefficient to filter a signal.
-    """
-    print("This function is deprecated and will be removed form the project")
-
-    # The "e" coefficients are used for the calculation of the bm coefficients for the Auditory Filtering Bank
-    e_coefficients_array = [0.0, 1.0, 11.0, 11.0, 1.0]
-    # Coefficients for the "Auditory Filtering Bank"
-    am_mod_coefficient_band_array = np.zeros(int(filter_order_k + 1), dtype=complex)
-    bm_mod_coefficient_band_array = np.zeros(int(filter_order_k + 1), dtype=complex)
-
-    # Coefficient calculation
-    bm_sum = 0.00
-
-    # "b" summation sequence, fraction calculation (F.12)
-    for j in range(filter_order_k - 1):
-        iterator_sum = j + 1
-        bm_sum = bm_sum + (
-            e_coefficients_array[iterator_sum] * (d_coefficients ** iterator_sum)
-        )
-
-    bm_numerator = (1 - d_coefficients) ** filter_order_k
-    bm_fraction = bm_numerator / bm_sum
-
-    # Implementation of the recursive formula (F.10)
-    for m in range(filter_order_k + 1):
-        exp_factor = (2 * math.pi * centre_freq * m) / fs
-        complex_exponential = 1j * exp_factor
-
-        # Binomial coefficient, "am" coefficient calculation (F.11)
-        if m == 0:
-            am_coefficient = 1.0
-        else:
-            binomial_coef_2 = float(sp.special.comb(filter_order_k, m, exact=True))
-
-            if (m % 2) == 0:
-                am_coefficient = (d_coefficients ** m) * binomial_coef_2
-            else:
-                am_coefficient = (-1) * (d_coefficients ** m) * binomial_coef_2
-
-        # "bm" coefficient calculation (F.12)
-        if m == filter_order_k:
-            bm_coefficient = 0.0
-        else:
-            bm_coefficient = (
-                bm_fraction * (d_coefficients ** m) * e_coefficients_array[m]
-            )
-
-        # Transformation from low-pass filter coefficients to band-pass.
-        # (*) In this step is where the error that is mentioned at the start of the Auditory Filtering Bank section
-        # was found
-        am_mod_coefficient_band = am_coefficient * cmath.exp(complex_exponential)
-        bm_mod_coefficient_band = bm_coefficient * cmath.exp(complex_exponential)
-
-        am_mod_coefficient_band_array[m] = am_mod_coefficient_band
-        bm_mod_coefficient_band_array[m] = bm_mod_coefficient_band
-
-    return am_mod_coefficient_band_array, bm_mod_coefficient_band_array
diff --git a/mosqito/functions/loudness_ecma_spain/afb_ppal_parameters.py b/mosqito/functions/loudness_ecma_spain/afb_ppal_parameters.py
deleted file mode 100644
index 61a90558..00000000
--- a/mosqito/functions/loudness_ecma_spain/afb_ppal_parameters.py
+++ /dev/null
@@ -1,93 +0,0 @@
-# -*- coding: utf-8 -*-
-"""
-@author: Daniel Jiménez-Caminero Costa
-"""
-import numpy as np
-import math
-import scipy as sp
-from scipy.special import comb
-
-
-def afb_ppal_parameters(fs, band_number, filter_order_k):
-    """Principal parameters for the Auditory Filtering Bank coefficients and other calculations in later sections:
-    central frequency, bandwidth, delay, d coefficient, block size and hop size.
-    This has been implemented as described in section 5.1.3 of ECMA-418-2.
-
-    Parameters
-    ----------
-    fs: float
-        'Hz', sampling frequency.
-
-    band_number: int
-        Band number in which the filter is going to be applied over the signal.
-
-    filter_order_k: int
-        Order of the auditory filtering bank.
-
-    Returns
-    -------
-    centre_freq: float
-        'Hz', central frequency of the filter.
-
-    f_bandwidth: float
-        Filter bandwidth.
-
-    t_delay: float
-        's'.
-
-    d_coefficients: float
-        "d coefficient", related with sampling rate and delay.
-
-    sb: int
-        Block size.
-
-    sh: int
-        Hop size.
-    """
-
-    print("This function is deprecated and will be removed form the project")
-
-    z_step_size = 0.50
-    af_f0 = 81.9289  # ECMA-418-2
-    c = 0.1618  # ECMA-418-2
-
-    # Exponent of the delay
-    exponent_1 = (2.0 * filter_order_k) - 1.0
-    # Binomial coefficient
-    n_binomial = (2.0 * filter_order_k) - 2.0
-    k_binomial = filter_order_k - 1.0
-    binomial_coef_1 = float(
-        sp.special.comb(int(n_binomial), int(k_binomial), exact=True)
-    )
-
-    # Critical band rate scale
-    z = (band_number + 1) * z_step_size
-    var = c * z
-
-    # Central frequency
-    centre_freq = (af_f0 / c) * math.sinh(var)
-
-    # Bandwidth
-    f_bandwidth = math.sqrt((af_f0 ** 2) + ((c * centre_freq) ** 2))
-
-    # Time constant, delay
-    t_delay = (1.0 / (2.0 ** exponent_1)) * binomial_coef_1 * (1.0 / f_bandwidth)
-
-    # "d" coefficient
-    d_coefficients = float(np.exp((-1.0) / (fs * t_delay)))
-
-    # Block length and hop size, for further calculations (Root-Mean-Square Values, section 5.1.6)
-    if z >= 13:
-        sb = 1024
-        sh = 256
-    elif 8.5 <= z < 13:
-        sb = 2048
-        sh = 512
-    elif 2 <= z < 8.5:
-        sb = 4096
-        sh = 1024
-    else:
-        sb = 8192
-        sh = 2048
-
-    return centre_freq, f_bandwidth, t_delay, d_coefficients, sb, sh
diff --git a/mosqito/functions/loudness_ecma_spain/comp_loudness.py b/mosqito/functions/loudness_ecma_spain/comp_loudness.py
deleted file mode 100644
index 6f37f2bd..00000000
--- a/mosqito/functions/loudness_ecma_spain/comp_loudness.py
+++ /dev/null
@@ -1,426 +0,0 @@
-# -*- coding: utf-8 -*-
-"""
-@author: Daniel Jiménez-Caminero Costa
-"""
-import numpy as np
-import scipy as sp
-import matplotlib.pyplot as plt
-from scipy.signal import welch
-
-# Project Imports
-from mosqito.functions.shared.load import load
-from mosqito.functions.hearing_model.sine_wave_generator import sine_wave_generator
-from mosqito.functions.hearing_model.ear_filter_design import ear_filter_design
-from mosqito.functions.hearing_model.gen_auditory_filters_centre_freq import (
-    gen_auditory_filters_centre_freq,
-)
-from mosqito.functions.hearing_model.gammatone import gammatone
-from mosqito.functions.hearing_model.afb_ppal_parameters import afb_ppal_parameters
-from mosqito.functions.hearing_model.segmentation_blocks import segmentation_blocks
-from mosqito.functions.hearing_model.segmentation_blocks import (
-    segmentation_blocks_test_a,
-)
-from mosqito.functions.hearing_model.segmentation_blocks import (
-    segmentation_blocks_test_b,
-)
-from mosqito.functions.hearing_model.nonlinear_common import nonlinear_common
-from mosqito.functions.hearing_model.nonlinear_repeat import nonlinear_repeat
-from mosqito.functions.hearing_model.sone2phone import sone2phone
-from mosqito.functions.hearing_model.sl_validation import annex_f_validation
-
-import sys
-
-sys.path.append("../../..")
-
-
-def comp_loudness(signal, validation=False):
-    """Calculation of specific loudness and total loudness of an input signal according to ECMA-418-2. It describes the
-    loudness excitation in a critical band per Bark. The output of the signal consists in a pair of array lists, one
-    for specific loudness and the other one for total loudness.
-
-    Parameters
-    ----------
-    validation: boolean
-        Default value is set to 'False'. It is used to mark whether complete validation is done or not.
-
-    signal: numpy.array
-        'Pa', time signal values. It can be, stereo (2 dimensions) or mono (1 dimension).The sampling frequency of the
-        signal must be 48000 Hz.
-
-    Returns
-    -------
-    n_array: numpy.array
-        'Sone per Bark'. It is a numpy array that is arranged as a matrix. Each matrix row represents a certain block of
-        audio (in chronological order) and each matrix column stores the value of the specific loudness for that block
-        and that specific band (column number). If the input signal is stereo, another matrix dimension is added.
-
-    t_array: numpy.array
-        'Sone per Bark'. As well as "n_array", it is a numpy array that is arranged as a matrix. Each matrix row
-        represents a certain block of audio (in chronological order). If the input signal is stereo, another matrix
-        dimension is added.
-    """
-    """ CONSTANTS AND VARIABLES
-
-    First, we are going to define some of the variables and constants that are needed in order to run our code, 
-    section by section.
-    """
-    # Signal dimensions
-    dim = signal.ndim
-
-    """ CONSIDERATION OF THRESHOLD IN QUIET """
-    # "The specific loudness in each band z is zero if it is at or below a critical-band-dependent specific loudness
-    # threshold LTQ(z)"
-    ltq_z = [
-        0.3310,
-        0.1625,
-        0.1051,
-        0.0757,
-        0.0576,
-        0.0453,
-        0.0365,
-        0.0298,
-        0.0247,
-        0.0207,
-        0.0176,
-        0.0151,
-        0.0131,
-        0.0115,
-        0.0103,
-        0.0093,
-        0.0086,
-        0.0081,
-        0.0077,
-        0.0074,
-        0.0073,
-        0.0072,
-        0.0071,
-        0.0072,
-        0.0073,
-        0.0074,
-        0.0076,
-        0.0079,
-        0.0082,
-        0.0086,
-        0.0092,
-        0.0100,
-        0.0109,
-        0.0122,
-        0.0138,
-        0.0157,
-        0.0172,
-        0.0180,
-        0.0180,
-        0.0177,
-        0.0176,
-        0.0177,
-        0.0182,
-        0.0190,
-        0.0202,
-        0.0217,
-        0.0237,
-        0.0263,
-        0.0296,
-        0.0339,
-        0.0398,
-        0.0485,
-        0.0622,
-    ]
-
-    """ AUDITORY FILTERING BANK """
-    # Order of the Outer and Middle ear filter
-    filter_order_k = 5
-    # Sampling frequency
-    fs = 48000.00
-    d_coefficients_array = []
-    # Step size
-    z_step_size = 0.50
-    # Centre frequencies F(z), 26.5/0.5 = 53
-    centre_freq_array = []
-    # F Bandwidth Af(z)
-    f_bandwidth_array = []
-    # Retardation
-    t_delay_array = []
-    # Block length and hop size lists
-    sb_array = []
-    sh_array = []
-    # Coefficients for the "Auditory Filtering Bank" part
-    am_mod_coefficient_array = np.zeros((53, int(filter_order_k + 1)), dtype=complex)
-    bm_mod_coefficient_array = np.zeros((53, int(filter_order_k + 1)), dtype=complex)
-    band_pass_signal_array = []
-
-    """ NON-LINEARITY. SOUND PRESSURE INTO SPECIFIC LOUDNESS """
-    p_0 = 2e-5
-    # c_N: In sones/bark
-    c_N = 0.0217406
-    alpha = 1.50
-    threshold_db_array = np.array([15.0, 25.0, 35.0, 45.0, 55.0, 65.0, 75.0, 85.0])
-    v_i_array = np.array(
-        [1.0, 0.6602, 0.0864, 0.6384, 0.0328, 0.4068, 0.2082, 0.3994, 0.6434]
-    )
-    M_exponents = len(threshold_db_array)
-    a_exponent_array, pt_threshold_array = nonlinear_common(
-        p_0, alpha, M_exponents, v_i_array, threshold_db_array
-    )
-
-    specific_loudness_all_bands_array = []
-
-    """ TOTAL LOUDNESS """
-    total_loudness_array = []
-
-    """ OUTER AND MIDDLE EAR FILTERING (5.1.2)
-
-    It is important to use the "filtfilt" version of the filter in order to reduce the lag in the results.
-    The Auditory Filtering Bank must be done by filtering with a non zero-phase filter in order to avoid filtering the
-    signal much more than it is intended to.
-    """
-    sos_ear = ear_filter_design()
-    # signal_filtered is pom(n). We have to assign the input signal to the filtered signal in order to implement
-    # the "sos" filtering
-    signal_filtered = signal
-    signal_filtered = sp.signal.sosfiltfilt(sos_ear, signal_filtered, axis=0)
-
-    """ AUDITORY FILTERING BANK (5.1.3)
-
-    In the following section, the signal is filtered by a series of 53 asymmetric and overlapping filters. They are 
-    supposed to replicate the activation process of the auditory hair cells, and their shape matches the one from the 
-    gammatone filters. First, we have to calculate the parameters (central frequency, bandwidth, delay, and "d" 
-    coefficient) that are going to define our filter and that are used in the recursive formula number 10. Second, 
-    we have to link the critical band rate scale "z", with a certain hop size "sh_array" and its band dependent block 
-    size "sb_array" (Table G.1). After that, it comes the coefficient calculation (Formula 11 and 12), in which we 
-    compute the values for the filter coefficients (*), and we store them in a list. Finally, the end of this section 
-    comes with the filtering of the signal.
-
-    (*) An error in the standard has been found in expressions 13 and 14. The prior paragraph says that in order to 
-    obtain the approximation of the band-pass filter, the low-pass filter coefficients shall be modified by adding a 
-    negative exponential with the transformation parameters "complex_exponential". If we try to develop the filter with 
-    these guidelines, we find that the resulting filter is a low-pass filter, not a band-pass. The actual way of 
-    obtaining a band-pass filter by transforming a low-pass filter is to multiply the low pass filter coefficients in 
-    time with a positive exponential. After the multiplication we will have complex filter coefficients, we have 2 
-    options in order to fix this. Either we make 2 transformations, one to the right and another one to the left, and 
-    we adapt the transformation parameter, or we take for granted the last transformation by discarding the complex part 
-    of the filter and multiplying by "2" (two transformations) and by a cosine with the original transformation 
-    parameter. We decided to implement the second option as they did on ECMA-418-2.
-    """
-    centre_freq = gen_auditory_filters_centre_freq()
-    for band_number in range(53):
-        (
-            bm_mod_coefficient_array[band_number],
-            am_mod_coefficient_array[band_number],
-        ) = gammatone(centre_freq[band_number], order=5, fs=fs)
-
-        """ 
-        "scipy.signal.lfilter" instead of "scipy.signal.filtfilt" in order to maintain consistency. That process 
-        makes possible to obtain a signal "band_pass_signal" that does not line up in time with the original signal 
-        because of the non zero-phase filtering of "lfilter", but it has a more appropriate slope than filtfilt. 
-        By using filtfilt the slope is that high that filters too much the signal. 
-        """
-        band_pass_signal = (
-            2.0
-            * (
-                sp.signal.lfilter(
-                    bm_mod_coefficient_array[band_number],
-                    am_mod_coefficient_array[band_number],
-                    signal_filtered,
-                    axis=0,
-                )
-            ).real
-        )
-
-        """ RECTIFICATION (5.1.5)
-
-        This part acts as the activation of the auditory nerves when the basilar membrane vibrates in a certain 
-        direction. In order to rectify the signal we are using "np.clip" which establish a minimum and a maximum value
-        for the signal. "a_min" is set to 0 float, while "a_max" is set to "None" in order to consider the positive 
-        value of the signal. 
-        """
-        band_pass_signal_hr = np.clip(band_pass_signal, a_min=0.00, a_max=None)
-        band_pass_signal_array.append(band_pass_signal_hr)
-
-        """ SEGMENTATION OF THE SIGNAL INTO BLOCKS (5.1.4)
-
-        The segmentation of the signal is done in order to obtain results for intervals of time, not for the whole 
-        duration of the signal. The reason behind this decision resides in the fact that processing the signal in its 
-        full length at one time could end up in imprecise results. By using a "for loop", we are able to decompose the 
-        signal array "band_pass_signal_hr" into blocks. "sb_array" is the block size which changes depending on the 
-        "band_number" in which we are processing the signal. "sh_array" is the step size, the time shift to the next 
-        block.
-        """
-        _, _, _, _, sb, sh = afb_ppal_parameters(fs, band_number, filter_order_k)
-        block_array = segmentation_blocks(band_pass_signal_hr, sb, sh, dim)
-
-        """ ROOT-MEAN-SQUARE (5.1.6)
-
-        After the segmentation of the signal into blocks, root-mean square values of each block are calculated in 
-        Formula 17). Some of the necessary steps that are general to all of the blocks have been priorly calculated 
-        in order to reach a better performance. 
-        We have to initialize here "a_array" (array of A values) and "specific_loudness_array" (array of specific 
-        loudness values), in order to treat each band-pass signal separately.
-        """
-        if signal.ndim == 2:
-            empty_loudness_array = np.zeros((1, 2), dtype=float)
-
-        else:
-            empty_loudness_array = np.zeros((1, 1), dtype=float)
-
-        # Auxiliary array list for loudness assignation before adding the first block
-        specific_loudness_aux = []
-
-        for l_number, l_block in enumerate(block_array):
-            # The next line makes possible to save the value of the first block that it is not going to be repeated
-            if l_number > 0:
-                # We have to specify the type of the list in order to avoid some problems reported with this library,
-                # where it results in having negative values as output
-                rms_factor = np.mean(np.power(l_block.astype(float), 2), axis=0)
-
-                # "rms_block_value" is the result in RMS-values of each block.
-                rms_block_value = np.sqrt(2.00 * rms_factor)
-
-                # NEXT PART REMAINS FOR TESTING PURPOSES - band number 17 is the one centered on 1 kHz (1027.0247 Hz)
-                # It makes possible to print the RMS value of the block number 10 at 1 kHz, as well as to plot its
-                # rectified signal.
-                """
-                if band_number == 17 and l_number == 10:
-                    print('\nTEST - RMS block value at 1027.0247 Hz\n\tRMS block value: ' + str(rms_block_value)
-                          + '\n\tCentral Frequency: ' + str(centre_freq)
-                          + '\n\tz value: ' + str(z_step_size * band_number))
-
-                    block_time = np.arange(len(l_block)) / fs
-                    block_x_limit = len(l_block) / fs
-
-                    plt.figure(figsize=(10, 5))
-                    plt.plot(block_time, l_block)
-                    plt.xlim(left=0, right=block_x_limit)
-                    plt.ylim(bottom=0)
-                    plt.xlabel('Time [s]')
-                    plt.ylabel('Amplitude [Pa]')
-                    plt.title('Rectified signal, block number 10, 1 kHz band')
-                    plt.grid(which='both', linestyle='-', color='grey')
-                    plt.show()
-                """
-
-                """ NON-LINEARITY (5.1.7)
-
-                This section covers the other part of the calculations needed to consider the non-linear transformation 
-                of sound pressure to specific loudness that does the the auditory system. After this point, the 
-                computation is done equally to every block in which we have divided our signal. The following loop makes 
-                possible the product of the sequence shown in Formula 18. "M_exponents" equals "M", which is the number 
-                of given exponents.
-                """
-                sequence_product = 1.0
-                for i_position in range(M_exponents):
-                    # "a_base" is used to store the value of the "i_position" element of the sequence
-                    a_base = (
-                        np.power(
-                            (rms_block_value / pt_threshold_array[i_position]), alpha
-                        )
-                        + 1.00
-                    )
-
-                    # The first loop we do not have a previous element in the sequence ("sequence_product"). On each
-                    # iteration the previous multiplication is multiplied by a new element
-                    sequence_product *= np.power(a_base, a_exponent_array[i_position])
-
-                # "a_value" is the variable used to store the resulting value after the appliance of the non-linearity.
-                # It is the same value as the "specific loudness" without consideration of the threshold in quiet
-                a_value = c_N * (rms_block_value / p_0) * sequence_product
-
-                """ SPECIFIC LOUDNESS CONSIDERING THE THRESHOLD IN QUIET (5.1.8)
-
-                The next calculation helps us obtain the result for the specific loudness - specific loudness with 
-                consideration of the lower threshold of hearing. After all the calculations that have been done, 
-                now its time to consider if the signal is stereo or mono. "signal.ndim" gives the number of dimensions 
-                that has the array in order to difference mono and stereo calculations. Also, we are using  one-line 
-                sentences in order to simplify the reading.
-
-                ** ATTENTION ** This part gives problems later on, if we leave this part as it is it converts negative 
-                results to zero. That conversion makes sense but it causes that the Figure 5 is not well represented,
-                it will give the same result for values that are lower than the threshold.
-                """
-                if signal.ndim == 2:
-                    left_value = (
-                        0.00
-                        if (a_value[0] < ltq_z[band_number])
-                        else a_value[0] - ltq_z[band_number]
-                    )
-                    right_value = (
-                        0.00
-                        if (a_value[1] < ltq_z[band_number])
-                        else a_value[1] - ltq_z[band_number]
-                    )
-                    specific_loudness_aux.append(np.array([left_value, right_value]))
-
-                else:
-                    mono_value = (
-                        0.00
-                        if (a_value < ltq_z[band_number])
-                        else a_value - ltq_z[band_number]
-                    )
-
-                    specific_loudness_aux.append(np.array([mono_value]))
-
-        # Next step line calls the function "nonlinear_repeat" that wil fill the necessary array slots in order to
-        # then rearrange all the bands.
-        # "np.insert" adds the empty array to the start of the blocks.
-        band_blocks_array = np.insert(
-            nonlinear_repeat(np.array(specific_loudness_aux), band_number),
-            0,
-            empty_loudness_array,
-            axis=0,
-        )
-        # print('Band: ' + str(band_number) + ', length: ' + str(len(band_blocks_array)))
-
-        # specific_loudness_all_bands_array.append(band_blocks_array)
-        specific_loudness_all_bands_array.append(specific_loudness_aux)
-
-    """
-    When we did the specific loudness calculation, we did it for all the blocks at a time, band per band. Therefore, 
-    we obtained "specific_loudness_allbands_array", which is an array list of specific loudness per bands. In order to 
-    obtain arrays in which the elements of every array are the results per band for that specific block of signal, 
-    we have to use the "zip" function. "zip" creates new arrays with elements that have the same index.
-    """
-    # n_array = np.array(
-    #     [
-    #         loudness_per_block
-    #         for loudness_per_block in zip(*specific_loudness_all_bands_array)
-    #     ]
-    # )
-    n_array = specific_loudness_all_bands_array
-
-    """ TOTAL LOUDNESS (5.1.8)
-
-    Inside the next "for loop", we are calculating the total loudness (Formula 21) for each block of audio in which 
-    we divided our signal. After that, we append that value (stereo or mono) to an array list "total_loudness_array". 
-    "n_block" is an array of all the N values of each band in that specific audio part (block). After the for loop and 
-    the creation of the array list "total_loudness_array", we have to convert that list to a numpy array in order to 
-    maintain stereo format when needed.
-    """
-    for n_block in n_array:
-        total_loudness_value = np.sum(n_block, axis=0) * z_step_size
-        total_loudness_array.append(total_loudness_value)
-
-    t_array = np.array(total_loudness_array)
-
-    """
-    On the one hand, if "validation" is set to default value (False) it does not calculate the complete validation 
-    steps. On the other hand, if "validation" is set to True, it runs the validation function.
-    """
-    if validation:
-        annex_f_validation(
-            fs,
-            signal,
-            signal_filtered,
-            sos_ear,
-            centre_freq_array,
-            f_bandwidth_array,
-            t_delay_array,
-            d_coefficients_array,
-            sb_array,
-            sh_array,
-            am_mod_coefficient_array,
-            bm_mod_coefficient_array,
-            band_pass_signal_array,
-            n_array,
-            t_array,
-        )
-
-    return n_array, t_array
diff --git a/mosqito/functions/loudness_ecma_spain/comp_loudness_alt.py b/mosqito/functions/loudness_ecma_spain/comp_loudness_alt.py
deleted file mode 100644
index 1922d638..00000000
--- a/mosqito/functions/loudness_ecma_spain/comp_loudness_alt.py
+++ /dev/null
@@ -1,95 +0,0 @@
-# -*- coding: utf-8 -*-
-"""
-@author: Daniel Jiménez-Caminero Costa
-"""
-import numpy as np
-import matplotlib.pyplot as plt
-
-# Project Imports
-from mosqito.functions.loudness_ecma_spain.rectified_band_pass_signals import (
-    rectified_band_pass_signals,
-)
-from mosqito.functions.loudness_ecma_spain.loudness_function import (
-    loudness_function as loudness_function,
-)
-
-# Data import
-# Threshold in quiet
-from mosqito.functions.loudness_ecma_spain.loudness_ecma_data import ltq_z
-
-
-def comp_loudness(signal, sb=2048, sh=1024):
-    """Calculation of specific loudness and total loudness of an input signal according to ECMA-418-2. It describes the
-    loudness excitation in a critical band per Bark. The output of the signal consists in a pair of array lists, one
-    for specific loudness and the other one for total loudness.
-
-    Parameters
-    ----------
-    signal: numpy.array
-        time signal valuesin 'Pa'. It can be, stereo (2 dimensions) or mono (1 dimension).The sampling frequency of the
-        signal must be 48000 Hz.
-    sb: int or list of int
-        block size.
-    sh: int or list of int
-        Hop size.
-
-    Returns
-    -------
-    n_array: numpy.array
-        'Sone per Bark'. It is a numpy array that is arranged as a matrix. Each matrix row represents a certain block of
-        audio (in chronological order) and each matrix column stores the value of the specific loudness for that block
-        and that specific band (column number). If the input signal is stereo, another matrix dimension is added.
-
-    t_array: numpy.array
-        'Sone per Bark'. As well as "n_array", it is a numpy array that is arranged as a matrix. Each matrix row
-        represents a certain block of audio (in chronological order). If the input signal is stereo, another matrix
-        dimension is added.
-    """
-
-    """Rectified band-pass signals (5.1.5)"""
-    # # sb and sh for Tonality
-    # z = np.linspace(0.5, 26.5, num=53, endpoint=True)
-    # sb = np.ones(53, dtype="int")
-    # sh = np.ones(53, dtype="int")
-    # sb[z <= 1.5] = 8192
-    # sh[z <= 1.5] = 2048
-    # sb[np.all([z >= 2, z <= 8], axis=0)] = 4096
-    # sh[np.all([z >= 2, z <= 8], axis=0)] = 1024
-    # sb[np.all([z >= 8.5, z <= 12.5], axis=0)] = 2048
-    # sh[np.all([z >= 8.5, z <= 12.5], axis=0)] = 512
-    # sb[z >= 13] = 1024
-    # sh[z >= 13] = 256
-
-    block_array_rect = rectified_band_pass_signals(signal, sb, sh)
-
-    n_array = []
-    for band_number in range(53):
-        """ROOT-MEAN-SQUARE (5.1.6)
-
-        After the segmentation of the signal into blocks, root-mean square values of each block are calculated in
-        Formula 17).
-        """
-
-        rms_block_value = np.sqrt(
-            2 * np.mean(np.array(block_array_rect[band_number]) ** 2, axis=1)
-        )
-
-        """NON-LINEARITY (5.1.7)
-
-        This section covers the other part of the calculations needed to consider the non-linear transformation
-        of sound pressure to specific loudness that does the the auditory system. After this point, the
-        computation is done equally to every block in which we have divided our signal.
-        """
-        a_prime = loudness_function(rms_block_value)
-
-        """SPECIFIC LOUDNESS CONSIDERING THE THRESHOLD IN QUIET (5.1.8)
-
-        The next calculation helps us obtain the result for the specific loudness - specific loudness with
-        consideration of the lower threshold of hearing. 
-        """
-
-        a_prime[a_prime < ltq_z[band_number]] = ltq_z[band_number]
-        N_prime = a_prime - ltq_z[band_number]
-        n_array.append(N_prime)
-
-    return n_array
diff --git a/mosqito/functions/loudness_ecma_spain/loudness_function.py b/mosqito/functions/loudness_ecma_spain/loudness_function.py
deleted file mode 100644
index 322feed0..00000000
--- a/mosqito/functions/loudness_ecma_spain/loudness_function.py
+++ /dev/null
@@ -1,93 +0,0 @@
-import numpy as np
-
-
-def loudness_function(p):
-
-    """ NON-LINEARITY. SOUND PRESSURE INTO SPECIFIC LOUDNESS """
-    p_0 = 2e-5
-    # c_N: In sones/bark
-    c_N = 0.0217406
-    alpha = 1.50
-    v_i = np.array(
-        [1.0, 0.6602, 0.0864, 0.6384, 0.0328, 0.4068, 0.2082, 0.3994, 0.6434]
-    )
-    thresh = np.array([15.0, 25.0, 35.0, 45.0, 55.0, 65.0, 75.0, 85.0])
-    p_ti = p_0 * 10 ** (thresh / 20)
-    # exp = (v_i[1:] - v_i[:-1]) / alpha
-    # prod = np.prod(
-    #     (1 + (p / p_ti[:, np.newaxis]) ** alpha)
-    #     ** exp[:, np.newaxis],
-    #     axis=0,
-    # )
-    # a_prime = c_N * p / p_0 * prod
-
-    a_prime = np.ones(p.shape)
-    for i in range(8):
-        a_prime *= (1 + (p / p_ti[i]) ** alpha) ** ((v_i[i + 1] - v_i[i]) / alpha)
-    a_prime *= c_N * p / p_0
-
-    return a_prime
-
-
-def loudness_function_alt(p):
-
-    p_ref = np.array(
-        [
-            6.598919,
-            8.838655,
-            11.250792,
-            14.17477,
-            17.785425,
-            20.533108,
-            23.620234,
-            26.878784,
-            31.16102,
-            35.784645,
-            40.063473,
-            44.85609,
-            48.621613,
-            53.92802,
-            58.548725,
-            64.02948,
-            70.02354,
-            76.3624,
-            82.87268,
-            87.32732,
-            92.63762,
-            97.607506,
-            101.3779,
-            104.46454,
-        ]
-    )
-
-    loud_ref = np.array(
-        [
-            0.002630106,
-            0.0051425043,
-            0.010632731,
-            0.021176357,
-            0.047157302,
-            0.076508306,
-            0.11306698,
-            0.17023043,
-            0.2377996,
-            0.32601807,
-            0.3997006,
-            0.50857294,
-            0.6120943,
-            0.8082804,
-            0.9908526,
-            1.4907696,
-            2.2847574,
-            3.9156787,
-            6.8367267,
-            10.09929,
-            15.481487,
-            25.099924,
-            36.40073,
-            52.800743,
-        ]
-    )
-
-    a_prime = np.interp(p, 2e-5 * 10 ** (p_ref / 20), loud_ref)
-    return a_prime
\ No newline at end of file
diff --git a/mosqito/functions/loudness_ecma_spain/nonlinear_common.py b/mosqito/functions/loudness_ecma_spain/nonlinear_common.py
deleted file mode 100644
index 25c936ec..00000000
--- a/mosqito/functions/loudness_ecma_spain/nonlinear_common.py
+++ /dev/null
@@ -1,54 +0,0 @@
-# -*- coding: utf-8 -*-
-"""
-@author: Daniel Jiménez-Caminero Costa
-"""
-
-import numpy as np
-import math
-
-
-def nonlinear_common(p_0, alpha, m_exponents, v_i_array, threshold_db_array):
-    """Array lists that are necessary for the calculation of the non-linearity and are general to each band. This has
-    been implemented as described in formula F.18 (section F.3.6) of Annex F (ECMA-74).
-
-    Parameters
-    ----------
-    p_0: numpy.array
-        '20 uPa'.
-
-    alpha: float
-        Constant for the exponent.
-
-    m_exponents: int
-        Max index of the multiplication sequence in Formula 18.
-
-    v_i_array: numpy.array
-        Exponents for the multiplication sequence in Formula 18.
-
-    threshold_db_array: numpy.array
-        Thresholds for their corresponding "vi" exponent.
-
-    Returns
-    -------
-    a_exponent_array: numpy.array
-
-    pt_threshold_array: numpy.array
-        'dB'
-    """
-
-    print("This function is deprecated and will be removed form the project")
-
-    pt_threshold_array = np.zeros(m_exponents, dtype=float)
-
-    # Numpy array for the exponent in the non-linearity function
-    a_exponent_array = np.array(np.diff(v_i_array) / alpha)
-
-    # COMMON CALCULATIONS
-    for i_position in range(m_exponents):
-        # "pt_threshold" changes to the value of each pt threshold (Table F.2)
-        th_exponent = threshold_db_array[i_position] / 20
-        pt_threshold = p_0 * math.pow(10, th_exponent)
-        # Numpy array for the threshold in the non-linearity function
-        pt_threshold_array[i_position] = pt_threshold
-
-    return a_exponent_array, pt_threshold_array
diff --git a/mosqito/functions/loudness_ecma_spain/nonlinear_repeat.py b/mosqito/functions/loudness_ecma_spain/nonlinear_repeat.py
deleted file mode 100644
index 391f6983..00000000
--- a/mosqito/functions/loudness_ecma_spain/nonlinear_repeat.py
+++ /dev/null
@@ -1,49 +0,0 @@
-# -*- coding: utf-8 -*-
-"""
-@author: Daniel Jiménez-Caminero Costa
-"""
-
-import numpy as np
-
-
-def nonlinear_repeat(specific_loudness_array, band_number):
-    """It repeats certain blocks depending on the band that you are in, in order pseudo-link bands on time, because
-    lower bands have larger blocks than the higher ones.
-
-    Parameters
-    ----------
-    specific_loudness_array: numpy.array
-        'sones hms' array of specific loudness values.
-
-    band_number: int
-        Band number in which the filter is going to be applied over the signal.
-
-    Parameters
-    ----------
-
-    Returns
-    -------
-    sl_array_rep: numpy.array
-        'sones hms'.
-    """
-
-    print("This function is deprecated and will be removed form the project")
-
-    # Critical band rate scale
-    z = (band_number + 1) * 0.50
-
-    # It makes possible to create arrays of repeated values based on their critical band rate scale, and consequently
-    # on their block size
-    if z >= 13:
-        sl_array_rep = specific_loudness_array
-
-    elif 8.5 <= z < 13:
-        sl_array_rep = np.repeat(specific_loudness_array, repeats=2, axis=0)
-
-    elif 2 <= z < 8.5:
-        sl_array_rep = np.repeat(specific_loudness_array, repeats=4, axis=0)
-
-    else:
-        sl_array_rep = np.repeat(specific_loudness_array, repeats=8, axis=0)
-
-    return sl_array_rep
diff --git a/mosqito/functions/loudness_ecma_spain/segmentation_blocks.py b/mosqito/functions/loudness_ecma_spain/segmentation_blocks.py
deleted file mode 100644
index 7d8dc77e..00000000
--- a/mosqito/functions/loudness_ecma_spain/segmentation_blocks.py
+++ /dev/null
@@ -1,166 +0,0 @@
-# -*- coding: utf-8 -*-
-"""
-@author: Daniel Jiménez-Caminero Costa
-"""
-
-import numpy as np
-
-
-def segmentation_blocks(band_pass_signal_hr, sb, sh, dim):
-    """Function used for the segmentation of the signal into smaller parts of audio (blocks). This has been implemented
-    as described in Formula 16 (section 5.1.4) of ECMA-418-2.
-
-    Parameters
-    ----------
-    band_pass_signal_hr: numpy.array
-        time signal values of the input signal
-
-    sb: int
-        block size
-
-    sh: int
-        Hop size.
-
-    dim: int
-        Signal dimensions.
-
-    Returns
-    -------
-    block_array: List[None]
-        Array list of blocks in which the signal has been segmented.
-    """
-
-    l = 0
-    block_array = []
-    while l * sh <= len(band_pass_signal_hr):
-        i_begin = l * sh - sb
-        i_end = l * sh
-        if i_begin < 0 and i_end <= 0:
-            block_array.append(np.zeros(sb))
-        elif i_begin < 0:
-            block_array.append(
-                np.concatenate((np.zeros(-i_begin), band_pass_signal_hr[0:i_end]))
-            )
-        else:
-            block_array.append(band_pass_signal_hr[i_begin:i_end])
-        l += 1
-
-    # # Creation of the first empty block
-    # if dim == 2:
-    #     first_block = np.zeros((sb, dim), dtype=float)
-    # else:
-    #     first_block = np.zeros(sb, dtype=float)
-
-    # # Joins the first empty block with the rest of the signal
-    # band_pass_signal_hr_complete = np.concatenate(
-    #     (first_block, band_pass_signal_hr), axis=0
-    # )
-
-    # # For loop that helps to cut the signal into different blocks
-    # block_array = [
-    #     band_pass_signal_hr_complete[block_s : block_s + sb]
-    #     for block_s in range(0, len(band_pass_signal_hr_complete), sh)
-    # ]
-
-    return block_array
-
-
-def segmentation_blocks_test_a(band_pass_signal_hr, sb, sh, dim):
-    """THIS SECTION REMAINS FOR TESTING PURPOSES
-
-    Function used for the segmentation of the signal into smaller parts of audio (blocks). This implementation
-    has been modified in order to make a simpler version of the original one. After the tests, it has been proven that
-    this algorithm is way faster than the original, it reduces some uncertainty and the number of blocks remains the
-    same in most cases. Also, it treats the void that leaves the standard in blocks that have samples that ar out
-    of bounds.
-
-    Parameters
-    ----------
-    band_pass_signal_hr: numpy.array
-        Time signal values of the input signal.
-
-    sb: int
-        Block size.
-
-    sh: int
-        Hop size.
-
-    dim: int
-        Signal dimensions.
-
-    Returns
-    -------
-    block_array: List[List[float]]
-        Array list of blocks in which the signal has been segmented.
-    """
-    # Creation of the Array list of blocks with a specific block size (sb) and hop size (sh)
-    block_array = [
-        band_pass_signal_hr[block_s : block_s + sb]
-        for block_s in range(0, len(band_pass_signal_hr), sh)
-    ]
-
-    return block_array
-
-
-def segmentation_blocks_test_b(band_pass_signal_hr, sb, sh, dim):
-    """THIS SECTION REMAINS FOR TESTING PURPOSES
-
-    Function used for the segmentation of the signal into smaller parts of audio (blocks). This has been implemented
-    as described in formula F.16 (section F.3.5) of Annex F (ECMA-74). It has some issues  when the input signal is
-    stereo.
-
-    Parameters
-    ----------
-    band_pass_signal_hr: numpy.array
-        Time signal values of the input signal.
-
-    sb: int
-        Block size.
-
-    sh: int
-        Hop size.
-
-    dim: int
-        Signal dimensions.
-
-    Returns
-    -------
-    block_array: List[List[float]]
-        Array list of blocks in which the signal has been segmented.
-    """
-    # Array list of blocks
-    block_array = []
-
-    # Assignation of the empty  value for the block cell
-    if band_pass_signal_hr.ndim == 2:
-        empty_value = [0.00, 0.00]
-    else:
-        empty_value = 0.00
-
-    """
-    Total number of blocks. The sequence goes from "block_number = 0" to "block_number = total_blocks - 1". The result 
-    of the number of blocks is rounded in order to prevent an error in the value assignment. Another option could be to 
-    limit the call to higher positions in the array. The standard does not treat as it should the segmentation.
-    """
-    total_blocks = round((len(band_pass_signal_hr) + 1) / sh)
-
-    for block_number in range(total_blocks):
-        # Creation of the block
-        block = np.zeros((sb, band_pass_signal_hr.ndim), dtype=float)
-
-        # Here is stated the range for the assignment
-        p_range = sb - (block_number * sh)
-        # Assignment position
-        p_segmented = (block_number * sh) - sb
-
-        for n_segmentation in range(sb):
-            block[n_segmentation] = (
-                band_pass_signal_hr[p_segmented + n_segmentation]
-                if (n_segmentation >= p_range)
-                else empty_value
-            )
-
-        # Block assignment
-        block_array.append(block)
-
-    return block_array
diff --git a/mosqito/functions/loudness_ecma_spain/sl_validation.py b/mosqito/functions/loudness_ecma_spain/sl_validation.py
deleted file mode 100644
index ac0ecd87..00000000
--- a/mosqito/functions/loudness_ecma_spain/sl_validation.py
+++ /dev/null
@@ -1,479 +0,0 @@
-# -*- coding: utf-8 -*-
-"""
-@author: Daniel Jiménez-Caminero Costa
-"""
-import numpy as np
-import math
-import scipy as sp
-import matplotlib.pyplot as plt
-from scipy.signal import welch
-import sys
-sys.path.append('../../..')
-
-
-def out_mid_validation(fs, signal, signal_filtered, sos_ear):
-    """ Function that serves for the validation of the outer and middle ear filtering  section.
-
-    Parameters
-    ----------
-    fs: float
-        'Hz', sampling frequency.
-
-    signal: numpy.array
-        Time signal values of the input signal.
-
-    signal_filtered: numpy.array
-        Time signal values of the filtered signal.
-
-    sos_ear: List[List[float]]
-        "sos coefficients".
-
-    Returns
-    -------
-
-    """
-    dim = signal.ndim
-    time = np.arange(len(signal)) / fs
-    x_limit = len(signal) / fs
-
-    print('Outer & Middle Ear Filter Coefficients:\n' + str(sos_ear) + '\n')
-
-    """ GRAPHIC REPRESENTATION """
-    # Frequency response 'sos'
-    plt.figure(figsize=(10, 5))
-    w_om, h_om = sp.signal.sosfreqz(sos_ear, worN=round(fs / 2))
-    w_om_hertz = w_om * (fs / (2.0 * math.pi))
-    h_om_db = 20.0 * np.log10(np.maximum(np.abs(h_om), 1e-7))
-    plt.semilogx(w_om_hertz, h_om_db)
-    # Legend created following the guidelines from:
-    #   https://stackoverflow.com/questions/4700614/how-to-put-the-legend-out-of-the-plot
-    plt.legend(['F. Response (sos filter)'], bbox_to_anchor=(1.04, 0.5), borderaxespad=0, loc='center left', ncol=1)
-    plt.subplots_adjust(right=0.75)
-
-    plt.xlim(20, 20000)
-    plt.ylim(-25, 10)
-    plt.xlabel('Frequency [Hz]')
-    plt.ylabel('Amplitude [dB]')
-    plt.title('Frequency Response of Outer and Middle Ear Filter')
-    plt.grid(which='both', linestyle='-', color='grey')
-    plt.xticks([20, 50, 100, 200, 500, 1000, 2000, 5000, 10000, 20000],
-               ["20", "50", "100", "200", "500", "1K", "2K", "5K", "10K", "20K"])
-    plt.show()
-
-    if dim == 2:
-        # Left channel
-        plt.figure(figsize=(10, 5))
-        plt.plot(time, signal_filtered[:, 0], 'r', label='Filtered signal after OMF (L channel)')
-        plt.plot(time, signal[:, 0], 'b', label='Input signal (L channel)')
-        plt.legend(bbox_to_anchor=(1.04, 0.5), borderaxespad=0, loc='center left', ncol=1)
-        plt.subplots_adjust(right=0.68)
-
-        plt.xlim(0, len(signal_filtered))
-        plt.xlabel('Time [s]')
-        plt.ylabel('Amplitude [Pa]')
-        plt.title('Comparison between filtered signal and input signal (Left channel)')
-        plt.grid(which='both', linestyle='-', color='grey')
-        plt.show()
-
-        # Right channel
-        plt.figure(figsize=(10, 5))
-        plt.plot(time, signal_filtered[:, 1], 'r', label='Filtered signal after OMF (R channel)')
-        plt.plot(time, signal[:, 1], 'b', label='Input signal (R channel)')
-        plt.legend(bbox_to_anchor=(1.04, 0.5), borderaxespad=0, loc='center left', ncol=1)
-        plt.subplots_adjust(right=0.68)
-
-        plt.xlim(0, len(signal_filtered))
-        plt.xlabel('Time [s]')
-        plt.ylabel('Amplitude [Pa]')
-        plt.title('Comparison between filtered signal and input signal (Right channel)')
-        plt.grid(which='both', linestyle='-', color='grey')
-        plt.show()
-
-    elif dim == 1:
-        # Mono channel
-        plt.figure(figsize=(10, 5))
-        plt.plot(time, signal_filtered, 'r', label='Filtered signal after OMF (Mono)')
-        plt.plot(time, signal, 'b', label='Input signal (Mono)')
-        plt.legend(bbox_to_anchor=(1.04, 0.5), borderaxespad=0, loc='center left', ncol=1)
-        plt.subplots_adjust(right=0.68)
-
-        plt.xlim(0, x_limit)
-        plt.xlabel('Time [s]')
-        plt.ylabel('Amplitude [Pa]')
-        plt.title('Comparison between filtered signal and input signal (Mono channel)')
-        plt.grid(which='both', linestyle='-', color='grey')
-        plt.show()
-
-
-def afb_validation(fs, dim, centre_freq_array, f_bandwidth_array, t_delay_array, d_coefficients_array, sb_array,
-                   sh_array, am_mod_coefficient_array, bm_mod_coefficient_array, band_pass_signal_array):
-    """ Function that serves for the validation of the Auditory Filtering Bank.
-
-    Parameters
-    ----------
-    fs: float
-        'Hz', sampling frequency.
-
-    dim: int
-        Number of dimensions that has the input signal.
-
-    centre_freq_array: List[float]
-        Central frequencies for the auditory filtering bank.
-
-    f_bandwidth_array: List[float]
-        Bandwidth for the auditory filtering bank.
-
-    t_delay_array: List[float]
-        Time constant (delay) for the auditory filtering bank.
-
-    d_coefficients_array: List[float]
-        "d" coefficients for the auditory filtering bank.
-
-    sb_array: List[int]
-        Block size.
-
-    sh_array: List[int]
-        Hop size.
-
-    am_mod_coefficient_array: numpy.array
-        Real part of the "am" coefficients for the auditory filtering bank.
-
-    bm_mod_coefficient_array: numpy.array
-        Real part of the "bm" coefficients for the auditory filtering bank.
-
-    band_pass_signal_array: numpy.array
-        Band-pass filtered signals.
-
-    Returns
-    -------
-
-    """
-    # Representative band-pass signals that we want to display
-    rep_bands = np.array([4, 9, 17, 27, 37, 47])
-    time = np.arange(len(band_pass_signal_array[0])) / fs
-    x_limit = len(band_pass_signal_array[0]) / fs
-
-    print('\n\n_____________AUDITORY FILTERING BANK_____________\n')
-    print('Central Frequencies:\n' + str(centre_freq_array) + '\n')
-    print('Filter Bandwidths:\n' + str(f_bandwidth_array) + '\n')
-    print('T. Delay:\n' + str(t_delay_array) + '\n')
-    print('D coefficients:\n' + str(d_coefficients_array) + '\n')
-    print('Sb values - block size:\n' + str(sb_array) + '\n')
-    print('Sh values - hop size:\n' + str(sh_array) + '\n')
-    print('Am coefficients - band-pass filter:\n' + str(am_mod_coefficient_array) + '\n')
-    print('Bm coefficients - band-pass filter:\n' + str(bm_mod_coefficient_array) + '\n')
-
-    """ GRAPHIC REPRESENTATION """
-    # Frequency response of the filters
-    plt.figure(figsize=(17.5, 7.5))
-    for afb_band in range(0, 53):
-        band_number = afb_band + 1
-        w_afb, h_afb = sp.signal.freqz(bm_mod_coefficient_array[afb_band],
-                                       am_mod_coefficient_array[afb_band], worN=round(fs/2))
-        w_afb_hertz = w_afb * (fs / (2.0 * math.pi))
-        h_afb_db = 20.0 * np.log10(np.abs(h_afb))
-        plt.semilogx(w_afb_hertz, h_afb_db,
-                     label='FR band: ' + str(band_number)
-                           + ' (' + str(round(float(centre_freq_array[afb_band]), 2)) + 'Hz filter)')
-
-    plt.legend(bbox_to_anchor=(1.04, 1), borderaxespad=0, loc='upper left', ncol=2)
-    plt.subplots_adjust(right=0.65)
-
-    plt.xlim(right=20000)
-    plt.ylim(bottom=-250)
-    plt.xlabel('Frequency [Hz]')
-    plt.ylabel('Amplitude [dB]')
-    plt.title('Auditory Filtering Bank')
-    plt.grid(which='both', linestyle='-', color='grey')
-    plt.xticks([20, 50, 100, 200, 500, 1000, 2000, 5000, 10000, 20000],
-               ["20", "50", "100", "200", "500", "1K", "2K", "5K", "10K", "20K"])
-    plt.show()
-
-    if dim == 2:
-        # Left channel
-        plt.figure(figsize=(10, 5))
-        [plt.plot(time, band_pass_signal_array[band_num, :, 0],
-                  label='Band-pass signal ' + str(band_num + 1) + ' - '
-                        + str(round(float(centre_freq_array[band_num]), 2)) + 'Hz)') for band_num in rep_bands]
-        plt.legend(bbox_to_anchor=(1.04, 0.5), borderaxespad=0, loc='center left', ncol=1)
-        plt.subplots_adjust(right=0.65)
-
-        plt.xlim(0, x_limit)
-        plt.ylim(bottom=0)
-        plt.xlabel('Time [s]')
-        plt.ylabel('Amplitude [Pa]')
-        plt.title('Band-pass signals after filtering with AFB (Left channel)')
-        plt.grid(which='both', linestyle='-', color='grey')
-        plt.show()
-
-        # Right channel
-        plt.figure(figsize=(10, 5))
-        [plt.plot(time, band_pass_signal_array[band_num, :, 1],
-                  label='Band-pass signal ' + str(band_num + 1) + ' - '
-                        + str(round(float(centre_freq_array[band_num]), 2)) + 'Hz)') for band_num in rep_bands]
-        plt.legend(bbox_to_anchor=(1.04, 0.5), borderaxespad=0, loc='center left', ncol=1)
-        plt.subplots_adjust(right=0.65)
-
-        plt.xlim(0, x_limit)
-        plt.ylim(bottom=0)
-        plt.xlabel('Time [s]')
-        plt.ylabel('Amplitude [Pa]')
-        plt.title('Band-pass signals after filtering with AFB (Right channel)')
-        plt.grid(which='both', linestyle='-', color='grey')
-        plt.show()
-
-    elif dim == 1:
-        # Mono channel
-        plt.figure(figsize=(10, 5))
-        [plt.plot(time, band_pass_signal_array[band_num],
-                  label='Band-pass signal ' + str(band_num + 1) + ' - '
-                        + str(round(float(centre_freq_array[band_num]), 2)) + 'Hz)') for band_num in rep_bands]
-        plt.legend(bbox_to_anchor=(1.04, 0.5), borderaxespad=0, loc='center left', ncol=1)
-        plt.subplots_adjust(right=0.65)
-
-        plt.xlim(0, x_limit)
-        plt.ylim(bottom=0)
-        plt.xlabel('Time [s]')
-        plt.ylabel('Amplitude [Pa]')
-        plt.title('Band-pass signals after filtering with AFB (Mono channel)')
-        plt.grid(which='both', linestyle='-', color='grey')
-        plt.show()
-
-
-def sl_validation(dim, centre_freq_array, n_array):
-    """ Function that serves for the validation of the specific loudness.
-
-    Parameters
-    ----------
-    dim: int
-        Number of dimensions that has the input signal.
-
-    centre_freq_array: List[float]
-        Central frequencies for the auditory filtering bank.
-
-    n_array: numpy.array
-        Specific loudness.
-
-    Returns
-    -------
-
-    """
-    # Representative blocks that we want to display
-    n_blocks = len(n_array)
-    rep_blocks = np.array([str(round(n_blocks / 12)), str(round(n_blocks / 8)), str(round(n_blocks / 4)),
-                           str(round(n_blocks / 2)), str(int(float((n_blocks * 3) / 4)))])
-
-    """ EXTRACTION OF PARAMETERS & GRAPHIC REPRESENTATION """
-    # If n_array has 2 dimensions (stereo), we have to extract the data dimension by dimension
-    if dim == 2:
-        # Left channel
-        print('Specific loudness values (Left channel, stereo signal):\n' + str(n_array[:, :, 0]) + '\n')
-        np.savetxt("specific_loudness_left.txt", n_array[:, :, 0], fmt='%f',
-                   delimiter='\t')
-
-        plt.figure(figsize=(10, 5))
-        [plt.semilogx(centre_freq_array, n_array[int(block_num), :, 0],
-                      label='SL values, block: ' + str(int(block_num + 1))) for block_num in rep_blocks]
-        plt.legend(bbox_to_anchor=(1.04, 0.5), borderaxespad=0, loc='center left', ncol=1)
-        plt.subplots_adjust(right=0.75)
-
-        plt.xlim(left=50, right=15000)
-        plt.ylim(bottom=0)
-        plt.xlabel('Frequency [Hz]')
-        plt.ylabel('Sound pressure [sones]')
-        plt.title('Specific loudness contours (Left channel), block nº: ')
-        plt.grid(which='both', linestyle='-', color='grey')
-        plt.xticks([50, 100, 200, 500, 1000, 2000, 5000, 10000, 15000],
-                   ["50", "100", "200", "500", "1K", "2K", "5K", "10K", "15K"])
-        plt.show()
-
-        # Right channel
-        print('Specific loudness values (Right channel, stereo signal):\n' + str(n_array[:, :, 1]) + '\n')
-        np.savetxt("specific_loudness_right.txt", n_array[:, :, 1], fmt='%f',
-                   delimiter='\t')
-
-        plt.figure(figsize=(10, 5))
-        [plt.semilogx(centre_freq_array, n_array[int(block_num), :, 1],
-                      label='SL values, block: ' + str(int(block_num + 1))) for block_num in rep_blocks]
-        plt.legend(bbox_to_anchor=(1.04, 0.5), borderaxespad=0, loc='center left', ncol=1)
-        plt.subplots_adjust(right=0.75)
-
-        plt.xlim(left=50, right=15000)
-        plt.ylim(bottom=0)
-        plt.xlabel('Frequency [Hz]')
-        plt.ylabel('Sound pressure [sones]')
-        plt.title('Specific loudness contours (Right channel), block nº: ')
-        plt.grid(which='both', linestyle='-', color='grey')
-        plt.xticks([50, 100, 200, 500, 1000, 2000, 5000, 10000, 15000],
-                   ["50", "100", "200", "500", "1K", "2K", "5K", "10K", "15K"])
-        plt.show()
-
-    elif dim == 1:
-        # Mono channel
-        print('Specific loudness values (Mono signal):\n' + str(n_array) + '\n')
-        np.savetxt("specific_loudness.txt", n_array[:, :, 0], fmt='%f',
-                   delimiter='\t')
-
-        plt.figure(figsize=(10, 5))
-        [plt.semilogx(centre_freq_array, n_array[int(block_num), :, 0],
-                      label='SL values, block: ' + str(int(block_num))) for block_num in rep_blocks]
-        plt.legend(bbox_to_anchor=(1.04, 0.5), borderaxespad=0, loc='center left', ncol=1)
-        plt.subplots_adjust(right=0.75)
-
-        plt.xlim(left=50, right=15000)
-        plt.ylim(bottom=0)
-        plt.xlabel('Frequency [Hz]')
-        plt.ylabel('Sound pressure [sones]')
-        plt.title('Specific loudness contours (Mono channel), block nº: ')
-        plt.grid(which='both', linestyle='-', color='grey')
-        plt.xticks([50, 100, 200, 500, 1000, 2000, 5000, 10000, 15000],
-                   ["50", "100", "200", "500", "1K", "2K", "5K", "10K", "15K"])
-        plt.show()
-
-
-def tl_validation(dim, t_array):
-    """ Function that serves for the validation of the total loudness.
-
-    Parameters
-    ----------
-    dim: int
-        Number of dimensions that has the input signal.
-
-    t_array: numpy.array
-        Total loudness.
-
-    Returns
-    -------
-
-    """
-    print('Total loudness values:\n' + str(t_array) + '\n')
-
-    """ EXTRACTION OF PARAMETERS & GRAPHIC REPRESENTATION """
-    # If n_array has 2 dimensions (stereo), we have to extract the data dimension by dimension
-    if dim == 2:
-        # Left channel
-        print('Total loudness values (Left channel, stereo signal):\n' + str(t_array[:, 0]) + '\n')
-        np.savetxt("total_loudness_left.txt", t_array[:, 0], fmt='%f',
-                   delimiter='\t')
-
-        plt.figure(figsize=(10, 5))
-        plt.plot(t_array[:, 0])
-        plt.legend(['Total loudness (signal)'], bbox_to_anchor=(1.04, 0.5), borderaxespad=0, loc='center left', ncol=1)
-        plt.subplots_adjust(right=0.75)
-
-        plt.xlim(left=0, right=len(t_array))
-        plt.ylim(bottom=0)
-        plt.xlabel('Block number')
-        plt.ylabel('Sound pressure [sones]')
-        plt.title('Total loudness over time (Left channel)')
-        plt.grid(which='both', linestyle='-', color='grey')
-        plt.show()
-
-        # Right channel
-        print('Total loudness values (Right channel, stereo signal):\n' + str(t_array[:, 0]) + '\n')
-        np.savetxt("total_loudness_right.txt", t_array[:, 1], fmt='%f',
-                   delimiter='\t')
-
-        plt.figure(figsize=(10, 5))
-        plt.plot(t_array[:, 0])
-        plt.legend(['Total loudness (signal)'], bbox_to_anchor=(1.04, 0.5), borderaxespad=0, loc='center left', ncol=1)
-        plt.subplots_adjust(right=0.75)
-
-        plt.xlim(left=0, right=len(t_array))
-        plt.ylim(bottom=0)
-        plt.xlabel('Block number')
-        plt.ylabel('Sound pressure [sones]')
-        plt.title('Total loudness over time (Right channel)')
-        plt.grid(which='both', linestyle='-', color='grey')
-        plt.show()
-
-    elif dim == 1:
-        print('Total loudness values (Mono signal):\n' + str(t_array) + '\n')
-        np.savetxt("total_loudness.txt", t_array, fmt='%f', delimiter='\t')
-
-        plt.figure(figsize=(10, 5))
-        plt.plot(t_array[:, 0])
-        plt.legend(['Total loudness (signal)'], bbox_to_anchor=(1.04, 0.5), borderaxespad=0, loc='center left', ncol=1)
-        plt.subplots_adjust(right=0.75)
-
-        plt.xlim(left=0, right=len(t_array))
-        plt.ylim(bottom=0)
-        plt.xlabel('Block number')
-        plt.ylabel('Sound pressure [sones]')
-        plt.title('Total loudness over time (Mono channel)')
-        plt.grid(which='both', linestyle='-', color='grey')
-        plt.show()
-
-
-def annex_f_validation(fs, signal, signal_filtered, sos_ear, centre_freq_array, f_bandwidth_array, t_delay_array,
-                       d_coefficients_array, sb_array, sh_array, am_mod_coefficient_array, bm_mod_coefficient_array,
-                       band_pass_signal_array, n_array, t_array):
-    """ Function that serves for the validation of the principal parameters presented in Annex F of ECMA-74. As well as
-    the general validation, here it is also verified the final result for specific and total loudness.
-
-    Parameters
-    ----------
-    fs: float
-        'Hz', sampling frequency.
-
-    signal: numpy.array
-        Time signal values of the input signal.
-
-    signal_filtered: numpy.array
-        Time signal values of the filtered signal.
-
-    sos_ear: List[List[float]]
-        "sos coefficients".
-
-    centre_freq_array: List[float]
-        Central frequencies for the auditory filtering bank.
-
-    f_bandwidth_array: List[float]
-        Bandwidth for the auditory filtering bank.
-
-    t_delay_array: List[float]
-        Time constant (delay) for the auditory filtering bank.
-
-    d_coefficients_array: List[float]
-        "d" coefficients for the auditory filtering bank.
-
-    sb_array: List[int]
-        Block size.
-
-    sh_array: List[int]
-        Hop size.
-
-    am_mod_coefficient_array: numpy.array
-        Real part of the "am" coefficients for the auditory filtering bank.
-
-    bm_mod_coefficient_array: numpy.array
-        Real part of the "bm" coefficients for the auditory filtering bank.
-
-    band_pass_signal_array: numpy.array
-        Band-pass filtered signals.
-
-    n_array: numpy.array
-        Specific loudness.
-
-    t_array: numpy.array
-        Total loudness.
-
-    Returns
-    -------
-
-    """
-    signal_dimension = signal.ndim
-    print("Mono signal (1) / Stereo signal (2): " + str(signal_dimension) + '\n')
-
-    """ OUTER AND MIDDLE EAR FILTERING """
-    out_mid_validation(fs, signal, signal_filtered, sos_ear)
-
-    """ AUDITORY FILTERING BANK """
-    afb_validation(fs, signal_dimension, centre_freq_array, f_bandwidth_array, t_delay_array, d_coefficients_array,
-                   sb_array, sh_array, am_mod_coefficient_array, bm_mod_coefficient_array, band_pass_signal_array)
-
-    """ SPECIFIC LOUDNESS WITH CONSIDERATION OF THRESHOLD IN QUIET """
-    sl_validation(signal_dimension, centre_freq_array, n_array)
-
-    """ TOTAL LOUDNESS """
-    tl_validation(signal_dimension, t_array)
diff --git a/mosqito/functions/loudness_ecma_spain/sine_wave_generator.py b/mosqito/functions/shared/sine_wave_generator.py
similarity index 78%
rename from mosqito/functions/loudness_ecma_spain/sine_wave_generator.py
rename to mosqito/functions/shared/sine_wave_generator.py
index 67af0d74..438f57a1 100644
--- a/mosqito/functions/loudness_ecma_spain/sine_wave_generator.py
+++ b/mosqito/functions/shared/sine_wave_generator.py
@@ -1,29 +1,21 @@
 # -*- coding: utf-8 -*-
-"""
-@author: Daniel Jiménez-Caminero Costa
-"""
+
 import numpy as np
 import math
-import sys
-
-sys.path.append('../../..')
 
 
 def sine_wave_generator(fs, t, spl_value, freq):
-    """ It creates a sine wave signal given some input parameters like frequency, duration, sampling rate or sound
+    """It creates a sine wave signal given some input parameters like frequency, duration, sampling rate or sound
     pressure level.
 
     Parameters
     ----------
     fs: int
         'Hz', sampling frequency.
-
     t: int
         's', signal duration.
-
     spl_value: int
         'dB SPL', signal sound pressure level.
-
     freq: int
         'Hz', sine wave frequency.
 
@@ -31,17 +23,17 @@ def sine_wave_generator(fs, t, spl_value, freq):
     -------
     signal: numpy.array
         'Pa', time signal values. For ECMA-418-2 the sampling frequency of the signal must be 48000 Hz.
-
     time: numpy.array
         Time scale arranged in a numpy array.
     """
+
     # "Peak" value in Pascals (amplitude)
     p_ref = 2e-5
     pressure_rms = p_ref * (10.00 ** (spl_value / 20.00))
 
     # Sample range
     # samples = np.linspace(0, t, int(fs * t), endpoint=False)
-    time = np.arange(0, t, 1/fs)
+    time = np.arange(0, t, 1 / fs)
 
     # Theta lets you specify the sine wave value at time 0
     theta = 0
diff --git a/validations/loudness_ecma/hearing_model_validation.py b/validations/loudness_ecma/hearing_model_validation.py
index 42398d03..92f41f18 100644
--- a/validations/loudness_ecma/hearing_model_validation.py
+++ b/validations/loudness_ecma/hearing_model_validation.py
@@ -1,215 +1,109 @@
-# -*- coding: utf-8 -*-
-"""
-@author: Daniel Jiménez-Caminero Costa
-"""
-import numpy as np
-import math
-import scipy as sp
-import matplotlib.pyplot as plt
-from scipy.signal import welch
-
-# Project Imports
-from mosqito.functions.loudness_ecma_spain.sone2phone import sone2phone
-from mosqito.functions.loudness_ecma_spain.phone2spl import phone2spl
-from mosqito.functions.loudness_ecma_spain.comp_loudness import comp_loudness
-from mosqito.functions.loudness_ecma_spain.equal_loudness_contours import (
-    equal_loudness_contours,
-)
-from mosqito.functions.loudness_ecma_spain.sine_wave_generator import (
-    sine_wave_generator,
-)
-
-import sys
-
-sys.path.append("../../..")
-
-
-def hearing_model_validation():
-    """Function that serves for the validation of the hearing model presented in Annex F of ECMA-74.
-
-    Parameters
-    ----------
-
-    Returns
-    -------
-
-    """
-    # Duration of the signal
-    duration = 1
-    # Sampling frequency
-    fs = 48000.0
-
-    # Array of frequencies for loudness contours
-    freq_array = np.array(
-        [
-            20.0,
-            25.0,
-            31.5,
-            40.0,
-            50.0,
-            63.0,
-            80.0,
-            100.0,
-            125.0,
-            160.0,
-            200.0,
-            250.0,
-            315.0,
-            400.0,
-            500.0,
-            630.0,
-            800.0,
-            1000.0,
-            1250.0,
-            1600.0,
-            2000.0,
-            2500.0,
-            3150.0,
-            4000.0,
-            5000.0,
-            6300.0,
-            8000.0,
-            10000.0,
-            12500.0,
-        ]
-    )
-
-    n_frequencies = len(freq_array)
-
-    phons = [20, 40, 60, 80]
-    # phons = [15.45, 40.1, 67.25, 94.9]
-
-    n_phons = len(phons)
-
-    phon_loudness_array = np.zeros((n_phons, n_frequencies), dtype=float)
-    phon_final_array = np.zeros((n_phons, n_frequencies), dtype=float)
-
-    for i_phone_contours in range(n_phons):
-        # The next sentence returns an array with the dB SPL values for an specific phon value
-        # spl_contours_array = equal_loudness_contours(phons[i_phone_contours])[0]
-        print("Phone cont: " + str(phons[i_phone_contours]))
-        phon_ref_value = phons[i_phone_contours]
-
-        for i_freq in range(n_frequencies):
-            signal, samples = sine_wave_generator(
-                fs, duration, phon_ref_value, freq_array[i_freq]
-            )
-            print("Freq: " + str(freq_array[i_freq]))
-
-            """
-            The next sentence returns the "t_array" from the function comp_loudness to then calculate the mean value 
-            of the array. That makes possible to retain the resulting loudness.
-            """
-            t_array = comp_loudness(signal, validation=False)[1]
-            # CHANGED TO MAX
-            mean_loudness_value = float(np.mean(t_array[:, 0]))
-            print("t_array: " + str(mean_loudness_value))
-
-            """
-            Here is calculated the phon value for the result that is in sones.
-            """
-            phon_loudness_value = sone2phone(mean_loudness_value)
-
-            phon_loudness_array[i_phone_contours][i_freq] = phon_loudness_value
-            print("phone_loudness_value: " + str(phon_loudness_value))
-
-            phon_diff = abs(phon_ref_value - phon_loudness_value)
-
-            phon_final_value = (
-                phon_ref_value - phon_diff
-                if (phon_ref_value <= phon_loudness_value)
-                else phon_ref_value + phon_diff
-            )
-
-            phon_final_array[i_phone_contours][i_freq] = phon_final_value
-            print("phone_loudness_value: " + str(phon_final_value))
-
-    spl_array_0_phons, frequencies = equal_loudness_contours(0)
-    spl_array_20_phons = equal_loudness_contours(20)[0]
-    spl_array_40_phons = equal_loudness_contours(40)[0]
-    spl_array_60_phons = equal_loudness_contours(60)[0]
-    spl_array_80_phons = equal_loudness_contours(80)[0]
-
-    plt.figure(figsize=(10, 5))
-    # plt.semilogx(frequencies, spl_array_0_phons, 'g:')
-    plt.semilogx(frequencies, spl_array_20_phons, "b:")
-    plt.semilogx(frequencies, spl_array_40_phons, "g:")
-    plt.semilogx(frequencies, spl_array_60_phons, "r:")
-    plt.semilogx(frequencies, spl_array_80_phons, "c:")
-
-    plt.semilogx(freq_array, phon_final_array[0, :], "b", label="20 phons")
-    plt.semilogx(freq_array, phon_final_array[1, :], "g", label="40 phons")
-    plt.semilogx(freq_array, phon_final_array[2, :], "r", label="60 phons")
-    plt.semilogx(freq_array, phon_final_array[3, :], "c", label="80 phons")
-    # [
-    #     plt.semilogx(
-    #         freq_array,
-    #         phon_final_array[ii_phone_contour, :],
-    #         label=str(phons[ii_phone_contour]) + " phons",
-    #     )
-    #     for ii_phone_contour in range(n_phons)
-    # ]
-
-    plt.xlim(left=100, right=10100)
-    plt.ylabel("Sound pressure level [dB SPL]")
-    plt.xlabel("Frequency [Hz]")
-    plt.title("Specific loudness contours")
-    plt.grid(which="both", linestyle="-", color="grey")
-    plt.xticks(
-        [20, 50, 100, 200, 500, 1000, 2000, 5000, 10000, 20000],
-        ["20", "50", "100", "200", "500", "1K", "2K", "5K", "10K", "20K"],
-    )
-    plt.legend()
-    plt.show()
-
-    print("Validation passed")
-
-
-if __name__ == "__main__":
-    # hearing_model_validation()
-    sig, samp = sine_wave_generator(48000, 5, 100, 1000)
-
-    # plt.figure(figsize=(10, 5))
-    # plt.plot(samp, sig)
-    # plt.xlim(left=0)
-    # plt.ylabel('Sound pressure level [dB SPL]')
-    # plt.xlabel('Time (s)')
-    # plt.title('prueba')
-    # plt.grid(which='both', linestyle='-', color='grey')
-    # plt.show()
-    #
-    # spl_array_0_phons, frequencies = equal_loudness_contours(0)
-    # spl_array_20_phons = equal_loudness_contours(20)[0]
-    # spl_array_40_phons = equal_loudness_contours(40)[0]
-    # spl_array_60_phons = equal_loudness_contours(60)[0]
-    # spl_array_80_phons = equal_loudness_contours(80)[0]
-    #
-    # plt.figure(figsize=(10, 5))
-    # plt.semilogx(frequencies, spl_array_0_phons, 'green')
-    # plt.semilogx(frequencies, spl_array_20_phons, 'purple')
-    # plt.semilogx(frequencies, spl_array_40_phons, 'yellow')
-    # plt.semilogx(frequencies, spl_array_60_phons, 'red')
-    # plt.semilogx(frequencies, spl_array_80_phons, 'blue')
-    # plt.xlim(right=20000)
-    # plt.ylabel('Sound pressure level [dB SPL]')
-    # plt.xlabel('Frequency [Hz]')
-    # plt.title('Specific loudness contours (Right channel)')
-    # plt.grid(which='both', linestyle='-', color='grey')
-    # plt.show()
-
-    # fs = 48000
-    # n = len(sig)
-    # plt.figure(figsize=(10, 5))
-    # freqs = np.linspace(0, fs / 2, n)
-    # spectrum = 20 * np.log10(abs(np.fft.fft(sig * np.blackman(n))))
-
-    # plt.plot(freqs, spectrum)
-
-    # plt.xscale('log')
-    # plt.xlim(10, 10000)
-    # plt.ylim(0, 200)
-    # plt.xlabel("Frequency [Hz]")
-    # plt.ylabel("Amplitude [dB]")
-    # plt.show()
-
-    hearing_model_validation()
+# -*- coding: utf-8 -*-
+
+import numpy as np
+import matplotlib.pyplot as plt
+from scipy.optimize import root_scalar
+
+# Project Imports
+from mosqito.functions.loudness_ecma.comp_loudness import comp_loudness
+from mosqito.functions.loudness_ecma.equal_loudness_contours import (
+    equal_loudness_contours,
+)
+from mosqito.functions.shared.sine_wave_generator import (
+    sine_wave_generator,
+)
+
+
+def comp_loudness_wrapper(spl, freq=1000, n_1kHz=0):
+    """Return the difference between the loudness of a tone
+    at <freq> Hz and <spl> dB and a tone at 1 KHz and <spl> dB
+
+    Parameters
+    ----------
+    spl: float
+        Test tone Sound Pressure Level [dB]
+    freq: float
+        Test tone frequency [Hz]
+    n_1kHz: float
+        Total loudness value of a 1 kHz / <spl> dB test tone
+
+    Returns
+    -------
+    n_diff: float
+        Difference between the test tone and the 1 kHz tone loudness
+
+    """
+
+    # test signal generation
+    duration = 0.25
+    fs = 48000.0
+    signal, _ = sine_wave_generator(fs, duration, spl, freq)
+
+    # compute loudness for the test signal
+    n_specific, _ = comp_loudness(signal)
+    n_specific = np.array(n_specific)
+    n_tot = np.sum(n_specific, axis=0)
+    n_tot_mean = np.mean(n_tot[5:])
+
+    # Return n_tot_mean - n_1kHz
+    n_diff = n_tot_mean - n_1kHz
+
+    return n_diff
+
+
+def hearing_model_validation():
+    """Validation of the ECMA loudness implementation according
+    to ECMA-418-2 annex A
+    """
+
+    phons = [80, 60, 40, 20]
+
+    col_vec = ["tab:blue", "tab:red", "tab:orange", "tab:purple", "tab:green"]
+
+    for phon, col in zip(phons, col_vec):
+        print(str(phon) + " Phon equal loudness contour")
+
+        # Get ISO 226 equal loudness contour (index 17 = 1000 Hz)
+        spl_iso_vec, freq_vec = equal_loudness_contours(phon)
+
+        # Compute ecma loudness of a 1 kHz tone
+        n_1kHz = comp_loudness_wrapper(spl=phon)
+        spl_ecma_vec = np.zeros(spl_iso_vec.shape)
+        spl_ecma_vec[17] = phon
+
+        # For each frequency...
+        for i in np.arange(5, freq_vec.shape[0]):
+            if i != 17:
+                spl = spl_iso_vec[i]
+                freq = freq_vec[i]
+                # ... find the spl of the tone that produces the
+                # same loudness as the 1 kHz tone
+                spl_ecma_vec[i] = root_scalar(
+                    comp_loudness_wrapper,
+                    x0=spl,
+                    args=(freq, n_1kHz),
+                    bracket=[-5, 120],
+                ).root
+        # Plot the ISO 226 equal loudness contour
+        plt.semilogx(freq_vec, spl_iso_vec, ":", color=col)
+        # Plot the ECMA equal contour loudness
+        if phon == 0:
+            label = "lower treshold of hearing"
+        else:
+            label = str(phon) + " phon"
+        plt.semilogx(freq_vec[5:], spl_ecma_vec[5:], color=col, label=label)
+
+    plt.xlim(left=80, right=11000)
+    plt.ylim((-15, 105))
+    plt.ylabel("Sound Pressure Level [dB]")
+    plt.xlabel("Frequency [Hz]")
+    plt.grid(which="both", linestyle="-", color="grey")
+    plt.legend()
+    plt.savefig(
+        "./validations/loudness_ecma/output/" + "ecma_hearing_model_validation.png",
+        format="png",
+    )
+
+
+if __name__ == "__main__":
+    hearing_model_validation()
diff --git a/validations/loudness_ecma/hearing_model_validation_alt.py b/validations/loudness_ecma/hearing_model_validation_alt.py
deleted file mode 100644
index 831655fe..00000000
--- a/validations/loudness_ecma/hearing_model_validation_alt.py
+++ /dev/null
@@ -1,82 +0,0 @@
-# -*- coding: utf-8 -*-
-"""
-@author: Daniel Jiménez-Caminero Costa
-"""
-import numpy as np
-import matplotlib.pyplot as plt
-
-
-# Project Imports
-from mosqito.functions.loudness_ecma_spain.sone2phone import sone2phone
-from mosqito.functions.loudness_ecma_spain.comp_loudness_alt import comp_loudness
-from mosqito.functions.loudness_ecma_spain.equal_loudness_contours import (
-    equal_loudness_contours,
-)
-from mosqito.functions.loudness_ecma_spain.sine_wave_generator import (
-    sine_wave_generator,
-)
-
-
-def hearing_model_validation():
-    """Function that serves for the validation of the hearing model presented in Annex F of ECMA-74.
-
-    Parameters
-    ----------
-
-    Returns
-    -------
-
-    """
-    # Duration of the signal
-    duration = 1
-    # Sampling frequency
-    fs = 48000.0
-
-    phons = [20, 40, 60, 80]
-
-    plt.figure(figsize=(10, 5))
-    col_vec = ["b", "g", "r", "c"]
-
-    for phon, col in zip(phons, col_vec):
-        spl_vec, freq_vec = equal_loudness_contours(phon)
-
-        # The next sentence returns an array with the dB SPL values for an specific phon value
-        # spl_contours_array = equal_loudness_contours(phons[i_phone_contours])[0]
-        print("Phone cont: " + str(phon))
-        # phon_ref_value = phons[i_phone_contours]
-
-        loudness = []
-
-        for freq, spl in zip(freq_vec, spl_vec):
-            signal, _ = sine_wave_generator(fs, duration, spl, freq)
-            print("Freq: " + str(freq))
-
-            """
-            The next sentence returns the "t_array" from the function comp_loudness to then calculate the mean value 
-            of the array. That makes possible to retain the resulting loudness.
-            """
-            n_array = comp_loudness(signal)
-            specific_loudness = np.array(n_array)
-            tot_loudness = np.sum(specific_loudness, axis=0)
-            mean_tot_loudness = np.mean(tot_loudness)
-            loudness.append(sone2phone(mean_tot_loudness))
-
-        plt.semilogx(freq_vec, phon * np.ones(freq_vec.size), col + ":")
-        plt.semilogx(freq_vec, loudness, col, label=str(phon) + " phons")
-
-    plt.xlim(left=100, right=10100)
-    plt.ylabel("Loudness [Phons]")
-    plt.xlabel("Frequency [Hz]")
-    plt.grid(which="both", linestyle="-", color="grey")
-    plt.xticks(
-        [20, 50, 100, 200, 500, 1000, 2000, 5000, 10000, 20000],
-        ["20", "50", "100", "200", "500", "1K", "2K", "5K", "10K", "20K"],
-    )
-    plt.legend()
-    plt.show()
-
-    print("Validation passed")
-
-
-if __name__ == "__main__":
-    hearing_model_validation()
diff --git a/validations/loudness_ecma/output/ecma_hearing_model_validation.png b/validations/loudness_ecma/output/ecma_hearing_model_validation.png
new file mode 100644
index 00000000..cd97b738
Binary files /dev/null and b/validations/loudness_ecma/output/ecma_hearing_model_validation.png differ
diff --git a/validations/loudness_ecma/valid_ear_filter_design.py b/validations/loudness_ecma/valid_ear_filter_design.py
index 6df0d071..7b726290 100644
--- a/validations/loudness_ecma/valid_ear_filter_design.py
+++ b/validations/loudness_ecma/valid_ear_filter_design.py
@@ -10,8 +10,8 @@
 )
 import matplotlib.pyplot as plt
 
-from mosqito.functions.loudness_ecma_spain.ear_filter_design import ear_filter_design
-from mosqito.functions.loudness_ecma_spain.sine_wave_generator import (
+from mosqito.functions.loudness_ecma.ear_filter_design import ear_filter_design
+from mosqito.functions.shared.sine_wave_generator import (
     sine_wave_generator,
 )
 
diff --git a/validations/loudness_ecma/valid_gammatone.py b/validations/loudness_ecma/valid_gammatone.py
index 55fd244c..5b43a21a 100644
--- a/validations/loudness_ecma/valid_gammatone.py
+++ b/validations/loudness_ecma/valid_gammatone.py
@@ -4,10 +4,10 @@
 from numpy import log10, abs as np_abs
 import matplotlib.pyplot as plt
 
-from mosqito.functions.loudness_ecma_spain.gammatone import (
+from mosqito.functions.loudness_ecma.gammatone import (
     gammatone as mosqito_gamma,
 )
-from mosqito.functions.loudness_ecma_spain.gen_auditory_filters_centre_freq import (
+from mosqito.functions.loudness_ecma.gen_auditory_filters_centre_freq import (
     gen_auditory_filters_centre_freq,
 )
 
diff --git a/validations/loudness_ecma/valid_loudness_function.py b/validations/loudness_ecma/valid_loudness_function.py
deleted file mode 100644
index 124b4c06..00000000
--- a/validations/loudness_ecma/valid_loudness_function.py
+++ /dev/null
@@ -1,128 +0,0 @@
-import numpy as np
-import matplotlib.pyplot as plt
-
-from mosqito.functions.loudness_ecma_spain.loudness_function import (
-    loudness_function as loudness_function,
-)
-
-p = np.linspace(10, 100, 100)
-loudness = [loudness_function(2e-5 * 10 ** (i / 20)) for i in p]
-plt.semilogy(p, loudness, label="ECMA")
-
-p_ref = [
-    13.987476,
-    15.9721985,
-    18.109982,
-    20.782942,
-    23.303232,
-    26.129642,
-    29.567616,
-    33.082165,
-    37.132286,
-    42.405235,
-    47.44942,
-    51.80556,
-    56.620693,
-    62.19986,
-    68.238014,
-    72.82331,
-    79.31967,
-    84.89865,
-    90.40075,
-    94.90907,
-    98.729164,
-    102.396385,
-    105.2226,
-    107.590416,
-    109.57592,
-]
-
-loud_ref = [
-    0.020996287,
-    0.029044019,
-    0.039844647,
-    0.0555822,
-    0.07562439,
-    0.10120133,
-    0.13886018,
-    0.18895733,
-    0.2529075,
-    0.36485842,
-    0.49659508,
-    0.6591871,
-    0.86782557,
-    1.2211897,
-    1.7043209,
-    2.3195395,
-    3.431377,
-    4.909524,
-    7.2618318,
-    10.302499,
-    14.493973,
-    20.221483,
-    27.514244,
-    35.91071,
-    46.47905,
-]
-
-plt.semilogy(p_ref, loud_ref, label="Buus et al.")
-
-p_ref1 = [
-    6.598919,
-    8.838655,
-    11.250792,
-    14.17477,
-    17.785425,
-    20.533108,
-    23.620234,
-    26.878784,
-    31.16102,
-    35.784645,
-    40.063473,
-    44.85609,
-    48.621613,
-    53.92802,
-    58.548725,
-    64.02948,
-    70.02354,
-    76.3624,
-    82.87268,
-    87.32732,
-    92.63762,
-    97.607506,
-    101.3779,
-    104.46454,
-]
-
-loud_ref1 = [
-    0.002630106,
-    0.0051425043,
-    0.010632731,
-    0.021176357,
-    0.047157302,
-    0.076508306,
-    0.11306698,
-    0.17023043,
-    0.2377996,
-    0.32601807,
-    0.3997006,
-    0.50857294,
-    0.6120943,
-    0.8082804,
-    0.9908526,
-    1.4907696,
-    2.2847574,
-    3.9156787,
-    6.8367267,
-    10.09929,
-    15.481487,
-    25.099924,
-    36.40073,
-    52.800743,
-]
-
-plt.semilogy(p_ref1, loud_ref1, label="Buus et al. 1")
-plt.legend()
-plt.show()
-
-pass
\ No newline at end of file
diff --git a/validations/loudness_ecma/valid_segmentation.py b/validations/loudness_ecma/valid_segmentation.py
index 580a8042..6bf49356 100644
--- a/validations/loudness_ecma/valid_segmentation.py
+++ b/validations/loudness_ecma/valid_segmentation.py
@@ -3,10 +3,10 @@
 import matplotlib.pyplot as plt
 
 # Local application imports
-from mosqito.functions.loudness_ecma_spain.sine_wave_generator import (
+from mosqito.functions.shared.sine_wave_generator import (
     sine_wave_generator,
 )
-from mosqito.functions.loudness_ecma_spain.segmentation_blocks import (
+from mosqito.functions.loudness_ecma.segmentation_blocks import (
     segmentation_blocks,
 )
 
@@ -16,7 +16,7 @@
     spl_value=60,
     freq=40,
 )
-blocks = segmentation_blocks(signal, 8192, 2048, 1)
+blocks = segmentation_blocks(signal, 8192, 2048)
 
 plt.subplot(211)
 plt.plot(signal)


Test patch:
diff --git a/pytest.ini b/pytest.ini
index 010f6372..1336284d 100644
--- a/pytest.ini
+++ b/pytest.ini
@@ -1,5 +1,6 @@
 [pytest]
 markers =
+    loudness_ecma: marks tests related to loudness ecma
     loudness_zwst: marks tests related to loudness zwicker stationary (deselect with '-m "not loudness_zwst"')
     loudness_zwtv: marks tests related to loudness zwicker time-varying
     sharpness_din: marks tests related to sharpness DIN45692
diff --git a/tests/Audio/test_Audio.py b/tests/Audio/test_Audio.py
index 38b5c2af..d6bcfb36 100644
--- a/tests/Audio/test_Audio.py
+++ b/tests/Audio/test_Audio.py
@@ -169,11 +169,17 @@ def test_compute_tnr_pr(fixture_import_signal):
 @pytest.mark.audio
 def test_compute_tnr_pr_time(fixture_import_signal_time):
     audio = fixture_import_signal_time
+
     audio.compute_tnr_pr(method="all")
     audio.tonality["tnr"].plot_3D_Data(
-        "time",
         "freqs",
+        "time",
         is_2D_view=True,
         is_show_fig=is_show_fig,
         save_path=out_path + "test_Audio_compute_tnr_pr_time.png",
-    )
\ No newline at end of file
+    )
+
+
+# test of the function
+if __name__ == "__main__":
+    test_compute_tnr_pr_time()
\ No newline at end of file
diff --git a/tests/loudness/test_loudness_ecma.py b/tests/loudness/test_loudness_ecma.py
new file mode 100644
index 00000000..aa0c9685
--- /dev/null
+++ b/tests/loudness/test_loudness_ecma.py
@@ -0,0 +1,48 @@
+# -*- coding: utf-8 -*-
+
+# Third party imports
+import pytest
+import numpy as np
+
+# Local application imports
+from mosqito.functions.shared.sine_wave_generator import (
+    sine_wave_generator,
+)
+from mosqito.functions.loudness_ecma.comp_loudness import comp_loudness
+
+
+@pytest.mark.loudness_ecma  # to skip or run only loudness ecma tests
+def test_loudness_ecma():
+    """Test function for the Loudness_ecma calculation"""
+
+    # Generate a 1kHz / 80 dB test tone and compute loudness
+    signal, _ = sine_wave_generator(
+        fs=48000,
+        t=0.25,
+        spl_value=80,
+        freq=1000,
+    )
+    n_specific, _ = comp_loudness(signal)
+    n_specific = np.array(n_specific)
+    n_tot = np.sum(n_specific, axis=0)
+    n_tot_mean_1kHz = np.mean(n_tot[5:])
+
+    # Generate a 5kHz / 78.7 dB test tone and compute loudness
+    signal, _ = sine_wave_generator(
+        fs=48000,
+        t=0.25,
+        spl_value=78.73977248964925,
+        freq=5000,
+    )
+    n_specific, _ = comp_loudness(signal)
+    n_specific = np.array(n_specific)
+    n_tot = np.sum(n_specific, axis=0)
+    n_tot_mean_5kHz = np.mean(n_tot[5:])
+
+    # Compare loudness
+    np.testing.assert_almost_equal(n_tot_mean_1kHz, n_tot_mean_5kHz)
+
+
+# test de la fonction
+if __name__ == "__main__":
+    test_loudness_ecma()
diff --git a/tests/loudness/tst_loudness_ecma.py b/tests/loudness/tst_loudness_ecma.py
deleted file mode 100644
index d656b56d..00000000
--- a/tests/loudness/tst_loudness_ecma.py
+++ /dev/null
@@ -1,46 +0,0 @@
-# -*- coding: utf-8 -*-
-
-# Third party imports
-import pytest
-import numpy as np
-
-# Local application imports
-from mosqito.functions.loudness_ecma_spain.sine_wave_generator import (
-    sine_wave_generator,
-)
-from mosqito.functions.loudness_ecma_spain.comp_loudness_alt import comp_loudness
-from mosqito.functions.loudness_ecma_spain.sone2phone import sone2phone
-
-
-# @pytest.mark.loudness_ecma  # to skip or run only loudness ecma tests
-def tst_loudness_ecma():
-    """Test function for the Loudness_ecma calculation
-
-    Parameters
-    ----------
-    None
-
-    Outputs
-    -------
-    None
-    """
-
-    signal, _ = sine_wave_generator(
-        fs=48000,
-        t=1,
-        spl_value=60,
-        freq=1000,
-    )
-    n_array = comp_loudness(signal)
-    specific_loudness = np.array(n_array)
-    tot_loudness = np.sum(specific_loudness, axis=0)
-    mean_tot_loudness = np.mean(tot_loudness)
-    phon_loudness_value = sone2phone(mean_tot_loudness)
-    assert phon_loudness_value < 60.1 and phon_loudness_value > 59.9
-
-    pass
-
-
-# test de la fonction
-if __name__ == "__main__":
-    tst_loudness_ecma()