aus 0.1.8

// File: spectral_analysis_tools.rs
// 
// This file contains functionality for computing spectral features.
// 
// Note that the "compute..." functions are for use with the "analyzer" function.
// These functions help to reduce the number of duplicate calculations needed 
// (such as computing the sum of the magnitude spectrum, etc.) There are separate 
// functions that allow these operations to be run individually. The compute functions
// are not exposed via the parent module.
//
// Many of the spectral features extracted here are based on the formulas provided in
// Florian Eyben, "Real-Time Speech and Music Classification by Large Audio Feature Space Extraction," Springer, 2016.

use core::f64;
use num::Complex;
use rustfft::FftPlanner;
use crate::spectrum;
use crate::spectrum::SpectrumError;
use crate::util::*;
use crate::analysis::computation::*;

/// Calculates the alpha ratio from provided magnitude spectrum.
/// The alpha ratio is calculated by summing the bins from 50Hz to 1kHz,
/// and dividing this by the sum of the bins from 1kHz to 2kHz.
/// This function suggests the use of the magnitude spectrum, although
/// you can choose to use another spectrum.
/// (Eyben, p. 38)
/// 
/// # Example
/// ```
/// use aus::{spectrum, analysis};
/// let fft_size = 2048;
/// let audio = aus::read("myfile.wav").unwrap();
/// let imaginary_spectrum = spectrum::rfft(&audio.samples[0][..fft_size], fft_size);
/// let (magnitude_spectrum, phase_spectrum) = spectrum::complex_to_polar_rfft(&imaginary_spectrum);
/// let freqs = spectrum::rfftfreq(fft_size, audio.sample_rate);
/// let a_ratio = analysis::alpha_ratio(&magnitude_spectrum, &freqs);
/// ```
pub fn alpha_ratio(magnitude_spectrum: &[f64], rfft_freqs: &[f64]) -> f64 {
    let mut idx_50: usize = 0;
    let mut idx_1k: usize = 0;
    let mut idx_2k: usize = 0;
    for i in 0..rfft_freqs.len() {
        if rfft_freqs[i] < 50.0 {
            idx_50 += 1;
        }
        if rfft_freqs[i] < 1000.0 {
            idx_1k += 1;
        }
        if rfft_freqs[i] > 2000.0 {
            break;
        }
        idx_2k += 1;
    }

    let lower_sum: f64 = magnitude_spectrum[idx_50..idx_1k].iter().sum();
    let upper_sum: f64 = magnitude_spectrum[idx_1k..idx_2k].iter().sum();
    lower_sum / upper_sum
}

/// Computes the energy autocorrelation of a signal of length `fft_size` using the FFT method described in Eyben, 45:
/// $$
/// \textrm{ACF}_e(x)=\textrm{FFT}^{-1}\left(\textrm{FFT}(x)\cdot\overline{\textrm{FFT}(x)}\right)
/// $$
/// You must zero-pad the audio before calling this function if the `audio` length does not match the `fft_size`.
/// This autocorrelation method performs `N/2` zero-padding to the left and right as described in Eyben, 45.
/// The resulting `f64` vector contains only the computed values for `tau >= 0` and has length `fft_size`.
/// 
/// # Example
/// 
/// ```
/// use aus::{read, analysis::autocorrelation};
/// let fft_size: usize = 2048;
/// let audio = read("myfile.wav").unwrap();
/// let auto = autocorrelation(&audio.samples[0][..fft_size], fft_size).unwrap();
/// ```
pub fn autocorrelation(audio: &[f64], fft_size: usize) -> Result<Vec<f64>, SpectrumError> {
    if audio.len() != fft_size {
        return Err(SpectrumError{ error_msg: String::from("The audio length does not match the FFT size.")});
    }

    // prepare fft
    let padded_fft_size = fft_size * 2;
    let mut auto: Vec<f64> = vec![0.0; fft_size];
    let mut planner = FftPlanner::new();
    let fft = planner.plan_fft_forward(padded_fft_size);
    let ifft = planner.plan_fft_inverse(padded_fft_size);

    // prepare input audio vector
    let mut spectrum: Vec<Complex<f64>> = Vec::with_capacity(padded_fft_size);
    for _ in 0..fft_size / 2 {
        spectrum.push(Complex{re: 0.0, im: 0.0})
    }
    for i in 0..audio.len() {
        spectrum.push(Complex{re: audio[i], im: 0.0});
    }
    for _ in 0..fft_size / 2 {
        spectrum.push(Complex{re: 0.0, im: 0.0})
    }

    // compute autocorrelation
    fft.process(&mut spectrum);

    for i in 0..padded_fft_size {
        spectrum[i] = spectrum[i] * spectrum[i].conj();
    }

    ifft.process(&mut spectrum);

    // copy result into output vector
    for i in 0..fft_size {
        auto[i] = spectrum[i + fft_size].re;
    }
    
    Ok(auto)
}


