test
This commit is contained in:
parent
570149c2ee
commit
8d3f1564f4
3 changed files with 182 additions and 38 deletions
24
xfer_tracer/flattop.py
Executable file
24
xfer_tracer/flattop.py
Executable file
|
|
@ -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()
|
||||
|
|
@ -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.])
|
||||
|
||||
|
|
|
|||
|
|
@ -11,6 +11,7 @@ pkgs.mkShell {
|
|||
ps.pyusb
|
||||
ps.numpy
|
||||
ps.pyaudio
|
||||
ps.matplotlib]))
|
||||
ps.matplotlib
|
||||
ps.scipy]))
|
||||
];
|
||||
}
|
||||
|
|
|
|||
Loading…
Reference in a new issue