我正在研究使用 Goertzel 滤波器来破译编码为单个频率的二进制消息相移是否可行。
这是我的结果：
我希望底线反映顶线中可见的不连续性。
这是生成第二行和第三行的代码：
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 ) )
似乎第三行只是跟踪第二行的阶段。在信号开始之前，相位实际上是随机的（我在信号中加入了非常低的微观噪声水平）。
当第二条线振幅下降到 0 时，相位会发生震颤，这是有道理的：当振幅接近零时，它会到处传播。
真正令人担忧的是，当波反转相位时，过滤器似乎无法记录这一点，就好像前向推进器关闭而后向推进器打开一样。
我看不出如何从这个过程中获得有用的结果。
这种方法是否存在根本性缺陷？我从 Goertzel 算法的 Wikipedia 页面复制了数学，并与我从各种文章中找到的几个实现进行了交叉检查。
我怀疑也许实现实际上是正确的，我只是有一个错误的期望。
这是完整的 Python 脚本：
#!/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( )
编辑：链接 估计噪声中音调的开始时间？ http://asp.eurasipjournals.com/content/2012/1/56 http://www.mstarlabs.com/dsp/goertzel/goertzel.html
