A collection of audio analysis algorithms and some helper functions

audio-analyzer.js

 /** 
 * A collection of audio analysis algorithms and some helper functions
 */

// p5.js map function
function map(x, min, max, targetMin, targetMax) {
  return (x - min) / (max - min) * (targetMax - targetMin) + targetMin;
}

function clamp(x, min, max) {
  return Math.min(Math.max(x, min), max);
}

function idxWrapOver(x, length) {
  return (x % length + length) % length;
}

function getClippedIdx(x, length, mode = 'clip') {
  switch (mode) {
    case 'clip':
      return clamp(x, 0, length-1);
    case 'periodic':
      return idxWrapOver(x, length);
    case 'mirror':
      const period = 2 * (length - 1);
      const i = idxWrapOver(x, period);
      return (i > length - 1) ? period - i : i;
    default:
      return x;
  }
}

function smooth(array, amount = 2, parameter, filterType, idxClipMode) {
  const kernel = amount + 1;
  let filter = new Array(2 * kernel + 1).fill(0);
  filter = filter.map((_, i, arr) => {
    let x = i * 2 / (arr.length-1) - 1;
    return applyWindow(x, filterType, parameter);
  });
  let norm = filter.reduce((acc, val) => acc + val);
  return array.map((_, idx, arr) => {
    let segmentStart = idx - kernel;
    let sum = 0;
    for (let i = 0; i < filter.length; i++) {
      sum += (segmentStart+i >= 0 && segmentStart+i < arr.length) || idxClipMode !== 'zero' ? arr[getClippedIdx(segmentStart+i, arr.length, idxClipMode)]*filter[i]/norm : 0;
    }
    return sum;
  });
}

// dB convertion
function linearTodB(x) {
  return 20*Math.log10(x);
}

function dBToLinear(x) {
  return 10 ** (x/20);
}

// Linear interpolation (used for FFT bin interpolation)
function lerp(x, y, z) {
  return x*(1-z)+y*z;
}

function cubicInterp(x, y, z, w, i, tension = 0) {
  const tangentFactor = 1 - tension,
        m1 = tangentFactor * (z - x) / 2,
        m2 = tangentFactor * (w - y) / 2,
        squared = i ** 2,
        cubed = i ** 3;
  return (2 * cubed - 3 * squared + 1) * y + (cubed - 2 * squared + i) * m1 + (-2 * cubed + 3 * squared) * z + (cubed - squared) * m2;
}

function interp(arr, x, interpMode, interpParameter = 0, nthRoot = 1, logarithmic = false, clipMode = 'clip') {
  const intX = Math.trunc(x);
  switch (interpMode) {
    case 'nearest neighbor':
      return arr[getClippedIdx(Math.round(x), arr.length, clipMode)];
    case 'cubic':
      return absInvAscale(cubicInterp(absAscale(arr[getClippedIdx(intX-1, arr.length, clipMode)], nthRoot, logarithmic),
                                      absAscale(arr[getClippedIdx(intX, arr.length, clipMode)], nthRoot, logarithmic),
                                      absAscale(arr[getClippedIdx(intX+1, arr.length, clipMode)], nthRoot, logarithmic),
                                      absAscale(arr[getClippedIdx(intX+2, arr.length, clipMode)], nthRoot, logarithmic), x - intX, interpParameter), nthRoot, logarithmic);
    default:
      return absInvAscale(lerp((absAscale(arr[getClippedIdx(intX, arr.length, clipMode)], nthRoot, logarithmic)),
                  (absAscale(arr[getClippedIdx(intX+1, arr.length, clipMode)], nthRoot, logarithmic)), x - intX), nthRoot, logarithmic);
  }
}

function calcSmoothingTimeConstant(targetArr, sourceArr, factor = 0.5) {
  for (let i = 0; i < targetArr.length; i++) {
    targetArr[i] = (isNaN(targetArr[i]) ? 0 : targetArr[i])*(factor)+(isNaN(sourceArr[i]) ? 0 : sourceArr[i])*(1-factor);
  }
}

function calcFreqTilt(x, centerFreq = 440, amount = 3) {
  if (Math.abs(amount) > 0)
    return 10 ** (Math.log2(x/centerFreq)*amount/20);
  else
    return 1;
}

function applyWeight(x, weightAmount = 1, weightType = 'a') {
  const f2 = x ** 2;
  switch (weightType) {
    case 'a':
      return (1.2588966 * 148840000 * (f2 ** 2) /
      ((f2 + 424.36) * Math.sqrt((f2 + 11599.29) * (f2 + 544496.41)) * (f2 + 148840000))) ** weightAmount;
    case 'b':
      return (1.019764760044717 * 148840000 * (x ** 3) /
      ((f2 + 424.36) * Math.sqrt(f2 + 25122.25) * (f2 + 148840000))) ** weightAmount;
    case 'c':
      return (1.0069316688518042 * 148840000 * f2 /
      ((f2 + 424.36) * (f2 + 148840000))) ** weightAmount;
    case 'd':
      return ((x / 6.8966888496476e-5) * Math.sqrt(
               (
                 ((1037918.48 - f2)*(1037918.48 - f2) + 1080768.16*f2) /
                 ((9837328 - f2)*(9837328 - f2) + 11723776*f2)
               ) / ((f2 + 79919.29) * (f2 + 1345600))
             )) ** weightAmount;
    case 'm':
      const h1 = -4.737338981378384e-24*(f2 ** 3) + 2.043828333606125e-15*(f2 ** 2) - 1.363894795463638e-7*f2 + 1,
            h2 = 1.306612257412824e-19*(x ** 5) - 2.118150887518656e-11*(x ** 3) + 5.559488023498642e-4*x;

      return (8.128305161640991 * 1.246332637532143e-4 * x / Math.hypot(h1, h2)) ** weightAmount;
    default:
      return 1;
  }
}

