Sine wave oscillator bank using "inverse / reverse" Goertzel algorithm

goertzelbank.nim

#nim.cfg
#-d:release
#-d:danger
#-d:blas=openblas 
#-d:lapack=openblas
#-d:samplerate=44100

# Goertzel oscillator bank:
# run time (seconds: 113, nanosecond: 104692500)
# duration : 4440 s samplerate : 44100  samples : 195804000
# 5.776424000531143e-07 seconds per sample @ 1000 oscillators


import std/[math]
import arraymancer

const samplerate {.intdefine.} = 44100

type

  OscBank* = object
    multiplier: Tensor[float]
    current*: Tensor[float]
    previous: Tensor[float]
    d0: Tensor[float]
    n: int


func initOscBank*(frequency: seq[float], phase: seq[float]): OscBank =
  {.noSideEffect.}: #?
    let 
      ts = Tau / samplerate.float
      frequency = frequency.toTensor()
      phaseincrement = ts * frequency
      phase = phase.toTensor() 
    OscBank(
      multiplier: (cos(phaseincrement)) * 2,
      current: sin(phase),
      previous: sin(phase - phaseIncrement),
      d0: zeros[float]([frequency.shape[0]]),
      n: 0
    )

func next*(oscbank: var OscBank) =
  {.noSideEffect.}: #?
    if oscbank.n > 0:
      oscbank.d0 = elwise_mul(oscbank.multiplier, oscbank.current) - oscbank.previous
      oscbank.previous = oscbank.current
      oscbank.current = oscbank.d0
    else:
      inc oscbank.n  


when isMainModule:
  import std/[monotimes]

  func lmap*(t, minin, maxin, minout, maxout: float): float {.inline.} =
    ((t - minin) / (maxin - minin)) * (maxout - minout) + minout


  let 
    duration = 4440 #seconds, is 74 min, is cd length.
    ticks = duration * samplerate
    nOsc = 1000

  var freq = newSeq[float](nOsc)
  for i in 0..<freq.len:
    freq[i] = 20.0 + i.float
  var phase = newSeq[float](nOsc)
  for i in 0..<freq.len:
    phase[i] = lmap(i.float, -PI, PI, 0.0, ticks.float)

  let t3 = getMonoTime()
  var oscbank = initOscBank(freq, phase)
  for i in 0..<ticks:
    oscbank.next
  let t4 = getMonoTime()
  let run = (ticks(t4) - ticks(t3)).float / 1000000000
  echo "\nGoertzel oscillator bank:"
  echo t4 - t3
  echo "duration : ", duration, "s samplerate : ", samplerate, " samples : ", ticks
  echo run / ticks.float , " seconds per sample @ ", nOsc , " oscillators\n"