From baf107946e3bd6b8dc5d34ec0ca937b955346dc1 Mon Sep 17 00:00:00 2001 From: Angela Jia <42385351+22angiej@users.noreply.github.com> Date: Sat, 12 Oct 2024 14:21:30 -0700 Subject: [PATCH] Separate sample rates and add invert spectra weighting functionality --- acquire_automatic.py | 163 ++++++++++++++++++++++++------------------- util.py | 102 +++++++++++++-------------- 2 files changed, 138 insertions(+), 127 deletions(-) diff --git a/acquire_automatic.py b/acquire_automatic.py index d0bc862..9440918 100644 --- a/acquire_automatic.py +++ b/acquire_automatic.py @@ -6,15 +6,18 @@ import matplotlib.pyplot as plt import redpitaya_scpi as scpi import numpy as np -from scipy.fftpack import fft +from numpy.fft import fft +# from scipy.fft import fft # use numpy import math import util +from timeit import timeit +from datetime import datetime IP = 'rp-f0c04a.local' rp_s = scpi.scpi(IP) print('Connected to ' + IP) -def run_one_shot(use_freq, max_freq, start_freq=1, end_freq=1000, ampl=0.1, decimation=8192, +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'): """Runs one shot of driving the speaker with a waveform and collecting the relevant data. @@ -30,15 +33,16 @@ def run_one_shot(use_freq, max_freq, start_freq=1, end_freq=1000, ampl=0.1, deci N = 16384 # Number of samples in buffer SMPL_RATE_DEC1 = 125e6 # sample rate for decimation=1 in Samples/s (Hz) - smpl_rate = SMPL_RATE_DEC1//decimation - burst_time = N / smpl_rate + gen_smpl_rate = SMPL_RATE_DEC1//gen_dec + acq_smpl_rate = SMPL_RATE_DEC1//acq_dec + burst_time = N / gen_smpl_rate wave_form = 'ARBITRARY' freq = 1 / burst_time - t, y = util.bounded_frequency_waveform(start_freq, end_freq, length=N, sample_rate=smpl_rate, - use_freq=use_freq, max_freq=max_freq) + t, y = util.bounded_frequency_waveform(start_freq, end_freq, length=N, sample_rate=gen_smpl_rate, invert=True) y = util.linear_convert(y) # convert range of waveform to [-1, 1] to properly set ampl + if plot_data: plt.plot(t, y) plt.show() @@ -54,89 +58,102 @@ def run_one_shot(use_freq, max_freq, start_freq=1, end_freq=1000, ampl=0.1, deci # Enable output rp_s.tx_txt('OUTPUT1:STATE ON') rp_s.tx_txt('SOUR1:TRig:INT') + # print("output enabled", datetime.now().time()) ##### Acqusition ##### + pd_data = [] + speaker_data = [] + vel_data = [] # Function for configuring Acquisition - rp_s.acq_set(dec=decimation, trig_delay=0) - rp_s.tx_txt('ACQ:START') - time.sleep(1) - rp_s.tx_txt('ACQ:TRig CH2_PE') - - # Wait for trigger - while 1: - rp_s.tx_txt('ACQ:TRig:STAT?') # Get Trigger Status - if rp_s.rx_txt() == 'TD': # Triggered? - break - ## ! OS 2.00 or higher only ! ## - while 1: - rp_s.tx_txt('ACQ:TRig:FILL?') - if rp_s.rx_txt() == '1': - break - - ##### Analaysis ##### - # Read data and plot function for Data Acquisition - pd_data = np.array(rp_s.acq_data(chan=1, convert=True)) # Volts - speaker_data = np.array(rp_s.acq_data(chan=2, convert=True)) # Volts - velocity_data, converted_signal, freq = util.velocity_waveform(speaker_data, smpl_rate, use_freq, max_freq) - displ_data, _, _ = util.displacement_waveform(speaker_data, smpl_rate, use_freq, max_freq) - y_vel, y_converted, _ = util.velocity_waveform(ampl*y, smpl_rate, use_freq, max_freq) - time_data = np.linspace(0, N-1, num=N) / smpl_rate - + 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 + + 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: - fig, ax = plt.subplots(nrows=3) - - ax[0].plot(time_data, pd_data, color='blue', label='Observed PD') - ax[0].plot(time_data, speaker_data, color='black', label='Observed Drive') - ax[0].plot(time_data, ampl*y, label='Drive Output', alpha=0.5) - ax[0].legend() + 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) + y_vel, y_converted, _ = util.velocity_waveform(ampl*y, gen_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(time_data, ampl*y, label='Drive Output', alpha=0.5) + # ax[0].legend() ax[0].set_ylabel('Amplitude (V)') - ax[0].set_xlabel('Time (s)') - - ax[1].plot(freq, np.abs(fft(speaker_data)), color='black', label='Observed Drive', marker='.') - ax[1].plot(freq, np.abs(converted_signal), color='green', label='Expected Observed Vel', marker='.') - ax[1].plot(freq, np.abs(fft(ampl*y)), color='blue', label='Expected Drive', marker='.') - ax[1].plot(freq, np.abs(y_converted), color='orange', label='Expected Ideal Vel', marker='.') - ax[1].loglog() - ax[1].set_xlabel('Frequency (Hz)') - ax[1].set_ylabel('$|\^{V}|$') - ax[1].legend() - - ax[2].plot(time_data, velocity_data, color='black', label='Observed Drive') - ax[2].plot(time_data, y_vel, label='Drive Output') + ax[0].set_xlabel('Samples') + + 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[2].legend() + # ax[1].legend() + + ax[3].plot(gen_time_data, ampl*y) + ax[3].set_xlabel('Time (s)') + ax[3].set_ylabel('Amplitude (V)') plt.tight_layout() plt.show() - - 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 = { - 'Time (s)': time_data, - 'Speaker (V)': speaker_data, - 'Speaker (Microns/s)': velocity_data, - 'PD (V)': pd_data, - 'Speaker (Microns)': displ_data - } - - util.write_data(file_path, entries) ##### Reset when closing program ##### rp_s.tx_txt('GEN:RST') rp_s.tx_txt('ACQ:RST') + num_shots = 2000 amplitude = 0 -end_freq = 1000 -max_freq = 10*end_freq +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(f"{i}: ampl = {amplitude}") - run_one_shot(True, max_freq, 30, end_freq, ampl=amplitude, decimation=256, store_data=True, plot_data=False, - filename='test_max10kHz_30to1kHz_2kshots_dec=256_randampl.h5py') - # print(i) \ No newline at end of file + 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()) \ No newline at end of file diff --git a/util.py b/util.py index fdd0613..256d5f4 100644 --- a/util.py +++ b/util.py @@ -1,8 +1,38 @@ import numpy as np -from scipy.fftpack import fft, ifft, fftfreq +from numpy.fft import fft, ifft, fftfreq import h5py -def bounded_frequency_waveform(start_frequency, end_frequency, length, sample_rate, use_freq, max_freq): +# Constants from calibration_rp using RPRPData.csv +f0 = 257.20857316296724 +Q = 15.804110908084784 +k = 33.42493417407945 +c = -3.208233068626455 + +def A(f): + """Calculates the expected displacement of the speaker at an inputted drive amplitude 'ampl' for a given frequency 'f', + based on the calibration fit at 0.2Vpp. + + Args: + f (1darr): frequencies at which to calculate expected displacement + + Returns: + 1darr: expected displacement/V_ampl in microns/V + """ + return (k * f0**2) / np.sqrt((f0**2 - f**2)**2 + f0**2*f**2/Q**2) + +def phase(f): + """Calculates the phase delay between the speaker voltage waveform and the photodiode response + at a given frequency 'f'. + + Args: + f (1darr): frequencies at which to calculate expected displacement + + Returns: + 1darr: phase in radians + """ + return np.arctan2(f0/Q*f, f**2 - f0**2) + c + +def bounded_frequency_waveform(start_frequency, end_frequency, length, sample_rate, invert=False): """Generates a random waveform within the given frequency range of a given length. Args: @@ -17,20 +47,20 @@ def bounded_frequency_waveform(start_frequency, end_frequency, length, sample_ra # Create an evenly spaced time array t = np.linspace(0, 1.0, length, False) # 1 second # Generate a random frequency spectrum between the start and end frequencies - if use_freq: - freq = np.linspace(0, max_freq, length//2, False) # replaced sample_rate/2 with end_frequency - else: - freq = np.linspace(0, sample_rate/2, length//2, False) - spectrum = np.random.uniform(0, 1, len(freq)) + freq = np.linspace(0, sample_rate/2, length//2, False) + spectrum = np.random.uniform(0.0, 1.0, len(freq)) spectrum = np.where((freq >= start_frequency) & (freq <= end_frequency), spectrum, 0) c = np.random.rayleigh(np.sqrt(spectrum*(freq[1]-freq[0]))) # See Phys. Rev. A 107, 042611 (2023) ref 28 for why we use the Rayleigh distribution here # Unless we use this distribution, the random noise will not be Gaussian distributed - phase = np.random.uniform(-np.pi, np.pi, len(freq)) + phi = np.random.uniform(-np.pi, np.pi, len(freq)) # Use the inverse Fourier transform to convert the frequency domain signal back to the time domain # Also include a zero phase component - spectrum = np.hstack([c*spectrum*np.exp(1j*phase), np.zeros_like(spectrum)]) - y = np.real(ifft(spectrum)) + spectrum = spectrum * c * np.exp(1j*phi) + if invert: + spectrum = np.divide(spectrum, A(freq) * np.exp(1j*phase(freq))) + spectrum = np.hstack([spectrum, np.zeros_like(spectrum)]) + y = np.real(ifft(spectrum, norm="ortho")) y = np.fft.fftshift(y) return t, y @@ -70,37 +100,7 @@ def write_data(file_path, entries): maxshape=(None, col_data.shape[0]), chunks=True) -# Constants from calibration_rp using RPRPData.csv -f0 = 257.20857316296724 -Q = 15.804110908084784 -k = 33.42493417407945 -c = -3.208233068626455 - -def A(f): - """Calculates the expected displacement of the speaker at an inputted drive amplitude 'ampl' for a given frequency 'f', - based on the calibration fit at 0.2Vpp. - - Args: - f (1darr): frequencies at which to calculate expected displacement - - Returns: - 1darr: expected displacement/V_ampl in microns/V - """ - return (k * f0**2) / np.sqrt((f0**2 - f**2)**2 + f0**2*f**2/Q**2) - -def phase(f): - """Calculates the phase delay between the speaker voltage waveform and the photodiode response - at a given frequency 'f'. - - Args: - f (1darr): frequencies at which to calculate expected displacement - - Returns: - 1darr: phase in radians - """ - return np.arctan2(f0/Q*f, f**2 - f0**2) + c - -def displacement_waveform(speaker_data, sample_rate, use_freq, max_freq): +def displacement_waveform(speaker_data, sample_rate): """Calculates the corresponding displacement waveform based on the given voltage waveform using calibration. @@ -113,21 +113,18 @@ def displacement_waveform(speaker_data, sample_rate, use_freq, max_freq): converted displacement waveform in frequency domain, frequency array (Hz) """ - speaker_spectrum = fft(speaker_data) + speaker_spectrum = fft(speaker_data, norm="ortho") n = speaker_data.size - if use_freq: - sample_spacing = 1/(2*max_freq) - else: - sample_spacing = 1/sample_rate + sample_spacing = 1/sample_rate freq = fftfreq(n, d=sample_spacing) # units: cycles/s = Hz # Multiply signal by transfer func in freq domain, then return to time domain converted_signal = speaker_spectrum * A(freq) * np.where(freq < 0, np.exp(-1j*phase(-freq)), np.exp(1j*phase(freq))) - y = np.real(ifft(converted_signal)) + y = np.real(ifft(converted_signal, norm="ortho")) return y, converted_signal, freq -def velocity_waveform(speaker_data, sample_rate, use_freq, max_freq): +def velocity_waveform(speaker_data, sample_rate): """Calculates the corresponding velocity waveform based on the given voltage waveform using calibration. @@ -140,16 +137,13 @@ def velocity_waveform(speaker_data, sample_rate, use_freq, max_freq): converted velocity waveform in frequency domain, frequency array (Hz) """ - speaker_spectrum = fft(speaker_data) + speaker_spectrum = fft(speaker_data, norm="ortho") n = speaker_data.size - if use_freq: - sample_spacing = 1/(2*max_freq) - else: - sample_spacing = 1/sample_rate + sample_spacing = 1/sample_rate freq = fftfreq(n, d=sample_spacing) # units: cycles/s = Hz # Multiply signal by transfer func in freq domain, then return to time domain converted_signal = 1j*freq * speaker_spectrum * A(freq) * np.where(freq < 0, np.exp(-1j*phase(-freq)), np.exp(1j*phase(freq))) - v = np.real(ifft(converted_signal)) + v = np.real(ifft(converted_signal, norm="ortho")) return v, converted_signal, freq \ No newline at end of file