/**
 * Apply window function
 * @param {number} x - The position of the window function from -1 to 1
 * @param {string} [windowType=Hann] - The specified window function to use
 * @param {number} [windowParameter=1] - The parameter of the window function (Adjustable window functions only)
 * @param {boolean} [truncate=true] - Zeroes out if x is more than 1 or less than -1
 * @param {number} [windowSkew=0] - Skew the window function to make symmetric windows asymmetric
 * @returns {number} The gain of window function
 */
function applyWindow(posX, windowType = 'Hann', windowParameter = 1, truncate = true, windowSkew = 0) {
  let x = windowSkew > 0 ? ((posX/2-0.5)/(1-(posX/2-0.5)*10*(windowSkew ** 2)))/(1/(1+10*(windowSkew ** 2)))*2+1 :
                           ((posX/2+0.5)/(1+(posX/2+0.5)*10*(windowSkew ** 2)))/(1/(1+10*(windowSkew ** 2)))*2-1;
  
  if (truncate && Math.abs(x) > 1)
    return 0;
  
  switch (windowType.toLowerCase()) {
    default:
      return 1;
    case 'hanning':
    case 'cosine squared':
    case 'hann':
      return Math.cos(x*Math.PI/2) ** 2;
    case 'raised cosine':
    case 'hamming':
      return 0.54 + 0.46 * Math.cos(x*Math.PI);
    case 'power of sine':
      return Math.cos(x*Math.PI/2) ** windowParameter;
    case 'tapered cosine':
    case 'tukey':
      return Math.abs(x) <= 1-windowParameter ? 1 : 
      (x > 0 ? 
       (-Math.sin((x-1)*Math.PI/windowParameter/2)) ** 2 :
       Math.sin((x+1)*Math.PI/windowParameter/2) ** 2);
    case 'blackman':
      return 0.42 + 0.5 * Math.cos(x*Math.PI) + 0.08 * Math.cos(x*Math.PI*2);
    case 'nuttall':
      return 0.355768 + 0.487396 * Math.cos(x*Math.PI) + 0.144232 * Math.cos(2*x*Math.PI) + 0.012604 * Math.cos(3*x*Math.PI);
    case 'flat top':
    case 'flattop':
      return 0.21557895 + 0.41663158 * Math.cos(x*Math.PI) + 0.277263158 * Math.cos(2*x*Math.PI) + 0.083578947 * Math.cos(3*x*Math.PI) + 0.006947368 * Math.cos(4*x*Math.PI);
    case 'kaiser':
      return Math.cosh(Math.sqrt(1-(x ** 2))*(windowParameter ** 2))/Math.cosh(windowParameter ** 2);
    case 'gauss':
    case 'gaussian':
      return Math.exp(-(windowParameter ** 2)*(x ** 2));
    case 'bartlett':
    case 'triangle':
    case 'triangular':
      return 1 - Math.abs(x);
    case 'poisson':
    case 'exponential':
      return Math.exp(-Math.abs(x * (windowParameter ** 2)));
    case 'hyperbolic secant':
    case 'sech':
      return 1/Math.cosh(x * (windowParameter ** 2));
    case 'quadratic spline':
      return Math.abs(x) <= 0.5 ? -((x*Math.sqrt(2)) ** 2)+1 : (Math.abs(x*Math.sqrt(2))-Math.sqrt(2)) ** 2;
    case 'parzen':
      return Math.abs(x) > 0.5 ? -2 * ((-1 + Math.abs(x)) ** 3) : 1 - 24 * (Math.abs(x/2) ** 2) + 48 * (Math.abs(x/2) ** 3);
    case 'welch':
      return 1 - (x ** 2);
  }
}

/**
 * Frequency scaling
 */

function fscale(x, freqScale = 'logarithmic', freqSkew = 0.5) {
  switch(freqScale.toLowerCase()) {
    default:
      return x;
    case 'log':
    case 'logarithmic':
      return Math.log2(x);
    case 'mel':
      return Math.log2(1+x/700);
    case 'critical bands':
    case 'bark':
      return (26.81*x)/(1960+x)-0.53;
    case 'equivalent rectangular bandwidth':
    case 'erb':
      return Math.log2(1+0.00437*x);
    case 'cam':
    case 'cams':
      return Math.log2((x/1000+0.312)/(x/1000+14.675));
    case 'sinh':
    case 'arcsinh':
    case 'asinh':
      return Math.asinh(x/(10 ** (freqSkew*4)));
    case 'shifted log':
    case 'shifted logarithmic':
      return Math.log2((10 ** (freqSkew*4))+x);
    case 'nth root':
      return x ** (1/(11-freqSkew*10));
    case 'negative exponential':
      return -(2 ** (-x/(2 ** (7+freqSkew*8))));
    case 'adjustable bark':
      return (26.81 * x)/((10 ** (freqSkew*4)) + x);
    case 'period':
      return 1/x;
  }
}

