analyzers.lib

//################################ analyzers.lib ##########################################
// Analyzers library. Its official prefix is `an`.
//########################################################################################

ma = library("maths.lib");
ba = library("basics.lib");
ro = library("routes.lib");
si = library("signals.lib");
fi = library("filters.lib");
an = library("analyzers.lib"); // for compatible copy/paste out of this file

declare name "Faust Analyzer Library";
declare version "0.1";

/************************************************************************
FAUST library file, jos section

Except where noted otherwise, The Faust functions below in this
section are Copyright (C) 2003-2017 by Julius O. Smith III <jos@ccrma.stanford.edu>
([jos](http://ccrma.stanford.edu/~jos/)), and released under the
(MIT-style) [STK-4.3](#stk-4.3-license) license.

All MarkDown comments in this section is Copyright 2016-2017 by Romain
Michon and Julius O. Smith III, and is released under the
[CCA4I](https://creativecommons.org/licenses/by/4.0/) license (TODO: if/when Romain agrees)

************************************************************************/

//==============================Amplitude Tracking========================================
//========================================================================================

//------------------`(an.)abs_envelope_rect`-----------------------------------
// Absolute value average with moving-average algorithm.
//
// #### Usage
//
// ```
// _ : abs_envelope_rect(period) : _;
// ```
//
//  Where:
//
// * `period`: sets the averaging frame in seconds
//-----------------------------------------------------------------------------
declare abs_envelope_rect author "Dario Sanfilippo and Julius O. Smith III";
declare abs_envelope_rect copyright "Copyright (C) 2020 Dario Sanfilippo 
      <sanfilippo.dario@gmail.com> and 
       2003-2020 by Julius O. Smith III <jos@ccrma.stanford.edu>";
declare abs_envelope_rect license "MIT-style STK-4.3 license";
abs_envelope_rect(period, x) = abs(x) : fi.avg_rect(period);


//------------------`(an.)abs_envelope_tau`------------------------------------
// Absolute value average with one-pole lowpass and tau response. 
// (See filters.lib.)
//
// #### Usage
//
// ```
// _ : abs_envelope_tau(period) : _;
// ```
//
//  Where:
//
// * `period`: (time to decay by 1/e) sets the averaging frame in secs
//-----------------------------------------------------------------------------
declare abs_envelope_tau author "Dario Sanfilippo and Julius O. Smith III";
declare abs_envelope_tau copyright "Copyright (C) 2020 Dario Sanfilippo 
      <sanfilippo.dario@gmail.com> and 
       2003-2020 by Julius O. Smith III <jos@ccrma.stanford.edu>";
declare abs_envelope_tau license "MIT-style STK-4.3 license";
abs_envelope_tau(period, x) = abs(x) : fi.avg_tau(period);


//------------------`(an.)abs_envelope_t60`------------------------------------
// Absolute value average with one-pole lowpass and t60 response.
// (See filters.lib.)
//
// #### Usage
//
// ```
// _ : abs_envelope_t60(period) : _;
// ```
//
//  Where:
//
// * `period`: (time to decay by 60 dB) sets the averaging frame in secs
//-----------------------------------------------------------------------------
declare abs_envelope_t60 author "Dario Sanfilippo and Julius O. Smith III";
declare abs_envelope_t60 copyright "Copyright (C) 2020 Dario Sanfilippo 
      <sanfilippo.dario@gmail.com> and 
       2003-2020 by Julius O. Smith III <jos@ccrma.stanford.edu>";
declare abs_envelope_t60 license "MIT-style STK-4.3 license";
abs_envelope_t60(period, x) = abs(x) : fi.avg_t60(period);


//------------------`(an.)abs_envelope_t19`------------------------------------
// Absolute value average with one-pole lowpass and t19 response. 
// (See filters.lib.)
//
// #### Usage
//
// ```
// _ : abs_envelope_t19(period) : _;
// ```
//
//  Where:
//
// * `period`: (time to decay by 1/e^2.2) sets the averaging frame in secs
//-----------------------------------------------------------------------------
declare abs_envelope_t19 author "Dario Sanfilippo and Julius O. Smith III";
declare abs_envelope_t19 copyright "Copyright (C) 2020 Dario Sanfilippo 
      <sanfilippo.dario@gmail.com> and 
       2003-2020 by Julius O. Smith III <jos@ccrma.stanford.edu>";
declare abs_envelope_t19 license "MIT-style STK-4.3 license";
abs_envelope_t19(period, x) = abs(x) : fi.avg_t19(period);