/// Computes the autocorrelation of a signal of length `fft_size` using the power spectrum. 
/// $$
/// \textrm{ACF}_e(x)=\textrm{FFT}^{-1}\left(\textrm{FFT}(x)\cdot|\textrm{FFT}(x)|\right)
/// $$
/// This function is based on the librosa `autocorrelate` function. See the librosa documentation at 
/// <https://librosa.org/doc/latest/generated/librosa.autocorrelate.html#librosa.autocorrelate>.
/// You will need to slice the audio down to an appropriate size and zero-pad it before
/// running this function. This function does not slice the autocorrelation output - if you wish to specify
/// a maximum size that is smaller than the length of the provided audio, you will need to slice the output
/// vector manually.
/// 
/// # Example
/// 
/// ```
/// use aus::{read, analysis::autocorrelation_power_spectrum};
/// let fft_size: usize = 2048;
/// let audio = read("myfile.wav").unwrap();
/// let auto = autocorrelation_power_spectrum(&audio.samples[0][..fft_size], fft_size).unwrap();
/// ```
pub fn autocorrelation_power_spectrum(audio: &[f64], fft_size: usize) -> Result<Vec<f64>, SpectrumError> {
    if audio.len() != fft_size {
        return Err(SpectrumError{error_msg: String::from(format!("The audio vector length ({}) was not the same as the FFT size ({}).", audio.len(), fft_size))});
    }

    let complex_spectrum = spectrum::rfft(audio, fft_size);
    let (mut magnitude_spectrum, phase_spectrum) = spectrum::complex_to_polar_rfft(&complex_spectrum);
    for i in 0..magnitude_spectrum.len() {
        magnitude_spectrum[i] = magnitude_spectrum[i] * magnitude_spectrum[i];
    }
    let complex_power_spectrum = match spectrum::polar_to_complex_rfft(&magnitude_spectrum, &phase_spectrum) {
        Ok(x) => x,
        Err(err) => return Err(err)
    };
    match spectrum::irfft(&complex_power_spectrum, fft_size) {
        Ok(x) => Ok(x),
        Err(err) => Err(err)
    }
}

/// Calculates the Hammarberg index from provided magnitude spectrum.
/// This function suggests the use of the magnitude spectrum, although
/// you can choose to use another spectrum.
/// (Eyben, p. 38)
/// 
/// # Example
/// ```
/// use aus::{spectrum, analysis};
/// let fft_size = 2048;
/// let audio = aus::read("myfile.wav").unwrap();
/// let imaginary_spectrum = spectrum::rfft(&audio.samples[0][..fft_size], fft_size);
/// let (magnitude_spectrum, phase_spectrum) = spectrum::complex_to_polar_rfft(&imaginary_spectrum);
/// let freqs = spectrum::rfftfreq(fft_size, audio.sample_rate);
/// let h_index = analysis::hammarberg_index(&magnitude_spectrum, &freqs);
/// ```
pub fn hammarberg_index(magnitude_spectrum: &[f64], rfft_freqs: &[f64]) -> f64 {
    let mut idx_2k: usize = 0;
    for i in 0..rfft_freqs.len() {
        if rfft_freqs[i] > 2000.0 {
            break;
        }
        idx_2k += 1;
    }

    let lower_max = match maxargmax(&magnitude_spectrum[..idx_2k]) {
        Some((val, _)) => val,
        None => 0.0
    };
    let upper_max = match maxargmax(&magnitude_spectrum[idx_2k..]) {
        Some((val, _)) => val,
        None => 0.0
    };

    lower_max / upper_max
}

/// Calculates the harmonicity of a magnitude pectrum (determines how harmonic a signal is).
/// (Eyben, 43-44)
/// 
/// # Example
/// ```
/// use aus::{spectrum, analysis};
/// let fft_size = 2048;
/// let audio = aus::read("myfile.wav").unwrap();
/// let imaginary_spectrum = spectrum::rfft(&audio.samples[0][..fft_size], fft_size);
/// let (magnitude_spectrum, phase_spectrum) = spectrum::complex_to_polar_rfft(&imaginary_spectrum);
/// let h_index = analysis::harmonicity(&magnitude_spectrum, true);
/// ```
pub fn harmonicity(magnitude_spectrum: &[f64], normalize_by_total_energy: bool) -> f64 {
    let mut argmaxmin: Vec<usize> = Vec::new();
    let mut maxmin: Vec<f64> = Vec::new();

    // Find spectral peaks and minima by looking at 2 neighbors on either side.
    // Note that argmaxmin will contain alternating maxima and minima because you
    // cannot have two adjacent maxima or minima with the definition we are using.
    for i in 2..magnitude_spectrum.len() - 2 {
        if magnitude_spectrum[i-2] < magnitude_spectrum[i] && 
            magnitude_spectrum[i-1] < magnitude_spectrum[i] && 
            magnitude_spectrum[i+1] < magnitude_spectrum[i] && 
            magnitude_spectrum[i+2] < magnitude_spectrum[i] {
            argmaxmin.push(i);
            maxmin.push(magnitude_spectrum[i]);
        }
        else if magnitude_spectrum[i-2] > magnitude_spectrum[i] && 
            magnitude_spectrum[i-1] > magnitude_spectrum[i] && 
            magnitude_spectrum[i+1] > magnitude_spectrum[i] && 
            magnitude_spectrum[i+2] > magnitude_spectrum[i] {
            argmaxmin.push(i);
            maxmin.push(magnitude_spectrum[i]);
        }
    }

    // Compute the distances of maxima to minima
    let mut distance_sum: f64 = 0.0;
    for i in 1..maxmin.len() {
        distance_sum += f64::abs(maxmin[i] - maxmin[i-1]);
    }

    // Normalize the harmonicity by either the number of bins or the energy
    if normalize_by_total_energy {
        distance_sum / magnitude_spectrum.iter().sum::<f64>()
    } else {
        distance_sum / magnitude_spectrum.len() as f64
    }
}