function invFscale(x, freqScale = 'logarithmic', freqSkew = 0.5) {
  switch(freqScale.toLowerCase()) {
    default:
      return x;
    case 'log':
    case 'logarithmic':
      return 2 ** x;
    case 'mel':
      return 700 * ((2 ** x) - 1);
    case 'critical bands':
    case 'bark':
      return 1960 / (26.81/(x+0.53)-1);
    case 'equivalent rectangular bandwidth':
    case 'erb':
      return (1/0.00437) * ((2 ** x) - 1);
    case 'cam':
    case 'cams':
      return (14.675 * (2 ** x) - 0.312)/(1-(2 ** x)) * 1000;
    case 'sinh':
    case 'arcsinh':
    case 'asinh':
      return Math.sinh(x)*(10 ** (freqSkew*4));
    case 'shifted log':
    case 'shifted logarithmic':
      return (2 ** x) - (10 ** (freqSkew*4));
    case 'nth root':
      return x ** ((11-freqSkew*10));
    case 'negative exponential':
      return -Math.log2(-x)*(2 ** (7+freqSkew*8));
    case 'adjustable bark':
      return (10 ** (freqSkew*4)) / (26.81 / x - 1);
    case 'period':
      return 1/x;
  }
}

// Amplitude scaling
function ascale(x, nthRoot = 1, logarithmic = false, dBRange = 70, useAbsoluteValue = true) {
  if (logarithmic)
    return map(20*Math.log10(x), -dBRange, 0, 0, 1);
  else
    return map(x ** (1/nthRoot), useAbsoluteValue ? 0 : dBToLinear(-dBRange) ** (1/nthRoot), 1, 0, 1);
}

function absAscale(x, nthRoot = 1, logarithmic = false) {
  if (logarithmic)
    return 20 * Math.log10(x);
  else 
    return x ** (1/nthRoot);
}

function absInvAscale(x, nthRoot = 1, logarithmic = false) {
  if (logarithmic)
    return 10 ** (x/20);
  else 
    return x ** (nthRoot);
}

// Hz and FFT bin conversion
function hertzToFFTBin(x, y = 'round', bufferSize = 4096, sampleRate = 44100) {
  const bin = x * bufferSize / sampleRate;
  let func = y;
  
  if (!['floor','ceil','trunc'].includes(func))
    func = 'round'; // always use round if you specify an invalid/undefined value
  
  return Math[func](bin);
}

function fftBinToHertz(x, bufferSize = 4096, sampleRate = 44100) {
  return x * sampleRate / bufferSize;
}

// Frequency bands generator
function generateFreqBands(N = 128, low = 20, high = 20000, freqScale, freqSkew, bandwidth = 0.5) {
  let freqArray = [];
  for (let i = 0; i < N; i++) {
    freqArray.push({
      lo: invFscale( map(i-bandwidth, 0, N-1, fscale(low, freqScale, freqSkew), fscale(high, freqScale, freqSkew)), freqScale, freqSkew),
      ctr: invFscale( map(i, 0, N-1, fscale(low, freqScale, freqSkew), fscale(high, freqScale, freqSkew)), freqScale, freqSkew),
      hi: invFscale( map(i+bandwidth, 0, N-1, fscale(low, freqScale, freqSkew), fscale(high, freqScale, freqSkew)), freqScale, freqSkew)
    });
  }
  return freqArray;
}

function generateOctaveBands(bandsPerOctave = 12, lowerNote = 4, higherNote = 124, detune = 0, bandwidth = 0.5) {
  const root24 = 2 ** ( 1 / 24 );
  const c0 = 440 * root24 ** -114; // ~16.35 Hz
  const groupNotes = 24/bandsPerOctave;
  let bands = [];
  for (let i = Math.round(lowerNote*2/groupNotes); i <= Math.round(higherNote*2/groupNotes); i++) {
    bands.push({
      lo: c0 * root24 ** ((i-bandwidth)*groupNotes+detune),
      ctr: c0 * root24 ** (i*groupNotes+detune),
      hi: c0 * root24 ** ((i+bandwidth)*groupNotes+detune)
    });
  }
  return bands;
}

// Calculate band frequencies, method derived from Avee Player
function generateAveePlayerFreqs(N, minFreq = 20, maxFreq = 20000, hzLinearFactor = 0, bandwidth = 0.5) {
  const freqBands = [];
  for (let i = 0; i < N; i++) {
    freqBands[i] = {
      lo: logSpace(minFreq, maxFreq, i-bandwidth, N-1, hzLinearFactor),
      ctr: logSpace(minFreq, maxFreq, i, N-1, hzLinearFactor),
      hi: logSpace(minFreq, maxFreq, i+bandwidth, N-1, hzLinearFactor)
    };
  }
  return freqBands;
}

