#!/usr/bin/env python
"""
Collection of functions to extract spectral features from audio signals
License: New BSD License
"""
# =============================================================================
# Load the modules
# =============================================================================
# Import external modules
import numpy as np
import pandas as pd
from scipy import interpolate
from numpy.lib.function_base import _quantile_is_valid
from scipy.optimize import root
# Import internal modules
from maad.util import moments, power2dB
from maad import sound
# define internal functions
def _interpolate_peak_location(pxx):
""" Estimate peak using quadratic interpolation """
idxmax = pxx.index.get_loc(pxx.idxmax())
if (idxmax==0) | (idxmax+1==len(pxx)):
# if maximum is at beginning or at the end os spectrum, take maximum
peak = pxx.idxmax()
amplitude = pxx[peak]
else:
# use interpolation
x = pxx.iloc[[idxmax-1, idxmax, idxmax+1]].index.values
y = pxx.iloc[[idxmax-1, idxmax, idxmax+1]].values
f = interpolate.interp1d(x,y, kind='quadratic')
xnew = np.arange(x[0], x[2], 1)
ynew = f(xnew)
peak = xnew[ynew.argmax()]
amplitude = ynew[ynew.argmax()]
return peak, amplitude
#%%
# =============================================================================
# public functions
# =============================================================================
#%%
[docs]
def spectral_moments (X, axis=None):
"""
Computes the first 4th moments of an amplitude spectrum (1d) or spectrogram (2d),
mean, variance, skewness, kurtosis.
Parameters
----------
X : ndarray of floats
Amplitude spectrum (1d) or spectrogram (2d).
axis : interger, optional, default is None
if spectrogram (2d), select the axis to estimate the moments.
Returns
-------
mean : float
mean of the audio
var : float
variance of the audio
skew : float
skewness of the audio
kurt : float
kurtosis of the audio
Examples
--------
>>> from maad import sound, features
>>> s, fs = sound.load('../data/spinetail.wav')
>>> Sxx_power,_,_,_ = sound.spectrogram(s, fs)
Compute spectral moments on the mean spectrum
>>> import numpy as np
>>> S_power = sound.avg_power_spectro(Sxx_power)
>>> sm, sv, ss, sk = features.spectral_moments (S_power)
>>> print('mean: %2.8f / var: %2.10f / skewness: %2.2f / kurtosis: %2.2f' % (sm, sv, ss, sk))
mean: 0.00000228 / var: 0.0000000001 / skewness: 5.84 / kurtosis: 40.49
Compute spectral moments of the spectrogram along the time axis
>>> from maad.features import spectral_moments
>>> sm_per_bin, sv_per_bin, ss_per_bin, sk_per_bin = spectral_moments(Sxx_power, axis=1)
>>> print('Length of sk_per_bin is : %2.0f' % len(sk_per_bin))
Length of sk_per_bin is : 512
"""
# force P to be ndarray
X = np.asarray(X)
return moments(X,axis)
#%%
[docs]
def peak_frequency(s, fs, method='best', nperseg=1024, roi=None, amp=False, as_pandas=False, **kwargs):
"""
Compute the peak frequency for audio signal. The peak frequency is the frequency of
maximum power. If a region of interest with time and spectral limits is provided,
the peak frequency is computed on the selection.
Parameters
----------
s : 1D array
Input audio signal
fs : float
Sampling frequency of audio signal
method : {'fast', 'best'}, optional
Method used to compute the peak frequency.
The default is 'best'.
nperseg : int, optional
Length of segment to compute the FFT. The default is 1024.
roi : pandas.Series, optional
Region of interest where peak frequency will be computed.
Series must have a valid input format with index: min_t, min_f, max_t, max_f.
The default is None.
as_pandas: bool
Return data as a pandas.Series or pandas.DataFrame, when amp
is False or True, respectively. Default is False.
kwargs : additional keyword arguments
Additional keyword arguments to pass to `sound.spectrum`.
Returns
-------
peak_freq : float
Peak frequency of audio segment.
amp_peak_freq : float
Amplitud of the peak frequency of audio segment.
Examples
--------
>>> from maad import features, sound
>>> s, fs = sound.load('../data/spinetail.wav')
Compute peak frequency
>>> peak_freq, peak_freq_amp = features.peak_frequency(s, fs, amp=True)
>>> print('Peak Frequency: {:.5f}, Amplitude: {:.5f}'.format(peak_freq, peak_freq_amp))
Peak Frequency: 6634.16016, Amplitude: 0.00013
"""
if roi is None:
pxx, fidx = sound.spectrum(s, fs, nperseg, **kwargs)
else:
pxx, fidx = sound.spectrum(s, fs, nperseg,
tlims=[roi.min_t, roi.max_t],
flims=[roi.min_f, roi.max_f],
**kwargs)
pxx = pd.Series(pxx, index=fidx)
if method=='fast':
# simplest form to get peak frequency, but less precise
peak_freq = pxx.idxmax()
amp_peak_freq = pxx[peak_freq]
elif method=='best':
# use interpolation get better estimate of peak frequency
peak_freq, amp_peak_freq = _interpolate_peak_location(pxx)
else:
raise Exception("Invalid method. Method should be 'fast' or 'best' ")
if amp:
if as_pandas:
out = pd.DataFrame({"freq":peak_freq, "amp":amp_peak_freq}, index=[0])
else:
out = np.array([peak_freq, amp_peak_freq])
else:
if as_pandas:
out = pd.Series([peak_freq], index=["peak_freq"])
else:
out = peak_freq
return out
#%%
[docs]
def spectral_quantile(s, fs, q=[0.05, 0.25, 0.5, 0.75, 0.95], nperseg=1024, roi=None,
as_pandas=False, amp=False, **kwargs):
"""
Compute the q-th quantile of the power spectrum. If a region of interest with time and
spectral limits is provided, the q-th quantile is computed on the selection.
Parameters
----------
s : 1D array
Input audio signal
fs : float
Sampling frequency of audio signal
q : array or float, optional
Quantile or sequence of quantiles to compute, which must be between 0 and 1
inclusive.
The defaul is [0.05, 0.25, 0.5, 0.75, 0.95].
nperseg : int, optional
Length of segment to compute the FFT. The default is 1024.
roi : pandas.Series, optional
Region of interest where peak frequency will be computed.
Series must have a valid input format with index: min_t, min_f, max_t, max_f.
The default is None.
as_pandas: bool
Return data as a pandas.Series or pandas.DataFrame, when amp
is False or True, respectively. Default is False.
amp: bool
Enable quantiles amplitude output. Default is False.
kwargs : additional keyword arguments
If `window='hann'`, additional keyword arguments to pass to
`sound.spectrum`.
Returns
-------
Pandas Series or Numpy array
Quantiles of power spectrum.
Examples
--------
>>> from maad import sound, features
>>> s, fs = sound.load('../data/spinetail.wav')
Compute the q-th quantile of the power spectrum
>>> qs = features.spectral_quantile(s, fs, [0.05, 0.25, 0.5, 0.75, 0.95], as_pandas=True)
>>> print(qs)
0.05 6029.296875
0.25 6416.894531
0.50 6632.226562
0.75 6890.625000
0.95 9216.210938
dtype: float64
Compute the q-th quantile of the power spectrum and its amplitude
>>> qs = features.spectral_quantile(s, fs, [0.05, 0.25, 0.5, 0.75, 0.95], amp=True, as_pandas=True)
>>> print(qs)
freq amp
0.05 6029.296875 0.000007
0.25 6416.894531 0.000067
0.50 6632.226562 0.000087
0.75 6890.625000 0.000022
0.95 9216.210938 0.000004
"""
q = np.asanyarray(q)
if not _quantile_is_valid(q):
raise ValueError("Percentiles must be in the range [0, 1]")
# Compute spectrum
# if roi is None:
# pxx, fidx = sound.spectrum(s, fs, nperseg, **kwargs)
# else:
# pxx, fidx = sound.spectrum(s, fs, nperseg,
# tlims=[roi.min_t, roi.max_t],
# flims=[roi.min_f, roi.max_f],
# **kwargs)
# pxx = pd.Series(pxx, index=fidx)
if roi is None:
Sxx,_,fn,_ = sound.spectrogram(s, fs, nperseg=nperseg, **kwargs)
else:
Sxx,_,fn,_ = sound.spectrogram(s, fs, nperseg=nperseg,
tlims=[roi.min_t, roi.max_t],
flims=[roi.min_f, roi.max_f],
**kwargs)
pxx = pd.Series(np.average(Sxx,axis=1), index=fn)
# Compute spectral q
norm_cumsum = pxx.cumsum()/pxx.sum()
spec_quantile = []
for quantile in q:
spec_quantile.append(pxx.index[np.where(norm_cumsum>=quantile)[0][0]])
if amp:
if as_pandas:
out = pd.DataFrame({"freq":spec_quantile, "amp":pxx[spec_quantile].values}, index=q)
else:
out = np.transpose(np.array([q, spec_quantile, pxx[spec_quantile]]))
else:
if as_pandas:
out = pd.Series(spec_quantile, index=q)
else:
out = np.array(spec_quantile)
return out
#%%
[docs]
def spectral_bandwidth(s, fs, nperseg=1024, roi=None, as_pandas=False, **kwargs):
"""
Compute the bandwith of the power spectrum. If a region of interest with time
and spectral limits is provided, the spectral bandwidth is computed on the selection.
Parameters
----------
s : 1D array
Input audio signal
fs : float
Sampling frequency of audio signal
nperseg : int, optional
Length of segment to compute the FFT. The default is 1024.
roi : pandas.Series, optional
Region of interest where peak frequency will be computed.
Series must have a valid input format with index: min_t, min_f, max_t, max_f.
The default is None.
as_pandas : bool
Return data as a pandas.Series. This is usefull when computing multiple features
over a signal. Default is False.
kwargs : additional keyword arguments
If `window='hann'`, additional keyword arguments to pass to
`sound.spectrum`.
Returns
-------
bandwidth_50 : float
Bandwidth 50% of the audio
bandwidth_90 : float
Bandwidth 90% of the audio
Examples
--------
>>> from maad import features, sound
>>> s, fs = sound.load('../data/spinetail.wav')
>>> bw_50, bw_90 = features.spectral_bandwidth(s, fs, nperseg=1024)
>>> print("Bandwidth 50% : {:.4f} / Bandwidth 90% : {:.4f}".format(bw_50, bw_90))
Bandwidth 50% : 473.7305 / Bandwidth 90% : 3186.9141
"""
# Compute spectral bandwith with the quantiles
q = spectral_quantile(
s, fs, [0.05, 0.25, 0.75, 0.95], nperseg, roi, False, True, **kwargs)
out = np.array([np.abs(q[2,1]-q[1,1]), np.abs(q[3,1]-q[0,1])])
if as_pandas:
out = pd.Series(out, index=["bw_50", "bw_90"])
return out
#%%
[docs]
def all_spectral_features(s, fs, nperseg=1024, roi=None, method='fast', **kwargs):
"""
Compute all the spectral features for a signal.
Parameters
----------
s : 1D array
Input audio signal
fs : float
Sampling frequency of audio signal
nperseg : int, optional
Length of segment to compute the FFT. The default is 1024.
roi : pandas.Series, optional
Region of interest where peak frequency will be computed.
Series must have a valid input format with index: min_t, min_f, max_t, max_f.
The default is None.
method : {'fast', 'best'}, optional
Method used to compute the peak frequency.
The default is 'fast'.
kwargs : additional keyword arguments
If `window='hann'`, additional keyword arguments to pass to
`sound.spectrum`.
Returns
-------
spectral_features : pandas DataFrame
DataFrame with all spectral features computed in the spectrum
Examples
--------
>>> from maad import features, sound
>>> s, fs = sound.load('../data/spinetail.wav')
Compute all the spectral features
>>> features.all_spectral_features(s, fs, nperseg=1024, roi=None)
sm 2.276330e-06
sv 8.118042e-11
ss 5.844664e+00
sk 4.048891e+01
freq_05 6.029297e+03
freq_25 6.416895e+03
freq_50 6.632227e+03
freq_75 6.890625e+03
freq_95 9.216211e+03
peak_freq 6.632227e+03
bw_50 4.737305e+02
bw_90 3.186914e+03
dtype: float64
"""
# Compute transformations
Sxx_power,_,_,_ = sound.spectrogram (s, fs, window='hann', nperseg=nperseg)
S_power = sound.avg_power_spectro(Sxx_power)
# Compute features
sm = spectral_moments(S_power)
qs = spectral_quantile(s, fs, [0.05, 0.25, 0.5, 0.75, 0.95], nperseg, roi, **kwargs)
bw_50, bw_90 = spectral_bandwidth(s, fs, nperseg, roi, **kwargs)
peak_freq = peak_frequency(s, fs, method, nperseg, roi, **kwargs)
# Organize data into a Dataframe
spectral_features = pd.Series({
"sm":sm[0],
"sv":sm[1],
"ss":sm[2],
"sk":sm[3],
"freq_05":qs[0],
"freq_25":qs[1],
"freq_50":qs[2],
"freq_75":qs[3],
"freq_95":qs[4],
"peak_freq":peak_freq,
"bw_50":bw_50,
"bw_90":bw_90})
return spectral_features
if __name__ == "__main__":
import doctest
doctest.testmod()