Separate sample rates and add invert spectra weighting functionality

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)
+        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())

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