informationEOC.lib

// =============================================================================
// ========== informationEOC.lib ===============================================
// =============================================================================
// 
// Information processing functions library including low-level and high-level
// algorithms both based on hard-coded and adaptive mechanisms. The
// low-level functions provide time-domain techniques for feature
// extraction that are normally based on FFT processing, such as
// spectral centroid and spectral flatness (noisiness). The high-level
// functions provide an analysis of the state space of low-level
// information signals to determine, based on notions of complexity theory and
// music perception, characteristics such as dynamicity, heterogeneity,
// and complexity of audio streams.
// 
// All functions below that use one-pole lowpass filters for
// accumulation are based on a time constant that is 2pi*tau.
// 
// The environment prefix is "ip".
// 
// List of functions:
// 
//    a_weighting, 
//    complexity, 
//    dynamicity, 
//    env_fol, 
//    highest_partial, 
//    heterogeneity, 
//    heterogeneity10, 
//    instant_amp, 
//    instant_freq, 
//    instant_ph, 
//    loudness, 
//    lowest_partial, 
//    mle,
//    noisiness, 
//    peaks, 
//    peak_env,
//    peak_env_AHR_cascade,
//    peak_env_AHR_switch,
//    peak_env_AR_cascade, 
//    peak_env_AR_switch, 
//    peak_hold,
//    peak_hold_H,
//    peak_hold_L,
//    peak_hold_LH,
//    recurrence, 
//    rms, 
//    rms4, 
//    roughness, 
//    spec_bal, 
//    spec_balN, 
//    spec_peakN, 
//    spec_ten, 
//    spec_tenN, 
//    spec_ten_lite, 
//    zc, 
//    zcr, 
//    zcr4.
// 
// Copyright (c) 2019-2020, Dario Sanfilippo <sanfilippo.dario at gmail dot com>
// All rights reserved.

declare name "Information Processing Library";
declare author "Dario Sanfilippo";
declare copyright "Copyright (c) 2019-2020, Dario Sanfilippo <sanfilippo.dario
      at gmail dot com>";
declare version "2.1.0";
declare license "GPLv2.0";

au = library("auxiliary.EOClib");
ba = library("basics.lib");
fi = library("filters.lib");
f2 = library("filtersEOC.lib");
ip = library("informationEOC.lib");
ma = library("maths.lib");
m2 = library("mathsEOC.lib");
ro = library("routes.lib");
si = library("signals.lib");

// ip.a_weighting(F[n]); -------------------------------------------------------
//
// Function for equal-loudness contour based on the A-weighting equation. 
// The function is calibrated to output unity (0 dB) at 1000 Hz.
//
// 1 inputs:
//    F[n], frequency in Hz;
//
// 1 outputs:
//    y[n], linear amplitude at F[n].
//
a_weighting(f) = (12194 ^ 2 * f ^ 4) / ((f ^ 2 + 20.6 ^ 2) * 
      sqrt((f ^ 2 + 107.7 ^ 2) * (f ^ 2 + 737.9 ^ 2)) * 
           (f ^ 2 + 12194 ^ 2)) + .206;
// -----------------------------------------------------------------------------

// ip.complexity(max_dt, dt[n], R[n], x[n]); -----------------------------------
//
// Complexity measurement as edge-of-chaos detection, that is,  accumulation 
// of the variations in the heterogeneity of a state space. The function
// outputs an index in the [0; 1] range. The input signal is assumed in
// the [0; 1] range, and the resolution is of 10 sub-regions in the state space.
//
// 3 inputs:
//    dt[n], differentiation period in seconds;
//    R[n], responsiveness of the system in Hz;
//    x[n].
//
// 1 outputs:
//    y[n], complexity index for x[n].
//
// 1 compile-time arguments:
//    max_dt, maximum differentiation period in seconds.
//
complexity(max_dt, dt, window, x) = ip.heterogeneity10(window, x) :
      abs(m2.delta(max_dt, dt)) : f2.lp1p(1 / max(ma.EPSILON, window));
// -----------------------------------------------------------------------------

