Implementation of the sliding windowed infinite Fourier transform with functionality to cascade sDFTs serially as well as the analog-style spectrum analyzer implementation

analog-style-analyzer.js

/**
 * Single file implementation of analog-style spectrum analyzer using cascaded biquad bandpass filter bank (works best on truly-logarithmic frequency scale and constant-Q resolution)
 *
 * The frequency bands data is formatted like:
 * {lo: lowerBound,
 *  ctr: center,
 *  hi: higherBound}
 *
 * where lo and hi are used for calculating the necessary bandwidth for variable-Q transform spectrum visualizations and ctr for center frequency. This is generated using functions like generateFreqBands
 * Note: This one uses its own biquad filter routine for calculating biquad bandpass filter instead of WebAudio API provided BiquadFilterNode, so it allows visualization of data that isn't WebAudio and no OfflineAudioContext is needed
 */

class AnalogStyleAnalyzer {
  constructor(...args) {
    // initialize the analog-style analyzer coefficients
    this.calcCoeffs(args);
    this.spectrumData = [];
  }
  
  calcCoeffs(freqBands, order = 4, timeRes = Infinity, bandwidth = 1, sampleRate = 44100, prewarpQ = true) {
    this._coeffs = freqBands.map(x => {
      // biquad bandpass filter (cascaded biquad bandpass is not Butterworth nor Bessel, rather it is something called "critically-damped" since each filter stage shares the same every biquad coefficients)
      const K = Math.tan(Math.PI * x.ctr/sampleRate),
            bw = Math.abs(x.hi-x.lo) * bandwidth + 1/(timeRes/500),
            qCompensationFactor = prewarpQ ? (Math.PI * x.ctr/sampleRate)/K : 1, // an option to compensate Q values for frequency warping associated with bilinear transform
            Q = x.ctr/bw * qCompensationFactor,
            norm = 1 / (1 + K / Q + K * K),
            a0 = K / Q * norm,
            a1 = 0,
            a2 = -a0,
            b1 = 2 * (K * K - 1) * norm,
            b2 = (1 - K / Q + K * K) * norm,
            zs = [];
      for (let i = 0; i < order; i++) {
        zs[i] = {
          z1: 0,
          z2: 0,
          out: 0
        }
      }
      return {
        a0: a0,
        a1: a1,
        a2: a2,
        b1: b1,
        b2: b2,
        zs: zs
      };
    });
  }

  analyze(samples) {
    const newSpectrumData = new Array(this._coeffs.length).fill(0);
    for (const x of samples) {
      for (let i = 0; i < this._coeffs.length; i++) {
        for (let j = 0; j < this._coeffs[i].zs.length; j++) {
          const input = j <= 0 ? x : this._coeffs[i].zs[j-1].out;
          this._coeffs[i].zs[j].out = input * this._coeffs[i].a0 + this._coeffs[i].zs[j].z1;
          this._coeffs[i].zs[j].z1 = input * this._coeffs[i].a1 + this._coeffs[i].zs[j].z2 - this._coeffs[i].b1 * this._coeffs[i].zs[j].out;
          this._coeffs[i].zs[j].z2 = input * this._coeffs[i].a2 - this._coeffs[i].b2 * this._coeffs[i].zs[j].out;
        }
        newSpectrumData[i] = Math.max(newSpectrumData[i], Math.abs(this._coeffs[i].zs[this._coeffs[i].zs.length-1].out));
      }
    }
    this.spectrumData = newSpectrumData.map(x => x/2);
  }
}

iir-sliding-dft.js

/**
 * Single file implementation of sliding windowed infinite Fourier transform (SWIFT)
 *
 * The frequency bands data is formatted like:
 * {lo: lowerBound,
 *  ctr: center,
 *  hi: higherBound}
 *
 * where lo and hi are used for calculating the necessary bandwidth for variable-Q transform spectrum visualizations and ctr for center frequency. This is generated using functions like generateFreqBands
 */
class SWIFT {
  constructor(...args) {
    // initialize the sDFT coefficients
    this.calcCoeffs(args);
    this.spectrumData = [];
  }
  
  calcCoeffs(freqBands, order = 4, timeRes = 600, bandwidth = 1, sampleRate = 44100) {
    // calcCoeffs() can be called anywhere else to re-initialize sliding DFT after changes in frequency band distributions and note that x and y are used instead of real and imaginary since vector rotation is the equivalent of the complex one
    this._coeffs = [];
    freqBands.map((x, i) => {
      // rX and rY are calculated in advance here since calculating sin and cos functions are pretty slow af
      this._coeffs[i] = {
        rX: Math.cos(x.ctr*Math.PI/sampleRate*2),
        rY: Math.sin(x.ctr*Math.PI/sampleRate*2),
        decay: Math.E ** (-Math.abs(x.hi-x.lo) * 4 * bandwidth / sampleRate - 1/(timeRes*sampleRate/2000)),
        coeffs: []
      };
      for (let j = 0; j < order; j++) {
        this._coeffs[i].coeffs[j] = {
          x: 0,
          y: 0
        };
      }
    });
  }
  
  analyze(dataArray) {
    const newSpectrumData = new Array(this._coeffs.length).fill(0);
    for (const x of dataArray) {
      for (let i = 0; i < this._coeffs.length; i++) {
        for (let j = 0; j < this._coeffs[i].coeffs.length; j++) {
          const input = j <= 0 ? {
            x: x,
            y: 0,
          } : this._coeffs[i].coeffs[j-1],
                outX = (this._coeffs[i].coeffs[j].x * this._coeffs[i].rX - this._coeffs[i].coeffs[j].y * this._coeffs[i].rY) * this._coeffs[i].decay + input.x * (1-this._coeffs[i].decay),
                outY = (this._coeffs[i].coeffs[j].x * this._coeffs[i].rY + this._coeffs[i].coeffs[j].y * this._coeffs[i].rX) * this._coeffs[i].decay + input.y * (1-this._coeffs[i].decay);
          
          this._coeffs[i].coeffs[j].x = outX;
          this._coeffs[i].coeffs[j].y = outY;
        }
        newSpectrumData[i] = Math.max(newSpectrumData[i],
                                      this._coeffs[i].coeffs[this._coeffs[i].coeffs.length-1].x ** 2 +
                                      this._coeffs[i].coeffs[this._coeffs[i].coeffs.length-1].y ** 2);
      }
    }
    this.spectrumData = newSpectrumData.map((x) => Math.sqrt(x));
  }
}