Extracting phase using Goertzel filter in Python

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