/// Creates a log spectrum based on a provided magnitude or power spectrum.
/// Since the dB spectrum is often created using either 10 or 20 as
/// the coefficient, you need to provide the coefficient here.
/// 
/// You need to provide a floor value to avoid large negative numbers.
/// If you provide a ceiling value, the log spectrum will be scaled so the
/// max value is the ceiling value.
/// 
/// # Example
/// 
/// ```
/// use aus::{spectrum, analysis};
/// let fft_size = 2048;
/// let audio = aus::read("myfile.wav").unwrap();
/// let imaginary_spectrum = spectrum::rfft(&audio.samples[0][..fft_size], fft_size);
/// let (magnitude_spectrum, phase_spectrum) = spectrum::complex_to_polar_rfft(&imaginary_spectrum);
/// let power_spectrum = analysis::make_power_spectrum(&magnitude_spectrum);
/// let log_spectrum = analysis::make_log_spectrum(&power_spectrum, 10.0, -10e-8, Some(0.0));
/// ```
#[inline]
pub fn make_log_spectrum(spectrum: &[f64], coef: f64, floor: f64, ceil: Option<f64>) -> Vec<f64> {
    let mut log_spec: Vec<f64> = Vec::with_capacity(spectrum.len());
    for i in 0..spectrum.len() {
        let logval = coef * spectrum[i].log10();
        if logval < floor {
            log_spec.push(floor);
        } else {
            log_spec.push(logval);
        }
    }

    //if the ceiling was provided, scale by that
    if let Some(ceil_val) = ceil {
        let maxval = match crate::util::max(&log_spec) {
            Some(x) => x,
            None => 0.0
        };
        let adj_offset = ceil_val - maxval;
        for i in 0..log_spec.len() {
            log_spec[i] += adj_offset;
            if log_spec[i] < floor {
                log_spec[i] = floor;
            }        
        }
    }

    log_spec
}

/// Creates a log spectrogram based on a provided magnitude spectrogram.
/// You need to provide a floor value to avoid large negative numbers.
/// 
/// # Example
/// 
/// ```
/// use aus::{spectrum, analysis};
/// let fft_size = 2048;
/// let audio = aus::read("myfile.wav").unwrap();
/// let imaginary_spectrogram = spectrum::rstft(&audio.samples[0], fft_size, fft_size / 2, aus::WindowType::Hanning);
/// let (magnitude_spectrogram, _) = spectrum::complex_to_polar_rstft(&imaginary_spectrogram);
/// let power_spectrogram = analysis::make_power_spectrogram(&magnitude_spectrogram);
/// let log_spectrogram = analysis::make_log_spectrogram(&power_spectrogram, 10.0, 10e-8, Some(-1.0));
/// ```
#[inline]
pub fn make_log_spectrogram(spectrogram: &[Vec<f64>], coef: f64, floor: f64, ceil: Option<f64>) -> Vec<Vec<f64>> {
    let mut log_spectrogram: Vec<Vec<f64>> = Vec::with_capacity(spectrogram.len());
    for i in 0..spectrogram.len() {
        log_spectrogram.push(make_log_spectrum(&spectrogram[i], coef, floor, ceil));
    }
    log_spectrogram
}

