Extracting phase using Goertzel filter in Python

Question

I'm investigating whether it's practical to use a Goertzel filter to decipher a binary message phase-shift encoded into a single frequency.

Here is my result:

I would like the bottom line to reflect the discontinuities visible in the top line.

Here is the code generating the second and third lines:

samps_per_sec = SAMP_RATE
cycles_per_sec = CARRIER_FREQ_Hz
rad_per_sec = two_pi * cycles_per_sec
rad_per_samp = rad_per_sec / samps_per_sec

import cmath

def Goertzel( signal ):

  w = rad_per_samp
  exp_w = np.exp( 1j*w )
  exp_nw = np.exp( -1j*w ) # or exp_w.conj()

  _s = 0.
  __s = 0.

  out_s= []
  out_phase= []
  for i in range( 0, len( signal ) ):
    s = signal[i] + 2.*exp_w.real * _s - __s
    out_s += [s]

    y = s - exp_nw * _s

    phase = cmath.phase( y ) - i * w # subtract expected phase increment

    phasewrap = ( phase + np.pi) % two_pi - np.pi
    out_phase += [ phasewrap / np.pi ] # range -1 to +1

    __s = _s
    _s = s

  export_WAV_mono( "out_s.WAV", 0.001 * np.array( out_s ) )
  export_WAV_mono( "out_phase.WAV", 0.1 * np.array( out_phase ) )

It seems as though the third line is just tracking the phase of the second. Before the signal starts the phase is effectively random (I put a very low microscopic noise level into the signal).

When the second line amplitude goes down to 0 there is a tremor in the phase, which makes sense: as the amplitude approaching zero it is going to go everywhere.

What is really worrying is that when the wave inverts its phase, the filter doesn't seem to be able to register this, it is as if the forward thrusters turn off and the backward thrusters turn on.

I can't see how I can get a useful result from this process.

Is the approach fundamentally flawed? I've copied the maths from the Wikipedia page on Goertzel algorithm, and cross checked it with several implementations I found from various articles.

I suspect that maybe the implementation is actually correct, and I just had a false expectation.

Here is the complete Python script:

#!/usr/local/bin/python

import numpy as np
from itertools import *
from array import array


# generator expression
# similar to list comprehensions, but create one value at a time
# def white_noise( amp=1. ):
#   return ( np.random.uniform(-amp, +amp) for _ in count(0) )

two_pi = 2. * np.pi

SAMP_RATE = 44100
SAMPS_PER_WAVE = 20
WAVES_PER_HALF_BIT = 10

CARRIER_FREQ_Hz = SAMP_RATE / float( SAMPS_PER_WAVE )

SAMPS_PER_HALF_BIT = SAMPS_PER_WAVE * WAVES_PER_HALF_BIT

import wave

def export_WAV_mono( filepath, samps ):
  # the largest possible signed 16-bit integer
  S16MAX = float( 2 ** 15 - 1 )
  samps_sint16 = S16MAX * samps.clip(-1.,+1.)
  data = array( 'h', samps_sint16.astype(int) )

  wv = wave.open( filepath, "w" )

  wv.setparams( (
    1,                  # nchannels
    2,                  # sampwidth in bytes
    44100,              # framerate
    len( data ),        # nframes
    'NONE',             # comptype
    'not compressed'    # compname
    ) )

  # this is the crucial step: the .tostring() method forces this into the string format that AIFC requires
  wv.writeframes( data.tostring() )

  wv.close()

# Fs -> 2 pi
# 1 -> 2 pi / Fs
# CARRIER_FREQ_Hz -> ?
#theta = ( two_pi / SAMP_RATE ) * CARRIER_FREQ_Hz

samps_per_sec = SAMP_RATE
cycles_per_sec = CARRIER_FREQ_Hz
rad_per_sec = two_pi * cycles_per_sec
rad_per_samp = rad_per_sec / samps_per_sec

import cmath