function logSpace(x, y, z, w, l) {
  const centerFreq = x * ((y/x) ** (z/w));
  return centerFreq * (1-l) + (x + ((y - x) * z * (1/w))) * l;
}

/**
 * Generate frequency bands with piecewise lin-log scales
 * This function behaves similarly to KeyboardAudioVisualizer by DarthAffe
 * Dependent on both sample rate and DFT buffer length
 */

function generatePiecewiseLinLogBands(N = 64, low = 20, high = 15000, freqScale = 'log', freqSkew = 0.7, bufferSize = 2048, sampleRate = 44100) {
  const fromIdx = hertzToFFTBin(low, 'round', bufferSize, sampleRate),
        toIdx = hertzToFFTBin(high, 'round', bufferSize, sampleRate),
        usableSrcData = Math.max(N, (toIdx - fromIdx) + 1),
        fftBinsPerBand = usableSrcData / N,
        ratio = usableSrcData ** (1/N);
  
  let bands = [],
      fftBinCounter = 0,
      logCalculation = 1,
      idx = fromIdx;
  for (let i = 0; i < N; i++) {
    let count;
    switch (freqScale) {
      case 'linear':
        fftBinCounter += fftBinsPerBand;
        count = fftBinCounter;
        bands.push({
          lo: fftBinToHertz(idx, bufferSize, sampleRate),
          ctr: fftBinToHertz(idx+count/2, bufferSize, sampleRate),
          hi: fftBinToHertz(idx+count, bufferSize, sampleRate)
        });
        
        idx += count;
        fftBinCounter -= count;
        break;
      case 'nth root':
        const gamma = 11-freqSkew*10;
        count = Math.max(1, ((((i + 1)/N) ** gamma)*usableSrcData) - idx);
        bands.push({
          lo: fftBinToHertz(idx, bufferSize, sampleRate),
          ctr: fftBinToHertz(idx+count/2, bufferSize, sampleRate),
          hi: fftBinToHertz(idx+count, bufferSize, sampleRate)
        });
        idx += count;
        break;
      default:
        logCalculation *= ratio;
        count = Math.max(1, Math.trunc(logCalculation) - idx);
        bands.push({
          lo: fftBinToHertz(idx, bufferSize, sampleRate),
          ctr: fftBinToHertz(idx+count/2, bufferSize, sampleRate),
          hi: fftBinToHertz(idx+count, bufferSize, sampleRate)
        });
        idx += count;
    }
  }
  return bands;
}

/**
 * Pseudo-logarithmic frequency bands, derived from Classic Spectrum Analyzer for Winamp
 * This function generates pseudo-logarithmic bands that changes its frequency scale depending on frequency range, the number of bands, the specified FFT size and the sample rate
 */
function generatePseudoLogBands(N = 128, low = 0, high = 16000, freqSkew = 0, bufferSize = 2048, sampleRate = 44100) {
  // LogBarValueTable part
  let pnBarTable = new Array(N).fill(1),
      numUnassigned = bufferSize - N,
      numLastAssign = N - 1,
      bands = [];
  
  const fDiv = numUnassigned ** (1/N*(1-freqSkew)),
        nyquist = sampleRate / 2,
        delta = (high - low) / 2; // Delta between specified high frequencies and low frequencies, to naively allows setting the starting frequency from the first FFT bin
  
  // make the last bucket stop at the specified frequency
  if (nyquist > delta) {
    const maxFreqDiv = nyquist / (nyquist - delta);
    if (1 + numUnassigned - numUnassigned / fDiv < bufferSize / maxFreqDiv) {
      pnBarTable[N-1] = Math.trunc(bufferSize / maxFreqDiv);
      numUnassigned = bufferSize - N + 1;
      numUnassigned -= pnBarTable[N-1];
      numLastAssign -= 1;
    }
  }
  
  while (numUnassigned > 0) {
    let nAssign = Math.max(Math.trunc(numUnassigned - numUnassigned / fDiv + 0.5), 1);

    // currently O(N) complexity, which is slow on higher number of bands settings
    for (let nW = numLastAssign; nW >= 0 && numUnassigned > 0; nW--) {
      pnBarTable[nW] += nAssign;
      numUnassigned -= nAssign;
      nAssign = Math.max(Math.trunc(numUnassigned - numUnassigned / fDiv + 0.5), 1);
    }
  }
  
  // done with LogBarValueTable thing
  let currentFreq = hertzToFFTBin(low, 'round', bufferSize, sampleRate); // starts with approximate frequency from specified low frequency to display
  for (let i = 0; i < N; i++) {
    const nextFreq = currentFreq + pnBarTable[i];
    bands.push({
      lo: fftBinToHertz(currentFreq, bufferSize, sampleRate),
      ctr: fftBinToHertz(currentFreq + (nextFreq - currentFreq)/2, bufferSize, sampleRate),
      hi: fftBinToHertz(nextFreq, bufferSize, sampleRate)
    });
    currentFreq = nextFreq;
  }
  return bands;
}

// Calculate the FFT
function calcFFT(input) {
  let fft = input.map(x => x);
  let fft2 = input.map(x => x);
  transform(fft, fft2);
  let output = new Array(Math.round(fft.length/2)).fill(0);
  for (let i = 0; i < output.length; i++) {
    output[i] = Math.hypot(fft[i], fft2[i])/(fft.length);
  }
  return output;
}