/// Creates a power spectrum based on a provided magnitude spectrum.
/// 
/// # Example
/// 
/// ```
/// use aus::{spectrum, analysis};
/// let fft_size = 2048;
/// let audio = aus::read("myfile.wav").unwrap();
/// let imaginary_spectrum = spectrum::rfft(&audio.samples[0][..fft_size], fft_size);
/// let (magnitude_spectrum, phase_spectrum) = spectrum::complex_to_polar_rfft(&imaginary_spectrum);
/// let power_spectrum = analysis::make_power_spectrum(&magnitude_spectrum);
/// ```
#[inline]
pub fn make_power_spectrum(magnitude_spectrum: &[f64]) -> Vec<f64> {
    let mut power_spec: Vec<f64> = vec![0.0; magnitude_spectrum.len()];
    for i in 0..magnitude_spectrum.len() {
        power_spec[i] = magnitude_spectrum[i].powf(2.0);
    }
    power_spec
}

/// Creates a power spectrum based on a provided magnitude spectrum.
/// 
/// # Example
/// 
/// ```
/// use aus::{spectrum, analysis};
/// let fft_size = 2048;
/// let audio = aus::read("myfile.wav").unwrap();
/// let imaginary_spectrogram = spectrum::rstft(&audio.samples[0], fft_size, fft_size / 2, aus::WindowType::Hanning);
/// let (magnitude_spectrogram, _) = spectrum::complex_to_polar_rstft(&imaginary_spectrogram);
/// let power_spectrogram = analysis::make_power_spectrogram(&magnitude_spectrogram);
/// ```
#[inline]
pub fn make_power_spectrogram(magnitude_spectrogram: &[Vec<f64>]) -> Vec<Vec<f64>> {
    let mut power_spectrogram: Vec<Vec<f64>> = Vec::with_capacity(magnitude_spectrogram.len());
    for i in 0..magnitude_spectrogram.len() {
        let mut power_spec: Vec<f64> = vec![0.0; magnitude_spectrogram[i].len()];
        for j in 0..magnitude_spectrogram[i].len() {
            power_spec[j] = magnitude_spectrogram[i][j].powf(2.0);
        }
        power_spectrogram.push(power_spec);
    }
    power_spectrogram
}

/// Generates the spectrum power mass function (PMF) based on provided power spectrum 
/// and sum of power spectrum.
/// (Eyben, p. 40)
#[inline]
pub fn make_spectrum_pmf(power_spectrum: &[f64], power_spectrum_sum: f64) -> Vec<f64> {
    let mut pmf_vector: Vec<f64> = vec![0.0; power_spectrum.len()];
    for i in 0..power_spectrum.len() {
        pmf_vector[i] = power_spectrum[i] / power_spectrum_sum;
    }
    pmf_vector
}

/// Scales a spectrogram by a scaling coefficient. If you provide a coefficient,
/// it will be used. Otherwise, the maximum value from the spectrogram will be used
/// as the scaling coefficient. This function will return the scaling coefficient
/// that was used, so you can reuse it later for consistency.
/// 
/// # Example
/// ```
/// use aus::spectrum;
/// use aus::analysis::normalize_spectrogram;
/// let file = aus::read("myfile.wav").unwrap();
/// let mut imaginary_spectrogram = spectrum::rstft(&file.samples[0], 2048, 1024, aus::WindowType::Hanning);
/// let (mut magnitude_spectrogram, mut phase_spectrogram) = spectrum::complex_to_polar_rstft(&imaginary_spectrogram);
/// let scaling_coef = normalize_spectrogram(&mut magnitude_spectrogram, None);
/// ```
#[inline]
pub fn normalize_spectrogram(spectrogram: &mut Vec<Vec<f64>>, scaling_coef: Option<f64>) -> f64 {
    let maxval: f64 = match scaling_coef {
        Some(val) => val,
        None => {
            let mut max = f64::NEG_INFINITY;
            for i in 0..spectrogram.len() {
                for j in 0..spectrogram[i].len() {
                    if spectrogram[i][j] > max {
                        max = spectrogram[i][j];
                    }
                }
            }
            max
        }
    };

    for i in 0..spectrogram.len() {
        for j in 0..spectrogram[i].len() {
            spectrogram[i][j] /= maxval;
        }
    }
    maxval
}

/// Scales a spectrum by a scaling coefficient. If you provide a coefficient,
/// it will be used. Otherwise, the maximum value from the spectrum will be used
/// as the scaling coefficient. This function will return the scaling coefficient
/// that was used, so you can reuse it later for consistency.
/// 
/// # Example
/// ```
/// use aus::spectrum::{rfft, complex_to_polar_rfft};
/// use aus::analysis::normalize_spectrum;
/// let fft_size: usize = 2048;
/// let mut pseudo_audio = vec![0.0, 0.1, 0.3, -0.4, 0.1, -0.51];
/// // zero-pad the audio
/// pseudo_audio.extend(vec![0.0; fft_size - pseudo_audio.len()]);
/// let spectrum = rfft(&pseudo_audio, fft_size);
/// let (mut magnitude_spectrum, _) = complex_to_polar_rfft(&spectrum);
/// let scaling_coef = normalize_spectrum(&mut magnitude_spectrum, None);
/// ```
#[inline]
pub fn normalize_spectrum(spectrum: &mut Vec<f64>, scaling_coef: Option<f64>) -> f64 {
    let maxval: f64 = match scaling_coef {
        Some(val) => val,
        None => {
            let mut max = f64::NEG_INFINITY;
            for i in 0..spectrum.len() {
                if spectrum[i] > max {
                    max = spectrum[i];
                }
            }
            max
        }
    };

    for i in 0..spectrum.len() {
        spectrum[i] /= maxval;
    }
    maxval
}

