Compare sine wave oscillators for speed and accuracy. Goertzel, Rotating Vector vs. sin()

goertzel_bench.nim

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

# Result (12th Gen Intel(R) Core(TM) i9-12900   2.40 GHz):
#
# GSinOsc  : (seconds: 0, nanosecond: 516204700)
# RVSinOsc : (seconds: 1, nanosecond: 672888300)
# STDSinOsc: (seconds: 5, nanosecond: 908334300)
# MSerror Goertzel vs sin()    : 2.108412096000499e-15
# MSerror Rotating vec vs sin(): 0.001964335532599123


import std/[complex, math]

const samplerate {.intdefine.} = 44100

type
  GSinOsc* = object  #inverse / reverse Goertzel
    multiplier: float
    current*: float
    previous: float
    n: int

  RvSinOsc = object   #rotating vector
    multiplier: Complex64
    value*: Complex64


func initGsinOsc(frequency: float, phase: float = 0.0): GSinOsc =
  let 
    phaseIncrement = frequency * Tau / samplerate.float
  GSinosc(
    multiplier: 2 * cos(phaseIncrement),
    current: sin(phase),
    previous: sin(phase - phaseIncrement),
    n: 0
  )

func next(osc: var GSinOsc){.inline.} =
  if osc.n > 0:
    let d0 = osc.multiplier * osc.current - osc.previous 
    osc.previous = osc.current
    osc.current = d0
  else:
    inc osc.n



func initRvSinOsc(frequency:float, samplerate:int): RvSinOsc =
  let phaseIncrement = frequency * Tau / samplerate.float
  RvSinOsc(
    multiplier: complex64(cos(phaseIncrement), sin(phaseIncrement)),
    value: complex64(1.0, 0.0)
  )

func next(osc: var RvSinOsc){.inline.} =
  osc.value = osc.value * osc.multiplier



when isMainModule:
  import std/[monotimes]

  func mse(s1, s2: seq[float]):float =
    var mseSum = 0.0
    for i in 0..<s1.len:
      let diff = s1[i] - s2[i]
      mseSum += float(diff * diff)
    mseSum / float(s1.len)


  let t = 4440 #seconds, is 74 min, is cd length.
  let ticks = t * samplerate
  let frequency = 440.0

  var 
    gsinOsc = initGsinOsc(frequency)
    gbuff = newSeq[float](ticks)
  let t0 = getMonoTime()
  for i in 0..<ticks:
    gsinOsc.next()
    gbuff[i] = gsinOsc.current
  echo "GSinOsc  : ", getMonoTime() - t0

  var 
    rvsinOsc = initRvSinOsc(frequency, samplerate)
    rvbuff = newSeq[float](ticks)
  let t1 = getMonoTime()
  for i in 0..<ticks:
    rvsinOsc.next()
    rvbuff[i] = rvsinOsc.value.im 
  echo "RVSinOsc : ", getMonoTime() - t1

  var 
    stdbuff = newSeq[float](ticks)
  let
    phaseIncrement = frequency * Tau / samplerate.float
    t2 = getMonoTime()
  for tick in 0..<ticks:
    stdbuff[tick] = sin(tick.float * phaseIncrement)
  echo "STDSinOsc: ", getMonoTime() - t2

  echo "MSerror Goertzel vs sin()    : ", mse(stdbuff, gbuff)
  echo "MSerror Rotating vec vs sin(): ", mse(stdbuff, rvbuff)