#
# Loudspeaker impedance meter
# Method: stepped sine
#
import sounddevice as sd
import numpy as np
from scipy.io.wavfile import write as write_wav
from scipy.signal.windows import tukey as window
import matplotlib.pyplot as plt
import time
# -----------------[Mode (calibration/measurement)]---------------------
CALIB_MODE = False
# ----------------------[Measurement settings]--------------------------
START_FREQUENCY = 5.0
END_FREQUENCY = 40000.0
FREQ_STEP_PER_OCT = 12
MIN_SIGNAL_PERIODS = 10
MIN_SIGNAL_TIME = 0.4
MIN_SILENCE_PERIODS = 4
MIN_SILENCE_TIME = 0.05
TIME_OFFSET = -0.04
VOLUME = 0.25
R_SENSE = 15.24
USE_WINDOWING = True
USE_CALIB_VALUES = True
# ------------------------[Debug settings]------------------------------
USE_SAVE_WAV = False
USE_REPLAY = False
# ------------------------[Device settings]-----------------------------
# device name:
#DEVICE_NAME = 'sysdefault' #linux
#DEVICE_NAME = 'E-MU 0202 | USB: Audio (hw:1,0)' #linux
#DEVICE_NAME = 'Yamaha Steinberg USB ASIO' #win10
DEVICE_NAME = 'Steinberg UR22: USB Audio (hw:1,0)' #linux
# sample rate:
FS = 192000
# ----------------------------------------------------------------------
# ----------------------------------------------------------------------
# ----------------------------------------------------------------------
# audio device selection and setting
print('Searching sound devices...')
device_number = -1
for i, val in enumerate(sd.query_devices()):
print ('\t' + str(i) + '\t' + val['name'])
if val['name'] == DEVICE_NAME:
device_number = i
if device_number == -1:
print ('ERROR! "' + DEVICE_NAME + '" not found!')
raise Exception('Device not found')
print ('Initialize device: #' + str(device_number) + ' ' + DEVICE_NAME)
sd.default.device = device_number
sd.default.samplerate = FS
sd.default.channels = 2
sd.default.dtype = 'float64'
#sd.default.dither_off = True
# ---------------------------- m a i n ---------------------------------
start_time = time.time()
calib_arr_C = []
zma_file_content = []
graph_freq = []
graph_abs = []
graph_deg = []
if CALIB_MODE:
print('Calibration mode')
else:
calib_arr_C = np.load('calibrate.npy')
print('Measurement mode')
#---------------------------meas loop-----------------------------------
print('| frequency |peak channel 1 |peak channel 2 |mod| magnitude | angle |')
print('|-------------|---------------|---------------|---|-------------|--------------|')
index = 0
freq = START_FREQUENCY
while freq <= END_FREQUENCY:
# ----- generate stimulant signal -----
# calculate number of sine periodes with aligned to the end of the last period
num_periodes = -int(-max(MIN_SIGNAL_TIME * freq, MIN_SIGNAL_PERIODS))
# signal length in samples
num_samples = int(num_periodes * FS / freq)
# silence before/after the sine signal (in sample)
silence_samples = int(max(MIN_SILENCE_TIME, MIN_SILENCE_PERIODS / freq) * FS)
# generate silence signal
signal = np.float64([0] * int(silence_samples * 2 + num_samples))
# sine wave position
position = silence_samples + int(TIME_OFFSET * FS)
# generated sine wave into signal
signal[position:position + num_samples] = np.float64(VOLUME * (window(num_samples, 0.5) if USE_WINDOWING else 1.0) * np.sin(2.0 * np.pi * freq * np.arange(0, num_samples) / FS))
# ----- play and record signals with sound device -----
# soundcard play and rec signal (full duplex)
sound_record = sd.playrec(signal, samplerate=FS, channels=2, dtype='float64', blocking=True)
# index swap: [sample][L/R]->[L/R][sample]
sound_record_t = np.transpose(sound_record)
# save generated and digitized signals to wav files (for debug)
if USE_SAVE_WAV:
write_wav('py_records/signal_'+str(freq)+'.wav', FS, np.array(signal * np.iinfo(np.int32).max).astype(np.int32))
write_wav('py_records/record_'+str(freq)+'.wav', FS, np.array(sound_record * np.iinfo(np.int32).max).astype(np.int32))
# replay digitized signal (for debug)
if USE_REPLAY:
sd.play(sound_record_t, blocking=True)
# calculate peak values of input channels
peak_ch1_dB = 20 * np.log10(np.max(np.abs(sound_record_t[0])))
peak_ch2_dB = 20 * np.log10(np.max(np.abs(sound_record_t[1])))
# ----- signal processing -----
# calculated sin and cos correlations
len_rec = len(sound_record_t[0])
base_sin = np.sin(2.0 * np.pi * freq * np.arange(0, len_rec) / FS)
base_cos = np.cos(2.0 * np.pi * freq * np.arange(0, len_rec) / FS)
# correlation scalars (sin and cos components) for both channels
corr_sin_ch1 = np.average(sound_record_t[0] * base_sin)
corr_cos_ch1 = np.average(sound_record_t[0] * base_cos)
corr_sin_ch2 = np.average(sound_record_t[1] * base_sin)
corr_cos_ch2 = np.average(sound_record_t[1] * base_cos)
# conversion sin/cos components to complex numbers
ch1_C = corr_sin_ch1 + corr_cos_ch1 * 1j
ch2_C = corr_sin_ch2 + corr_cos_ch2 * 1j
# transfer function calculation
meas_C = ch1_C / ch2_C # ch1/ch2 relative of complex number
# ----- calculate physical values and writing results -----
if CALIB_MODE:
# make calibration
calib_arr_C.append(meas_C)
#print ('C:\t%3.2f\t%3.2f' % (np.abs(meas_C), np.rad2deg(np.angle(meas_C))))
print ('|{:>10.2f} Hz|{:>10.2f} dBfs|{:>10.2f} dBfs|CAL|{:>10.2f} dB|{:>10.2f} deg|'.format(freq,peak_ch1_dB,peak_ch2_dB,np.abs(meas_C),np.angle(meas_C, deg=True)))
else:
# correction of the measurement value with the calibration values
if USE_CALIB_VALUES:
meas_C = meas_C / calib_arr_C[index]
# calculate complex impedance
Z_C = R_SENSE / ((1 / meas_C) - 1)
# conversion of complex number to magnitude and phase
Z_abs = np.abs(Z_C)
Z_deg = np.angle(Z_C, deg=True)
print ('|{:>10.2f} Hz|{:>10.2f} dBfs|{:>10.2f} dBfs|Z |{:>9.2f} Ohm|{:>10.2f} deg|'.format(freq,peak_ch1_dB,peak_ch2_dB,Z_abs,Z_deg))
# make zma
zma_file_content.append(str(freq) + '\t' + str(Z_abs) + '\t' + str(Z_deg))
# make graph arrays
graph_freq.append(freq)
graph_abs.append(Z_abs)
graph_deg.append(Z_deg)
# calculate next frequency
index += 1
freq = START_FREQUENCY * np.power(2, index / FREQ_STEP_PER_OCT)
#------------------------end of meas loop-------------------------------
execution_time = int(time.time() - start_time)
print ('Execution time: '+str(int(execution_time / 60))+':'+str(execution_time % 60))
# -----------------witing result files and make grap--------------------
if CALIB_MODE:
# writing calibration values to file
print('writing calibrate file ...')
np.save('calibrate.npy', calib_arr_C)
else:
# writing zma file
print('writing impedance.zma file ...')
with open('impedance.zma', 'w') as f:
f.write('\n'.join(zma_file_content))
# show impedance curves into graph
print('plotting impedance graph ...')
figure, axis = plt.subplots(2, 1)
axis[0].semilogx(graph_freq, graph_abs)
axis[0].grid()
axis[0].grid(which='minor', axis='x', linestyle=':')
axis[1].semilogx(graph_freq, graph_deg)
axis[1].grid()
axis[1].grid(which='minor', axis='x', linestyle=':')
plt.show()
print('Done')