Voices from harmonics

import matplotlib.pyplot as plt
import numpy as np
from scipy.interpolate import RegularGridInterpolator

import sonore as so

plt.rcParams.update({"font.size": 9, "axes.titlesize": 10, "figure.dpi": 100})


def finish(snd):
    """How every sound in the gallery is played: 5 ms ramps, RMS 0.1, peak at most 0.95."""
    snd = snd.ramp(5e-3).normalize(rms=0.1)
    return snd.normalize(peak=0.95) if snd.peak > 0.95 else snd


# The sentence at its native 16 kHz, and its F0 track (WORLD Harvest, 0 where unvoiced).
# Sources: docs/speech/SOURCES.md.
sentence = finish(so.load("docs/speech/bdl_arctic_a0131.flac"))
fs, duration = sentence.fs, sentence.duration
f0_times, f0_harvest = np.loadtxt("docs/speech/bdl_arctic_a0131_f0.csv", delimiter=",", skiprows=2).T

stft = so.STFT(sentence, win_dur=0.040, hop_dur=0.005)
cep = so.Cepstrum(stft)
t_cep, f0_cep, _ = cep.f0(f_lo=75, f_hi=400)
track = so.f0_track(sentence)

# Three contours, each a (times, values) pair with 0 where unvoiced; an F0Track also works as is.
contours = {
    "cepstral": (t_cep, f0_cep[0]),
    "so.f0_track": (track.t, track.f0[0]),
    "Harvest": (f0_times, f0_harvest),
}
FMAX = 4000


def show(snd, title, contour=None):
    """Waveform and narrowband spectrogram (Hann 33 ms, harmonics resolved), on one time axis."""
    fig = plt.figure(figsize=(10, 4.2), layout="constrained")
    ax0, ax1 = fig.subplots(2, 1, sharex=True, height_ratios=[0.45, 1])
    snd.plot(ax0, color="k", lw=0.4)
    ax0.set(title=f"Waveform: {title}", xlabel="")
    so.STFT(snd, win_dur=0.0333, hop_dur=0.002).plot(ax1, db_range=60, colorbar=False, fmax=FMAX)
    ax1.set_title("Spectrogram (Hann 33 ms)")
    if contour is not None:
        t, f = contour
        ax1.plot(t, np.where(f > 0, f, np.nan), color="c", lw=1.2)
    for ax in (ax0, ax1):
        ax.set_xlim(0, duration)
    return fig, [ax0, ax1]

for name, (_, f) in contours.items():
    print(f"{name:12s} voiced on {np.mean(f > 0):.0%} of its time windows")

voiced_t, voiced_f0 = track.t[track.voiced[0]], track.f0[0][track.voiced[0]]
cutoffs = 0.5 / np.exp(np.interp(cep.t, voiced_t, np.log(voiced_f0)))  # half a period, per time window
envelope = cep.lifter(cutoffs).envelope()[0]  # (frequencies, time windows)
log_envelope = RegularGridInterpolator(
    (cep.t, stft.f), np.log(envelope.T), bounds_error=False, fill_value=None
)


def formants(t, f):
    """The envelope's gain at times t and frequencies f, interpolated in dB between time windows and bins."""
    return np.exp(log_envelope(np.stack([t, f], axis=-1)))

noise = so.gaussian_noise(duration, fs, rng=2)
coloured = so.STFT(noise, win_dur=0.040, hop_dur=0.005)
coloured.data = envelope[None] * np.exp(1j * np.angle(coloured.data))
breath = coloured.to_sound().data[: len(sentence), 0]

# Voicing between time windows, 1 where voiced; the harmonics are switched at the same time windows.
t_samples = np.arange(len(sentence)) / fs
voicing = np.interp(t_samples, track.t, track.voiced[0].astype(float))


def rms(x):
    return np.sqrt(np.mean(x**2))


# How loud the unvoiced stretches of the sentence are, relative to the voiced ones.
x = sentence.data[:, 0]
balance = rms((1 - voicing) * x) / rms(voicing * x)
print(f"unvoiced stretches: {20 * np.log10(balance):.1f} dB relative to the voiced ones")


def with_breath(voice, noise=breath):
    """The voiced sound, plus envelope-shaped noise where unvoiced, at the sentence's own balance."""
    gap = (1 - voicing) * noise
    return so.Sound(voice.data[:, 0] + gap * balance * rms(voice.data[:, 0]) / rms(gap), fs)

voiced = track.voiced[0]
median = np.median(track.f0[0][voiced])
print(f"median voiced F0: {median:.0f} Hz")


def on(values):
    """The resynthesis on another contour, at so.f0_track's window times."""
    contour = (track.t, np.where(voiced, values, 0.0))
    return finish(with_breath(so.harmonic_complex(duration, fs, contour, amplitudes=formants))), contour

def resynthesize(snd, seed=2):
    """A pitch track and a cepstral envelope from snd, put back together as above:
    harmonics on so.f0_track shaped by the envelope, and envelope-shaped noise where unvoiced."""
    snd_stft = so.STFT(snd, win_dur=0.040, hop_dur=0.005)
    snd_cep = so.Cepstrum(snd_stft)
    snd_track = so.f0_track(snd)
    is_voiced = snd_track.voiced[0]
    voiced_times, voiced_f0s = snd_track.t[is_voiced], snd_track.f0[0][is_voiced]
    lifter_cutoffs = 0.5 / np.exp(np.interp(snd_cep.t, voiced_times, np.log(voiced_f0s)))
    snd_envelope = snd_cep.lifter(lifter_cutoffs).envelope()[0]
    snd_log_envelope = RegularGridInterpolator(
        (snd_cep.t, snd_stft.f), np.log(snd_envelope.T), bounds_error=False, fill_value=None
    )

    def snd_formants(t, f):
        return np.exp(snd_log_envelope(np.stack([t, f], axis=-1)))

    snd_duration = snd.duration
    harmonics = so.harmonic_complex(snd_duration, fs, snd_track, amplitudes=snd_formants).data[:, 0]
    noise_stft = so.STFT(so.gaussian_noise(snd_duration, fs, rng=seed), win_dur=0.040, hop_dur=0.005)
    noise_stft.data = snd_envelope[None] * np.exp(1j * np.angle(noise_stft.data))
    snd_breath = noise_stft.to_sound().data[: len(snd), 0]
    sample_times = np.arange(len(snd)) / fs
    snd_voicing = np.interp(sample_times, snd_track.t, is_voiced.astype(float))
    original = snd.data[:, 0]
    snd_balance = rms((1 - snd_voicing) * original) / rms(snd_voicing * original)
    gap = (1 - snd_voicing) * snd_breath
    out = harmonics + gap * snd_balance * rms(harmonics) / rms(gap)
    return finish(so.Sound(out, fs)), (snd_track.t, snd_track.f0[0])


sentence_female = finish(so.load("docs/speech/slt_arctic_a0131.flac"))
resynthesis_female, contour_female = resynthesize(sentence_female)
f0_female = contour_female[1]
print(f"median voiced F0: {np.median(f0_female[f0_female > 0]):.0f} Hz (slt), {median:.0f} Hz (bdl)")