function calcComplexFFT(input) {
  let fft = input.map(x => x);
  let fft2 = input.map(x => x);
  transform(fft, fft2);
  return input.map((_, i, arr) => {
    return {
      re: fft[i]/(arr.length/2),
      im: fft2[i]/(arr.length/2),
      magnitude: Math.hypot(fft[i], fft2[i])/(arr.length/2),
      phase: Math.atan2(fft2[i], fft[i])
    };
  });
}

function calcGoertzel(waveform, coeff) {
  let f1 = 0, f2 = 0, sine;
  for (const x of waveform) {
    sine = x + coeff * f1 - f2;
    f2 = f1;
    f1 = sine;
  }
  return Math.sqrt(f1 ** 2 + f2 ** 2 - coeff * f1 * f2) / waveform.length;
}

function calcBandpower(fftCoeffs, dataIdx, endIdx, average = true, useRMS = true, sum = false, interpolate = false, interpMode = 'linear', interpParam = 0, interpNthRoot = 1, interpLogScale = false) {
  let amp = 0;
  const low = clamp(dataIdx, 0, fftCoeffs.length-1),
        high = clamp(endIdx, 0, fftCoeffs.length-1),
        diff = endIdx - dataIdx + 1;
  if (interpolate) {
    const xFloat = low,
          xInt = Math.trunc(xFloat);
    if (xFloat - xInt <= 0)
      amp = fftCoeffs[dataIdx];
    else
      amp = interp(fftCoeffs, xFloat, interpMode, interpParam, interpNthRoot, interpLogScale);
  }
  else {
    for (let i = low; i <= high; i++) {
      if (average)
        amp += fftCoeffs[i] ** (useRMS ? 2 : 1);
      else
        amp = Math.max(amp, fftCoeffs[i]);
    }
    if (average)
      amp = useRMS ? Math.sqrt(amp/(sum ? 1 : diff)) : amp/(sum ? 1 : diff);
  }
  return amp;
}

/**
 * Calculate frequency bands power spectrum from FFT
 *
 * @param {number[]} fftCoeffs - The array of FFT magnitude coefficients
 * @param {string} [mode=interpolate] - FFT to frequency bands spectrum mode, setting to interpolate to do FFT bin interpolation, otherwise to insert zeroes, or bandpower to get average amount of amplitude between 2 frequency boundaries
 * @param {boolean} [average=true] - Average to sum FFT coefficients into bandpower coefficient
 * @param {FrequencyBands[]} [hzArray=generateFreqBands()] - The array of frequency bands in Hz
 * @returns {number[]} Output of frequency bands spectrum from linearly-spaced FFT coefficients
 */
function calcSpectrum(fftCoeffs, mode = 'interpolate', average = true, useRMS = true, sum = false, fftBinInterpolation, interpParam, interpNthRoot, interpLogScale, hzArray = generateFreqBands(), bufferSize = fftCoeffs.length*2, sampleRate = 44100) {
  let boundArr = [];
  let nBars = 0;
  for (let i = 0; i < hzArray.length; i++) {
    let minIdx = hertzToFFTBin(hzArray[i].lo, mode === 'bandpower' || mode === 'gaps' ? 'round' : 'ceil', bufferSize, sampleRate);
    let maxIdx = hertzToFFTBin(hzArray[i].hi, mode === 'bandpower' || mode === 'gaps' ? 'round' : 'floor', bufferSize, sampleRate);
    let minIdx2 = hertzToFFTBin(hzArray[i].lo, 'floor', bufferSize, sampleRate);
    if (mode === 'interpolate' && fftBinInterpolation !== 'zero insertion') {
      if (minIdx > maxIdx)
        nBars++;
      else {
        if (nBars > 1) {
          for (let j = 1; j <= nBars; j++)
            boundArr[boundArr.length - j].factor = (nBars - j) / nBars;
        }
        nBars = 1;
      }
    }
    boundArr.push({
      dataIdx: minIdx <= maxIdx || mode !== 'interpolate' || fftBinInterpolation === 'zero insertion' ? minIdx : minIdx2,
      endIdx: maxIdx - (mode === 'gaps'),
      factor: 0,
      interpolate: minIdx > maxIdx && mode === 'interpolate' && fftBinInterpolation !== 'zero insertion'
    });
  }
  return boundArr.map((x) => {
    return calcBandpower(fftCoeffs, x.dataIdx + x.factor, x.endIdx, average, useRMS, sum, x.interpolate, fftBinInterpolation, interpParam, interpNthRoot, interpLogScale);
  });
}

/**
 * Calculates spectrum from frequency bands by computing DTFT
 * The algorithm for calculating full spectrum using this have a O(N^2) complexity
 * @param {number[]} waveform - An array of audio data
 * @param {number} [octavesResolution=12] - Constant-Q transform frequency resolution
 * @param {string} [cqtMode=cbw] - Mode when comes to calculating per-bin bandwidth. Setting it to 'cbw' makes it constant-bandwidth, resembling Spectrum2 from Avee Player
 * @param (number) [cqtAlignment=1] - Temporal alignment of CQT/VQT bins, setting it to default value makes it causal and negative values makes it anti-causal
 * @param (number) [timeResolution=1000] - Desired time resolution at the center frequency specified by desiredTresAtFreq
 * @param (number) [desiredTresAtFreq=50] - The center frequency where VQT time resolution matches to timeResolution parameter in milliseconds
 * @param (string) frequencyScale - The specified frequency scale to use when comes to calculating VQT bandwidth
 * @param {FrequencyBands[]} [hzArray=generateFreqBands()] - The array of center frequency bands
 * @returns {number[]} Spectrum2 output, either DTFT or CQT
 */