// ip.dynamicity(max_dt, dt[n], R[t], x[n]); -----------------------------------
//
// Dynamicity index as accumulation of the magnitude of the delta at a 
// specified dt (secs) and accumulation rate (Hz).
// It is supposed to be used with infrasonic low-level information signals such 
// as RMS, spectral tendency, noisiness.
//
// 3 inputs:
//    dt[n], differentiation period in seconds;
//    R[n], responsiveness of the system in Hz;
//    x[n].
//
// 1 outputs:
//    y[n], dynamicity index for x[n].
//
// 1 compile-time arguments:
//    max_dt, maximum differentiation period in seconds.
//
dynamicity(max_dt, dt, window, in) = abs(m2.delta(max_dt, dt, in)) :
      f2.lp1p(window);
// -----------------------------------------------------------------------------

// ip.env_fol(R[n], x[n]); -----------------------------------------------------
//
// Envelope following (Dodge and Jerse).
//
// 2 inputs:
//    R[n], responsiveness in Hz;
//    x[n].
//
// 1 outputs:
//    y[n], amplitude envelope of x[n].
//
env_fol(window, in) = f2.lp1p(window, abs(in));
// -----------------------------------------------------------------------------

// ip.highest_partial(R[n], x[n]); ---------------------------------------------
//
// It detects the highest partial in a signal. This function is an
// extension of the spectral tendency algorithm by including a positive
// feedback loop at the top of the chain: the system recursively removes 
// low-frequency componenets through a highpass until no components are
// left on the upper side of the spectrum except the last partial.
// The SNR and sensitivity of the filter, approximately, could be adjusted
// by changing the order of the filters.
// 
// 2 inputs:
//    R[n], responsiveness in Hz;
//    x[n].
// 
// 1 outputs:
//    y[n], frequency of the highest partial in x[n].
// 
highest_partial(window, in) =  
      (_ <:    _ ,
               (   _ , 
                   in : hp1p1z2 <: _ ,
                                   _) : (m2.div(specbal, ip.rms(window)) / 
      ma.SR * window : f2.int_clip(0, 1) ^ 2)) 
      ~ * (m2.ny)
      with {
           specbal = f2.xover1p1z :    ip.rms(window) + .000001 , 
                                       ip.rms(window) 
               : ro.cross(2) : -;
           hp1p1z2(cf, x) = x : seq(i, 2, f2.hp1p1z(cf));
      };
// -----------------------------------------------------------------------------

// ip.heterogeneity(N, L[n], H[n], T[n], x[n]); --------------------------------
//
// Provides an index of how heterogeneous the state space of a signal is, which
// correlates to the unpredictability of the signal. The algorithm
// follows principles of recurrence quantification analysis by dividing the
// state space into sub-regions of the same size. If the signal is in a
// sub-region, then an integrator is triggered to keep track of the
// recurrence of that state. When the signal moves to a different
// sub-region, then the recurrence value in the non-active sub-region is
// held for a specified time before being reset to 0. The output of the
// peak holder is smoothed out by a one-pole lowpass with the same period
// as the peak holder for attack and decay smoothing. The final
// heterogeneity index is then computed by processing the outputs of the
// sub-regions recurrences through normalised average absolute deviation (AAD). 
// AAD is used instead of standard deviation to have an output in the range 
// [0; 1].
//
// 4 inputs:
//    L[n], lower edge of the state space;
//    H[n], upper edge of the state space;
//    T[n], analysis period, it sets the memory of the system, in
//         seconds, determining hold and decay times;
//    x[n].
//
// 1 outputs:
//    y[n], heterogeneity index for x[n].
//
// 1 compile-time arguments:
//    N, state space resolution, that is, (integer) number of
//         equally-spaced regions in the state space.
//
heterogeneity(N, l, h, t, in) = 
      par(i, N,    region(l, h, N, i) ,
                   t ,
                   in : ip.recurrence : ip.peak_hold(t) : 
                       f2.lp1p(m2.div(1, t))) : m2.aad(N)
      with {
           region(lower, upper, N, i) = 
               lower + ((upper - lower) / N) * i ,
               lower + ((upper - lower) / N) * (i + 1);
      };
// -----------------------------------------------------------------------------

// ip.heterogeneity10(T[n], x[n]); ---------------------------------------------
//
// Shortcut for a heterogeneity index for a 10-region state space in a [0; 1] 
// range.
//
// 2 inputs:
//    T[n], analysis period, it sets the memory of the system, in
//         seconds, determining hold and decay times;
//    x[n].
//
// 1 outputs:
//    y[n], heterogeneity index for x[n] with a 10-region resolution in the 
//    state space, assuming that its range is [0; 1].
//
heterogeneity10(t, in) = ip.heterogeneity(10, 0, 1, t, in);
// -----------------------------------------------------------------------------

