This notebook investigates the difference between two spectral envelope parametrization / recovering methods below:
code_spectral_envelope
and decode_spectral_envelope
I measuered normalized mean square error (NMSE) between spectral envelope and reconstructed spectral envelope. Synthesized audio examples are available in the notebook.
%pylab inline
Populating the interactive namespace from numpy and matplotlib
import matplotlib
import seaborn
seaborn.set(style="dark")
rcParams['figure.figsize'] = (16, 4)
import pysptk
import pyworld
import librosa
import librosa.display
from IPython.display import Audio
def decompose(x, fs, period=5.0):
"""Decompose speech signal into f0, spectral envelope and aperiodicity using WORLD
"""
f0, timeaxis = pyworld.harvest(x, fs, frame_period=period, f0_floor=71.0, f0_ceil=800.0)
sp = pyworld.cheaptrick(x, f0, timeaxis, fs)
ap = pyworld.d4c(x, f0, timeaxis, fs)
return f0, timeaxis, sp, ap
def vis(sp_param, approx_sp, x, y, top_title):
"""Visualize compressed spectral parameter, recovered spectral envelope and reconstructed speech wavevform
"""
figure(figsize=(16,14))
subplots_adjust(hspace=0.4)
subplot(3,1,1)
librosa.display.specshow(sp_param.T, sr=fs, hop_length=hop_length, x_axis="time", y_axis="linear")
title(top_title)
subplot(3,1,2)
librosa.display.specshow(10*log(approx_sp).T, sr=fs, hop_length=hop_length, x_axis="time", y_axis="linear")
title("20log|H(w)|")
subplot(3,1,3)
plot(y, "r-+", label="reconstructed speech signal")
plot(x, label="original speech signal")
xlabel("sample")
legend(fontsize=14)
def mcepalpha(fs):
"""Determine frequency warping parameter from fs (simplest implementation)
See https://bitbucket.org/happyalu/mcep_alpha_calc/ for details.
"""
if fs == 16000:
return 0.41
elif fs == 22050:
return 0.455
elif fs == 44100:
return 0.544
elif fs == 48000:
return 0.554
else:
raise NotImplementedError
path = pysptk.util.example_audio_file() # from cmu_arctic
x, fs = librosa.load(path, sr=None)
x = x.astype(np.float64)
librosa.display.waveplot(x, sr=fs)
Audio(x, rate=fs)
frame_period = 5.0
hop_length = int(fs * frame_period * 0.001)
fftlen = pyworld.get_cheaptrick_fft_size(fs)
f0, timeaxis, sp, ap = decompose(x, fs, frame_period)
print(hop_length)
print(sp.shape)
print(fftlen)
80 (801, 513) 1024
plot(timeaxis, f0, linewidth=2, label="F0 contour estimated by Harvest")
xlabel("Time [sec]")
ylabel("Frequency [Hz]")
legend(fontsize=18);
librosa.display.specshow(10*log(sp).T, sr=fs, hop_length=hop_length, x_axis="time", y_axis="linear")
colorbar()
<matplotlib.colorbar.Colorbar at 0x7f9131cd92e8>
librosa.display.specshow(20*log(ap).T, sr=fs, hop_length=hop_length, x_axis="time", y_axis="linear")
colorbar()
<matplotlib.colorbar.Colorbar at 0x7f9130361550>
dim = 60
coded_sp = pyworld.code_spectral_envelope(sp, fs, dim)
decoded_sp = pyworld.decode_spectral_envelope(coded_sp, fs, fftlen)
y = pyworld.synthesize(f0, decoded_sp, ap, fs, frame_period)
nmse = np.linalg.norm(np.log(sp) - np.log(decoded_sp)) / np.linalg.norm(np.log(sp))
print("NMSE: {}".format(nmse))
vis(coded_sp, decoded_sp, x, y, "Mel-cepstrum")
Audio(y, rate=fs)
NMSE: 0.25376685579810726
dim = 60
alpha = mcepalpha(fs)
mc = np.apply_along_axis(pysptk.sp2mc, 1, sp, dim-1, alpha)
approximate_sp = np.apply_along_axis(pysptk.mc2sp, 1, mc, alpha, fftlen)
y = pyworld.synthesize(f0, approximate_sp, ap, fs, frame_period)
nmse = np.linalg.norm(np.log(sp) - np.log(approximate_sp)) / np.linalg.norm(np.log(sp))
print("NMSE: {}".format(nmse))
vis(mc, approximate_sp, x, y, "Coded spectral envelope by WORLD")
Audio(y, rate=fs)
NMSE: 0.027851414093334522