function calcSpectrum2(waveform, octavesResolution = 12, cqtMode = 'cbw', cqtAlignment = 1, timeResolution = 1000, desiredTresAtFreq = 50, freqResolutionScale, freqResolutionLinearFactor, windowFunction, windowParameter, windowSkew, hzArray = generateFreqBands(), sampleRate = 44100) {
  return hzArray.map(x => {
    const coeff = 2 * Math.cos(2*Math.PI*x.ctr/sampleRate);
    if (cqtMode === 'cbw') {
      return calcGoertzel(waveform, coeff);
    }
    else {
      const vqtBandwidth = 1/(2*timeResolution/1000),
            slo = fscale(desiredTresAtFreq - vqtBandwidth, freqResolutionScale, freqResolutionLinearFactor),
            shi = fscale(desiredTresAtFreq + vqtBandwidth, freqResolutionScale, freqResolutionLinearFactor),
            lowerBound = invFscale(map(map(fscale(x.ctr, freqResolutionScale, freqResolutionLinearFactor), slo, shi, 0, 1)-0.5, 0, 1, slo, shi), freqResolutionScale, freqResolutionLinearFactor),
            higherBound = invFscale(map(map(fscale(x.ctr, freqResolutionScale, freqResolutionLinearFactor), slo, shi, 0, 1)+0.5, 0, 1, slo, shi), freqResolutionScale, freqResolutionLinearFactor);
      let f1 = 0,
          f2 = 0,
          sine,
          tlen = cqtMode === 'vqt' ? 1/(higherBound - lowerBound) : 2 / (x.ctr / octavesResolution + (sampleRate/waveform.length) * 2 * (cqtMode === 'pcqt')),
          scaledTlen = Math.min(tlen*sampleRate/waveform.length, 1),
          //lowerIdx = Math.ceil(map(-scaledTlen, -1, 1, 0, waveform.length-1)),
          //higherIdx = Math.floor(map(scaledTlen, -1, 1, 0, waveform.length-1)),
	  offset = Math.trunc((waveform.length-scaledTlen*waveform.length)*(0.5+cqtAlignment/2)),
          lowerIdx = offset,
          higherIdx = Math.trunc(scaledTlen*waveform.length)+offset-1,
          norm = 0;
      for (let i = lowerIdx; i <= higherIdx; i++) {
        let posX = map(i, lowerIdx, higherIdx, -1, 1);//(i * 2 / (waveform.length-1) - 1) / scaledTlen;
        let w = applyWindow(posX, windowFunction, windowParameter, true, windowSkew);
        norm += w;
        
        // Goertzel transform
        sine = waveform[i]*w + coeff * f1 - f2;
        f2 = f1;
        f1 = sine;
      }
      return Math.sqrt(f1 ** 2 + f2 ** 2 - coeff * f1 * f2) / norm;
    }
  });
}

/**
 * Brown-Puckette constant-Q transform algorithm, frequency domain filterbank to turn FFT into CQT
 * 
 * @param {ComplexFFT[]} complexFFTCoeffs - The complex FFT coefficients to process
 * @param {number} [octaveSpacing=12] - Fractional octave, higher values provides frequency higher resolution
 * @param {boolean} [constantBandwidth=false] - If set to true, the algorithm will be constant-bandwidth instead
 * @param {FrequencyBands[]} [hzArray=generateFreqBands()] - Array of center frequency bands
 * @returns {number[]} Brown-Puckette's CQT output
 */