// ip.instant_amp(x[n]); -------------------------------------------------------
//
// Instantaneous amplitude.
//
// 1 inputs:
//    x[n].
//
// 1 outputs:
//    y[n], instantaneous amplitude of x[n], sqrt(re^2 + im^2).
//
instant_amp(x) = f2.analytic(x) : sqrt(pow(2) + pow(2));
// -----------------------------------------------------------------------------

// ip.instant_freq(x[n]); ------------------------------------------------------
//
// Instantaneous frequency.
//
// 1 inputs:
//    x[n].
//
// 1 outputs:
//    y[n], instantaneous frequency of x[n] as derivative of the instantaneous
//    phase.
//
instant_freq(x) = m2.diff(instant_ph(x)) / (2 * ma.PI) * ma.SR;
// -----------------------------------------------------------------------------

// ip.instant_ph(x[n]); --------------------------------------------------------
//
// Instantaneous phase.
//
// 1 inputs:
//    x[n].
//
// 1 outputs:
//    y[n], instantaneous phase of x[n] as the arcotangent of the ratio
//         between its imaginary and real parts.
//
instant_ph(x) = f2.analytic(x) : ro.cross(2) : atan2;
// -----------------------------------------------------------------------------

// ip.loudness(R[n], x[n]); ----------------------------------------------------
//
// Loudness measurement through A-weighting function and spectral centroid.
//
// 1 inputs:
//    R[n], responsiveness of the system in Hz;
//    x[n].
//
// 1 outputs:
//    y[n], loudness of x[n] as linear amplitude.
//
loudness(window, in) = 
      peak_env(1 / max(ma.EPSILON, window), in) * 
           (spec_ten(window, in) * m2.ny : a_weighting);
// -----------------------------------------------------------------------------

// ip.lowest_partial(R[n], x[n]); ----------------------------------------------
//
// It detects the lowest partial in a signal. This function is an
// extension of the spectral tendency algorithm by including a positive
// feedback loop at the top of the chain: the system recursively removes 
// high-frequency componenets through a lowpass until no components are
// left on the lower side of the spectrum except the last partial.
// The SNR and sensitivity of the filter, approximately, could be adjusted
// by changing the order of the filters.
// 
// 2 inputs:
//    R[n], responsiveness in Hz;
//    x[n].
// 
// 1 outputs:
//    y[n], frequency of the lowhest partial in x[n].
// 
lowest_partial(window, in) =   
      (_ <:    _ ,
               (   _ , 
                   in : lp1p1z2 <: _ ,
                                   _) : 
      (specbal / max(ma.EPSILON, rms(window)) / ma.SR * window : 
           f2.int_clip(0, 1) ^ 2)) 
      ~ (* (m2.ny) : max(10))
      with {
           specbal = f2.xover1p1z :    rms(window) , 
                                       rms(window) + .000001 
               : ro.cross(2) : -;
           lp1p1z2(cf, x) = x : seq(i, 2, f2.lp1p1z(cf));
      };
// -----------------------------------------------------------------------------

// ip.lyapunov_exp(T[n], x[n]); ------------------------------------------------
//
// Lyapunov Exponent as measure of sensitivity to initial conditions
// and unpredictability.
//
// 2 inputs:
//      T[n], analysis period in seconds;
//      x[n], input signal.
//
// 1 outputs:
//      y[n], MLE of x[n].
//
lyapunov_exp(window, x) = 
    f2.lp1p(1.0 / max(ma.EPSILON, window), 
        log(max(ma.EPSILON, abs(m2.diff(x)))));
// -----------------------------------------------------------------------------


// ip.noisiness(R[n], x[n]); ---------------------------------------------------
//
// Noisiness index normalised for a 96 kHz samplerate, that is, the
// output of the function is 1 for white noise, and 0 (or very close to 0)
// for sinusoids. Unlike noisiness drectly calculated through zero-crossing
// rate (ZCR), this algorithm measures noisiness based on the derivative of the
// ZCR, as non-periodicity characterises noisy signals, whereas
// high-frequency sinusoidal signals can still have a high ZCR.
//
// 2 inputs:
//    R[n], responsiveness in Hz;
//    x[n].
// 
// 1 outputs:
//    y[n], noisiness index of x[n] in the [0; 1] range.
// 
noisiness(window, in) = spec_ten(window, in) *  m2.ny : max(1) <:  
      m2.inv ,
      (    _ , 
           in : zcr4 : m2.unit_log(10)) : abs(m2.delta(1)) : f2.lp1p(window) / 
               .266 <: * : min(1);
// -----------------------------------------------------------------------------

// ip.peak_env(R[n], x[n]); ----------------------------------------------------
//
// Peak envelope function: infinitely fast attack and adjustable release.
// This function outputs the absolute value of the input if the
// magnitude of the signal is greater or equal than the output, 
// or it performs an exponential decay of the last detected peak 
// when the input is smaller than the output.
//
// 2 inputs:
//    R[n], release time in seconds;
//    x[n].
//
// 1 outputs:
//    y[n], peak envelope of x[n].
//
peak_env(RT, x) =  max(abs(x))
                   ~ * (m2.rt55(RT));
// -----------------------------------------------------------------------------

// ip.peak_env_AHR_cascade(AT[n], HT[n], RT[n], x[n]); -------------------------
//
// Peak envelope function with attack, hold, and release times in
// seconds. The effective attack, hold, and release times are
// interrelated as this design is based on cascaded filters, particularly,
// a peak holder (peak_hold) feeding into a peak_envelope (peak_env) feeding 
// into a one-pole lowpass (smooth). Considering the tau*2π time constant
// where most of the final value is reached after the attack time, we will
// only have a hold segment in the resulting envelope function if the hold
// time is greater than the attack time, and the hold segment will be
// approximately the hold time minus the attacktime. Alternatively, only 
// attack and release segments will result from this cascaded design. 
// Furthermore, the attack time should still be << than the release time to 
// minimally affect the decay.
//
// 4 inputs:
//    AT[n], attack time in seconds;
//    HT[n], hold time in seconds;
//    RT[n], release time in seconds;
//    x[n].
//
// 1 outputs:
//    y[n], envelope profile of x[n].
//
peak_env_AHR_cascade(AT, HT, RT, x) = 
      ip.peak_hold(HT, x) : ip.peak_env_AR_cascade(AT, RT);
// -----------------------------------------------------------------------------

// ip.peak_env_AHR_switch(AT[n], HT[n], RT[n], x[n]); --------------------------
//
// Envelope function with attack, hold, and release times in seconds.
// This design is based on switching filter sections depending on the
// attack or release phases of the input signal, which is determined by
// comparison between the output and the input. When the input is greater
// than the output, the system switches to a one-pole lowpass. Otherwise,
// the system switches to a peak holder (peak_hold) section, and finally
// into a peak envelope (peak_env) section when the hold time has passed and 
// no new peaks have been detected. This way, the resulting envelope curve
// will always include attack, hold, and release segments at the specified
// times, where the hold value is dependent on the value that has been reached 
// in the attack section.
//
// 4 inputs:
//    AT[n], attack time in seconds;
//    HT[n], hold time in seconds;
//    RT[n], release time in seconds;
//    x[n].
//
// 1 outputs:
//    y[n], envelope profile of x[n].
//
peak_env_AHR_switch(AT, HT, RT, x) =   loop
                                       ~ _
      with {
           loop(fb) = m2.if(cond1, 
                            attack, 
                            m2.if(  cond2,
                                    hold, 
                                    release))
               with {
                   cond1 = abs(x) >= fb;
                   cond2 = fi.pole(notcond1, notcond1) <= rint(HT * ma.SR);
                   notcond1 = 1 - cond1;
                   attack = abs(x) * (1 - AT_coeff) + fb * AT_coeff;
                   hold = fb;
                   release = RT_coeff * fb;
                   AT_coeff = m2.rt55(AT);
                   RT_coeff = m2.rt55(RT);
               };
      };
// -----------------------------------------------------------------------------

// ip.peak_env_AR_cascade(AT[n], RT[n], x[n]); ---------------------------------
//
// Peak envelope function with dependent attack and release times. The
// function applies two cascaded filter sections regardless of the
// attack or release phases of the input signal: a peak envelope filter 
// feeding into a one-pole lowpass filter. The attack time is assumed to 
// be << than the release time for best behaviour. Since the sections are
// cascaded, the effective release time will always be greater than the
// desired one. The difference will be less noticeable when attack and
// release times are far apart. This function provides a smoother
// transition between attack and release phases.
//
// 3 inputs:
//    AT[n], attack time in seconds;
//    RT[n], release time in seconds;
//    x[n].
//
// 1 outputs:
//    y[n], peak envelope of x[n].
//
peak_env_AR_cascade(AT, RT, x) = ip.peak_env(RT, x) : si.smooth(m2.rt55(AT));
// -----------------------------------------------------------------------------

