466 lab4 -- Sage

CES466 Lab 4

Exploring the effectiveness of the Goertzel Algorithm for different frequencies and sampling rates.

Here we first implement the algorithm itself and then make a little wrapper function for simulating it's response.

Notice that k and \omega are defined a bit differently than in the lab reference. This way gives improved SNR and matching to the target frequency.

X = [] data = [] def goertzel(samples, ft = 10000, fs = 30000, Ns = 24): w, k = var('w,k') #k = floor(1/2 + ft/fs*Ns) # this #w = 2*pi*k/Ns k = ft/fs w = 2*pi*k coef = 2*cos(w) Q1 = 0 Q2 = 0 for s in samples: Q0 = N(s + coef*Q1 - Q2) Q2 = Q1 Q1 = Q0 return N(Q2^2 + Q1^2 - coef*Q1*Q2) f_sin(t,fg,phi) = sin(2*pi*fg*(t + phi)) def f_real(t, fg=None, phi=0): t += phi period = 100e-6 while t > period: t -= period while t < 0: t += period #a = -100783.31852998595 #b = 0.1961888003229936 a,b = (-100783.3185299859, 0.549278242519127) f(u) = b*e^(a*u) if t < 1.76e-6: T = t/1.76e-6 val = (1-T)*0 + T*0.46 elif t < period/2: val = f(t) else: val = -f_real(t - period/2) return val def gen_data(fg = 10000, ft = 10000, fs = 30000, Ns = 24, phi = 0, fn = f_real): global X, data start = 0 stop = Ns/fs nw = 0 sw = 217 X = srange(start, stop, (stop-start)/Ns) data = [sw*fn(x,fg,phi) + nw*random() for x in X] return data def goertzel_sim(fg = 10000, ft = 10000, fs = 30000, Ns = 24, fn = f_real): data = gen_data(fg, ft, fs, Ns, fn=fn) return goertzel(data, ft, fs, Ns)

In this next section we look at how the values given by the Goertzel across a range of frequencies. Notice that not only does 20 samples per Goertzel result in a lower response than 24 samples, but its peak response is not on the target frequency.

Here we approximate fairly closely the actual signal received by the ADC on the MSP430F2012 from a square wave transmission, as opposed to a perfect sine wave. The maximum and minimum Goertzel results can be calculated. This process takes a long time to do and so isn't shown in this published version of the worksheet.

FS = 29950 real_sim = [goertzel(gen_data(fs=FS, phi=i*25.0/100e-6), fs=FS) for i in range(25)]