function calcCQT(complexFFTCoeffs, octaveSpacing = 12, cqtMode = 'cbw', bandwidth = 4, timeResolution = 1000, desiredTresAtFreq = 50, freqResolutionScale, freqResolutionLinearFactor, windowFunction, windowParameter, windowSkew, hzArray = generateFreqBands(), bufferSize = complexFFTCoeffs.length, sampleRate = 44100) {
  return hzArray.map((band) => {
    const freq = band.ctr;
    let result = 0,
        l1 = 0,
        l2 = 0,
        r1 = 0,
        r2 = 0,
        a1 = 0,
        a2 = 0,
        b1 = 0,
        b2 = 0,
        real, imag;
    const vqtBandwidth = 1/(2*timeResolution/1000),
          slo = fscale(desiredTresAtFreq - vqtBandwidth, freqResolutionScale, freqResolutionLinearFactor),
          shi = fscale(desiredTresAtFreq + vqtBandwidth, freqResolutionScale, freqResolutionLinearFactor),
          lowerBound = invFscale(map(map(fscale(freq, freqResolutionScale, freqResolutionLinearFactor), slo, shi, 0, 1)-0.5, 0, 1, slo, shi), freqResolutionScale, freqResolutionLinearFactor),
          higherBound = invFscale(map(map(fscale(freq, freqResolutionScale, freqResolutionLinearFactor), slo, shi, 0, 1)+0.5, 0, 1, slo, shi), freqResolutionScale, freqResolutionLinearFactor);
    
    let tlen = Math.min(cqtMode === 'cbw' ? Infinity : (cqtMode === 'vqt' ? 1 / (higherBound - lowerBound) : 2/(freq/octaveSpacing + (sampleRate/bufferSize) * 2 * (cqtMode === 'pcqt'))), bufferSize / sampleRate);
    let flen = Math.min(bandwidth * bufferSize / (tlen * sampleRate), bufferSize);
    let center = freq * bufferSize / sampleRate;
    let start = Math.ceil(center - flen/2);
    let end = Math.floor(center + flen/2);
    /*let len = end - start + 1;
    
    let coeffs = new Array(len).fill(0);
    for (let i = start; i <= end; i++) {
      let sign = (i & 1) ? (-1) : (1);
      let x = 2 * (i - center) / flen;
      let w = applyWindow(x, windowFunction, windowParameter, true, windowSkew);
      coeffs[i - start] = w*sign;
    }
    for (let i = 0; i < coeffs.length; i++) {
      let u = coeffs[i];
      let x = start + i;
      let y = bufferSize - x;
      a1 += u * complexFFTCoeffs[(x % complexFFTCoeffs.length + complexFFTCoeffs.length) % complexFFTCoeffs.length].re;
      a2 += u * complexFFTCoeffs[(x % complexFFTCoeffs.length + complexFFTCoeffs.length) % complexFFTCoeffs.length].im;
      b1 += u * complexFFTCoeffs[(y % complexFFTCoeffs.length + complexFFTCoeffs.length) % complexFFTCoeffs.length].re;
      b2 += u * complexFFTCoeffs[(y % complexFFTCoeffs.length + complexFFTCoeffs.length) % complexFFTCoeffs.length].im;
    }*/
    for (let i = start; i <= end; i++) {
      let sign = (i & 1) ? (-1) : (1);
      let x = 2 * (i - center) / flen;
      let w = applyWindow(x, windowFunction, windowParameter, true, windowSkew);
      let u = w*sign;
      let y = bufferSize - i;
      a1 += u * complexFFTCoeffs[(i % complexFFTCoeffs.length + complexFFTCoeffs.length) % complexFFTCoeffs.length].re;
      a2 += u * complexFFTCoeffs[(i % complexFFTCoeffs.length + complexFFTCoeffs.length) % complexFFTCoeffs.length].im;
      b1 += u * complexFFTCoeffs[(y % complexFFTCoeffs.length + complexFFTCoeffs.length) % complexFFTCoeffs.length].re;
      b2 += u * complexFFTCoeffs[(y % complexFFTCoeffs.length + complexFFTCoeffs.length) % complexFFTCoeffs.length].im;
    }
    l1 = a1 + b1;
    l2 = a2 - b2;
    r1 = b2 + a2;
    r2 = b1 - a1;
    real = Math.hypot(l1, l2);
    imag = Math.hypot(r1, r2);
    return Math.hypot(real, imag)/8;
  });
}

function calcAutocorrelation(waveform) {
  return waveform.map((_, lag, arr) => {
    let sum = 0;
    for (let i = 0; i < arr.length-lag; i++) {
      let lagged = i+lag;
      let x = waveform[i];
      let y = waveform[lagged];
      let product = x * y;
      sum += product;
    }
    
    return sum / arr.length;
  });
}

/**
 * FFT and convolution (JavaScript)
 * 
 * Copyright (c) 2017 Project Nayuki. (MIT License)
 * https://www.nayuki.io/page/free-small-fft-in-multiple-languages
 */

/* 
 * Computes the discrete Fourier transform (DFT) of the given complex vector, storing the result back into the vector.
 * The vector can have any length. This is a wrapper function.
 */
function transform(real, imag) {
	const n = real.length;
	if (n != imag.length)
		throw "Mismatched lengths";
	if (n <= 0)
		return;
	else if ((2 ** Math.trunc(Math.log2(n))) === n)  // Is power of 2
		transformRadix2(real, imag);
	else  // More complicated algorithm for arbitrary sizes
		transformBluestein(real, imag);
}


/* 
 * Computes the inverse discrete Fourier transform (IDFT) of the given complex vector, storing the result back into the vector.
 * The vector can have any length. This is a wrapper function. This transform does not perform scaling, so the inverse is not a true inverse.
 */
function inverseTransform(real, imag) {
	transform(imag, real);
}


/* 
 * Computes the discrete Fourier transform (DFT) of the given complex vector, storing the result back into the vector.
 * The vector's length must be a power of 2. Uses the Cooley-Tukey decimation-in-time radix-2 algorithm.
 */