/// Calculates the spectral centroid from provided magnitude spectrum.
/// (Eyben, pp. 39-40)
///
/// $$
/// S_{centroid}=\frac{\sum_m F(m)X(m)}{\sum_m X(m)}
/// $$
/// where $F(m)$ is the frequency of bin $m$ in Hz.
/// 
/// # Example
///
/// ```
/// use aus::{spectrum, analysis};
/// let fft_size = 2048;
/// let audio = aus::read("myfile.wav").unwrap();
/// let imaginary_spectrum = spectrum::rfft(&audio.samples[0][..fft_size], fft_size);
/// let (magnitude_spectrum, phase_spectrum) = spectrum::complex_to_polar_rfft(&imaginary_spectrum);
/// let freqs = spectrum::rfftfreq(fft_size, audio.sample_rate);
/// let centroid = analysis::spectral_centroid(&magnitude_spectrum, &freqs);
/// ```
pub fn spectral_centroid(magnitude_spectrum: &[f64], rfft_freqs: &[f64]) -> f64 {
    let magnitude_spectrum_sum: f64 = magnitude_spectrum.iter().sum();
    compute_spectral_centroid(magnitude_spectrum, rfft_freqs, magnitude_spectrum_sum)
}

/// Computes the spectral difference between two STFT frames using the L2 norm.
/// You can optionally choose to only consider positive spectral differences in this calculation.
/// (Eyben, 42)
/// 
/// $$
/// SD^{(k)}=\sqrt{\sum_m \left(X^{(k)}(m)-X^{(k-1)}(m)\right)^2}
/// $$
/// 
/// # Example
///
/// ```
/// use aus::{spectrum, analysis, WindowType};
/// let fft_size = 2048;
/// let audio = aus::read("myfile.wav").unwrap();
/// let imaginary_spectrogram = spectrum::rstft(&audio.samples[0], fft_size, fft_size / 2, WindowType::Hanning);
/// let (magnitude_spectrogram, phase_spectrogram) = spectrum::complex_to_polar_rstft(&imaginary_spectrogram);
/// let difference = analysis::spectral_difference(&magnitude_spectrogram[0], &magnitude_spectrogram[1], false);
/// ```
pub fn spectral_difference(magnitude_spectrum1: &[f64], magnitude_spectrum2: &[f64], positive_only: bool) -> Result<f64, SpectrumError> {
    if magnitude_spectrum1.len() != magnitude_spectrum2.len() {
        return Err(SpectrumError{error_msg: String::from(format!("The provided spectra have different lengths: magnitude_spectrum1 has length {} but magnitude_spectrum2 has length {}.", magnitude_spectrum1.len(), magnitude_spectrum2.len()))});
    }
    let mut difference: f64 = 0.0;
    for i in 0..magnitude_spectrum1.len() {
        let local_difference = magnitude_spectrum2[i] - magnitude_spectrum1[i];
        if !positive_only || local_difference > 0.0 {
            difference += local_difference * local_difference;
        }      
    }
    Ok(f64::sqrt(difference))
}

/// Computes the spectral flux between two STFT frames.
/// You can choose which normalization coefficients will be used
/// (None, which corresponds to 1, or the $L^1$ or $L^2$ norm.)
/// (Eyben, 42-43)
/// 
/// $$
/// S_{flux}^{(k)}=\sum_m \left(\frac{X^{(k)}(m)}{\mu_k}-\frac{X^{(k-1)}(m)}{\mu_{k-1}}\right)^2
/// $$
/// where $\mu_k$ is a normalization coefficient.
/// 
/// # Example
///
/// ```
/// use aus::{spectrum, analysis, util, WindowType};
/// let fft_size = 2048;
/// let audio = aus::read("myfile.wav").unwrap();
/// let imaginary_spectrogram = spectrum::rstft(&audio.samples[0], fft_size, fft_size / 2, WindowType::Hanning);
/// let (magnitude_spectrogram, phase_spectrogram) = spectrum::complex_to_polar_rstft(&imaginary_spectrogram);
/// let flux = analysis::spectral_flux(&magnitude_spectrogram[0], &magnitude_spectrogram[1], Some(util::Norm::L2));
/// ```
pub fn spectral_flux(magnitude_spectrum1: &[f64], magnitude_spectrum2: &[f64], normalization_type: Option<Norm>) -> Result<f64, SpectrumError> {
    if magnitude_spectrum1.len() != magnitude_spectrum2.len() {
        return Err(SpectrumError{error_msg: String::from(format!("The provided spectra have different lengths: magnitude_spectrum1 has length {} but magnitude_spectrum2 has length {}.", magnitude_spectrum1.len(), magnitude_spectrum2.len()))});
    }
    let mut flux = 0.0;
    let mut norm1 = 1.0;
    let mut norm2 = 1.0;
    match normalization_type {
        Some(x) => {
            norm1 = lnorm(magnitude_spectrum1, &x);
            norm2 = lnorm(magnitude_spectrum2, &x);
        },
        None => ()
    };
    for i in 0..magnitude_spectrum1.len() {
        let local_difference = magnitude_spectrum2[i] / norm2 - magnitude_spectrum1[i] / norm1;
        flux += local_difference * local_difference;      
    }
    Ok(flux)
}

