Skip to content

Commit

Permalink
Merge branch 'ml' of https://github.com/npeard/self-interferometry in…
Browse files Browse the repository at this point in the history 
…to ml
  • Loading branch information
@ange1a-j14
ange1a-j14 committed Dec 11, 2024
2 parents afad4b0 + ad0782a commit a97eef4
Show file tree
Hide file tree
Showing 2 changed files with 84 additions and 65 deletions.
145 changes: 82 additions & 63 deletions acquire_automatic.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6,16 +6,20 @@
import matplotlib.pyplot as plt
import redpitaya_scpi as scpi
import numpy as np
from numpy.fft import fft
from numpy.fft import fft, fftfreq
# from scipy.fft import fft # use numpy
import math
import util
from timeit import timeit
from datetime import datetime
from scipy.signal import butter, filtfilt

IP = 'rp-f0c04a.local'
rp_s = scpi.scpi(IP)
print('Connected to ' + IP)
# plt.rcParams.update({
# "text.usetex": True
# })

def run_one_shot(start_freq=1, end_freq=1000, ampl=0.1, gen_dec=8192, acq_dec=256, num_acq=1,
store_data=False, plot_data=False, filename='data.h5py'):
Expand Down Expand Up @@ -53,7 +57,7 @@ def run_one_shot(start_freq=1, end_freq=1000, ampl=0.1, gen_dec=8192, acq_dec=25

##### Generation #####
# Function for configuring Source
rp_s.sour_set(1, wave_form, ampl, freq, data=y)
rp_s.sour_set(1, wave_form, 1, freq, data=ampl*y)

# Enable output
rp_s.tx_txt('OUTPUT1:STATE ON')
Expand All @@ -68,75 +72,85 @@ def run_one_shot(start_freq=1, end_freq=1000, ampl=0.1, gen_dec=8192, acq_dec=25
rp_s.tx_txt('ACQ:RST')
rp_s.acq_set(dec=acq_dec, trig_delay=8192)

for i in range(num_acq):
rp_s.tx_txt('ACQ:START')
# ! OS 2.00 or higher only ! ##
time.sleep(0.1)
rp_s.tx_txt('ACQ:TRig NOW') # CH2_PE
# Wait for trigger
while 1:
rp_s.tx_txt('ACQ:TRig:STAT?') # Get Trigger Status
if rp_s.rx_txt() == 'TD': # Triggered?
# print("td", datetime.now().time())
break
while 1:
rp_s.tx_txt('ACQ:TRig:FILL?')
if rp_s.rx_txt() == '1':
# print("filled", datetime.now().time())
break
##### Analysis #####
# Read data and plot function for Data Acquisition
pds = np.array(rp_s.acq_data(chan=1, convert=True))
speaker = np.array(rp_s.acq_data(chan=2, convert=True))
vels, _, _= util.velocity_waveform(speaker, acq_smpl_rate)
acq_time_data = np.linspace(0, N-1, num=N) / acq_smpl_rate
rp_s.tx_txt('ACQ:START')
# ! OS 2.00 or higher only ! ##
time.sleep(0.4)
rp_s.tx_txt('ACQ:TRig NOW') # CH2_PE
time.sleep(0.4)
# Wait for trigger
while 1:
rp_s.tx_txt('ACQ:TRig:STAT?') # Get Trigger Status
if rp_s.rx_txt() == 'TD': # Triggered?
# print("td", datetime.now().time())
break
while 1:
rp_s.tx_txt('ACQ:TRig:FILL?')
if rp_s.rx_txt() == '1':
# print("filled", datetime.now().time())
break
##### Analysis #####
# Read data and plot function for Data Acquisition
pds = np.array(rp_s.acq_data(chan=1, convert=True))
speaker = np.array(rp_s.acq_data(chan=2, convert=True))
acq_time_data = np.linspace(0, N-1, num=N) / acq_smpl_rate

if plot_data:
pd_data.append(pds) # Volts
speaker_data.append(speaker) # Volts
vel_data.append(vels)

if store_data:
# Store data in h5py file
path = "/Users/angelajia/Code/College/SMI/data/"
# filename = "training_data.h5py"
file_path = os.path.join(path, filename)

entries = {
'Speaker (V)': speaker,
'Speaker (Microns/s)': vels,
'PD (V)': pds
}
util.write_data(file_path, entries)
if plot_data:
pd_data = pds
speaker_data = speaker
# pd_data.append(pds) # Volts
# speaker_data.append(speaker) # Volts
# vel_data.append(vels)

if store_data:
# Store data in h5py file
path = "/Users/angelajia/Code/College/SMI/data/"
# filename = "training_data.h5py"
file_path = os.path.join(path, filename)
entries = {
'Speaker (V)': speaker,
'Speaker (Microns/s)': vels,
'PD (V)': pds
}
util.write_data(file_path, entries)

if plot_data:
gen_time_data = np.linspace(0, N-1, num=N) / gen_smpl_rate
pd_data = np.concatenate(pd_data)
speaker_data = np.concatenate(speaker_data)
vel_data = np.concatenate(vel_data)
# pd_data = np.concatenate(pd_data)
# speaker_data = np.concatenate(speaker_data)
# vel_data = np.concatenate(vel_data)
y_vel, y_converted, _ = util.velocity_waveform(ampl*y, gen_smpl_rate)
avg_speaker = np.mean(np.reshape(speaker, (-1, 32)), axis=1)
print(avg_speaker.shape)
vel_data, vel_converted, freqs= util.velocity_waveform(avg_speaker, acq_smpl_rate)

fig, ax = plt.subplots(nrows=4)

ax[0].plot(pd_data, color='blue', label='Observed PD')
ax[0].plot(speaker_data, color='black', label='Observed Drive')
ax[0].plot(pd_data, color='blue', label='Observed PD')# , marker='.')
ax[0].plot(speaker_data, color='black', label='Observed Drive') #, marker='.')
# ax[0].plot(time_data, ampl*y, label='Drive Output', alpha=0.5)
# ax[0].legend()
ax[0].legend()
ax[0].set_ylabel('Amplitude (V)')
ax[0].set_xlabel('Samples')

#ax[1].set_title(r'$\tilde{F}^{-1}(VelTransferFunc(f) * \tilde{F}(Voltage Time Series)(f) )$')
ax[1].set_title("Vel for acq_dec=256")
ax[1].plot(vel_data)
ax[1].set_xlabel('Samples')
ax[1].set_ylabel('Expected Vel (Microns/s)')

ax[2].plot(gen_time_data, y_vel, label='Drive Output', alpha=0.7)
ax[2].set_ylabel('Expected Vel (Microns/s)')
ax[2].set_xlabel('Time (s)')

# ax[3].set_title("Vel for gen_dec=8192")
# ax[3].plot(gen_time_data, y_vel, label='Drive Output', alpha=0.7)
# ax[3].set_ylabel('Expected Vel (Microns/s)')
# ax[3].set_xlabel('Time (s)')
# ax[1].legend()

ax[3].plot(gen_time_data, ampl*y)
ax[3].set_xlabel('Time (s)')
ax[3].set_ylabel('Amplitude (V)')
# ax[2].set_title("Generated Speaker Voltage (gen_dec=8192)")
# ax[2].plot(gen_time_data, ampl*y)
# ax[2].set_xlabel('Time (s)')
# ax[2].set_ylabel('Amplitude (V)')

ax[2].plot(fftfreq(speaker.shape[0], d=1/acq_smpl_rate), np.abs(fft(speaker, norm='ortho')), marker='.')
ax[3].plot(freqs, np.abs(vel_converted), marker='.')
plt.tight_layout()
plt.show()

Expand All @@ -145,15 +159,20 @@ def run_one_shot(start_freq=1, end_freq=1000, ampl=0.1, gen_dec=8192, acq_dec=25
rp_s.tx_txt('ACQ:RST')


num_shots = 2000
filenames = ['test_1to1kHz_misaligned_invertspectra_trigdelay8192_sleep100ms_2kx1shots_randampl.h5py',
'val_1to1kHz_misaligned_invertspectra_trigdelay8192_sleep100ms_2kx1shots_randampl.h5py']
shots = [1, 0]
amplitude = 0
acq_num = 1
print("start", datetime.now().time())
for i in range(num_shots):
amplitude = np.random.uniform(0.1, 0.6)
if i % 400 == 0:
print("\t", datetime.now().time())
print(f"\t{i*acq_num}: ampl = {amplitude}")
run_one_shot(1, 1000, ampl=amplitude, gen_dec=8192, acq_dec=256, num_acq=acq_num, store_data=True, plot_data=False,
filename='test_1to1kHz_invertspectra_trigdelay8192_sleep100ms_2kx1shots_dec=256_8192_randampl.h5py')
print("end", datetime.now().time())
for (file, num_shots) in zip(filenames, shots):
print(file)
print(num_shots)
for i in range(num_shots):
amplitude = 0.1 # np.random.uniform(0.1, 0.6)
if i % 500 == 0:
print("\t", datetime.now().time())
print(f"\t{i*acq_num}: ampl = {amplitude}")
run_one_shot(100, 101, ampl=amplitude, gen_dec=8192, acq_dec=256, num_acq=acq_num, store_data=False, plot_data=True,
filename=file)
print("end", datetime.now().time())
4 changes: 2 additions & 2 deletions acquire_continuous.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -12,12 +12,12 @@
print('Connected to ' + IP)

wave_form = "SINE"
freq = 1000
freq = 120
ampl = 0.1

N = 16384 # Number of samples in buffer
SMPL_RATE_DEC1 = 125e6 # sample rate for decimation=1 in Samples/s (Hz)
decimation = 32
decimation = 256# 32
smpl_rate = SMPL_RATE_DEC1//decimation

# Reset Generation and Acquisition
Expand Down

0 comments on commit a97eef4

Please sign in to comment.