# 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( )

• I haven't parsed through your code, but I'm wondering if the phase issue is similar to the one found in this question? Your question is different, however, because the signal's phase will change with time --- and with each new phase transition. I'm not sure you can find an analytic expression for what the phase should be for the standard Goertzel. You may have to introduce a forgetting factor, so that only two or so "chips" need to be accounted for. – Peter K. Mar 15 '14 at 1:37

## 2 Answers

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');


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.

• Thankyou HotPaw! I get it! And if I make my bin 1/3 of a bit-length, then cut my input signal into chunks of this size, then when the bit transitions somewhere between chunk[k] and chunk[k+1], the chunk either side i.e. k-1 and k+2 will be pure, so I can just compare chunk[u] with chuck[u-3]. – P i Mar 15 '14 at 19:14