// ip.peak_env_AR_switch(AT[n], RT[n], x[n]); ----------------------------------
//
// Peak envelope function with independent attack and release times. The
// function applies two independent filter sections based on whether 
// the input is in the attack or release phase, which is determined by 
// comparing the absolute value of the input signal with the output. 
// The attack-smoothing filter is a one-pole lowpass: a leaky integrator
// whose input is scaled down with the complement of the pole position.
// The release-smoothing filter, on the other hand, is simply an
// exponential decay of the detected peak.
//
// Ref: (2008) Udo Zölzer – Digital Audio Signal Processing 2nd Edition.
//
// 3 inputs:
//    AT[n], attack time in seconds;
//    RT[n], release time in seconds;
//    x[n].
//
// 1 outputs:
//    y[n], peak envelope of x[n].
//
peak_env_AR_switch(AT, RT, x) =    loop
                                   ~ _
      with {
           loop(fb) = m2.if(abs(x) > fb, attack, release)
               with {
                   attack = abs(x) * (1 - AT_coeff) + fb * AT_coeff;
                   release = fb * RT_coeff;
                   AT_coeff = m2.rt55(AT);
                   RT_coeff = m2.rt55(RT);
               };
      };
// -----------------------------------------------------------------------------

// ip.peak_hold(HT[n], x[n]); --------------------------------------------------
//
// Peak holder: it holds the peak of the absolute value of the input
// signal for a time specified in seconds. If no new peak occurs before
// the specified time, the peak will reset to whatever the absolute value
// of the input is. A new peak is detected if the input is greater or
// equal to the output.
//
// Note that the countdown for the hold time starts after a peak
// disappears. For example, if the function is tested with a step response
// of 20 samples, if the hold time is set to 48 samples (.001 seconds at
// SR = 48 kHz), we will have 1s from sample  n = 0 to n = 19 from the step
// function, 1s from n = 20 to n = 67 from the hold time, and the output will
// then be 0 from sample n = 68.
//
// 2 inputs:
//    HT[n], hold time in seconds (resulting in the closest integer
//         number of samples representing that length);
//    x[n].
//
// 1 outputs:
//    y[n], high peak of |x[n]|.
//
peak_hold(HT, x) = ip.peak_hold_H(HT, abs(x));
// -----------------------------------------------------------------------------

// ip.peak_hold_H(HT[n], x[n]); ------------------------------------------------
//
// High peak holder: it holds the highest peak of the input signal for a time 
// specified in seconds. If no new peak occurs before the specified time, 
// the peak will reset to whatever the value of the input is. A new peak is 
// detected if the input is greater or equal to the output.
//
// Note that the countdown for the hold time starts after a peak
// disappears. For example, if the function is tested with a step response
// of 20 samples, if the hold time is set to 48 samples (.001 seconds at
// SR = 48 kHz), we will have 1s from sample  n = 0 to n = 19 from the step
// function, 1s from n = 20 to n = 67 from the hold time, and the output will
// then be 0 from sample n = 68.
//
// 2 inputs:
//    HT[n], hold time in seconds (resulting in the closest integer 
//         number of samples representing that length);
//    x[n].
//
// 1 outputs:
//    y[n], high peak of x[n].
//
peak_hold_H(HT, x) =   loop_peak 
                       ~ _
      with {
           l = rint(HT * ma.SR);
           loop_peak(fb_p) = m2.if(cond1 | notcond2, x, fb_p)
               with {
                   cond1 = x >= fb_p;
                   cond2 = loop_timer <= l 
                           ~ _
                       with {
                           loop_timer(fb_t) = notcond1 & fb_t <: fi.pole;
                       };
                   notcond1 = 1 - cond1;
                   notcond2 = 1 - cond2;

               };
      };
// -----------------------------------------------------------------------------