/// Calculates the spectral entropy from provided magnitude spectrum.
/// (Eyben, pp. 23, 40, 41)
///
/// $$
/// S_{entropy}=-\sum_m p_X(m)\log_2 p_X(m)
/// $$
/// where $p_X(m)$ is the spectrum power mass function
/// $$
/// p_X(m)=\frac{X(m)}{\sum_m X(m)}
/// $$
/// 
/// # Example
///
/// ```
/// use aus::{spectrum, analysis};
/// let fft_size = 2048;
/// let audio = aus::read("myfile.wav").unwrap();
/// let imaginary_spectrum = spectrum::rfft(&audio.samples[0][..fft_size], fft_size);
/// let (magnitude_spectrum, phase_spectrum) = spectrum::complex_to_polar_rfft(&imaginary_spectrum);
/// let entropy = analysis::spectral_entropy(&magnitude_spectrum);
/// ```
pub fn spectral_entropy(magnitude_spectrum: &[f64]) -> f64 {
    let power_spectrum = make_power_spectrum(magnitude_spectrum);
    let spectrum_pmf = make_spectrum_pmf(&power_spectrum, power_spectrum.iter().sum());
    compute_spectral_entropy(&spectrum_pmf)
}

/// Calculates the spectral flatness from provided magnitude spectrum.
/// (Eyben, p. 39, <https://en.wikipedia.org/wiki/Spectral_flatness>)
///
/// $$
/// S_{flatness}=\frac{m\sqrt[m]{\prod_m X(m)}}{\sum_m X(m)}
/// $$
/// 
/// # Example
///
/// ```
/// use aus::{spectrum, analysis};
/// let fft_size = 2048;
/// let audio = aus::read("myfile.wav").unwrap();
/// let imaginary_spectrum = spectrum::rfft(&audio.samples[0][..fft_size], fft_size);
/// let (magnitude_spectrum, phase_spectrum) = spectrum::complex_to_polar_rfft(&imaginary_spectrum);
/// let flatness = analysis::spectral_flatness(&magnitude_spectrum);
/// ```
pub fn spectral_flatness(magnitude_spectrum: &[f64]) -> f64 {
    let magnitude_spectrum_sum: f64 = magnitude_spectrum.iter().sum();
    compute_spectral_flatness(magnitude_spectrum, magnitude_spectrum_sum)
}

/// Calculates the spectral kurtosis from provided magnitude spectrum and real FFT frequency list.
/// (Eyben, pp. 23, 39-40)
///
/// $$
/// S_{kurtosis} = \frac{1}{S_{variance}^2} \sum_m \left(F(m)-S_{centroid}\right)^4 p_X(m)
/// $$
/// where $F(m)$ is the frequency corresponding to bin $m$, $S_{centroid}$ is the spectral centroid, $S_{variance}$ is the spectral variance, and $p_X(m)$ is the spectral power mass function
/// $$
/// p_X(m)=\frac{X(m)}{\sum_m X(m)}
/// $$
/// 
/// # Example
///
/// ```
/// use aus::{spectrum, analysis};
/// let fft_size = 2048;
/// let audio = aus::read("myfile.wav").unwrap();
/// let imaginary_spectrum = spectrum::rfft(&audio.samples[0][..fft_size], fft_size);
/// let (magnitude_spectrum, phase_spectrum) = spectrum::complex_to_polar_rfft(&imaginary_spectrum);
/// let freqs = spectrum::rfftfreq(fft_size, audio.sample_rate);
/// let kurtosis = analysis::spectral_kurtosis(&magnitude_spectrum, &freqs);
/// ```
pub fn spectral_kurtosis(magnitude_spectrum: &[f64], rfft_freqs: &[f64]) -> f64 {
    let power_spectrum = make_power_spectrum(magnitude_spectrum);
    let spectrum_pmf = make_spectrum_pmf(&power_spectrum, power_spectrum.iter().sum());
    let spectral_centroid = compute_spectral_centroid(magnitude_spectrum, rfft_freqs, magnitude_spectrum.iter().sum());
    let spectral_variance = compute_spectral_variance(&spectrum_pmf, rfft_freqs, spectral_centroid);
    compute_spectral_kurtosis(&spectrum_pmf, rfft_freqs, spectral_centroid, spectral_variance)
}