//---------------------------`(an.)amp_follower`---------------------------
// Classic analog audio envelope follower with infinitely fast rise and
// exponential decay.  The amplitude envelope instantaneously follows
// the absolute value going up, but then floats down exponentially.
// 
// `amp_follower` is a standard Faust function.
//
// #### Usage
//
// ```
// _ : amp_follower(rel) : _
// ```
//
// Where:
//
// * `rel`: release time = amplitude-envelope time-constant (sec) going down
//
// #### Reference
//
// * Musical Engineer's Handbook, Bernie Hutchins, Ithaca NY, 1975 Electronotes
// Newsletter, Bernie Hutchins
//------------------------------------------------------------
amp_follower(rel) = abs : env with {
 	p = ba.tau2pole(rel);
 	env(x) = x * (1.0 - p) : (+ : max(x,_)) ~ *(p);
};

peak_envelope = amp_follower; // Synonym for more standard naming


//---------------------------`(an.)amp_follower_ud`---------------------------
// Envelope follower with different up and down time-constants
// (also called a "peak detector").
//
// #### Usage
//
// ```
//    _ : amp_follower_ud(att,rel) : _
// ```
//
// Where:
//
// * `att`: attack time = amplitude-envelope time constant (sec) going up
// * `rel`: release time = amplitude-envelope time constant (sec) going down
//
// #### Note
//
// We assume rel >> att.  Otherwise, consider rel ~ max(rel,att).
// For audio, att is normally faster (smaller) than rel (e.g., 0.001 and 0.01).
// Use `amp_follower_ar` below to remove this restriction.
//
// #### Reference
//
// * "Digital Dynamic Range Compressor Design --- A Tutorial and Analysis", by
//   Dimitrios Giannoulis, Michael Massberg, and Joshua D. Reiss
//   <https://www.eecs.qmul.ac.uk/~josh/documents/2012/GiannoulisMassbergReiss-dynamicrangecompression-JAES2012.pdf>
//------------------------------------------------------------
amp_follower_ud(att,rel) = amp_follower(rel) : si.smooth(ba.tau2pole(att));


//---------------`(an.)amp_follower_ar`----------------
// Envelope follower with independent attack and release times. The
// release can be shorter than the attack (unlike in `amp_follower_ud`
// above).
//
// #### Usage
//
// ```
// _ : amp_follower_ar(att,rel) : _;
// ```
// * Author Jonatan Liljedahl, revised by RM
//---------------------------------------------------------
amp_follower_ar(att,rel) = abs : si.lag_ud(att,rel);


//------------------`(an.)ms_envelope_rect`------------------------------------
// Mean square with moving-average algorithm.
//
// #### Usage
//
// ```
// _ : ms_envelope_rect(period) : _;
// ```
//
//  Where:
//
// * `period`: sets the averaging frame in secs
//-----------------------------------------------------------------------------
declare ms_envelope_rect author "Dario Sanfilippo and Julius O. Smith III";
declare ms_envelope_rect copyright "Copyright (C) 2020 Dario Sanfilippo 
      <sanfilippo.dario@gmail.com> and 
       2003-2020 by Julius O. Smith III <jos@ccrma.stanford.edu>";
declare ms_envelope_rect license "MIT-style STK-4.3 license";
ms_envelope_rect(period, x) = x * x : fi.avg_rect(period);


//------------------`(an.)ms_envelope_tau`-------------------------------------
// Mean square average with one-pole lowpass and tau response. 
// (See filters.lib.)
//
// #### Usage
//
// ```
// _ : ms_envelope_tau(period) : _;
// ```
//
//  Where:
//
// * `period`: (time to decay by 1/e) sets the averaging frame in secs
//-----------------------------------------------------------------------------
declare ms_envelope_tau author "Dario Sanfilippo and Julius O. Smith III";
declare ms_envelope_tau copyright "Copyright (C) 2020 Dario Sanfilippo 
      <sanfilippo.dario@gmail.com> and 
       2003-2020 by Julius O. Smith III <jos@ccrma.stanford.edu>";
declare ms_envelope_tau license "MIT-style STK-4.3 license";
ms_envelope_tau(period, x) = x * x : fi.avg_tau(period);


//------------------`(an.)ms_envelope_t60`-------------------------------------
// Mean square with one-pole lowpass and t60 response. 
// (See filters.lib.)
//
// #### Usage
//
// ```
// _ : ms_envelope_t60(period) : _;
// ```
//
//  Where:
//
// * `period`: (time to decay by 60 dB) sets the averaging frame in secs
//-----------------------------------------------------------------------------
declare ms_envelope_t60 author "Dario Sanfilippo and Julius O. Smith III";
declare ms_envelope_t60 copyright "Copyright (C) 2020 Dario Sanfilippo 
      <sanfilippo.dario@gmail.com> and 
       2003-2020 by Julius O. Smith III <jos@ccrma.stanford.edu>";
declare ms_envelope_t60 license "MIT-style STK-4.3 license";
ms_envelope_t60(period, x) = x * x : fi.avg_t60(period);


//------------------`(an.)ms_envelope_t19`-------------------------------------
// Mean square with one-pole lowpass and t19 response. 
// (See filters.lib.)
//
// #### Usage
//
// ```
// _ : ms_envelope_t19(period) : _;
// ```
//
//  Where:
//
// * `period`: (time to decay by 1/e^2.2) sets the averaging frame in secs
//-----------------------------------------------------------------------------
declare ms_envelope_t19 author "Dario Sanfilippo and Julius O. Smith III";
declare ms_envelope_t19 copyright "Copyright (C) 2020 Dario Sanfilippo 
      <sanfilippo.dario@gmail.com> and 
       2003-2020 by Julius O. Smith III <jos@ccrma.stanford.edu>";
declare ms_envelope_t19 license "MIT-style STK-4.3 license";
ms_envelope_t19(period, x) = x * x : fi.avg_t19(period);


//------------------`(an.)rms_envelope_rect`-----------------------------------
// Root mean square with moving-average algorithm.
//
// #### Usage
//
// ```
// _ : rms_envelope_rect(period) : _;
// ```
//
//  Where:
//
// * `period`: sets the averaging frame in secs
//-----------------------------------------------------------------------------
declare rms_envelope_rect author "Dario Sanfilippo and Julius O. Smith III";
declare rms_envelope_rect copyright "Copyright (C) 2020 Dario Sanfilippo 
      <sanfilippo.dario@gmail.com> and 
       2003-2020 by Julius O. Smith III <jos@ccrma.stanford.edu>";
declare rms_envelope_rect license "MIT-style STK-4.3 license";
rms_envelope_rect(period, x) = ms_envelope_rect(period, x) : sqrt;


//------------------`(an.)rms_envelope_tau`------------------------------------
// Root mean square with one-pole lowpass and tau response.
// (See filters.lib.)
//
// #### Usage
//
// ```
// _ : rms_envelope_tau(period) : _;
// ```
//
//  Where:
//
// * `period`: (time to decay by 1/e) sets the averaging frame in secs
//-----------------------------------------------------------------------------
declare rms_envelope_tau author "Dario Sanfilippo and Julius O. Smith III";
declare rms_envelope_tau copyright "Copyright (C) 2020 Dario Sanfilippo 
      <sanfilippo.dario@gmail.com> and 
       2003-2020 by Julius O. Smith III <jos@ccrma.stanford.edu>";
declare rms_envelope_tau license "MIT-style STK-4.3 license";
rms_envelope_tau(period, x) = ms_envelope_tau(period, x) : sqrt;


//------------------`(an.)rms_envelope_t60`------------------------------------
// Root mean square with one-pole lowpass and t60 response.
// (See filters.lib.)
//
// #### Usage
//
// ```
// _ : rms_envelope_t60(period) : _;
// ```
//
//  Where:
//
// * `period`: (time to decay by 60 dB) sets the averaging frame in secs
//-----------------------------------------------------------------------------
declare rms_envelope_t60 author "Dario Sanfilippo and Julius O. Smith III";
declare rms_envelope_t60 copyright "Copyright (C) 2020 Dario Sanfilippo 
      <sanfilippo.dario@gmail.com> and 
       2003-2020 by Julius O. Smith III <jos@ccrma.stanford.edu>";
declare rms_envelope_t60 license "MIT-style STK-4.3 license";
rms_envelope_t60(period, x) = ms_envelope_t60(period, x) : sqrt;


//------------------`(an.)rms_envelope_t19`------------------------------------
// Root mean square with one-pole lowpass and t19 response.
// (See filters.lib.)
//
// #### Usage
//
// ```
// _ : rms_envelope_t19(period) : _;
// ```
//
//  Where:
//
// * `period`: (time to decay by 1/e^2.2) sets the averaging frame in secs
//-----------------------------------------------------------------------------
declare rms_envelope_t19 author "Dario Sanfilippo and Julius O. Smith III";
declare rms_envelope_t19 copyright "Copyright (C) 2020 Dario Sanfilippo 
      <sanfilippo.dario@gmail.com> and 
       2003-2020 by Julius O. Smith III <jos@ccrma.stanford.edu>";
declare rms_envelope_t19 license "MIT-style STK-4.3 license";
rms_envelope_t19(period, x) = ms_envelope_t19(period, x) : sqrt;


//-----------------------`(an.)zcr`--------------------------------------------
// Zero-crossing rate (ZCR) with one-pole lowpass averaging based on the tau 
// constant. It outputs an index between 0 and 1 at a desired analysis frame. 
// The ZCR of a signal correlates with the noisiness [Gouyon et al. 2000] and 
// the spectral centroid [Herrera-Boyer et al. 2006] of a signal. 
// For sinusoidal signals, the ZCR can be multiplied by ma.SR/2 and used
// as a frequency detector. For example, it can be deployed as a
// computationally efficient adaptive mechanism for automatic Larsen
// suppression.
//
// #### Usage
//
// ```
// _ : zcr(tau) : _;
// ```
//
// Where:
//
// * `tau`: (time to decay by e^-1) sets the averaging frame in seconds.
declare zcr author "Dario Sanfilippo";
declare zcr copyright "Copyright (C) 2020 Dario Sanfilippo
      <sanfilippo.dario@gmail.com>";
declare zcr license "MIT-style STK-4.3 license";
zcr(period, x) = ma.zc(x) : fi.lptau(period);


//=============================Spectrum-Analyzers=========================================
// Spectrum-analyzers split the input signal into a bank of parallel signals, one for
// each spectral band. They are related to the Mth-Octave Filter-Banks in `filters.lib`.
// The documentation of this library contains more details about the implementation.
// The parameters are:
//
// * `M`: number of band-slices per octave (>1)
// * `N`: total number of bands (>2)
// * `ftop` = upper bandlimit of the Mth-octave bands (<SR/2)
//
// In addition to the Mth-octave output signals, there is a highpass signal
// containing frequencies from ftop to SR/2, and a "dc band" lowpass signal
// containing frequencies from 0 (dc) up to the start of the Mth-octave bands.
// Thus, the N output signals are:
// ```
// highpass(ftop), MthOctaveBands(M,N-2,ftop), dcBand(ftop*2^(-M*(N-1)))
// ```
//
// A Spectrum-Analyzer is defined here as any band-split whose bands span
// the relevant spectrum, but whose band-signals do not
// necessarily sum to the original signal, either exactly or to within an
// allpass filtering. Spectrum analyzer outputs are normally at least nearly
// "power complementary", i.e., the power spectra of the individual bands
// sum to the original power spectrum (to within some negligible tolerance).
//
// #### Increasing Channel Isolation
//
// Go to higher filter orders - see Regalia et al. or Vaidyanathan (cited
// below) regarding the construction of more aggressive recursive
// filter-banks using elliptic or Chebyshev prototype filters.
//
// #### References
//
// * "Tree-structured complementary filter banks using all-pass sections",
//   Regalia et al., IEEE Trans. Circuits & Systems, CAS-34:1470-1484, Dec. 1987
// * "Multirate Systems and Filter Banks", P. Vaidyanathan, Prentice-Hall, 1993
// * Elementary filter theory: <https://ccrma.stanford.edu/~jos/filters/>
//========================================================================================


//-------------------------`(an.)mth_octave_analyzer`----------------------------
// Octave analyzer.
// `mth_octave_analyzer[N]` are standard Faust functions.
//
// #### Usage
// ```
// _ : mth_octave_analyzer(O,M,ftop,N) : par(i,N,_); // Oth-order Butterworth
// _ : mth_octave_analyzer6e(M,ftop,N) : par(i,N,_); // 6th-order elliptic
// ```
//
// Also for convenience:
//
// ```
// _ : mth_octave_analyzer3(M,ftop,N) : par(i,N,_); // 3d-order Butterworth
// _ : mth_octave_analyzer5(M,ftop,N) : par(i,N,_); // 5th-roder Butterworth
// mth_octave_analyzer_default = mth_octave_analyzer6e;
// ```
//
// Where:
//
// * `O`: order of filter used to split each frequency band into two
// * `M`: number of band-slices per octave
// * `ftop`: highest band-split crossover frequency (e.g., 20 kHz)
// * `N`: total number of bands (including dc and Nyquist)
//------------------------------------------------------------
mth_octave_analyzer6e(M,ftop,N) = _ <: bsplit(N-1) with {
  fc(n) = ftop * 2^(float(n-N+1)/float(M)); // -3dB crossover frequencies
  lp(n) = fi.lowpass6e(fc(n));  // 6th-order elliptic - see other choices above
  hp(n) = fi.highpass6e(fc(n)); //   (search for lowpass* and highpass*)
  bsplit(0)  = _;
  bsplit(i) = hp(i), (lp(i) <: bsplit(i-1));
};

