介紹
從本文開始，將會開始一系列的語音特徵方面的介紹。第一個介紹的語音特徵就是基音週期，當發濁音的時候，此時會引起聲帶的振動，而這個振動會呈現一定的週期性，即基音週期，基音週期的倒數叫基頻(fundamental frequency,f0)。簡言之，就是去求一個非完全周期函式的近似週期。基音週期用途非常廣泛，可以用來檢測語音噪聲、特殊聲音檢測、男女判別，說話人識別等。針對於基音週期估計系列，主要介紹兩個經典的在時域估計基音週期的方法，Yin以及pYin。本文著重介紹Yin，下一篇介紹pYin。
Yin
時域內的基音週期估計方法主要是去找語音波形的最小正週期，主要是藉助於自相關函式。Yin方法求解基音週期主要有如下幾個步驟，
首先是計算自相關函式(ACF),Yin採用差函式的形式
d_{t}(\tau)=\sum_{j=t}^{t+W-1}(x_{j}-x_{j+\tau})^2
其中W表示窗的點數。對於差函式而言，峰谷所對應的時延\tau可以表示週期。做減法的優點在於：如果訊號幅度在逐漸增大，那麼訊號往後平移的越多，與原訊號相乘的時候結果就會越大。這樣就會導致所估計得到的基音週期過大，對應的基頻過小。而減法則避免了這個問題。但是基於減法得ACF同樣存在一個問題，就是總會在\tau = 0處取得最小值，將會對估計基音週期產生干擾。因此Yin方法提出了累積均值歸一化差值函式，將時延值\tau = 0這種情況單獨進行了考慮。
但是累積均值歸一化差值函式也會存在一個問題，可能出現下圖所示的情況
即峰谷不僅僅出現在上圖中的最左邊，並且右邊的峰谷值比最左邊的峰谷值還要小，此時如果直接求累積均值歸一化差值函式的最小值，那麼很有可能挑選到右邊的峰谷，從而估計的基音週期比原基音週期大兩倍，即產生了半頻錯誤。所以就引出了Yin論文中的Step 4, 通過設立一個閾值，如果第一個峰谷值比閾值還要大，那麼就選擇該峰谷值對應的時延作為基音週期。如果沒有峰谷值比閾值更大，此時有兩種處理方法，要麼選相對更小的峰谷要麼認為該幀不存在基音週期。最後為了從累積均值歸一化差值函式中得到準確的時延，可以先進行插值，然後再搜尋得到峰谷對應的時延值。
程式碼實現
程式碼主要是參考librosa開源框架庫.
import soundfile as sf
import numpy as np
from numpy.lib.stride_tricks import as_strided
def speech_frame(x, frame_length, hop_length, axis=-1):
x = np.ascontiguousarray(x)
n_frames = 1 + (x.shape[axis] - frame_length) // hop_length
strides = np.asarray(x.strides)
new_stride = np.prod(strides[strides > 0] // x.itemsize) * x.itemsize
if axis == -1:
shape = list(x.shape)[:-1] + [frame_length, n_frames]
strides = list(strides) + [hop_length * new_stride]
elif axis == 0:
shape = [n_frames, frame_length] + list(x.shape)[1:]
strides = [hop_length * new_stride] + list(strides)
return as_strided(x, shape=shape, strides=strides)
def cumulative_mean_normalized_difference(y_frames, frame_length, win_length, min_period, max_period):
# 自相關函式和功率譜互為傅立葉變換對.該種方式計算自相關函式可以節省計算量.
a = np.fft.rfft(y_frames, frame_length, axis=0)
# 疑問點1:y_frames[win_length:: -1, :],不知道為啥這樣取.
b = np.fft.rfft(y_frames[win_length:: -1, :], frame_length, axis=0)
acf_frames = np.fft.irfft(a * b, frame_length, axis=0)[win_length:]
acf_frames[np.abs(acf_frames) < 1e-6] = 0
energy_frames = np.cumsum(y_frames ** 2, axis=0)
#疑問點2: 不知道為啥這樣就能表示r_{t+\tau}
energy_frames = energy_frames[win_length:, :] - energy_frames[:-win_length, :]
energy_frames[np.abs(energy_frames) < 1e-6] = 0
y_yin_frames = energy_frames[0, :] + energy_frames - 2 * acf_frames
yin_numerator = y_yin_frames[min_period: max_period + 1, :]
tau_range = np.arange(1, max_period + 1)[:, None]
cumulative_mean = np.cumsum(y_yin_frames[1 : max_period + 1, :], axis=0) / tau_range
yin_denominator = cumulative_mean[min_period - 1: max_period, :]
yin_frames = yin_numerator / (yin_denominator + tiny(yin_denominator))
return yin_frames
def parabolic_interpolation(y_frames):
parabolic_shifts = np.zeros_like(y_frames)
parabola_a = (y_frames[:-2, :] + y_frames[2:, :] - 2 * y_frames[1:-1, :]) / 2
parabola_b = (y_frames[2:, :] - y_frames[:-2, :]) / 2
parabolic_shifts[1:-1, :] = -parabola_b / (2 * parabola_a + tiny(parabola_a))
parabolic_shifts[np.abs(parabolic_shifts) > 1] = 0
return parabolic_shifts
def localmin(x,axis = 0):
paddings = [(0, 0)] * x.ndim
paddings[axis] = (1, 1)
x_pad = np.pad(x, paddings, mode="edge")
inds1 = [slice(None)] * x.ndim
inds1[axis] = slice(0, -2)
inds2 = [slice(None)] * x.ndim
inds2[axis] = slice(2, x_pad.shape[axis])
return (x < x_pad[tuple(inds1)]) & (x <= x_pad[tuple(inds2)])
def tiny(x):
x = np.asarray(x)
if np.issubdtype(x.dtype, np.floating) or np.issubdtype(
x.dtype, np.complexfloating
):
dtype = x.dtype
else:
dtype = np.float32
return np.finfo(dtype).tiny
y, sr = sf.read('q1.wav')
frame_length = 2048
win_length = frame_length // 2
hop_length = frame_length // 4
y_frame = speech_frame(y, frame_length=frame_length, hop_length=hop_length)
fmax = 2000
fmin = 100
trough_threshold = 0.2
min_period = max(int(np.floor(sr / fmax)), 1)
max_period = min(int(np.ceil(sr / fmin)), frame_length - win_length - 1)
yin_frames = cumulative_mean_normalized_difference(y_frame, frame_length,
win_length, min_period, max_period)
parabolic_shifts = parabolic_interpolation(yin_frames)
is_tough = localmin(yin_frames, axis=0)
is_tough[0, :] = yin_frames[0, :] < yin_frames[1, :]
is_threshold_trough = np.logical_and(is_tough, yin_frames < trough_threshold)
global_min = np.argmin(yin_frames, axis=0)
yin_period = np.argmax(is_threshold_trough, axis=0)
no_trough_below_threshold = np.all(~is_threshold_trough, axis=0)
yin_period[no_trough_below_threshold] = global_min[no_trough_below_threshold]
yin_period = (
min_period
+ yin_period
+ parabolic_shifts[yin_period, range(yin_frames.shape[1])]
)
f0 = sr / yin_period
print(f0)其中q1.wav的波形圖,語譜圖以及頻率響應如下圖
在紅線處,其基頻大約在237hz,和估計的基本吻合.
本作品採用《CC 協議》，轉載必須註明作者和本文連結