/// Calculates the spectral roll off frequency from provided magnitude spectrum, real FFT frequency list, and roll-off point.
/// The parameter `n` (0.0 <= n <= 1.00) indicates the roll-off point we wish to calculate.
/// (Eyben, p. 41)
/// 
/// $$
/// \sum_{m=0}^{r-1} X_P(m) \leq \frac{n}{100}\sum_m X_P(m)
/// $$
/// where $X_P(m)$ is the power spectrum.
///
/// # Example
///
/// ```
/// use aus::{spectrum, analysis};
/// let fft_size = 2048;
/// let audio = aus::read("myfile.wav").unwrap();
/// let imaginary_spectrum = spectrum::rfft(&audio.samples[0][..fft_size], fft_size);
/// let (magnitude_spectrum, phase_spectrum) = spectrum::complex_to_polar_rfft(&imaginary_spectrum);
/// let freqs = spectrum::rfftfreq(fft_size, audio.sample_rate);
/// let roll_off = analysis::spectral_roll_off_point(&magnitude_spectrum, &freqs, 0.75);
/// ```
pub fn spectral_roll_off_point(magnitude_spectrum: &[f64], rfft_freqs: &[f64], n: f64) -> f64 {
    let power_spectrum = make_power_spectrum(magnitude_spectrum);
    let power_spectrum_sum: f64 = power_spectrum.iter().sum();
    compute_spectral_roll_off_point(&power_spectrum, rfft_freqs, power_spectrum_sum, n)
}

/// Calculates the spectral skewness from provided magnitude spectrum and real FFT frequency list.
/// (Eyben, pp. 23, 39-40)
///
/// $$
/// S_{skewness} = \frac{1}{S_{variance}^{\frac{3}{2}}} \sum_m \left(F(m)-S_{centroid}\right)^3 p_X(m)
/// $$
/// where $F(m)$ is the frequency corresponding to bin $m$, $S_{centroid}$ is the spectral centroid, $S_{variance}$ is the spectral variance, and $p_X(m)$ is the spectral power mass function
/// $$
/// p_X(m)=\frac{X(m)}{\sum_m X(m)}
/// $$
/// 
/// # Example
///
/// ```
/// use aus::{spectrum, analysis};
/// let fft_size = 2048;
/// let audio = aus::read("myfile.wav").unwrap();
/// let imaginary_spectrum = spectrum::rfft(&audio.samples[0][..fft_size], fft_size);
/// let (magnitude_spectrum, phase_spectrum) = spectrum::complex_to_polar_rfft(&imaginary_spectrum);
/// let freqs = spectrum::rfftfreq(fft_size, audio.sample_rate);
/// let skewness = analysis::spectral_skewness(&magnitude_spectrum, &freqs);
/// ```
pub fn spectral_skewness(magnitude_spectrum: &[f64], rfft_freqs: &[f64]) -> f64 {
    let power_spectrum = make_power_spectrum(magnitude_spectrum);
    let spectrum_pmf = make_spectrum_pmf(&power_spectrum, power_spectrum.iter().sum());
    let spectral_centroid = compute_spectral_centroid(magnitude_spectrum, rfft_freqs, magnitude_spectrum.iter().sum());
    let spectral_variance = compute_spectral_variance(&spectrum_pmf, rfft_freqs, spectral_centroid);
    compute_spectral_skewness(&spectrum_pmf, rfft_freqs, spectral_centroid, spectral_variance)
}

/// Calculates the spectral slope from provided magnitude spectrum.
/// (Eyben, pp. 35-38)
///
/// # Example
///
/// ```
/// use aus::{spectrum, analysis};
/// let fft_size = 2048;
/// let audio = aus::read("myfile.wav").unwrap();
/// let imaginary_spectrum = spectrum::rfft(&audio.samples[0][..fft_size], fft_size);
/// let (magnitude_spectrum, phase_spectrum) = spectrum::complex_to_polar_rfft(&imaginary_spectrum);
/// let slope = analysis::spectral_slope(&magnitude_spectrum);
/// ```
pub fn spectral_slope(magnitude_spectrum: &[f64]) -> f64 {
    let power_spectrum = make_power_spectrum(magnitude_spectrum);
    let power_spectrum_sum: f64 = power_spectrum.iter().sum();
    compute_spectral_slope(&power_spectrum, power_spectrum_sum)
}

