From c39e3cb87e9e152127991b9eb91c6df72808d1ef Mon Sep 17 00:00:00 2001
From: Angela Jia <42385351+22angiej@users.noreply.github.com>
Date: Tue, 16 Jul 2024 16:19:57 -0700
Subject: [PATCH] Get arbitrary waveform generation working

---
 acquire_automatic.py  | 97 +++++++++++++++++++++++++++++++++++++++++++
 acquire_continuous.py |  6 +--
 util.py               | 38 +++++++++++++++++
 3 files changed, 138 insertions(+), 3 deletions(-)
 create mode 100644 acquire_automatic.py
 create mode 100644 util.py

diff --git a/acquire_automatic.py b/acquire_automatic.py
new file mode 100644
index 0000000..e5917b4
--- /dev/null
+++ b/acquire_automatic.py
@@ -0,0 +1,97 @@
+#!/usr/bin/env python3
+
+import os
+import sys
+import time
+import matplotlib.pyplot as plt
+import redpitaya_scpi as scpi
+import numpy as np
+import math
+import util
+from datetime import datetime
+
+IP = 'rp-f0c04a.local'
+rp_s = scpi.scpi(IP)
+print('Connected to ' + IP)
+
+# Set up waveform
+wave_form = "ARBITRARY"
+freq = 100 # good range 10-200Hz
+ampl = 0.3 # good range 0-0.6V
+
+N = 16384 # Number of samples in buffer
+decimation = 85
+smpl_rate = 125e6//decimation
+# t, y from exampled RP arbitrary wavegen:
+# t = np.linspace(0, 1, N)*2*math.pi
+# y = np.sin(t) + 1/3*np.sin(3*t) # same overall period as regular sin wave
+t, y = util.bounded_frequency_waveform(20, 1000, length=N, sample_rate=smpl_rate)
+y = util.linear_convert(y) # convert range of waveform to [-1, 1] to properly set ampl
+plt.plot(t, y)
+plt.show()
+
+# Reset Generation and Acquisition
+rp_s.tx_txt('GEN:RST')
+rp_s.tx_txt('ACQ:RST')
+
+##### Generation #####
+# Function for configuring Source
+rp_s.sour_set(1, wave_form, ampl, freq, data=y)
+
+# Enable output
+rp_s.tx_txt('OUTPUT1:STATE ON')
+rp_s.tx_txt('SOUR1:TRig:INT')
+
+##### Acqusition #####
+# Function for configuring Acquisition
+rp_s.acq_set(dec=decimation, trig_delay=0)
+# print(rp_s.get_settings())
+rp_s.tx_txt('ACQ:START')
+time.sleep(1)
+rp_s.tx_txt('ACQ:TRig NOW')
+# print(rp_s.get_settings())
+time.sleep(1)
+
+# 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
+
+# Read data and plot
+# function for Data Acquisition
+print(rp_s.get_settings())
+pd_data = rp_s.acq_data(chan=1, convert=True) # Volts
+speaker_data = rp_s.acq_data(chan=2, convert=True) # Volts
+print("data shape:", np.array(pd_data).shape)
+print(rp_s.get_settings())
+
+plt.plot(speaker_data, color="black", label="Speaker")        
+plt.plot(pd_data, color="blue", label="PD")
+plt.legend()
+plt.ylabel('Amplitude [V]')
+plt.xlabel('Samples')
+plt.show()
+
+# Store data in txt file
+path = "/Users/angelajia/Code/College/SMI/data/"
+filename = f"{datetime.now()}.txt"
+file_path = os.path.join(path, filename)
+
+with open(file_path, 'x') as f:
+    # f.write(np.array2string(t, threshold=N+1))
+    # f.write("\n" + np.array2string(y, threshold=N+1))
+    f.write("\n")
+    for x in pd_data:
+        f.write(str(x))
+        f.write("\n")
+    f.write("STARTING SPEAKER DATA\n")
+    for y in speaker_data:
+        f.write(str(y))
+        f.write("\n")
\ No newline at end of file
diff --git a/acquire_continuous.py b/acquire_continuous.py
index 1943fe0..27a554e 100644
--- a/acquire_continuous.py
+++ b/acquire_continuous.py
@@ -4,6 +4,7 @@
 import time
 import matplotlib.pyplot as plt
 import redpitaya_scpi as scpi
+import numpy as np
 
 IP = 'rp-f0c04a.local'
 rp_s = scpi.scpi(IP)
@@ -27,8 +28,7 @@
 
 ##### Acqusition #####
 # Function for configuring Acquisition
-rp_s.acq_set(dec=1)
-
+rp_s.acq_set(dec=32, trig_delay=0)
 rp_s.tx_txt('ACQ:START')
 time.sleep(1)
 rp_s.tx_txt('ACQ:TRig AWG_PE')
@@ -50,7 +50,7 @@
 # function for Data Acquisition
 data = rp_s.acq_data(chan=1, convert=True)
 
-plt.plot(10*data)
+plt.plot(data)
 plt.ylabel('Amplitude [V]')
 plt.xlabel('Samples')
 plt.show()
\ No newline at end of file
diff --git a/util.py b/util.py
new file mode 100644
index 0000000..66ccca5
--- /dev/null
+++ b/util.py
@@ -0,0 +1,38 @@
+import numpy as np
+from scipy.fftpack import fft, ifft, fftfreq
+
+'''
+Generates a random waveform within the given frequency range of a given length. 
+'''
+def bounded_frequency_waveform(start_frequency, end_frequency, length=1000, sample_rate=1000):
+    # 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
+    freq = np.linspace(0, sample_rate/2, length//2, False)
+    spectrum = np.random.uniform(0, 1, len(freq))
+    spectrum = np.where((freq >= start_frequency) & (freq <= end_frequency), spectrum, 0)
+    c = np.random.rayleigh(np.sqrt(4*spectrum*(freq[1]-freq[0])))
+    # See Jiang 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))
+
+    # 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([spectrum*np.exp(1j*phase), np.zeros_like(spectrum)])
+    y = np.real(ifft(spectrum))
+    y = np.fft.fftshift(y)
+
+    return t, y
+
+'''
+Linearly scales data to a new range. 
+
+@param data: assumed to be 1D array 
+'''
+def linear_convert(data, new_min=-1, new_max=1):
+    old_min = np.min(data)
+    old_max = np.max(data)
+    old_range = old_max - old_min
+    new_range = new_max - new_min
+    return new_min + new_range * (data - old_min) / old_range
\ No newline at end of file