# Add MOSQITO to the Python path
import sys
sys.path.append('..')

# Import numpy
import numpy as np
# Import plot function
import matplotlib.pyplot as plt
# Import spectrum computation tool
from scipy.fft import fft, fftfreq
# Import mosqito functions
from mosqito.utils import load
from mosqito.sq_metrics import loudness_zwst
from mosqito.sq_metrics import loudness_zwst_perseg
from mosqito.sq_metrics import loudness_zwst_freq

# Import MOSQITO color sheme [Optional]
from mosqito import COLORS

# To get inline plots (specific to Jupyter notebook)
%matplotlib notebook

"""
Explanation: <h1>Table of Contents<span class="tocSkip"></span></h1>
<div class="toc"><ul class="toc-item"><li><span><a href="#Load-signal" data-toc-modified-id="Load-signal-1"><span class="toc-item-num">1&nbsp;&nbsp;</span>Load signal</a></span></li><li><span><a href="#Compute-loudness-of-the-whole-signal" data-toc-modified-id="Compute-loudness-of-the-whole-signal-2"><span class="toc-item-num">2&nbsp;&nbsp;</span>Compute loudness of the whole signal</a></span></li><li><span><a href="#Compute-loudness-per-signal-segments" data-toc-modified-id="Compute-loudness-per-signal-segments-3"><span class="toc-item-num">3&nbsp;&nbsp;</span>Compute loudness per signal segments</a></span></li><li><span><a href="#Compute-tonality-from-spectrum" data-toc-modified-id="Compute-tonality-from-spectrum-4"><span class="toc-item-num">4&nbsp;&nbsp;</span>Compute tonality from spectrum</a></span></li></ul></div>

How to compute acoustic Loudness according to Zwicker method for stationary signals
This tutorial explains how to use MOSQITO to compute the acoustic loudness of a signal according to the Zwicker method for stationary signals (as per ISO 532-1 section 5). For more information on the implementation and validation of the metric, you can refer to the documentation.
The following commands are used to import the necessary functions.
End of explanation
"""


# Define path to the .wav file
# To be replaced by your own path
path = "../validations/sq_metrics/loudness_zwst/input/ISO_532_1/Test signal 5 (pinknoise 60 dB).wav"
# load signal
sig, fs = load(path, wav_calib=2 * 2 **0.5)
# plot signal
t = np.linspace(0, (len(sig) - 1) / fs, len(sig))
plt.figure(1)
plt.plot(t, sig, color=COLORS[0])
plt.xlabel('Time [s]')
plt.ylabel('Acoustic pressure [Pa]')

"""
Explanation: Load signal
In this tutorial, the signal is imported from a .wav file. The tutorial Audio signal basic operations gives more information about the syntax of the import and the other supported file types. You can use any .wav file to perform the tutorial or you can download the pink noise signal from MOSQITO that is used in the following.
End of explanation
"""


N, N_specific, bark_axis = loudness_zwst(sig, fs, field_type="free")

"""
Explanation: Compute loudness of the whole signal
The acoustic loudness is computed by using the following command line. In addition to the signal (as ndarray) and the sampling frequency, the function takes 1 input argument: "field_type" that can be set to "free" or "diffuse" depending on the environment of the audio signal recording.
End of explanation
"""


print("N_zwst = {:.1f} sone".format(N) )

"""
Explanation: The function return the Loudness of the signal :
End of explanation
"""


plt.figure(2)
plt.plot(bark_axis, N_specific, color=COLORS[0])
plt.xlabel("Critical band rate [Bark]")
plt.ylabel("N'_zwst [sone/Bark]")

"""
Explanation: The specific loudness over critical band rate is also computed :
End of explanation
"""


N, N_specific, bark_axis, time_axis = loudness_zwst_perseg(
    sig, fs, nperseg=8192 * 2, noverlap=4096
)
plt.figure(3)
plt.plot(time_axis, N, color=COLORS[0])
plt.xlabel("Time [s]")
plt.ylabel("N_zwst [sone]")
plt.ylim((0, 15))

"""
Explanation: Compute loudness per signal segments
To compute the loudnes for successive, possibly overlaping, time segments, you can use the loudness_zwst_perseg function. It accepts two more input paramters:
- nperseg: to define the length of each segment
- noverlap: to define the number of points to overlap between segments
End of explanation
"""


# Compute spectrum
n = len(sig)
spec = np.abs(2 / np.sqrt(2) / n * fft(sig)[0:n//2])
freqs = fftfreq(n, 1/fs)[0:n//2]
# Compute Loudness
N, N_specific, bark_axis = loudness_zwst_freq(spec, freqs)

print("N_zwst = {:.1f} sone".format(N) )

"""
Explanation: Compute tonality from spectrum
The commands below shows how to compute the stationary loudness from a frequency spectrum either in complex values or amplitude values using the functions from MOSQITO. One should note that only stationary values can be computed from a frequency input. 
The input spectrum can be either 1D with size (Nfrequency) or 2D with size (fNrequency x Ntime). The corresponding time axis can be either the same for all the spectra, with size (Nfrequency) or different for each spectrum with size (Nfrequency x Ntime).
One should pay attention that the input spectrum must be in RMS values !
End of explanation
"""


from datetime import date
print("Tutorial generation date:", date.today().strftime("%B %d, %Y"))

"""
Explanation: 
End of explanation
"""