// ip.peak_hold_L(HT[n], x[n]); ------------------------------------------------
//
// Low peak holder: it holds the lowest peak of the input signal for a time 
// specified in seconds. If no new peak occurs before the specified time, 
// the peak will reset to whatever the value of the input is. A new peak is 
// detected if the input is less or equal to the output.
//
// Note that the countdown for the hold time starts after a peak
// disappears. For example, if the function is tested with a step response
// of 20 samples, if the hold time is set to 48 samples (.001 seconds at
// SR = 48 kHz), we will have 1s from sample  n = 0 to n = 19 from the step
// function, 1s from n = 20 to n = 67 from the hold time, and the output will
// then be 0 from sample n = 68.
//
// 2 inputs:
//    HT[n], hold time in seconds (resulting in the closest integer
//         number of samples representing that length);
//    x[n].
//
// 1 outputs:
//    y[n], low peak of x[n].
//
peak_hold_L(HT, x) =   loop_peak
                       ~ _
      with {
           l = rint(HT * ma.SR);
           loop_peak(fb_p) = m2.if(cond1 | notcond2, x, fb_p)
               with {
                   cond1 = x <= fb_p;
                   cond2 = loop_timer <= l
                           ~ _
                       with {
                           loop_timer(fb_t) = notcond1 & fb_t <: fi.pole;
                       };
                   notcond1 = 1 - cond1;
                   notcond2 = 1 - cond2;

               };
      };

// -----------------------------------------------------------------------------

// ip.peak_hold_LH(HT[n], x[n]); -----------------------------------------------
//
// It detects lower and upper peaks and, if no new peaks are detected,
// it holds them for a specified time in seconds before resetting the peaks
// to the new ones given by the incoming input.
//
// 2 inputs:
//    HT[n], hold time in seconds (resulting in the closest integer
//         number of samples representing that length);
//    x[n].
//
// 2 outputs:
//    y1[n], lower peak in x[n],
//    y2[n], upper peak in x[n].
//
peak_hold_LH(HT, x) =  peak_hold_L(HT, x) ,
                       peak_hold_H(HT, x);
// -----------------------------------------------------------------------------

// ip.recurrence(L[n], H[n], T[n], x[n]); --------------------------------------
//
// It detects and integrates occurrences of signals within a specified
// range using a leaky integrator with a specified decay time in seconds.
//
// 4 inputs:
//    L[n], lower bound (not including the edge for the detection);
//    H[n], upper bound (including the edge for the detection);
//    T[n], decay period in seconds;
//    x[n].
//
// 1 outputs:
//    y[n], recurrence of the signal x[n] within the specified range
//         ]lower; upper].
//
recurrence(lower, upper, t, in) =  in > lower ,
                                   in <= upper : & : fi.pole(m2.rt55(t));
// -----------------------------------------------------------------------------

// ip.rms(R[n], x[n]); ---------------------------------------------------------
//
// Root mean square measurement using one-pole lowpass for averaging.
//
// 2 inputs:
//    R[n], analysis rate in Hz;
//    x[n].
//
// 1 outputs:
//    y[n], RMS of x[n].
//
rms(window, in) = in <: * : f2.lp1p(window) : sqrt;
// -----------------------------------------------------------------------------

// ip.rms4(R[n], x[n]); --------------------------------------------------------
//
// RMS with four cascaded one-pole lowpass for averaging.
//
// 2 inouts:
//    R[n], analysis rate in Hz;
//    x[n].
//
// 1 outputs:
//    y[n], RMS of x[n].
//
rms4(window, in) = in <: * : seq(i, 4, f2.lp1p(window)) : sqrt;
// -----------------------------------------------------------------------------

// ip.roughness(R[n], x[n]); ---------------------------------------------------
//
// The roughness measurement is based on amplitude transients in the 
// 15-75 Hz range. Theoretically, roughness decreases below 15Hz, 
// peaks around 15-75 Hz, and starts descresing above 75Hz to 
// disappear around 150 Hz and above. The analysis rate parameter may affect
// the behaviour in the 0-15Hz range. An analysis rate of 1 Hz seems acceptable.
// References: [Helmholtz 2013; Vassilakis and Kendall 2010; Terhardt
// 1974; Molla and Torrésani 2004; Moore 1995].
//
// 2 inputs:
//    R[n], responsiveness in Hz;
//    x[n].
//
// 1 outputs:
//    y[n], roughness index of x[n] in the range [0; 1].
//
roughness(window, in) = 
      instant_amp(in) <: m2.div(m2.delta(.001, .001), rms(1)) : 
           rms4(75) : m2.delta(1 / 150, 1 / 150) : rms4(window) / .807242;
// -----------------------------------------------------------------------------

// ip.spec_bal(CF[n], R[n], x[n]); ---------------------------------------------
//
// Spectral power (RMS) difference at a splitting point, in Hz, of the
// spectrum of a signal. The input signal is filtered through a one-pole
// crossover; the difference of the RMS of the high and low spectra is the
// output.
//
// 3 inputs:
//    CF[n], crossover cut-off frequency in Hz;
//    R[n], RMS analysis rate in Hz;
//    x[n].
//
// 1 outputs:
//    y[n], RMS difference of the crossover outputs.
//
spec_bal(cf, window, in) = f2.xover1p1z(cf, in) :  rms(window) ,
                                                   rms(window) : 
      ro.cross(2) : -;
// -----------------------------------------------------------------------------

// ip.spec_balN(N, CF[n], R[n], x[n]); -----------------------------------------
//
// Spectral power (RMS) difference at a splitting point, in Hz, of the
// spectrum of a signal. The input signal is filtered through an
// Nth-order Butterworth crossover; the difference of the RMS of the high and 
// low spectra is the output.
//
// 3 inputs:
//    CF[n], crossover cut-off frequency in Hz;
//    R[n], RMS analysis rate in Hz;
//    x[n].
//
// 1 outputs:
//    y[n], RMS difference of the crossover outputs.
//
// 1 compile-time arguments:
//    N, (integer) order of the crossover.
//
spec_balN(N, cf, window, in) = f2.xover_butt(N, cf, in) :  rms(window) ,
                                                           rms(window) : 
      ro.cross(2) : -;
// -----------------------------------------------------------------------------

// ip.spec_peakN(N, R[n], x[n]); -----------------------------------------------
//
// Spectral peak with Nth-order spectral tendency calculation and
// 2nd-order normalised BP filter (BLTI). It works by recursively
// removing components on the sides of the spectral tendency frequency.
// Currently a prototype as it still requires a mechanism to detect new
// appearing signals that happen to fall within the filtered band.
// One possibility is to drive the Q through the complement of the
// derivative of a spectral tendency measure at the top of the chain. Still
// a work in progress but it works well to isolate the frequency of pure
// tones in a noisy background.
//
// 2 inputs:
//    R[n], responsiveness in Hz;
//    x[n].
//
// 1 outputs:
//    y[n], spectral peak in Hz.
//
// 1 compile-time arguments:
//    N, (integer) order of the crossover in the spectral centroid
//         analysis.
spec_peakN(N, window, in) = ((max(10), 1, in) : f2.bp2blti : 
      spec_tenN(N, window)) 
                            ~ * (m2.ny);
// -----------------------------------------------------------------------------

// ip.spec_sprN(N, Q[n], R[n], x[n]); ------------------------------------------
//
// Spectral spread: complement of the ratio between the RMS of the 
// spectral-centroid-centered bandpassed input and the input RMS. Nth-order
// spectral tendency, arbitrary Q-factor for the bandpass to tune-in
// specific applications.
//
// 3 inputs:
//    Q[n], Q-factor in the bandpass filter;
//    R[n], RMS analysis rate in Hz;
//    x[n].
//
// 1 outputs:
//    y[n], spectral spread index for x[n] in the range [0; 1] (0 no
//         spread; 1 maximum spread).
//
// 1 compile-time arguments:
//    N, (integer) order of the Butterworth crossover in the spectral
//         centroid analysis.
//
spec_sprN(N, q, window, in) =  (   spec_tenN(N, window, in) * m2.ny ,
                                   q ,
                                   in : f2.bpbi :  ip.rms(window)) ,
                               ip.rms(window, in) : m2.div : m2.complement;
// -----------------------------------------------------------------------------

// ip.spec_ten(R[n], x[n]); ----------------------------------------------------
//
// Spectral tendency (we can call it centroid as it finds a balancing
// point) equal-power spectral split-point. The adaptive algorithm finds
// the cut-off frequency of a crossover where the RMS difference of its
// outputs is 0 through a negative feedback mechanism. To be precise, the
// algorithm minimally oscillates, hence it is in dynamical equilibrium,
// around the equilibrium point. The input is filtered through a
// one-pole-one-zero crossover; the difference of the RMS of the outputs of the 
// filter is integrated; the result is squared and clipped for stability, mapped 
// over Nyquist, and fed back to drive the frequency of the crossover.
//
// 2 inputs:
//    R[n], RMS analysis rate in Hz (responsiveness);
//    x[n].
//
// 1 outputs:
//    y[n], spectral centroid of x[n] as an index in the [0; 1] range.
//
spec_ten(window, in) =  (   (   _ ,
                                window ,
                                in : spec_bal) ,
                            (   window ,
                                in : rms) : m2.div * window : 
      f2.int_clip(0, 1) ^ 2) 
                        ~ * (m2.ny);
// -----------------------------------------------------------------------------

// ip.spec_tenN(N, R[n], x[n]); ------------------------------------------------
//
// Spectral tendency (we can call it centroid as it finds a balancing
// point) equal-power spectral split-point. The adaptive algorithm finds
// the cut-off frequency of a crossover where the RMS difference of its
// outputs is 0 through a negative feedback mechanism. To be precise, the
// algorithm minimally oscillates, hence it is in dynamical equilibrium,
// around the equilibrium point. The input is filtered through an Nth-order 
// butterworth crossover; the difference of the RMS of the outputs of the
// filter is integrated; the result is squared and clipped for stability, mapped
// over Nyquist, and fed back to drive the frequency of the crossover.
// Spectral tendency (centroid) as equal-power spectral split-point with nth-order
// This design is improved as it implements a cubic nonlinear function
// that increases stability and accuracy.
//
// 2 inputs:
//    R[n], RMS analysis rate in Hz (responsiveness);
//    x[n].
//
// 1 outputs
//    y[n], spectral centroid of x[n] as an index in the [0; 1] range.
//
// 1 compile-time arguments:
//    N, (integer) order of the crossover.
//
spec_tenN(N, window, in) = ((  _ ,
                               window ,
                               in : spec_balN(N)) ,
                            (  window ,
                               in : rms) : m2.div ^ 3 : 
      f2.int_eu_clip(0, 1, window) ^ 2)
                           ~ * (m2.ny);
// -----------------------------------------------------------------------------

// ip.spec_ten_lite(x[n]); -----------------------------------------------------
//
// Lite spectral tendency algorithm: fixed responsiveness parameter; less
// accurate (especially at low frequencies).
//
// 1 inputs:
//    x[n].
//
// 1 outputs:
//    y[n], spectral centroid of x[n] in the range [0; 1].
//
spec_ten_lite(in) = (  _,
                       in : f2.xover1p1zraw : ro.cross(2) : balance : 
      (_ ,
       (in : rms) : m2.div) : f2.integrator <: *)
                    ~ _
      with {
           a0 = 100 / ma.SR;
           b1 = 1 - a0;
           rms(in) = in <: * : * (a0) :    +
                                           ~ * (b1) : sqrt;
           balance =   rms ,
                       rms : -;
      };
// -----------------------------------------------------------------------------

// ip.zc(x[n]); ----------------------------------------------------------------
//
// Zero-crossing (ZC) indicator function: it returns 1 if a ZC occurs, 
// 0 otherwise.
//
// 1 inputs:
//    x[n].
//
// 1 outputs:
//    y[n], ZC indication.
//
zc(x) = x * x' < 0;
// -----------------------------------------------------------------------------
                       
// ip.zcr(R[n], x[n]); ---------------------------------------------------------
//
// Zero-crossing rate (ZCR) using a one-pole lowpass filter for averaging.
// The ZCR correlates with the noisiness of signals, although it is
// effective only in specific cases, for example when comparing voiced and
// unvoiced sounds, or percussive and non-percussive ones. The ZCR also
// correlates with the spectral centroid of a signal. For sinusoidal
// signals, the ZCR output of this function can be mapped over Nyquist and
// can be used effectively as a frequency detector.
// References: Gouyon et al. 2000; Herrera-Boyer et al. 2006; 
// Peeters et al. 2011.
//
// 2 inputs:
//    R[n], analysis rate in Hz;
//    x[n].
//
// 1 outputs:
//    y[n], ZCR as an index in the range [0; 1].
//
zcr(window, x) = f2.lp1p(window, zc(x));
// -----------------------------------------------------------------------------

// ip.zcr4(R[n], x[n]); --------------------------------------------------------
//
// Zero-crossing rate with four cascaded one-pole lowpass filters for averaging.
//
// 2 inputs:
//    R[n], analysis rate in Hz;
//    x[n].
//
// 1 outputs:
//    y[n], ZCR as an index in the range [0; 1].
//
zcr4(window, x) = zc(x) : seq(i, 4, f2.lp1p(window));
// -----------------------------------------------------------------------------