# 101. TCLab in the frequency domain¶

[1]:

import numpy

[2]:

%matplotlib inline

[3]:

tau_p = 150
K_p = 0.33
theta = 15

[4]:

omega = numpy.logspace(-3, -2)
s = omega*1j

[5]:

G = K_p/(tau_p*s + 1)*numpy.exp(-theta*s)

[6]:

import matplotlib.pyplot as plt


Let’s choose 5 logarithmically spaced points around the corner frequency

[7]:

freqs = numpy.logspace(-2.8, -2, 5)

[8]:

def plotfreqs(ax):
for freq in freqs:
ax.axvline(freq)

[9]:

def bode(omega, G, gainax=None, phaseax=None, phasecorr=0):
if gainax is None:
fig, (gainax, phaseax) = plt.subplots(2, 1, sharex=True)
gainax.loglog(omega, numpy.abs(G))
angle = numpy.angle(G)
phaseax.semilogx(omega, numpy.unwrap(angle) + phasecorr)
return gainax, phaseax

[10]:

gainax, phaseax = bode(omega, G)
plotfreqs(gainax)


## 101.1. Direct frequency domain tests¶

[11]:

Qbar = 50
deltaQ = A = 10


How long does one sine wave take to repeat?

$P = 2\pi/\omega$
[12]:

P = 2*numpy.pi/freqs


We’ll do a couple of repeats.

[13]:

nperiods = 3

[14]:

switch_times = numpy.concatenate([[0], numpy.cumsum(P*nperiods)])

[15]:

switch_times

[15]:

array([    0.        , 11893.26574888, 19397.409123  , 24132.20349893,
27119.65678502, 29004.61237718])

[16]:

from tclab import runexperiment, labtime

[17]:

def send_sine_wave(t, lab):
print(f'\rTime: {t} Last sleep: {labtime.lastsleep} ', end='')
step = numpy.max(numpy.nonzero(switch_times <= t))
lab.Q1(Qbar + A*numpy.sin(freqs[step]*(t - switch_times[step])))

[18]:

totaltime = sum(nperiods*P)
totaltime

[18]:

29004.612377178186


Experiment will take this many hours

[19]:

totaltime/60/60

[19]:

8.056836771438386


Looks like we’ll be running it overnight.

[57]:

%%time
experiment = runexperiment(send_sine_wave, connected=True,
plot=False, twindow=1000,
time=int(totaltime),
speedup=1,
dbfile='sinetest.db')

TCLab version 0.4.6dev
NHduino connected on port /dev/cu.wchusbserial1410 at 115200 baud.
TCLab Firmware 1.3.0 Arduino Uno.
Time: 4872.0 Last sleep: 0.8694100379943848 TCLab disconnected successfully.
CPU times: user 15.5 s, sys: 28 s, total: 43.6 s
Wall time: 1h 21min 16s

[37]:

import tclab

[38]:

h = tclab.Historian((('Q1', lambda: [0, 0, 0, 0]),
('Q2', None),
('T1', None),
('T2', None)), dbfile='sinetest.db')

[39]:

# h = experiment.historian

[40]:

h.get_sessions()

[40]:

[(2, '2018-03-06 18:45:27', 13710),
(12, '2018-03-07 14:55:43', 2001),
(15, '2018-03-07 18:43:39', 7526),
(25, '2018-03-08 05:34:09', 5523),
(27, '2018-03-08 07:10:23', 4873),
(28, '2018-03-08 12:59:29', 55),
(29, '2018-03-08 13:00:31', 116),
(30, '2018-03-08 13:02:35', 1001),
(31, '2018-03-08 13:25:17', 2001),
(32, '2018-03-08 14:30:19', 0),
(33, '2018-03-08 14:30:32', 891),
(34, '2018-03-08 14:46:08', 536),
(35, '2018-03-08 14:55:18', 132),
(36, '2018-03-08 15:02:28', 2001),
(37, '2018-03-09 04:37:56', 0),
(38, '2018-03-09 04:39:17', 0),
(39, '2018-03-09 04:45:49', 0),
(40, '2018-03-09 04:46:44', 0),
(41, '2018-03-09 08:22:48', 0),
(42, '2018-03-09 08:23:20', 0)]

[41]:

import pandas

[42]:

%matplotlib inline

[43]:

# h = experiment.historian

[44]:

sine_sessions = [2, 15, 25, 27]

[45]:

