I wrote some complicated signal analysis in python to perform the conversion.
First I transform the two audio channels into two components that are orthogonal to each other. In threephase there are really three channels (left, right and center) that are 120° offsets from each other, the analysis is easier when working with two channels at 90° (orthogonal) angle.
I use an Hilbert transform to extract the phase envelope of the two channels. There is some relationship between the phase envelope and the alpha and beta coordinates, I derived an explicit expression for that relationship and use that to find the alpha, beta and volume at audio sample resolution. This is the most complicated part of the process, this method only works well for supercollider files, for other files a lookuptable method with some pre-filtering works better.
Code: Select all
t, (L, R), samplerate = load_wav('MistressAndControlBox-ESTIM-v3.wav', ...)
def algorithm_4(L, R):
"""
Calculate alpha, beta and volume from audio channels
Method:
Convert L, R into orthogonal X, Y
use hilbert transform on X, Y
use trig to obtain explicit expressions for position and volume.
"""
[x_in, y_in] = (ab_transform_inv @ potential_to_channel_matrix_inv)[:2, :2] @ (L, R)
print('hilbert')
h_x = scipy.signal.hilbert(x_in)
h_y = scipy.signal.hilbert(y_in)
print('other compute')
# basic equations
# x_in = a * cos(theta) + a * sin(theta) * i
# y_in = b * cos(theta + phi) + b * sin(theta + phi) * i
a = np.abs(h_x)
b = np.abs(h_y)
phi = np.angle(h_x / h_y)
# obtain expression for radius of <x_in, y_in>
# r**2 = a**2 * cos^2(theta) + b**2 * cos^2(theta + phi)
# r**2 * 2 = a**2 + b**2 + c * cos(2 * theta + psi)
# c ** 2 = a ** 4 + b ** 4 + 2 * a ** 2 * b ** 2 * cos(2 * phi)
# c ** 2 = a ** 4 + b ** 4 + 4 * a_b_cos(phi) ** 2 - 2 * a ** 2 * b ** 2
a_b_cos = np.real(h_x) * np.real(h_y) + np.imag(h_x) * np.imag(h_y) # = a * b * cos(phi)
c = np.sqrt(a ** 4 + b ** 4 + 2 * (2 * a_b_cos ** 2 - a ** 2 * b ** 2))
# find radius of (alpha, beta).
r_max = np.sqrt(0.5 * (a ** 2 + b ** 2 + c))
r_min = np.sqrt(0.5 * (a ** 2 + b ** 2 - c))
r = 1 - r_min / r_max
# calculate angle of max radius
a0 = np.pi * (1 / 2)
psi = np.arctan2(
(a ** 2 * np.sin(a0) + b ** 2 * np.sin(a0 + 2 * phi)),
(a ** 2 * np.cos(a0) + b ** 2 * np.cos(a0 + 2 * phi)),
)
psi = psi + 0.5 * np.pi + np.pi
theta = -psi / 2
alpha_max = a * np.cos(theta)
beta_max = b * np.cos(theta + phi)
r_max = np.linalg.norm((alpha_max, beta_max), axis=0)
alpha_min = a * np.sin(theta)
beta_min = b * np.sin(theta + phi)
r_min = np.linalg.norm((alpha_min, beta_min), axis=0)
r = 1 - r_min / r_max
max_angle = np.arctan2(beta_max, alpha_max)
theta = -max_angle * 2
alpha_out = r * np.cos(theta)
beta_out = r * np.sin(theta)
print('compute done')
return alpha_out, beta_out, r_max
a, b, vol = analysis.algorithm_4(L, R)
Generating the alpha and beta funscripts from this information is easy, I apply a low-pass filter at 30hz, throw away any samples with low volume, and then run Ramer-Douglas-Peucker to reduce file size of the funscript down to manageable levels.
Code: Select all
vol = vol * 2
a = np.nan_to_num(a)
b = np.nan_to_num(b)
a = low_pass_filter_2(a, 30, samplerate)
b = low_pass_filter_2(b, 30, samplerate)
mask = vol > 0.05
t_alpha, a_rdp = pybind11_rdp.rdp(np.vstack((t, a)).T[mask][::50], epsilon=0.05).T
t_beta, b_rdp = pybind11_rdp.rdp(np.vstack((t, b)).T[mask][::50], epsilon=0.05).T
export_to_funscript("MistressAndControlBox.alpha.funscript", t_alpha, remap(a_rdp, -1, 1))
export_to_funscript("MistressAndControlBox.beta.funscript", t_beta, remap(b_rdp, -1, 1))
For the volume funscript I employ scipy.signal.find_peaks to find the maximum volume at 0.1s resolution. I call this the volume envelope.
Code: Select all
peaks, peaks_info = scipy.signal.find_peaks(np.nan_to_num(vol), distance=int(0.1 * samplerate))
t_envelope, v_envelope = t[peaks], vol[peaks]
t_envelope, v_envelope = pybind11_rdp.rdp(np.vstack((t_envelope, v_envelope)).T, epsilon=0.005).T
export_to_funscript("MistressAndControlBox.volume.funscript", t_envelope, remap(v_envelope, 0, 1))
To calculate the pulse frequency and pulse width, for each sample I determine whether signal is 'high' or 'low' based on a threshold of 50% of the volume envelope.
Code: Select all
peaks_interp = np.interp(t, t_envelope, v_envelope)
is_high = ((vol > (peaks_interp * 0.5)) * (peaks_interp > .05)).astype(float)
is_low = 1 - is_high
I calculate the time the pulse is 'high' and the time between two subsequent 'low to high' transitions. These are used to directly calculate the pulse width and pulse frequency.
Code: Select all
low_to_high_indices = np.arange(1, len(is_high), 1)[(is_low[:-1] * (is_high[1:])).astype(bool)]
high_to_low_indices = np.arange(1, len(is_high), 1)[(is_high[:-1] * (is_low[1:])).astype(bool)]
if high_to_low_indices[0] < low_to_high_indices[0]:
high_to_low_indices = high_to_low_indices[1:]
if low_to_high_indices[-1] > high_to_low_indices[-1]:
low_to_high_indices = low_to_high_indices[:-1]
tt = t[low_to_high_indices][:-1]
high_span = high_to_low_indices - low_to_high_indices
high_span = high_span[:-1]
high_span = scipy.signal.convolve(high_span, np.ones(5), 'same') / 5 # moving average filter
carrier = 800
pw = high_span / float(samplerate) * 800
interval_span = low_to_high_indices[1:] - low_to_high_indices[:-1]
interval_span = scipy.signal.convolve(interval_span, np.ones(5), 'same') / 5 # moving average filter
pf = float(samplerate) / interval_span
continuous_indices = pw > 50
pw[continuous_indices] = 3
pf[continuous_indices] = 200
pw[np.roll(continuous_indices, 1)] = 3
pf[np.roll(continuous_indices, 1)] = 200
t_pw, pw_rdp = pybind11_rdp.rdp(np.vstack((tt, pw)).T, epsilon=0.2).T
pf = np.log(pf)
t_pf, pf_rdp = pybind11_rdp.rdp(np.vstack((tt, pf)).T, epsilon=0.02).T
pf_rdp = np.exp(pf_rdp)
export_to_funscript("MistressAndControlBox.pulse_frequency.funscript", t_pf, remap(pf_rdp, 0, 200))
export_to_funscript("MistressAndControlBox.pulse_width.funscript", t_pw, remap(pw_rdp, 3, 20))
This method for extracting the pulse width and pulse frequency is unlikely to work with other files.
This code isn't in the repository.
)
