Demonstartion of a Goertzel Algorithm in Swift

Goertzel.swift

import Foundation


/*
  Demo implementation of the famous Goertzel algorithm in Swift
  NB it is not currently possible to make this generic without
  also pulling in the Swift Numerics package.
 
  Goertzel is well described here : https://en.wikipedia.org/wiki/Goertzel_algorithm

  *in theory* this is less cumputationally expensive than running a DFT,
  in practice however, you are going to run this code in a hot loop and DFT
  operations can be vectorised, while this cannot, so YMMV

 */

public struct Goertzel {
  
  public struct Constants {
    
    let coeff  : Float
    let sine   : Float
    let cosine : Float
    let block  : Int
    
    public init ( block: Int, sampleRate: Int ) {
      
      let N = Float(block)
      let S = Float(sampleRate)
      
      let k = Float( Int( 0.5 + ((N * 1200) / S) ) )
      let w = ((2 * Float.pi) / N) * k // yeah yeah, omega, w/e
      
      cosine = cos(w)
      sine   = sin(w)
      coeff  = 2 * cosine
      
      self.block = block
    }
    
    
    
    public func compute( samples:[Float], using consts: Goertzel.Constants ) -> Float {
      
      assert(samples.count == consts.block)
      
      var q0 : Float = 0
      var q1 : Float = 0
      var q2 : Float = 0
      
      // classic algo that you'll find all over the place
      
      for sample in samples {                // this avoids pre filling two values, and is
        q0 = consts.coeff * q1 - q2 + sample // equivalent to coeff * samples[i-2] + samples[i-1]
        q2 = q1                              // so that's a loop variant and we can't vectorise it.
        q1 = q0                              // bugger.
      }
      
      // should probably be a scaling factor in here as well.
      // likely we should divide these by samples.count / 2
      // verify empirically that this matters though
      // because those / flops are expensive in a hot loop

      
      let real = (q1 - q2 * consts.cosine)
      let imag = (q2 * consts.sine)
      
      // compute absolute magnitude of complex
      return sqrt ( (real * real) + (imag * imag) )
    }
    
  }
  
  
}