def sinefit(session, freq, Tbar=40, inphase=0, gain=0.4, phase=0):
t = numpy.array(h.t)
Q1 = numpy.array(h.logdict['Q1'])
Q1_sine = A*numpy.sin(t*freq + inphase)
T1 = numpy.array(h.logdict['T1'])
T1_sine = A*gain*numpy.sin(t*freq + inphase + phase)
plt.plot(t, Q1 - Qbar, color='blue', alpha=0.4)
plt.plot(t, Q1_sine, color='blue')
plt.plot(t, T1 - Tbar, color='red', alpha=0.4)
plt.plot(t, T1_sine, color='red')
plt.ylim(-A, A)
print(f'Gain={gain}, Phase={phase-inphase}')

[46]:

sinefit(2, freqs[0], 40, 0, 0.3, 0, )

Gain=0.3, Phase=0

[47]:

from ipywidgets import interact

[48]:

interact(sinefit,
session=sine_sessions,
freq=freqs,
Tbar=(30., 50.),
inphase=(0, 2*numpy.pi),
gain=(0., 1., 0.01),
phase=(-2*numpy.pi, 0),)

[48]:

<function __main__.sinefit>

[124]:

gains = [0.3, 0.32, 0.25, 0.24, 0.19]
phases = [-0.28, -0.38, -0.58, -0.88, -1.28]

[158]:

gainax, phaseax = bode(omega, G)
gainax.scatter(freqs, gains, color='red')
phaseax.scatter(freqs, phases, color='red')

[158]:

<matplotlib.collections.PathCollection at 0x11a430ba8>


## 101.2. FFT based bode diagram¶

[20]:

next_values = {'time': 0, 'value': +A}
def random_noise(t, lab):
if t > next_values['time']:
next_values['time'] += numpy.random.uniform(50, 300)
next_values['value'] = -A if next_values['value'] == +A else +A
lab.Q1(Qbar + next_values['value'])

[21]:

next_values

[21]:

{'time': 0, 'value': 10}

[22]:

%matplotlib notebook

[23]:

%%time
experiment = runexperiment(random_noise, connected=True,
plot=True, twindow=2000,
time=8*60*60,
speedup=1,
dbfile='sinetest.db',
)

TCLab version 0.4.6dev
NHduino connected on port /dev/cu.wchusbserial1410 at 115200 baud.
TCLab Firmware 1.3.0 Arduino Uno.

TCLab disconnected successfully.

---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
<timed exec> in <module>()

~/Documents/Development/TCLab/tclab/experiment.py in runexperiment(function, connected, plot, twindow, time, dbfile, speedup, synced)
99     """
100     with Experiment(connected, plot, twindow, time, dbfile, speedup, synced) as experiment:
--> 101         for t in experiment.clock():
102             function(t, experiment.lab)
103     return experiment

~/Documents/Development/TCLab/tclab/experiment.py in clock(self)
78         else:
79             times = range(self.time)
---> 80         for t in times:
81             yield t
82             if self.plot:

~/Documents/Development/TCLab/tclab/labtime.py in clock(period, step, tol, adaptive)
98                            'Step size was {} s, but {:.2f} s elapsed '
99                            '({:.2f} too long). Consider increasing step.')
--> 100                 raise RuntimeError(message.format(step, elapsed, elapsed-step))
101         labtime.sleep(step - (labtime.time() - start) % step)
102         now = labtime.time() - start

RuntimeError: Labtime clock lost synchronization with real time. Step size was 1 s, but 3.22 s elapsed (2.22 too long). Consider increasing step.

[181]:

h = experiment.historian

[30]:

h.load_session(36)

[31]:

%matplotlib inline

[32]:

plt.plot(h.t, h.logdict['T1'])

[32]:

[<matplotlib.lines.Line2D at 0x118ee7a20>]

[33]:

import numpy

[72]:

startcut = 900

[73]:

Q1 = numpy.array(h.logdict['Q1'][startcut:])
T1 = numpy.array(h.logdict['T1'][startcut:])

[74]:

plt.plot(Q1)
plt.plot(T1)

[74]:

[<matplotlib.lines.Line2D at 0x1176f1278>]

[75]:

u = numpy.fft.rfft(Q1 - Qbar)
y = numpy.fft.rfft(T1 - numpy.mean(T1))

[76]:

omegafft = numpy.fft.rfftfreq(len(Q1), 1)
Gfft = y/u

[160]:

ga, pa = bode(omega, G)
bode(omegafft, Gfft, ga, pa, numpy.pi/4)
plt.xlim(1e-3, 1e-2)
ga.set_ylim(0.01, 1)
pa.set_ylim(-numpy.pi/2, 0)

[160]:

(-1.5707963267948966, 0)

[ ]:



[ ]: