Skip to content

Commit bf775da

Browse files
ange1a-j14ange1a-j14
committed
Add velocity conversion and automatic data collection
1 parent 7529b21 commit bf775da

File tree

2 files changed

+164
-117
lines changed

2 files changed

+164
-117
lines changed

‎acquire_automatic.py

Lines changed: 114 additions & 105 deletions
Original file line numberDiff line numberDiff line change
@@ -9,113 +9,122 @@
99
from scipy.fftpack import fft
1010
import math
1111
import util
12-
from datetime import datetime
13-
import csv
1412

1513
IP = 'rp-f0c04a.local'
1614
rp_s = scpi.scpi(IP)
1715
print('Connected to ' + IP)
1816

19-
store_data = False # whether or not to save data to CSV
20-
21-
##### Create Waveform #####
22-
23-
N = 16384 # Number of samples in buffer
24-
SMPL_RATE_DEC1 = 125e6 # sample rate for decimation=1 in Samples/s (Hz)
25-
decimation = 8192
26-
smpl_rate = SMPL_RATE_DEC1//decimation
27-
burst_time = N / smpl_rate
28-
29-
wave_form = 'ARBITRARY'
30-
freq = 1 / burst_time
31-
ampl = 0.3 # good range 0-0.6V
32-
33-
# t, y from exampled RP arbitrary wavegen:
34-
# t = np.linspace(0, 1, N)*2*math.pi
35-
# y = np.sin(t) + 1/3*np.sin(3*t) # same overall period as regular sin wave
36-
t, y = util.bounded_frequency_waveform(50, 55, length=N, sample_rate=smpl_rate)
37-
y = util.linear_convert(y) # convert range of waveform to [-1, 1] to properly set ampl
38-
plt.plot(t, y)
39-
plt.show()
40-
41-
##### Reset Generation and Acquisition ######
42-
rp_s.tx_txt('GEN:RST')
43-
rp_s.tx_txt('ACQ:RST')
44-
45-
##### Generation #####
46-
# Function for configuring Source
47-
rp_s.sour_set(1, wave_form, ampl, freq, data=y)
48-
49-
# Enable output
50-
rp_s.tx_txt('OUTPUT1:STATE ON')
51-
rp_s.tx_txt('SOUR1:TRig:INT')
52-
53-
##### Acqusition #####
54-
# Function for configuring Acquisition
55-
rp_s.acq_set(dec=decimation, trig_delay=0)
56-
rp_s.tx_txt('ACQ:START')
57-
time.sleep(1)
58-
rp_s.tx_txt('ACQ:TRig NOW')
59-
time.sleep(1)
60-
61-
# Wait for trigger
62-
while 1:
63-
rp_s.tx_txt('ACQ:TRig:STAT?') # Get Trigger Status
64-
print('got status')
65-
if rp_s.rx_txt() == 'TD': # Triggered?
66-
break
67-
print('triggered')
68-
## ! OS 2.00 or higher only ! ##
69-
while 1:
70-
rp_s.tx_txt('ACQ:TRig:FILL?')
71-
if rp_s.rx_txt() == '1':
72-
break
73-
74-
##### Analaysis #####
75-
# Read data and plot function for Data Acquisition
76-
pd_data = rp_s.acq_data(chan=1, convert=True) # Volts
77-
speaker_data = rp_s.acq_data(chan=2, convert=True) # Volts
78-
displacement_data, converted_signal, freq = util.displacement_waveform(np.array(speaker_data), smpl_rate, ampl)
79-
y_displ, _, _ = util.displacement_waveform(ampl*y, smpl_rate, ampl)
80-
81-
fig, ax = plt.subplots(nrows=3)
82-
83-
time_data = np.linspace(0, N-1, num=N) / smpl_rate
84-
ax[0].plot(time_data, pd_data, color='blue', label='PD')
85-
ax[0].plot(time_data, speaker_data, color='black', label='Speaker')
86-
ax[0].plot(time_data, ampl*y, label='Original signal')
87-
ax[0].legend()
88-
ax[0].set_ylabel('Amplitude (V)')
89-
ax[0].set_xlabel('Time (s)')
90-
91-
ax[1].plot(freq, np.abs(fft(speaker_data)), color='black', label='Speaker')
92-
ax[1].plot(freq, np.abs(converted_signal), color='green', label='Displacement')
93-
ax[1].loglog()
94-
ax[1].set_xlabel('Frequency (Hz)')
95-
ax[1].set_ylabel('$|\^{V}|$')
96-
ax[1].legend()
97-
98-
ax[2].plot(time_data, displacement_data, color='black', label='Speaker')
99-
ax[2].plot(time_data, y_displ, label='Original signal')
100-
ax[2].set_ylabel('Expected Displacement (Microns)')
101-
ax[2].set_xlabel('Time (s)')
102-
ax[2].legend()
103-
plt.tight_layout()
104-
plt.show()
105-
106-
if store_data:
107-
# Store data in csv file
108-
path = "/Users/angelajia/Code/College/SMI/data/"
109-
filename = f"{datetime.now()}.csv"
110-
file_path = os.path.join(path, filename)
111-
112-
with open(file_path, 'w', newline='') as f:
113-
writer = csv.writer(f)
114-
writer.writerow(['Time (s)', 'Speaker (V)', 'Speaker (Microns)', 'PD (V)'])
115-
all_data = np.vstack((time_data, speaker_data, displacement_data, pd_data)).T
116-
for row in all_data:
117-
writer.writerow(row)
118-
119-
##### Reset when closing program #####
120-
rp_s.tx_txt('GEN:RST')
121-
rp_s.tx_txt('ACQ:RST')
17+
'''
18+
Runs one shot of driving the speaker with a waveform and collecting the relevant data.
19+
20+
@param store_data: whether to save time, measured speaker & PD voltage,
21+
and expected speaker velocity data to h5py
22+
@param plot_data: whether to plot data after acquisition
23+
'''
24+
def run_one_shot(start_freq=1, end_freq=1000, decimation=8192, store_data=False, plot_data=False):
25+
##### Create Waveform #####
26+
27+
N = 16384 # Number of samples in buffer
28+
SMPL_RATE_DEC1 = 125e6 # sample rate for decimation=1 in Samples/s (Hz)
29+
smpl_rate = SMPL_RATE_DEC1//decimation
30+
burst_time = N / smpl_rate
31+
32+
wave_form = 'ARBITRARY'
33+
freq = 1 / burst_time
34+
ampl = 0.1 # good range 0-0.6V
35+
36+
t, y = util.bounded_frequency_waveform(start_freq, end_freq, length=N, sample_rate=smpl_rate)
37+
y = util.linear_convert(y) # convert range of waveform to [-1, 1] to properly set ampl
38+
if plot_data:
39+
plt.plot(t, y)
40+
plt.show()
41+
42+
##### Reset Generation and Acquisition ######
43+
rp_s.tx_txt('GEN:RST')
44+
rp_s.tx_txt('ACQ:RST')
45+
46+
##### Generation #####
47+
# Function for configuring Source
48+
rp_s.sour_set(1, wave_form, ampl, freq, data=y)
49+
50+
# Enable output
51+
rp_s.tx_txt('OUTPUT1:STATE ON')
52+
rp_s.tx_txt('SOUR1:TRig:INT')
53+
54+
##### Acqusition #####
55+
# Function for configuring Acquisition
56+
rp_s.acq_set(dec=decimation, trig_delay=0)
57+
rp_s.tx_txt('ACQ:START')
58+
time.sleep(1)
59+
rp_s.tx_txt('ACQ:TRig CH2_PE')
60+
time.sleep(1)
61+
62+
# Wait for trigger
63+
while 1:
64+
rp_s.tx_txt('ACQ:TRig:STAT?') # Get Trigger Status
65+
print('got status')
66+
if rp_s.rx_txt() == 'TD': # Triggered?
67+
break
68+
print('triggered')
69+
## ! OS 2.00 or higher only ! ##
70+
while 1:
71+
rp_s.tx_txt('ACQ:TRig:FILL?')
72+
if rp_s.rx_txt() == '1':
73+
break
74+
75+
##### Analaysis #####
76+
# Read data and plot function for Data Acquisition
77+
pd_data = np.array(rp_s.acq_data(chan=1, convert=True)) # Volts
78+
speaker_data = np.array(rp_s.acq_data(chan=2, convert=True)) # Volts
79+
velocity_data, converted_signal, freq = util.velocity_waveform(speaker_data, smpl_rate)
80+
y_vel, _, _ = util.velocity_waveform(ampl*y, smpl_rate)
81+
time_data = np.linspace(0, N-1, num=N) / smpl_rate
82+
83+
if plot_data:
84+
fig, ax = plt.subplots(nrows=3)
85+
86+
ax[0].plot(time_data, pd_data, color='blue', label='Observed PD')
87+
ax[0].plot(time_data, speaker_data, color='black', label='Observed Drive')
88+
ax[0].plot(time_data, ampl*y, label='Drive Output')
89+
ax[0].legend()
90+
ax[0].set_ylabel('Amplitude (V)')
91+
ax[0].set_xlabel('Time (s)')
92+
93+
ax[1].plot(freq, np.abs(fft(speaker_data)), color='black', label='Observed Drive')
94+
ax[1].plot(freq, np.abs(converted_signal), color='green', label='Expected Observed Vel')
95+
ax[1].loglog()
96+
ax[1].set_xlabel('Frequency (Hz)')
97+
ax[1].set_ylabel('$|\^{V}|$')
98+
ax[1].legend()
99+
100+
ax[2].plot(time_data, velocity_data, color='black', label='Observed Drive')
101+
ax[2].plot(time_data, y_vel, label='Drive Output')
102+
ax[2].set_ylabel('Expected Velocity (Microns/s)')
103+
ax[2].set_xlabel('Time (s)')
104+
ax[2].legend()
105+
plt.tight_layout()
106+
plt.show()
107+
108+
if store_data:
109+
# Store data in h5py file
110+
path = "/Users/angelajia/Code/College/SMI/data/"
111+
filename = "data.h5py"
112+
file_path = os.path.join(path, filename)
113+
114+
entries = {
115+
'Time (s)': time_data,
116+
'Speaker (V)': speaker_data,
117+
'Speaker (Microns/s)': velocity_data,
118+
'PD (V)': pd_data
119+
}
120+
121+
util.write_data(file_path, entries)
122+
123+
##### Reset when closing program #####
124+
rp_s.tx_txt('GEN:RST')
125+
rp_s.tx_txt('ACQ:RST')
126+
127+
128+
num_shots = 1
129+
for i in range(num_shots):
130+
run_one_shot(5, 6, plot_data=True)

‎util.py

Lines changed: 50 additions & 12 deletions
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,6 @@
11
import numpy as np
22
from scipy.fftpack import fft, ifft, fftfreq
3+
import h5py
34

45
'''
56
Generates a random waveform within the given frequency range of a given length.
@@ -34,21 +35,41 @@ def linear_convert(data, new_min=-1, new_max=1):
3435
new_range = new_max - new_min
3536
return new_min + new_range * (data - old_min) / old_range
3637

37-
f0 = 257.1722062654443
38-
Q = 15.680894756413974
39-
k = 32.638601262518705
40-
c = -3.219944358492212
38+
def write_data(file_path, entries):
39+
'''Add data to a given dataset in 'file'. Creates dataset if it doesn't exist; otherwise,
40+
appends.
41+
42+
Keyword arguments:
43+
file_path (string) - the name of the output HDF5 file to which to append data
44+
entries (dict) - dictionary of column name & corresponding data
45+
'''
46+
with h5py.File(file_path, 'a') as f:
47+
for col_name, col_data in entries.items():
48+
if col_name in f.keys():
49+
f[col_name].resize((f[col_name].shape[0] + 1), axis=0)
50+
new_data = np.expand_dims(col_data, axis=0)
51+
f[col_name][-1:] = new_data
52+
else:
53+
f.create_dataset(col_name,
54+
data=np.expand_dims(col_data, axis=0),
55+
maxshape=(None, col_data.shape[0]),
56+
chunks=True)
57+
58+
# Constants from calibration_rp using RPRPData.csv
59+
f0 = 257.20857316296724
60+
Q = 15.804110908084784
61+
k = 33.42493417407945
62+
c = -3.208233068626455
4163

4264
'''
4365
Calculates the expected displacement of the speaker at an inputted drive amplitude 'ampl'
4466
for a given frequency 'f', based on the calibration fit at 0.2Vpp.
4567
46-
@param f: optimal range 20Hz-1kHz
47-
@param ampl: optimal range 0-0.6V
68+
@param f: frequencies at which to calculate expected displacement
4869
@return: expected displacement in microns
4970
'''
50-
def A(f, ampl):
51-
return ampl * (k * f0**2) / np.sqrt((f0**2 - f**2)**2 + f0**2*f**2/Q**2)
71+
def A(f):
72+
return (k * f0**2) / np.sqrt((f0**2 - f**2)**2 + f0**2*f**2/Q**2)
5273

5374
'''
5475
Calculates the phase delay between the speaker voltage waveform and the photodiode response
@@ -64,14 +85,31 @@ def phase(f):
6485
Calculates the corresponding displacement waveform based on the given voltage waveform
6586
using calibration.
6687
'''
67-
def displacement_waveform(speaker_data, sample_rate, ampl, right=True):
68-
fourier_signal = fft(speaker_data)
88+
def displacement_waveform(speaker_data, sample_rate):
89+
speaker_spectrum = fft(speaker_data)
6990
n = speaker_data.size
7091
sample_spacing = 1/sample_rate
7192
freq = fftfreq(n, d=sample_spacing) # units: cycles/s = Hz
7293

7394
# Multiply signal by transfer func in freq domain, then return to time domain
74-
converted_signal = fourier_signal * A(freq, ampl) * np.where(freq < 0, np.exp(-1j*phase(-freq)), np.exp(1j*phase(freq)))
95+
converted_signal = speaker_spectrum * A(freq) * np.where(freq < 0, np.exp(-1j*phase(-freq)), np.exp(1j*phase(freq)))
7596
y = np.real(ifft(converted_signal))
7697

77-
return y, converted_signal, freq
98+
return y, converted_signal, freq
99+
100+
'''
101+
Calculates the corresponding velocity waveform based on the given voltage waveform
102+
using calibration.
103+
'''
104+
def velocity_waveform(speaker_data, sample_rate):
105+
speaker_spectrum = fft(speaker_data)
106+
n = speaker_data.size
107+
sample_spacing = 1/sample_rate
108+
freq = fftfreq(n, d=sample_spacing) # units: cycles/s = Hz
109+
110+
# Multiply signal by transfer func in freq domain, then return to time domain
111+
converted_signal = 1j*freq * speaker_spectrum * A(freq) * np.where(freq < 0, np.exp(-1j*phase(-freq)), np.exp(1j*phase(freq)))
112+
v = np.real(ifft(converted_signal))
113+
# y = np.fft.fftshift(y)
114+
115+
return v, converted_signal, freq

0 commit comments

Comments
 (0)