function transformRadix2(real, imag) {
	// Length variables
	const n = real.length;
	if (n != imag.length)
		throw "Mismatched lengths";
	if (n <= 1)  // Trivial transform
		return;
	const logN = Math.log2(n);
	if ((2 ** Math.trunc(logN)) !== n)
		throw "Length is not a power of 2";
	
	// Trigonometric tables
	let cosTable = new Array(n / 2);
	let sinTable = new Array(n / 2);
	for (let i = 0; i < n / 2; i++) {
		cosTable[i] = Math.cos(2 * Math.PI * i / n);
		sinTable[i] = Math.sin(2 * Math.PI * i / n);
	}
	
	// Bit-reversed addressing permutation
	for (let i = 0; i < n; i++) {
		let j = reverseBits(i, logN);
		if (j > i) {
			let temp = real[i];
			real[i] = real[j];
			real[j] = temp;
			temp = imag[i];
			imag[i] = imag[j];
			imag[j] = temp;
		}
	}
	
	// Cooley-Tukey decimation-in-time radix-2 FFT
	for (let size = 2; size <= n; size *= 2) {
		let halfsize = size / 2;
		let tablestep = n / size;
		for (let i = 0; i < n; i += size) {
			for (let j = i, k = 0; j < i + halfsize; j++, k += tablestep) {
				const l = j + halfsize;
				const tpre =  real[l] * cosTable[k] + imag[l] * sinTable[k];
				const tpim = -real[l] * sinTable[k] + imag[l] * cosTable[k];
				real[l] = real[j] - tpre;
				imag[l] = imag[j] - tpim;
				real[j] += tpre;
				imag[j] += tpim;
			}
		}
	}
	
	// Returns the integer whose value is the reverse of the lowest 'bits' bits of the integer 'x'.
	function reverseBits(x, bits) {
		let y = 0;
		for (let i = 0; i < bits; i++) {
			y = (y << 1) | (x & 1);
			x >>>= 1;
		}
		return y;
	}
}


/* 
 * Computes the discrete Fourier transform (DFT) of the given complex vector, storing the result back into the vector.
 * The vector can have any length. This requires the convolution function, which in turn requires the radix-2 FFT function.
 * Uses Bluestein's chirp z-transform algorithm.
 */
function transformBluestein(real, imag) {
	// Find a power-of-2 convolution length m such that m >= n * 2 + 1
	const n = real.length;
	if (n != imag.length)
		throw "Mismatched lengths";
	const m = 2 ** Math.trunc(Math.log2(n*2)+1);
	
	// Trignometric tables
	let cosTable = new Array(n);
	let sinTable = new Array(n);
	for (let i = 0; i < n; i++) {
		let j = i * i % (n * 2);  // This is more accurate than j = i * i
		cosTable[i] = Math.cos(Math.PI * j / n);
		sinTable[i] = Math.sin(Math.PI * j / n);
	}
	
	// Temporary vectors and preprocessing
	let areal = newArrayOfZeros(m);
	let aimag = newArrayOfZeros(m);
	for (let i = 0; i < n; i++) {
		areal[i] =  real[i] * cosTable[i] + imag[i] * sinTable[i];
		aimag[i] = -real[i] * sinTable[i] + imag[i] * cosTable[i];
	}
	let breal = newArrayOfZeros(m);
	let bimag = newArrayOfZeros(m);
	breal[0] = cosTable[0];
	bimag[0] = sinTable[0];
	for (let i = 1; i < n; i++) {
		breal[i] = breal[m - i] = cosTable[i];
		bimag[i] = bimag[m - i] = sinTable[i];
	}
	
	// Convolution
	let creal = new Array(m);
	let cimag = new Array(m);
	convolveComplex(areal, aimag, breal, bimag, creal, cimag);
	
	// Postprocessing
	for (let i = 0; i < n; i++) {
		real[i] =  creal[i] * cosTable[i] + cimag[i] * sinTable[i];
		imag[i] = -creal[i] * sinTable[i] + cimag[i] * cosTable[i];
	}
}


/* 
 * Computes the circular convolution of the given real vectors. Each vector's length must be the same.
 */
function convolveReal(x, y, out) {
	const n = x.length;
	if (n != y.length || n != out.length)
		throw "Mismatched lengths";
	convolveComplex(x, newArrayOfZeros(n), y, newArrayOfZeros(n), out, newArrayOfZeros(n));
}


/* 
 * Computes the circular convolution of the given complex vectors. Each vector's length must be the same.
 */
function convolveComplex(xreal, ximag, yreal, yimag, outreal, outimag) {
	const n = xreal.length;
	if (n != ximag.length || n != yreal.length || n != yimag.length
			|| n != outreal.length || n != outimag.length)
		throw "Mismatched lengths";
	
	xreal = xreal.slice();
	ximag = ximag.slice();
	yreal = yreal.slice();
	yimag = yimag.slice();
	transform(xreal, ximag);
	transform(yreal, yimag);
	
	for (let i = 0; i < n; i++) {
		const temp = xreal[i] * yreal[i] - ximag[i] * yimag[i];
		ximag[i] = ximag[i] * yreal[i] + xreal[i] * yimag[i];
		xreal[i] = temp;
	}
	inverseTransform(xreal, ximag);
	
	for (let i = 0; i < n; i++) {  // Scaling (because this FFT implementation omits it)
		outreal[i] = xreal[i] / n;
		outimag[i] = ximag[i] / n;
	}
}


function newArrayOfZeros(n) {
	let result = new Array(n).fill(0);
	return result;
}