goertzel.jl

This unit implements three Goertzel's algorithms for estimating DFT complex coefficients at a specific frequency.

The result of Goertzel's algorithm for a given frequency is practically identical to the output of an FFT for that frequency, but requieres only $t$ operations, where $t$ is the window length in samples. Since the complexity of FFT algirthms is $t*log(t)$, Goertzels' algorithm are advantageous when only one or a few coefficients are needed. Furthermore, the estimation of Goertzel's algorithms is not restricted to a Fourier discrete frequency, but can be in the whole positive real line up to the Nyquist frequency. Smart block on-line implementations are also possible.

References

Goertzel G (1958) An algorithm for the evaluation of finite trigonometric series. The American Mathematical Monthly, 65(1), 34-35.

See also here.

FourierAnalysis.goertzelFunction
function goertzel( x   :: Vector{T},
                   f   :: IntOrReal,
                   sr  :: Int,
                   wl  :: Int,
               startAT :: Int = 1) where T<:Union{Real, Complex}

Given a time series as input vector x sampled at sr sampling rate, return the DFT complex coefficient at the discrete Fourier frequency which is the closest to f. The coefficient is computed for the time window starting at startAt (one-based) and of lenght wl.

If startAT=1 (default) the first wl samples of the vector are considered.

wl does not need to be power of 2.

See also: goertzel_fast, goertzel2.

Examples:

using FourierAnalysis
sr, t, f, a = 128, 128, 5, 10
v=sinusoidal(a, f, sr, t, 0)
c=goertzel(v, f, sr, t) # should output 0+aim
source
FourierAnalysis.goertzel_fastFunction
function goertzel_fast( x  :: Vector{T},
                        wl :: Int,
                        a  :: Real,
                        c  :: Real,
                        s  :: Real,
                        d  :: Real,
                    startAT:: Int = 1) where T<:Union{Real, Complex}

Fast version of the goertzel function to be preferred if the function is invoked repeatedly and speed is of concern. The user provides as arguments:

  • a time series as input vector x,
  • wl, the number of samples in x (or the window length),
  • a = $2/wl$,
  • c = $cos(p)$ where $p$=2*π*round(UInt16, (f*wl)/sr)/wl, $f$ is the desired frequency and $sr$ the sampling rate,
  • s = $sin(p)$
  • d = 2*c
  • startAt as in the goertzel function.

See also: goertzel, goertzel2.

source
FourierAnalysis.goertzel2Function
function goertzel2( x  :: Vector{T},
                    f  :: IntOrReal,
                    sr :: Int,
                    wl :: Int,
                startAT:: Int = 1) where T<:Union{Real, Complex}

Like the goertzel function, but allows estimating the DFT coefficient in the whole positive real line up to the Nyquist frequency. This is useful when the DFT coefficient is sought for a frequency which is not a discrete Fourier Frequency.

Using an appropriate taper this function allows almost exact recovery of both amplitude and phase at all frequencies in case of one sinusoid, whereas the goertzel function can do so only at exact Fourier discrete frequencies.

See also: goertzel_fast, goertzel.

Examples:

using FourierAnalysis
sr, t, f, a = 128, 128, 5, 10
ftrue=f+0.15
v=sinusoidal(a, ftrue, sr, t, 0 )
c=goertzel(v, f, sr, t) # should be 0+aim, but clearly fails
d=goertzel2(v, ftrue, sr, t) # this get closer
source