// Butterworth analyzers may be cascaded with allpass
// delay-equalizers to make (allpass-complementary) filter banks:

mth_octave_analyzer(O,M,ftop,N) = _ <: bsplit(N-1) with {
  fc(n) = ftop * 2^(float(n-N+1)/float(M));
  lp(n) = fi.lowpass(O,fc(n)); // Order O Butterworth
  hp(n) = fi.highpass(O,fc(n));
  bsplit(0)  = _;
  bsplit(i) = hp(i), (lp(i) <: bsplit(i-1));
};

mth_octave_analyzer3(M,ftop,N) = mth_octave_analyzer(3,M,ftop,N);
mth_octave_analyzer5(M,ftop,N) = mth_octave_analyzer(5,M,ftop,N);
mth_octave_analyzer_default = mth_octave_analyzer6e; // default analyzer


//============================Mth-Octave Spectral Level===================================
// Spectral Level: Display (in bar graphs) the average signal level in each spectral band.
//========================================================================================


//------------------------`(an.)mth_octave_spectral_level6e`-------------------------
// Spectral level display.
//
// #### Usage:
//
// ```
// _ : mth_octave_spectral_level6e(M,ftop,NBands,tau,dB_offset) : _;
// ```
//
// Where:
//
// * `M`: bands per octave
// * `ftop`: lower edge frequency of top band
// * `NBands`: number of passbands (including highpass and dc bands),
// * `tau`: spectral display averaging-time (time constant) in seconds,
// * `dB_offset`: constant dB offset in all band level meters.
//
// Also for convenience:
//
// ```
// mth_octave_spectral_level_default = mth_octave_spectral_level6e;
// spectral_level = mth_octave_spectral_level(2,10000,20);
// ```
//------------------------------------------------------------
mth_octave_spectral_level6e(M,ftop,N,tau,dB_offset) = _<:
    _,mth_octave_analyzer6e(M,ftop,N) :
    _,(display:>_):attach with {
  display = par(i,N,dbmeter(i));
  dbmeter(i) = abs : si.smooth(ba.tau2pole(tau)) : max(ma.EPSILON) : ba.linear2db : +(dB_offset) :
     meter(N-i-1);
    meter(i) = speclevel_group(vbargraph("[%2i] [unit:dB]
     [tooltip: Spectral Band Level in dB]", -50, 10));
  O = int(((N-2)/M)+0.4999);
  speclevel_group(x)  = hgroup("[0] CONSTANT-Q SPECTRUM ANALYZER (6E), %N bands spanning
  	LP, %O octaves below %ftop Hz, HP
     [tooltip: See Faust's filters.lib for documentation and references]", x);
};

mth_octave_spectral_level_default = mth_octave_spectral_level6e;
spectral_level = mth_octave_spectral_level(2,10000,20);  // simple default


//---------------`(an.)[third|half]_octave_[analyzer|filterbank]`----------------
// A bunch of special cases based on the different analyzer functions described above:
//
// ```
// third_octave_analyzer(N) = mth_octave_analyzer_default(3,10000,N);
// third_octave_filterbank(N) = mth_octave_filterbank_default(3,10000,N);
// half_octave_analyzer(N) = mth_octave_analyzer_default(2,10000,N);
// half_octave_filterbank(N) = mth_octave_filterbank_default(2,10000,N);
// octave_filterbank(N) = mth_octave_filterbank_default(1,10000,N);
// octave_analyzer(N) = mth_octave_analyzer_default(1,10000,N);
// ```
//
// #### Usage
//
// See `mth_octave_spectral_level_demo` in `demos.lib`.
//------------------------------------------------------------
third_octave_analyzer(N) = mth_octave_analyzer_default(3,10000,N);
third_octave_filterbank(N) = mth_octave_filterbank_default(3,10000,N);
// Third-Octave Filter-Banks have been used in audio for over a century.
// See, e.g.,
//   Acoustics [the book], by L. L. Beranek
//   Amer. Inst. Physics for the Acoustical Soc. America,
//   <http://asa.aip.org/publications.html, 1986 (1st ed.1954)>

// Third-octave bands across the audio spectrum are too wide for current
// typical computer screens, so half-octave bands are the default:
half_octave_analyzer(N) = mth_octave_analyzer_default(2,10000,N);
half_octave_filterbank(N) = mth_octave_filterbank_default(2,10000,N);

octave_filterbank(N) = mth_octave_filterbank_default(1,10000,N);
octave_analyzer(N) = mth_octave_analyzer_default(1,10000,N);


//===============Arbritary-Crossover Filter-Banks and Spectrum Analyzers==================
// These are similar to the Mth-octave analyzers above, except that the
// band-split frequencies are passed explicitly as arguments.
//========================================================================================

// ACKNOWLEDGMENT
// Technique for processing a variable number of signal arguments due
// to Yann Orlarey (as is the entire Faust framework!)

//---------------`(an.)analyzer`--------------------------
// Analyzer.
//
// #### Usage
//
// ```
// _ : analyzer(O,freqs) : par(i,N,_); // No delay equalizer
// ```
//
// Where:
//
// * `O`: band-split filter order (ODD integer required for filterbank[i])
// * `freqs`: (fc1,fc2,...,fcNs) [in numerically ascending order], where
//           Ns=N-1 is the number of octave band-splits
//           (total number of bands N=Ns+1).
//
// If frequencies are listed explicitly as arguments, enclose them in parens:
//
// ```
// _ : analyzer(3,(fc1,fc2)) : _,_,_
// ```
//---------------------------------------------------
analyzer(O,lfreqs) = _ <: bsplit(nb) with
{
   nb = ba.count(lfreqs);
   fc(n) = ba.take(n, lfreqs);
   lp(n) = fi.lowpass(O,fc(n));
   hp(n) = fi.highpass(O,fc(n));
   bsplit(0) = _;
   bsplit(i) = hp(i), (lp(i) <: bsplit(i-1));
};

//================ Fast Fourier Transform (fft) and its Inverse (ifft) ===================
// Sliding FFTs that compute a rectangularly windowed FFT each sample.
//========================================================================================

//---------------`(an.)gortzelOpt` --------------------------
// Optimized Goertzel filter. 
//
// #### Usage
//
// ```
// _ : goertzelOpt(freq,N) : _;
// ```
//
// Where:
//
// * `freq`: frequency to be analyzed
// * `N`: the Goertzel block size
//
// #### Reference
//
// * <https://en.wikipedia.org/wiki/Goertzel_algorithm>
//---------------------------------------------------
goertzelOpt(freq,N,x) = mg 
with{
  mg = sqrt(eq^2 + eq'^2-eq*eq'*c) : ba.sAndH(reset0); // magnitude
  cnt = ba.time%N; // counter for windowing
  reset0 = cnt == (N-1); // reset when end of window
  reset1 = 1-(cnt == 0); // reset when beginning of window
  k = 0.5 + N*freq/ma.SR;
  w = (2*ma.PI/N)*k;
  c = 2*cos(w);
  eq = s // equation
  letrec{
    's = c*s*reset1 - s'*reset1*reset1' + x;
  };
};

//---------------`(an.)gortzelComp` --------------------------
// Complex Goertzel filter. 
//
// #### Usage
//
// ```
// _ : goertzelComp(freq,N) : _;
// ```
//
// Where:
//
// * `freq`: frequency to be analyzed
// * `N`: the Goertzel block size
//
// #### Reference
//
// * <https://en.wikipedia.org/wiki/Goertzel_algorithm>
//---------------------------------------------------
goertzelComp(freq,N,x) = mg 
with{
  mg = sqrt(real^2 + imag^2); // magnitude
  cnt = ba.time%N; // counter for windowing
  reset0 = cnt == (N-1); // reset when end of window
  reset1 = 1-(cnt == 0); // reset when beginning of window
  k = 0.5 + N*freq/ma.SR;
  w = (2*ma.PI/N)*k;
  sine = sin(w);
  cosine = cos(w);
  c = 2*cosine;
  eq = s
  letrec{
    's = c*s*reset1 - s'*reset1*reset1' + x;
  };
  real = eq - eq'*cosine;
  imag = eq'*sine;
};

//---------------`(an.)goertzel` --------------------------
// Same as [`goertzelOpt`](#goertzelopt). 
//
// #### Usage
//
// ```
// _ : goertzel(freq,N) : _;
// ```
//
// Where:
//
// * `freq`: frequency to be analyzed
// * `N`: the Goertzel block size
//
// #### Reference
//
// * <https://en.wikipedia.org/wiki/Goertzel_algorithm>
//---------------------------------------------------
goertzel = goertzelOpt;

// Undocumented utility functions used by fft and ifft:
c_magsq(N) = si.cbus(N) : par(i,N,(par(j,2,abs<:_*_):>_)) :> si.bus(N);
c_magdb(N) = si.cbus(N) : an.c_magsq(N) : par(i,N,(max(ma.EPSILON):log10:*(10.0)));
c_select_pos_freqs(2) = (_,_), (_,_); // both dc and SR/2 included with "positive frequencies"
c_select_pos_freqs(N) = si.cbus(N) : par(i,N/2+1,(_,_)),par(i,N/2-1,(!,!)) : si.cbus(N/2+1); // for complex spectra
  select_pos_freqs(2) = _,_; // both dc and SR/2 included
  select_pos_freqs(N) = si.bus(N) : par(i,N/2+1, _), par(i,N/2-1, !) : si. bus(N/2+1); // real power spectra etc.

rtorv(N,x) = par(i,N,x@i); // convert real scalar signal to length N real vector
rtocv(N,x) = par(i,N,(x@i,0)); // convert real scalar signal to length N complex vector with 0 imag part
rvtocv(N) = si.bus(N), par(i,N,0) : ro.interleave(N,2); // convert real N-vector to complex with 0 imag part

bit_reverse_selector(N,0) = 0;
bit_reverse_selector(N,i) = int(int(N)>>1)*(i&1) + bit_reverse_selector(int(N)>>1,(i>>1));

// decimation in time does this to the input:
bit_reverse_shuffle(N) = si.bus(N) <: par(i,N,bit_reverse_permuter(N,i)) with {
    bit_reverse_permuter(N,k) = ba.selector(bit_reverse_selector(N,k),N);
};

c_bit_reverse_shuffle(N) = si.cbus(N) <: par(i,N,c_bit_reverse_permuter(N,i)) with {
    c_bit_reverse_permuter(N,k) = ba.cselector(bit_reverse_selector(N,k),N);
};


//---------------`(an.)fft` --------------------------
// Fast Fourier Transform (FFT).
//
// #### Usage
//
// ```
// si.cbus(N) : fft(N) : si.cbus(N);
// ```
//
// Where:
//
// * `si.cbus(N)` is a bus of N complex signals, each specified by real and imaginary parts:
//   (r0,i0), (r1,i1), (r2,i2), ...
// * `N` is the FFT size (must be a power of 2: 2,4,8,16,...)
// * `fft(N)` performs a length `N` FFT for complex signals (radix 2)
// * The output is a bank of N complex signals containing the complex spectrum over time:
//   (R0, I0), (R1,I1), ...
//   - The dc component is (R0,I0), where I0=0 for real input signals.
//
// FFTs of Real Signals:
//
// * To perform a sliding FFT over a real input signal, you can say
// ```
// process = signal : an.rtocv(N) : an.fft(N);
// ```
// where `an.rtocv` converts a real (scalar) signal to a complex vector signal having a zero imaginary part.
//
//   * See `an.rfft_analyzer_c` (in `analyzers.lib`) and related functions for more detailed usage examples.
//
//   * Use `an.rfft_spectral_level(N,tau,dB_offset)` to display the power spectrum of a real signal.
//
//   * See `dm.fft_spectral_level_demo(N)` in `demos.lib` for an example GUI driving `an.rfft_spectral_level()`.
//
// #### Reference
//
// * [Decimation-in-time (DIT) Radix-2 FFT](https://cnx.org/contents/zmcmahhR@7/Decimation-in-time-DIT-Radix-2)
//
//---------------------------------------------------
fft(N) = si.cbus(N) : an.c_bit_reverse_shuffle(N) : an.fftb(N); // shuffle off to the butterflies:
fftb(1) = _,_; // each complex number is represented as (real,imag)
fftb(N) = si.cbus(N) : (fftb(No2) <: (si.cbus(No2), si.cbus(No2))), (fftb(No2)
	  <: (si.cbus(N):twiddleOdd(N))) :> si.cbus(N)
with {
  No2 = int(N)>>1;
  // Half of these multiplies can go away since w(k) = - w(k+N/2), and others as well (1,j,-j,-1,...)
  twiddleOdd(N) = par(k,N,si.cmul(cos(w(k)),0-sin(w(k))))
  with {
    w(k) = 2.0*ma.PI*float(k)/float(N);
  };
};

// `rfft`
// Slow to compile: rfft(N) = si.bus(N) : an.bit_reverse_shuffle(N) : an.rvtocv(N) : an.fftb(N);
// Order of magnitude faster to compile but takes a scalar input, so too different from fft:
// rfft(N) = an.rtocv(N) : an.fft(N);

//---------------`(an.)ifft`--------------------------
// Inverse Fast Fourier Transform (IFFT).
//
// #### Usage
//
// ```
// si.cbus(N) : ifft(N) : si.cbus(N);
// ```
//
// Where:
//
// * N is the IFFT size (power of 2)
// * Input is a complex spectrum represented as interleaved real and imaginary parts:
//   (R0, I0), (R1,I1), (R2,I2), ...
// * Output is a bank of N complex signals giving the complex signal in the time domain:
//   (r0, i0), (r1,i1), (r2,i2), ...
//---------------------------------------------------
ifft(N) = si.cbus(N) : an.c_bit_reverse_shuffle(N) : an.ifftb(N); // input is shuffled off to the butterflies:
ifftb(1) = _,_;
ifftb(N) = si.cbus(N) : (ifftb(No2) <: (si.cbus(No2), si.cbus(No2))), (ifftb(No2)
	<: (si.cbus(N):twiddleOddConj(N))) :> si.cbus(N) : par(i,2*N,/(2.0))
with {
  No2 = int(N)>>1;
  // Half of these multiplies can go away since w(k) = - w(k+N/2), and others as well (1,j,-j,-1,...)
  twiddleOddConj(N) = par(k,N,si.cmul(cos(w(k)),sin(w(k))))
  with {
    w(k) = 2.0*ma.PI*float(k)/float(N);
  };
};


// ========== FFT Analyzers ==========
rfft_analyzer_c(N) = an.rtocv(N) : an.fft(N) : an.c_select_pos_freqs(N); // complex spectral bins 0 to N/2
rfft_analyzer_db(N) = an.rfft_analyzer_c(N) : an.c_magdb(N/2+1); // assumes real input
rfft_analyzer_magsq(N) = an.rfft_analyzer_c(N) : an.c_magsq(N/2+1); // assumes real input

rfft_spectral_level(N,tau,dB_offset) = _<: _, an.rfft_analyzer_db(N) : _,(display:>_):attach with {
  display = par(i,N/2+1,dbmeter(i));
  dbmeter(i) = si.smooth(ba.tau2pole(tau)) : +(dB_offset) : meter(i);
    meter(i) = speclevel_group(vbargraph("[%2i] [unit:dB]
	       [tooltip: FFT Spectral Band Level in dB]", -50, 10));
  speclevel_group(x)  = hgroup("[0] FFT SPECTRUM ANALYZER, %N bands
			[tooltip: fft_spectral_level in Faust's analyzers.lib]", x);
};

// end jos section
/************************************************************************
************************************************************************
FAUST library file, GRAME section

Except where noted otherwise, Copyright (C) 2003-2017 by GRAME,
Centre National de Creation Musicale.
----------------------------------------------------------------------
GRAME LICENSE

This program is free software; you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License as
published by the Free Software Foundation; either version 2.1 of the
License, or (at your option) any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU Lesser General Public License for more details.

You should have received a copy of the GNU Lesser General Public
License along with the GNU C Library; if not, write to the Free
Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
02111-1307 USA.

EXCEPTION TO THE LGPL LICENSE : As a special exception, you may create a
larger FAUST program which directly or indirectly imports this library
file and still distribute the compiled code generated by the FAUST
compiler, or a modified version of this compiled code, under your own
copyright and license. This EXCEPTION TO THE LGPL LICENSE explicitly
grants you the right to freely choose the license for the resulting
compiled code. In particular the resulting compiled code has no obligation
to be LGPL or GPL. For example you are free to choose a commercial or
closed source license or any other license if you decide so.
************************************************************************
************************************************************************/

// TODO: Add GRAME functions here

//########################################################################################
/************************************************************************
FAUST library file, further contributions section

All contributions below should indicate both the contributor and terms
of license.  If no such indication is found, "git blame" will say who
last edited each line, and that person can be emailed to inquire about
license disposition, if their license choice is not already indicated
elsewhere among the libraries.  It is expected that all software will be
released under LGPL, STK-4.3, MIT, BSD, or a similar FOSS license.
************************************************************************/