def Goertzel( signal ):
  # https://sites.google.com/site/hobbydebraj/home/goertzel-algorithm-dtmf-detection
  # https://dsp.stackexchange.com/questions/145/estimating-onset-time-of-a-tone-burst-in-noise/151#151

  w = rad_per_samp
  exp_w = np.exp( 1j*w )
  exp_nw = np.exp( -1j*w ) # or exp_w.conj()

  _s = 0.
  __s = 0.

  out_s= []
  out_phase= []
  for i in range( 0, len( signal ) ):
    s = signal[i] + 2.*exp_w.real * _s - __s
    out_s += [s]

    y = s - exp_nw * _s

    phase = cmath.phase( y ) - i * w # subtract expected phase increment

    phasewrap = ( phase + np.pi) % two_pi - np.pi
    out_phase += [ phasewrap / np.pi ] # range -1 to +1

    __s = _s
    _s = s

  export_WAV_mono( "out_s.WAV", 0.001 * np.array( out_s ) )
  export_WAV_mono( "out_phase.WAV", 0.1 * np.array( out_phase ) )


def main( ):
  binary_signal = [1,0,1] # ,1,0,1,0,1]

  BEFORE = 100
  AFTER = 100

  signal = [ ]

  PHASE_SHIFT = np.pi

  phase = 0.
  for b in binary_signal:
    phase += PHASE_SHIFT

    counter = 0
    while True:
      for i in range( 0, SAMPS_PER_HALF_BIT ):
        s = np.sin( phase )
        signal += [s]
        phase += rad_per_samp

      counter += 1
      if counter == 2:
        break

      if b:
        phase += PHASE_SHIFT



  print len( signal)
  noise_len = BEFORE+len(signal)+AFTER
  amp = 0.001
  sig = amp * ( 2. * np.random.random( noise_len ) - 1. )

  # SIGNAL_SAMP_OFFSET = 50
  for i in range( 0, len(signal) ):
    sig[BEFORE + i] += signal[ i ];

  export_WAV_mono( "signal.WAV", 0.25 * sig )


  Goertzel( sig )

if __name__ == '__main__':
  main( )

EDIT: Links Estimating onset time of a tone burst in noise? http://asp.eurasipjournals.com/content/2012/1/56 http://www.mstarlabs.com/dsp/goertzel/goertzel.html

Answer 1

I think you may be using a Phillips head screwdriver to drive a torx head screw (you're using the wrong tool for the problem).

What you're trying to detect is the phase discontinuities in the original data signal. These are local discontinuities: they happen once, and then we don't care about past ones any more.

What the Goertzel algorithm is doing is different: It's going to calculate the Fourier coefficient (amplitude and phase) of the entire data sequence --- including all past phase transitions.

The only way you can know what that should look like is if you knew the data sequence modulating your sinusoid in the first place.

To try to understand this myself, I wrote the scilab code below (sorry, python is too much of a stretch on a Friday night!). What I'm trying to calculate is the "true" phase of the signal at every time instant... as you can see, the phase is a little... unpredictable --- because the data is unknowable (well, random).

In the figure: the red circles are the "true" phase at the end of each bit period, the blue line is the "true" phase for every time instant, the red and green lines are the phase (and the negative of the phase) of the carrier.

The bit sequence that was generated for this particular realization was:

1.  1. -1.  1. -1. -1. -1.  1. -1. -1.  

scialab script only below.

bits = bool2s(rand(1,10,'normal')>0)*2-1;    
N = 32;
f0 = 2/N;
ONE = cos(2*%pi*f0*[0:N-1] + %pi/3);
data = kron(bits,ONE);

clear phi

for k=1:length(data)
    ft = sum(data(1:k).*exp(-%i*2*%pi*[0:k-1]*f0));    
    phi(k) = atan(imag(ft),real(ft));
    disp(idx)
end

clf
plot(phi)
boundries = [N:N:length(phi)];
plot(boundries,phi(boundries),'ro');
plot([0 N*10], [%pi/3 %pi/3],'r');
plot([0 N*10], -[%pi/3 %pi/3],'g');

Answer 2

Instead of using a changing length Goertzel filter, use two fixed length Goertzel filters. Your current filter is getting one point longer at each output point, which continually changes the frequency response. A sliding Goertzel filter, if you can solve stability issues, may require less computation than recomputing a fixed length filter at every input sample.

With two fixed length filters, each one bit-time long, the pair offset by one bit-time, you can check for maxima in the phase difference between the two complex Goertzel (or 1-bin DFT) filters to mark a phase shift in the input at the point right between the two filters.