Traceback (click to the left of this block for traceback)
...
__SAGE__

Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "_sage_input_3.py", line 10, in <module>
    exec compile(u'open("___code___.py","w").write("# -*- coding: utf-8 -*-\\n" + _support_.preparse_worksheet_cell(base64.b64decode("RlMgPSAyOTk1MApyZWFsX3NpbSA9IFtnb2VydHplbChnZW5fZGF0YShmcz1GUywgcGhpPWkqMjUuMC8xMDBlLTYpLCBmcz1GUykgZm9yIGkgaW4gcmFuZ2UoMjUpXQ=="),globals())+"\\n"); execfile(os.path.abspath("___code___.py"))
  File "", line 1, in <module>
    
  File "/tmp/tmpcKFo4T/___code___.py", line 4, in <module>
    exec compile(u'real_sim = [goertzel(gen_data(fs=FS, phi=i*_sage_const_25p0 /_sage_const_100en6 ), fs=FS) for i in range(_sage_const_25 )]
  File "", line 1, in <module>
    
  File "/tmp/tmpDu2rbv/___code___.py", line 53, in gen_data
    data = [sw*fn(x,fg,phi) + nw*random() for x in X]
  File "/tmp/tmpDu2rbv/___code___.py", line 27, in f_real
    while t > period: t -= period
KeyboardInterrupt
__SAGE__

min(real_sim), max(real_sim)

Traceback (click to the left of this block for traceback)
...
NameError: name 'real_sim' is not defined

Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "_sage_input_4.py", line 10, in <module>
    exec compile(u'open("___code___.py","w").write("# -*- coding: utf-8 -*-\\n" + _support_.preparse_worksheet_cell(base64.b64decode("bWluKHJlYWxfc2ltKSwgbWF4KHJlYWxfc2ltKQ=="),globals())+"\\n"); execfile(os.path.abspath("___code___.py"))
  File "", line 1, in <module>
    
  File "/tmp/tmpymbfNL/___code___.py", line 2, in <module>
    exec compile(u'min(real_sim), max(real_sim)
  File "", line 1, in <module>
    
NameError: name 'real_sim' is not defined

In this next section we look at how the values given by the Goertzel across a range of frequencies.

%hide T = srange(5000,14000,50) target = 10000 def make_line(Ns = 24, color=(0,0,1)): r1 = [goertzel_sim(fg = f, ft = target, fs = 3*target, Ns=Ns, fn=f_sin) for f in T] m1 = max(r1) mi1 = r1.index(m1) html("<p>N=%d Max value: $%.4f$, At freq: $%d$Hz"%(Ns, m1, T[mi1])) l1 = line(zip(T,r1), color=color) return m1, l1 html("Frequency response for $F_T = %d$Hz and $F_S = %d$Hz"%(target,3*target)) m1, l1 = make_line(20, color=(0,0,1)) m2, l2 = make_line(24, color=(0,.8,0)) t = line(((target,0),(target,max(m1,m2))), color=(1,0,0)) l1+l2+t

Frequency response for F_T = 10000Hz and F_S = 30000Hz
N=20 Max value: 4599208.6887, At freq: 9900Hz
N=24 Max value: 6780816.0000, At freq: 10000Hz

Frequency response for F_T = 10000Hz and F_S = 30000Hz
N=20 Max value: 4599208.6887, At freq: 9900Hz
N=24 Max value: 6780816.0000, At freq: 10000Hz

Next we examine the minimum relative signal response as a function of how long the true signal is active during the sampling period. Not pictures, but on average, the signal needs to be active over half of the time for a 50% response. However, the minimum represents a worst case scenario. Interestingly, the best case requires still just 2 fewer live samples than the worst case. Worst is shown in blue and best in gray.

%hide import random as rand import numpy as np S = 24 # N x = range(S) # controls dead samples runs = 100 # number of random runs to make worst = [10000000.0] * S best = [0.0] * S for i in range(runs): partial = gen_data(Ns = S, fs = 30000, phi = rand.random()) next = [goertzel([0]*(S-u) + partial[u:]) for u in x] worst = [min(next[i], worst[i]) for i in range(S)] best = [max(next[i], best[i]) for i in range(S)] html("N = $%d$"%S) best = np.array(best) best = best/max(worst) worst = np.array(worst) worst = worst/max(worst) l_worst = line(zip(x,list(worst))) l_best = line(zip(x,list(best)), color=(.7,.7,.7)) l_fifty = line(((0,0.5), (24,0.5)), color = (.8,.5,.2)) failure = list(worst<0.5).index(True) l_min = line(((failure - 1, 0), (failure - 1, 1)), color=(0,.8,0)) html("After $%d$ experiments, the signal can be dead for up to $%d$ of $%d$ samples"%(runs, failure - 1, S)) graph = l_min + l_worst + l_fifty + l_best graph.axes_labels(['# dead samples', 'relative strength']) graph.show()

Traceback (click to the left of this block for traceback)
...
__SAGE__

Traceback (most recent call last):    best = [0.0] * S
  File "", line 1, in <module>
    
  File "/tmp/tmpMb5bFn/___code___.py", line 12, in <module>
    next = [goertzel([_sage_const_0 ]*(S-u) + partial[u:]) for u in x]
  File "/tmp/tmpuLIKXO/___code___.py", line 17, in goertzel
    Q0 = N(s + coef*Q1 - Q2)
KeyboardInterrupt
__SAGE__

Here is some data generated across a range of input frequencies and number of samples for the target frequency of 10kHz and sampling at 30ksps.

I = srange(7500,11000,50) # Received frequency J = srange(16,40,1) # number of samples processed data = [(i, j, goertzel_sim(fg = i, Ns = j)) for i in I for j in J]

And then it is plotted. As you can see, a higher value of N has a great effect. But also notice that some increases in N have much greater effect than others (plateaus and slopes along the rising ridge).

%hide lp = list_plot3d(data, point_list=True, interpolation_type='linear') lp.show(viewer='tachyon')

466 lab4

191 days ago by jay

CES466 Lab 4