From 8d3f1564f41563ac8b3eee4efcfd1772a8de32e4 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?S=C3=A9l=C3=A8ne=20Corbineau?= <selene.corbineau@ens.fr>
Date: Tue, 17 Dec 2024 17:55:11 +0100
Subject: [PATCH] test

---
 xfer_tracer/flattop.py |  24 +++++
 xfer_tracer/script.py  | 193 +++++++++++++++++++++++++++++++++--------
 xfer_tracer/shell.nix  |   3 +-
 3 files changed, 182 insertions(+), 38 deletions(-)
 create mode 100755 xfer_tracer/flattop.py

diff --git a/xfer_tracer/flattop.py b/xfer_tracer/flattop.py
new file mode 100755
index 0000000..ef06607
--- /dev/null
+++ b/xfer_tracer/flattop.py
@@ -0,0 +1,24 @@
+#!/usr/bin/env python3
+
+import numpy as np
+from scipy import signal
+from scipy.fft import fft, fftshift
+import matplotlib.pyplot as plt
+
+window = signal.windows.flattop(2500)
+plt.plot(window)
+plt.title("Flat top window")
+plt.ylabel("Amplitude")
+plt.xlabel("Sample")
+
+plt.figure()
+A = fft(window,65536) / (len(window)/2.0)
+freq = np.linspace(-0.5*2500, 0.5*2500, len(A))
+response = 20 * np.log10(np.abs(fftshift(A / abs(A).max())))
+plt.plot(freq, response)
+#plt.axis([-0.5, 0.5, -120, 0])
+plt.title("Frequency response of the flat top window")
+plt.ylabel("Normalized magnitude [dB]")
+plt.xlabel("Normalized frequency [cycles per window]")
+
+plt.show()
diff --git a/xfer_tracer/script.py b/xfer_tracer/script.py
index f0e63fc..72fb8a2 100755
--- a/xfer_tracer/script.py
+++ b/xfer_tracer/script.py
@@ -5,52 +5,171 @@ import time
 
 import pyaudio as pa
 import numpy as np
+import scipy as sc
+
+
+current_phase = 0.0
+normalized_f  = 0.0 # in cycles per sample : f = normalized_f * f_sampling
+
+trans_len = 2000
+intro = 0.5 * (1 - np.cos(2 * np.linspace(0, np.pi/2, trans_len)))
+outro = 0.5 * (1 + np.cos(2 * np.linspace(0, np.pi/2, trans_len)))
+
+intro = np.append(intro, np.ones(1024))
+outro = np.append(outro, np.zeros(1024))
+
+change_f = True
+trans_done = False
+new_f = 0.0
+fade_time = 0
 
-phase = 0.0 # in cycles
-scaled_frequency = 0.0
-# actual frequency = scaled_frequency * sample_rate
-# For repeatability, scaled_frequency should ideally fit in
-# much fewer than 53 bits.
-# More precisely, if scaled_frequency * frame_count is always exact, we get
-# negligible phase errors over 1hr.
 
 def audio_callback(in_data, frame_count, time_info, status):
-	global phase
-	L = np.sin(2 * np.pi * (phase + np.arange(frame_count) * scaled_frequency)).astype(np.float32) # cast to float
-	delta_phase  = scaled_frequency * frame_count # exact
-	delta_phase -= np.floor(delta_phase) # exact
-	phase += delta_phase # error at most ulp(1) = 2^-52
-	if(phase >= 1):
-		phase -= int(phase)
-	outbytes = L.tobytes()
-	#print(L, frame_count)
-	return (outbytes, pa.paContinue)
+	global fade_time
+	global change_f
+	global new_f
+	global current_phase
+	global normalized_f
+	global trans_done
+
+	if change_f == True:
+		change_f = False
+		trans_done = False
+		fade_time = 0
+
+	if(fade_time != 2*trans_len):
+		if(fade_time < trans_len):
+			last_fade = fade_time + frame_count
+			envelope = outro[fade_time:last_fade]
+			fade_time = min(last_fade, trans_len)
+		elif(fade_time >= trans_len and fade_time < 2*trans_len):
+			normalized_f = new_f
+			if(new_f == 0):
+				current_phase = 0. # Avoid weird sub-audio behavior
+			last_fade = fade_time + frame_count
+			envelope = intro[fade_time - trans_len:last_fade - trans_len]
+			fade_time = min(last_fade, 2*trans_len)
+
+	if(fade_time == 2*trans_len):
+		trans_done = True
+
+	sine = np.sin(2 * np.pi * (current_phase + np.arange(frame_count)*normalized_f))
+	delta = frame_count*normalized_f
+	delta -= np.floor(delta)
+
+	if(delta + current_phase > 1):
+		current_phase += delta - 1
+	else:
+		current_phase += delta
+
+	if(fade_time != 2*trans_len):
+		sine *= envelope
+
+	return (sine.astype(np.float32).tobytes(), pa.paContinue)
 
 paudio = pa.PyAudio()
