diff --git a/xfer_tracer/flattop.py b/xfer_tracer/flattop.py deleted file mode 100755 index ef06607..0000000 --- a/xfer_tracer/flattop.py +++ /dev/null @@ -1,24 +0,0 @@ -#!/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 72fb8a2..56a2588 100755 --- a/xfer_tracer/script.py +++ b/xfer_tracer/script.py @@ -6,7 +6,7 @@ import time import pyaudio as pa import numpy as np import scipy as sc - +import matplotlib.pyplot as plt current_phase = 0.0 normalized_f = 0.0 # in cycles per sample : f = normalized_f * f_sampling @@ -71,11 +71,13 @@ paudio = pa.PyAudio() audio_strm = paudio.open(format=pa.paFloat32, channels=1, - rate=48000, + rate=44100, output=True, input=False, stream_callback=audio_callback, start=True) +print(paudio.get_default_output_device_info()) +print(audio_strm.is_active()) def acquire_samples(sc_frequency, oscillo): global new_f @@ -89,16 +91,15 @@ def acquire_samples(sc_frequency, oscillo): oscillo.write("acquire:state run") oscillo.write("data:source CH1") - input_ch = oscillo.query_binary("curve?", container=np.array) + input_ch = oscillo.query_binary_values("curve?",datatype="b",container=np.array).astype(np.float64) oscillo.write("data:source CH2") - output_ch = oscillo.query_binary("curve?", container=np.array) + output_ch = oscillo.query_binary_values("curve?",datatype="b",container=np.array).astype(np.float64) 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] + hscales = [0.25, 100e-3, 50e-3, 25e-3, 10e-3, 5e-3, 2.5e-3, 1e-3, 500e-6, 250e-6, 100e-6, 50e-6, 25e-6, 10e-6] + hscales_str = ["2.5e-1", "1.0e-1", "5.0e-2", "2.5e-2", "1.0e-2", "5.0e-3", "2.5e-3", "1.0e-3", "5.0e-4", "2.5e-4", "1.0e-4", "5.0e-5", "2.5e-5", "1.0e-5"] hscale_index = 0 current_hscale = hscales[hscale_index] @@ -112,13 +113,14 @@ def manage_measures(f_list, oscillo): # These are NOT normalized frequencies W = sc.signal.windows.flattop(2500) for f in f_list: - sc_frequency = np.round(2**20/48000 * f)/2**20 + sc_frequency = np.round(2**20/44100 * 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)) + print("horizontal:main:scale " + hscales_str[hscale_index]) + oscillo.write("horizontal:main:scale " + hscales_str[hscale_index]) cal_done = False print(f, " : ", current_hscale,"s/div, ", num_periods, " periods") @@ -126,14 +128,15 @@ def manage_measures(f_list, oscillo): # These are NOT normalized frequencies if not cal_done: cal_A, cal_B = acquire_samples(0., oscillo) cal_done = True - + print(A) + print(B) 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) + h = np.sum(B*V)/np.sum(A*V) print("H(", f, ") = ", h) H.append(h) return H @@ -147,18 +150,20 @@ def find_oscillo(): # exit(1) return rm.open_resource(R[0]) + osci = find_oscillo() +osci.timeout = 10000 print(osci.query("*IDN?")) osci.write("*RST") -time.sleep(5) +time.sleep(10) -osci.write("ch1:scale 1.0e-1") +osci.write("ch1:scale 2.5e-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:scale 2.5e-1") osci.write("ch2:bandwidth ON") osci.write("ch2:coupling AC") osci.write("select:ch2 1") @@ -170,6 +175,16 @@ 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.]) +osci.write("data:encdg RIBinary") +decade = [1.0, 1.2, 1.5, 1.8, 2.2, 2.7, 3.3, 3.9, 4.7, 5.6, 6.8, 8.2] +X = [68., 82] + list(np.array(decade)*([100.]*len(decade))) + list(np.array(decade) * ([1000.]*len(decade))) + [1e4, 1.2e4, 1.5e4, 1.8e4] +H = manage_measures(X, osci) +Y = 20 * np.log10(np.abs(H)) + +fig, axes = plt.subplots(nrows = 2, ncols=1, sharex = True) +axes[0].semilogx(X,Y) +axes[0].grid() +axes[1].semilogx(X, np.angle(H, deg=True)) +axes[1].grid() +plt.show()