/// Calculates the spectral slope from provided magnitude spectrum, between the frequencies
/// specified. The frequencies specified do not have to align with the frequencies in the real FFT frequency list.
/// (Eyben, pp. 35-38)
///
/// # Example
///
/// ```
/// use aus::{spectrum, analysis};
/// let fft_size = 2048;
/// let audio = aus::read("myfile.wav").unwrap();
/// let imaginary_spectrum = spectrum::rfft(&audio.samples[0][..fft_size], fft_size);
/// let (magnitude_spectrum, phase_spectrum) = spectrum::complex_to_polar_rfft(&imaginary_spectrum);
/// let freqs = spectrum::rfftfreq(fft_size, audio.sample_rate);
/// let slope = analysis::spectral_slope_region(&magnitude_spectrum, &freqs, 105.0, 852.0, audio.sample_rate);
/// ```
pub fn spectral_slope_region(magnitude_spectrum: &[f64], rfft_freqs: &[f64], f_lower: f64, f_upper: f64, sample_rate: u32) -> f64 {
    let power_spectrum = make_power_spectrum(magnitude_spectrum);
    compute_spectral_slope_region(&power_spectrum, rfft_freqs, f_lower, f_upper, sample_rate)
}

/// Calculates the spectral variance from provided magnitude spectrum and real FFT frequency list.
/// (Eyben, pp. 23, 39-40)
///
/// $$
/// S_{variance} = \sum_m \left(F(m)-S_{centroid}\right)^2 p_X(m)
/// $$
/// where $F(m)$ is the frequency corresponding to bin $m$, $S_{centroid}$ is the spectral centroid, and $p_X(m)$ is the spectral power mass function
/// $$
/// p_X(m)=\frac{X(m)}{\sum_m X(m)}
/// $$
/// 
/// # Example
///
/// ```
/// use aus::{spectrum, analysis};
/// let fft_size = 2048;
/// let audio = aus::read("myfile.wav").unwrap();
/// let imaginary_spectrum = spectrum::rfft(&audio.samples[0][..fft_size], fft_size);
/// let (magnitude_spectrum, phase_spectrum) = spectrum::complex_to_polar_rfft(&imaginary_spectrum);
/// let freqs = spectrum::rfftfreq(fft_size, audio.sample_rate);
/// let variance = analysis::spectral_variance(&magnitude_spectrum, &freqs);
/// ```
pub fn spectral_variance(magnitude_spectrum: &[f64], rfft_freqs: &[f64]) -> f64 {
    let power_spectrum = make_power_spectrum(magnitude_spectrum);
    let spectrum_pmf = make_spectrum_pmf(&power_spectrum, power_spectrum.iter().sum());
    let spectral_centroid = compute_spectral_centroid(magnitude_spectrum, rfft_freqs, magnitude_spectrum.iter().sum());
    compute_spectral_variance(&spectrum_pmf, rfft_freqs, spectral_centroid)
}

#[cfg(test)]
mod test {
    use super::*;
    use std::io;
    use std::io::Write;

    const AUDIO: &str = "myfile.wav";
    const FFT_SIZE: usize = 2048;

    #[test]
    fn test_log_spec() {
        const EPSILON: f64 = 1e-6;
        let in_vec: Vec<f64> = vec![6.76895304e-06, 1.36028451e-06, 2.77367273e-06, 1.89267587e-05,
            4.07858842e-04, 3.19873457e-04, 5.04253261e-06, 8.52845397e-06,
            4.99160938e-06, 5.40207921e-06, 3.29421796e-04, 7.69145971e-05,
            6.61126227e-06, 2.21605018e-04, 4.42958155e-05, 1.83195179e-04,
            9.42967421e-05, 1.23944028e-04, 1.34227360e-04, 1.39157067e-04];

        // this vector generated by Python code
        let expected_output = vec![
            -97.7998838, -104.7688013, -101.67454666, -93.33433635, -80.0,
            -81.05531678, -99.07861167, -96.79639572, -99.1226929, -98.77948934,
            -80.92757551, -87.24501112, -97.90225495, -82.64930292, -89.6414718,
            -83.47595841, -86.36013193, -85.17284276, -84.82668833, -84.67004614];
        let output = make_log_spectrum(&in_vec, 10.0, -10e8, Some(-80.0));

        for i in 0..output.len() {
            assert!(f64::abs(expected_output[i] - output[i]) < EPSILON);
        }
    }
    
    #[test]
    fn test_autocorrelation() {
        let audio = crate::read(AUDIO).unwrap();
        let auto = autocorrelation(&audio.samples[0][44100..44100+FFT_SIZE], FFT_SIZE).unwrap();
        
        for val in auto.iter() {
            print!("{} ", val);
            io::stdout().flush().unwrap();
        }

        let mut max = auto[0];
        let mut max_idx: usize = 0;
        for i in 0..auto.len() {
            if auto[i] > max {
                max = auto[i];
                max_idx = i;
            }
        }

        println!("\nMax: {}, max_idx: {}", max, max_idx);
    }
}