-output_strm = paudio.open(format=pa.paFloat32,
+
+audio_strm = paudio.open(format=pa.paFloat32,
 	channels=1,
 	rate=48000,
 	output=True,
-	stream_callback=audio_callback)
-output_strm.start_stream()
+	input=False,
+	stream_callback=audio_callback,
+	start=True)
 
-input_strm = paudio.open(format=pa.paFloat32,
-    channels=1,
-    rate=48000,
-    input=True)
-input_strm.stop_stream()
+def acquire_samples(sc_frequency, oscillo):
+	global new_f
+	global change_f
+	global trans_done
+	new_f = sc_frequency
+	change_f = True
+	while (change_f == True) or (trans_done == False):
+		time.sleep(0.05)
+	time.sleep(0.2)
 
-def do_frequency(scaled_f):
-    scaled_frequency = scaled_f
-    time.sleep(0.2) # Stall for latency/transitory behavior
-    input_strm.start_stream()
-    L = input_strm.read(int(48000 * 0.1)) # We get 1/sqrt(n) noise rejection and 1/n spurious tone rejection
-    input_strm.stop_stream()
-    return np.frombuffer(L, dtype=np.float32)
+	oscillo.write("acquire:state run")
+	oscillo.write("data:source CH1")
+	input_ch = oscillo.query_binary("curve?", container=np.array)
+	oscillo.write("data:source CH2")
+	output_ch = oscillo.query_binary("curve?", container=np.array)
 
-import matplotlib.pyplot as plt
-L = do_frequency(440.0/48000)
-print(L)
-plt.plot(L)
-plt.show()
+	return input_ch, output_ch
+	#return (input_ch, output_ch)
+
+def manage_measures(f_list, oscillo): # These are NOT normalized frequencies
+	hscales = [25e-3, 10e-3, 5e-3, 2.5e-3, 1e-3,
+	           500e-6, 250e-6, 100e-6, 50e-6, 25e-6, 10e-6]
+
+	hscale_index = 0
+	current_hscale = hscales[hscale_index]
+
+	cal_A = []
+	cal_B = []
+	cal_done = False
+
+	H = []
+
+	W = sc.signal.windows.flattop(2500)
+
+	for f in f_list:
+		sc_frequency = np.round(2**20/48000 * f)/2**20
+		num_periods = f * 10 * current_hscale
+		while(num_periods >= 5 * 2.5):
+			hscale_index += 1
+			current_hscale = hscales[hscale_index]
+			num_periods = f * 10 * current_hscale
+			oscillo.write("horizontal:main:scale " + str(hscales))
+			cal_done = False
+
+		print(f, " : ", current_hscale,"s/div, ", num_periods, " periods")
+		A, B = acquire_samples(sc_frequency, oscillo)
+		if not cal_done:
+			cal_A, cal_B = acquire_samples(0., oscillo)
+			cal_done = True
+
+		A, B = (A - cal_A)*W, (B - cal_B)*W
+		f_acq = 250./current_hscale
+		rel_f = f / f_acq
+
+		# Obtain the Fourier transforms of A, B @ rel_f, with flattop windowing 
+		V = np.exp(-2.0j * np.pi * np.arange(2500) * rel_f)
+		h = np.sum(B*V)/np.sum(A.V)
+		print("H(", f, ") = ", h)
+		H.append(h)
+	return H
+
+def find_oscillo():
+	rm = pyvisa.ResourceManager()
+	R = rm.list_resources()
+	print(R)
+#	if R.empty():
+#		print("No device found")
+#		exit(1)
+	return rm.open_resource(R[0])
+
+osci = find_oscillo()
+print(osci.query("*IDN?"))
+
+osci.write("*RST")
+time.sleep(5)
+
+osci.write("ch1:scale 1.0e-1")
+osci.write("ch1:bandwidth ON")
+osci.write("ch1:coupling AC")
+osci.write("select:ch1 1")
+
+osci.write("ch2:scale 5.0e-2")
+osci.write("ch2:bandwidth ON")
+osci.write("ch2:coupling AC")
+osci.write("select:ch2 1")
+
+
+osci.write("trigger:main:edge:source LINE")
+osci.write("trigger:main:mode NORMAL")
+osci.write("trigger:main:type EDGE")
+
+osci.write("acquire:mode sample")
+osci.write("acquire:stopafter sequence")
+osci.write("data:encdg ascii")
+manage_measures([330., 470., 680., 820., 1000.])
 
diff --git a/xfer_tracer/shell.nix b/xfer_tracer/shell.nix
index b504312..e0221a0 100644
--- a/xfer_tracer/shell.nix
+++ b/xfer_tracer/shell.nix
@@ -11,6 +11,7 @@ pkgs.mkShell {
 			ps.pyusb
 			ps.numpy
 			ps.pyaudio
-			ps.matplotlib]))
+			ps.matplotlib
+			ps.scipy]))
 	];
 }