母音[a], [i], [u], [e], [o]とフォルマント周波数には次のような関係があることが知られている[1].

フォルマント周波数と母音の関係

・第一フォルマントは口腔の広さに比例する。
母音の第一フォルマントを比べると, 「あ」が一番高い. そして「い、う」が一番低く, 「え、お」が中間になる. 「あ」＝広母音「え、お」＝半狭母音「い、う」＝狭母音と呼ばれる. よって第一フォルマントは口腔の狭めに反比例する。
・第二フォルマントは舌が前に来る「い」は高くなり、奥にいく「お」は低くなる。

本記事では母音を録音し, ケプストラム解析によりフォルマント周波数を導出し実際にその傾向になることを確認する.

音声分析の流れ
前処理
プリエンファシス処理
ハニング窓
FFT解析(ケプストラム解析とは別枠)
ケプストラム解析

音声分析の流れ

本記事で行う音声処理の流れを下図に示す通りシンプルである.
しかしその前に音声データのサンプリング周波数やデータ数等を確認する.

前処理

まず,次のプログラムを用いて母音音声情報の基礎的情報を確認する.

for i in range(len(wav_files)):
    wav_file = wav_files[i]

    wav = wave.open(wav_file)
    num_samples = wav.getnframes() # 総フレーム数
    print("num_samples = ",num_samples)
    framerate = wav.getframerate()      # サンプリング周波数 Hz
    print("framerate = ",framerate)
    waveform = wav.readframes(num_samples)
    waveform = np.frombuffer(waveform,dtype = np.int16)
    print("len(waveform)",len(waveform))
    wav.close()
    time_axis = np.arange(num_samples)/sampling_frequency

    plt.figure(figsize = (8,1))
    plt.plot(time_axis,waveform)
    plt.xlim([0, num_samples / sampling_frequency])
    plt.xlabel('time, s')
    plt.grid()
    plt.ylabel('Amplitude')
    plt.show()

上からそれぞれ[a], [i], [u], [e], [o]の音声データである.

5つの音声すべて48khzであり, また[a]だけ録音された音声が短いことがわかる.

プリエンファシス処理

プリエンファシス処理は高周波成分を強調するために音声分析で用いられる.
音声分析で利用される周波数は数千Hz以上と高い。そのため1-1000hzくらいの低周波成分のゲインを弱めて, 相対的に高周波を強調する.

pre_emphasis = 0.97
for i in range(len(wav_files)):
    wav_file = wav_files[i]
    wav = wave.open(wav_file)
    num_samples = wav.getnframes()
    waveform = wav.readframes(num_samples)
    waveform = np.frombuffer(waveform,dtype = np.int16)
    wav.close()
    time_axis = np.arange(num_samples)/sampling_frequency
    emphasized_audio = np.append(waveform[0], waveform[1:] - pre_emphasis * waveform[:-1])

    plt.figure(figsize=(8,1))
    plt.plot(time_axis, waveform,label="Original waveform ")
    plt.plot(time_axis, emphasized_audio,alpha=0.6,label="Pre-emphasized waveform")
    plt.xlim([0, num_samples / sampling_frequency])
    plt.legend(fontsize=7,loc="upper left")
    plt.xlabel('time, s')
    plt.grid()
    plt.ylabel('Amplitude')
    # plt.savefig("result/compare_preem.pdf")
    plt.show()

周波数特性を確認するためのプログラムおよび結果を次に示す.
ここでFFTをしているが詳細は「FFT解析(ケプストラム解析とは別枠)」で述べる.

for i in range(len(wav_files)):
    wav_file = wav_files[i]
    wav = wave.open(wav_file)
    num_samples = wav.getnframes()
    waveform = wav.readframes(num_samples)
    waveform = np.frombuffer(waveform,dtype = np.int16)
    
    N_signal = len(waveform)
    emphasized_audio = np.append(waveform[0], waveform[1:] - pre_emphasis * waveform[:-1])
    fft_original = np.fft.fft(waveform)
    fft_emphasized = np.fft.fft(emphasized_audio)

    fft_original = abs(fft_original / (N_signal ))
    fft_emphasized = abs(fft_emphasized / (N_signal ))

    fft_original = 20*np.log10(fft_original+1e-10)
    fft_emphasized = 20*np.log10(fft_emphasized+1e-10)


    fft_freq = np.fft.fftfreq(N_signal, d=1/sampling_frequency)
    valid_idx = fft_freq > 0 #valid_idxは0より大きければTrue,0以下だったらFlseを返す。
    fft_freq = fft_freq[valid_idx]#周波数軸で０より大きい範囲で絞る
    Spectrum_original=fft_original[valid_idx]
    Spectrum_emphasized=fft_emphasized[valid_idx]

    plt.figure(figsize=(8,1))
    plt.plot(fft_freq, Spectrum_original, label='Original Spectrum')#常用
    plt.plot(fft_freq, Spectrum_emphasized, label='emphasized Spectrum')#常用
    plt.xlabel("Frequency, Hz")
    plt.xscale("log")
    plt.ylabel("spectrum, dB")
    plt.legend(fontsize=8,loc="upper left")
    plt.grid(True)
    plt.show()

ハニング窓

プリエンファシス処理した音声はハニング窓にかけられる.
これは, フーリエ変換の前準備であり, フーリエ変換を行うためには, 信号に周期性を持たせる必要があるためである.
そしてしゃべり始めやしゃべり終わりは定常特性とは別の特性を示すことが予想される.
そのためFFTの際には, 信号を切りその信号に対して窓をかける処理が行われる.
次に示すプログラムがそれを確認するためのものである.[a]だけ始まりが遅いため「extra」変数で後ろのほうに調整してある.中央から$\pm0.2$秒切り抜いている.

pre_emphasis = 0.97
window_sec = 0.2
extra = 5000

for i in range(len(wav_files)):
    wav_file = wav_files[i]
    wav = wave.open(wav_file)
    num_samples = wav.getnframes()
    waveform = wav.readframes(num_samples)
    waveform = np.frombuffer(waveform,dtype = np.int16)
    wav.close()
    time_axis = np.arange(num_samples)/sampling_frequency
    emphasized_audio = np.append(waveform[0], waveform[1:] - pre_emphasis * waveform[:-1])
    
    window_samples = int(window_sec*sampling_frequency)
    center_idx =num_samples//2
    # print("center_idx = ",center_idx)
    # print("window_samples = ",window_samples)
    start = max(center_idx-window_samples+extra,0)
    end = min(center_idx+window_samples+extra,num_samples)
    cuted_signal = emphasized_audio[start:end]
    # ==== cuted_signal の時間軸 ====
    cuted_length = len(cuted_signal)
    time_axis_cut = np.arange(start, end) / sampling_frequency
    


    plt.figure(figsize=(8, 1))
    plt.plot(time_axis,emphasized_audio, label="pre-emphasized waveform")
    plt.plot(time_axis_cut, cuted_signal,alpha = 0.6, label="extracted pre-emphasized waveform")
    plt.xlim([0, num_samples / sampling_frequency])
    plt.legend(fontsize=8,loc="upper left")
    plt.xlabel('time, s')
    plt.grid()
    plt.ylabel('Amplitude')
    # plt.savefig("result/compare_preem.pdf")
    plt.show()

結果を次に示す.
青色がプリエンファシス処理した信号で, オレンジが切り抜いた信号である.

切り抜いた信号は次のように, ハニング窓をかけられる.



for i in range(len(wav_files)):
    wav_file = wav_files[i]
    wav = wave.open(wav_file)
    num_samples = wav.getnframes()
    waveform = wav.readframes(num_samples)
    waveform = np.frombuffer(waveform,dtype = np.int16)
    wav.close()
    time_axis = np.arange(num_samples)/sampling_frequency
    emphasized_audio = np.append(waveform[0], waveform[1:] - pre_emphasis * waveform[:-1])

    window_samples = int(window_sec*sampling_frequency)
    center_idx =num_samples//2
    print("center_idx = ",center_idx)
    print("window_samples = ",window_samples)
    start = max(center_idx-window_samples+extra,0)
    end = min(center_idx+window_samples+extra,num_samples)
    cuted_signal = emphasized_audio[start:end]
    # ==== cuted_signal の時間軸 ====
    cuted_length = len(cuted_signal)
    time_axis_cut = np.arange(start, end) / sampling_frequency

    window = np.hanning(len(cuted_signal))  # 信号の長さに基づくハニング窓
    windowed_signal = cuted_signal * window

    


    plt.figure(figsize=(8,1))
    plt.plot(time_axis_cut, cuted_signal,alpha=0.6,label="extracted waveform")
    plt.plot(time_axis_cut, windowed_signal,label="windowed waveform ")
    plt.xlim([0, num_samples / sampling_frequency])
    plt.legend()
    plt.xlabel('time, s')
    plt.grid()
    plt.xlim(2,5)
    plt.ylabel('Amplitude')
    # plt.savefig("result/compare_preem.pdf")
    plt.show()

FFT解析(ケプストラム解析とは別枠)

FFTは用途によって見せ方を変える必要がある.
ここでは対数振幅スペクトルを示す.
対数振幅スペクトルは, 時系列データをを$x[n]$を表すと, 次のように表現される.
\begin{gather}
\mathcal F[x[n]]=X[k]\\
\rm{対数振幅スペクトル}=20\log_{10}(\frac{|X[k]|}{N})
\end{gather}
ここで$N$はFFT長, $k$は周波数ビンであり、$n$はサンプル数である。
$X[k]$は複素数であり,
\begin{gather}
X[k] = Re{X[k]}+jIm{X[k]}
\end{gather}であり,
振幅は
\begin{gather}
|X[k]| = \sqrt{Re(X[k])^2+Im(X[k])^2}
\end{gather}であらわされる.

for i in range(len(wav_files)):
    wav_file = wav_files[i]
    wav = wave.open(wav_file)
    num_samples = wav.getnframes()
    waveform = wav.readframes(num_samples)
    waveform = np.frombuffer(waveform,dtype = np.int16)
    wav.close()
    time_axis = np.arange(num_samples)/sampling_frequency
    
    emphasized_audio = np.append(waveform[0], waveform[1:] - pre_emphasis * waveform[:-1])
    window_samples = int(window_sec*sampling_frequency)
    center_idx =num_samples//2
    start = max(center_idx-window_samples+extra,0)
    end = min(center_idx+window_samples+extra,num_samples)
    cuted_signal = emphasized_audio[start:end]
    # ==== cuted_signal の時間軸 ====
    cuted_length = len(cuted_signal)
    time_axis_cut = np.arange(start, end) / sampling_frequency

    window = np.hanning(len(cuted_signal))  # 信号の長さに基づくハニング窓
    windowed_signal = cuted_signal * window

    fft_result = np.fft.fft(windowed_signal)
    N = len(windowed_signal)
    fft_freq = np.fft.fftfreq(N, d=1/sampling_frequency)
    Spectrum = abs(fft_result/(N))#振幅スペクトル
    # 対数振幅スペクトル
    Spectrum_db = 20 * np.log10(Spectrum + 1e-10)#常用対数

    


    plt.figure(figsize=(8,1))
    plt.plot(fft_freq[1:int(N/2)],Spectrum_db[1:int(N/2)])
    # plt.xlim([0, num_samples / sampling_frequency])
    # plt.legend()
    plt.xscale('log')
    plt.grid()
    plt.xlabel('Frequency [Hz]')
    plt.ylabel('Amplitude [dB]')
    plt.show()

上に示した対数振幅スペクトルから
・数百Hzあたりで山がある.
・数千Hzあたりで対数振幅スペクトルが下がる谷がある
ことがわかる.

ケプストラム解析

リフタリング

ケプストラムは、信号をフーリエ変換して得られる振幅スペクトルの対数（対数振幅スペクトル）を逆フーリエ変換することにより得られる。
(対数振幅スペクトルの2/Nまでと、2/Nからが対象であるため, IFFTしてもFFTしても同じ結果が得られる)

時系列データを$x[n]$とすると,
\begin{gather}
\rm{ケプストラム}=\mathcal F^{-1}[\log_{10}(\frac{|\mathcal F\{x[n]\}|}{N})]
\end{gather}
により得られる.

声道特性は低周波特性であるため, ここではケフレンシー＝2ms以内を取り出す.
そのため, 2ms以降を０と置き、マスクしている.
この処理はリフタリングと呼ばれる.

for i in range(len(wav_files)):
    wav_file = wav_files[i]
    wav = wave.open(wav_file)
    num_samples = wav.getnframes()
    waveform = wav.readframes(num_samples)
    waveform = np.frombuffer(waveform,dtype = np.int16)
    wav.close()
    time_axis = np.arange(num_samples)/sampling_frequency
    emphasized_audio = np.append(waveform[0], waveform[1:] - pre_emphasis * waveform[:-1])

    window_samples = int(window_sec*sampling_frequency)
    
    center_idx =num_samples//2
    # print("center_idx = ",center_idx)
    # print("window_samples = ",window_samples)
    start = max(center_idx-window_samples+extra,0)
    end = min(center_idx+window_samples+extra,num_samples)
    cuted_signal = emphasized_audio[start:end]
    # ==== cuted_signal の時間軸 ====
    cuted_length = len(cuted_signal)
    time_axis_cut = np.arange(start, end) / sampling_frequency

    window = np.hanning(len(cuted_signal))  # 信号の長さに基づくハニング窓
    windowed_signal = cuted_signal * window

    fft_result = np.fft.fft(windowed_signal)
    N = len(windowed_signal)
    fft_freq = np.fft.fftfreq(N, d=1/sampling_frequency)
    Spectrum = abs(fft_result/(N))#振幅スペクトル
    # 対数振幅スペクトル
    Spectrum_db = 20 * np.log10(Spectrum + 1e-10)#常用対数
    #ケプストラムを計算
    cepstrum = np.fft.ifft(Spectrum_db).real
    #ケフレンシー軸
    quefrency_axis = np.arange(len(cepstrum)) / sampling_frequency

    lifter_threshold = int(quefrency * sampling_frequency) 
    threshold_quefrency = lifter_threshold / sampling_frequency
  
    
    
    plt.figure(figsize=(8,1))
    plt.grid()
    # plt.plot(fft_freq[1:int(N/2)],Spectrum_db[1:int(N/2)])
    # plt.xlim([0, num_samples / sampling_frequency])
    # plt.legend()
    # plt.xscale('log')
    plt.plot(quefrency_axis[:len(cepstrum) // 2], cepstrum[:len(cepstrum) // 2])
    plt.axvline(
    x=threshold_quefrency,
    color='red',
    linestyle='--',
    linewidth=2,
    label=f'Lifter threshold = {threshold_quefrency:.4f} s')
    plt.xlabel('Quefrency, s')
    plt.ylabel('Amplitude')
    plt.show()

スペクトル包絡

リフタリングしたケプストラムを再度フーリエ変換することで周波数領域に戻す.
こうすることで音声に含まれる低周波成分だけを残ることができる.




for i in range(len(wav_files)):
    wav_file = wav_files[i]
    wav = wave.open(wav_file)
    num_samples = wav.getnframes()
    waveform = wav.readframes(num_samples)
    waveform = np.frombuffer(waveform,dtype = np.int16)
    wav.close()
    time_axis = np.arange(num_samples)/sampling_frequency
    emphasized_audio = np.append(waveform[0], waveform[1:] - pre_emphasis * waveform[:-1])

    window_samples = int(window_sec*sampling_frequency)
    center_idx =num_samples//2
    start = max(center_idx-window_samples+extra,0)
    end = min(center_idx+window_samples+extra,num_samples)
    cuted_signal = emphasized_audio[start:end]
    cuted_length = len(cuted_signal)
    time_axis_cut = np.arange(start, end) / sampling_frequency

    window = np.hanning(len(cuted_signal))  # 信号の長さに基づくハニング窓
    windowed_signal = cuted_signal * window

    fft_result = np.fft.fft(windowed_signal)
    N = len(windowed_signal)
    fft_freq = np.fft.fftfreq(N, d=1/sampling_frequency)
    Spectrum = abs(fft_result/(N))#振幅スペクトル
    # 対数振幅スペクトル
    Spectrum_db = 20 * np.log10(Spectrum + 1e-10)#常用対数
    #ケプストラムを計算
    cepstrum = np.fft.ifft(Spectrum_db).real
    #ケフレンシー軸
    quefrency_axis = np.arange(len(cepstrum)) / sampling_frequency
    lifter_threshold = int(quefrency * sampling_frequency) 
    """
    #ケプストラムの表示
    plt.figure(figsize=(10, 4))
    plt.grid()
    threshold_quefrency = lifter_threshold / sampling_frequency
    plt.plot(quefrency_axis[:len(cepstrum) // 2], cepstrum[:len(cepstrum) // 2])
    plt.axvline(
    x=threshold_quefrency,
    color='red',
    linestyle='--',
    linewidth=2,
    label=f'Lifter threshold = {threshold_quefrency:.4f} s')
    plt.xlabel('Quefrency, s')
    plt.ylabel('Cepstrum')
    plt.show()
    """
    
     # 低ケフレンシー成分の抽出
    low_quefrency_cepstrum = np.copy(cepstrum)
    low_quefrency_cepstrum[lifter_threshold:] = 0

    # スペクトル包絡の計算
    envelope_log_spectrum = np.fft.fft(low_quefrency_cepstrum).real#a+jbのaを取り出す
    
    valid_idx = fft_freq > 0 #valid_idxは0より大きければTrue,0以下だったらFlseを返す。
    fft_freq = fft_freq[valid_idx]#周波数軸で０より大きい範囲で絞る
    
    Spectrum_db=Spectrum_db[valid_idx]
    envelope_log_spectrum=envelope_log_spectrum[valid_idx]

    # グラフ表示
    plt.figure(figsize=(8, 1))
    plt.plot(fft_freq, Spectrum_db, label='Original Spectrum')#常用
    plt.plot(fft_freq,envelope_log_spectrum, label='Spectral Envelope', linewidth=2)#常用対数
    plt.legend()
    plt.xscale("log")
    plt.xlabel("Frequency, Hz")
    plt.ylabel("spectrum, dB")
    plt.grid(True)
    plt.show()

フォルマント周波数

前節で述べたスペクトル包絡からピークとなる周波数(フォルマント周波数)を探す.


from scipy.signal import find_peaks


for i in range(len(wav_files)):
    wav_file = wav_files[i]
    wav = wave.open(wav_file)
    num_samples = wav.getnframes()
    waveform = wav.readframes(num_samples)
    waveform = np.frombuffer(waveform,dtype = np.int16)
    wav.close()
    time_axis = np.arange(num_samples)/sampling_frequency
    emphasized_audio = np.append(waveform[0], waveform[1:] - pre_emphasis  * waveform[:-1])


    window_samples = int(window_sec*sampling_frequency)
    
    center_idx =num_samples//2
    start = max(center_idx-window_samples+extra,0)
    end = min(center_idx+window_samples+extra,num_samples)
    cuted_signal = emphasized_audio[start:end]
    # ==== cuted_signal の時間軸 ====
    cuted_length = len(cuted_signal)
    time_axis_cut = np.arange(start, end) / sampling_frequency

    window = np.hanning(len(cuted_signal))  # 信号の長さに基づくハニング窓
    windowed_signal = cuted_signal * window

    fft_result = np.fft.fft(windowed_signal)
    N = len(windowed_signal)
    fft_freq = np.fft.fftfreq(N, d=1/sampling_frequency)
    Spectrum = abs(fft_result/(N))#振幅スペクトル
    # 対数振幅スペクトル
    Spectrum_db = 20 * np.log10(Spectrum + 1e-10)#常用対数
    #ケプストラムを計算
    cepstrum = np.fft.ifft(Spectrum_db).real
    #ケフレンシー軸
    quefrency_axis = np.arange(len(cepstrum)) / sampling_frequency
    lifter_threshold = int(quefrency * sampling_frequency)
    threshold_quefrency = lifter_threshold / sampling_frequency
     # 低ケフレンシー成分の抽出
    low_quefrency_cepstrum = np.copy(cepstrum)
    low_quefrency_cepstrum[lifter_threshold:] = 0
    # スペクトル包絡の計算
    envelope_log_spectrum = np.fft.fft(low_quefrency_cepstrum).real#a+jbのaを取り出す    
    valid_idx = fft_freq > 0 #valid_idxは0より大きければTrue,0以下だったらFlseを返す。
    fft_freq = fft_freq[valid_idx]#周波数軸で０より大きい範囲で絞る
    
    Spectrum_db=Spectrum_db[valid_idx]
    envelope_log_spectrum=envelope_log_spectrum[valid_idx]
    
    # F1,F2を探す周波数範囲
    f1_min = 0
    f1_max = 4000
    
    f1_range = (fft_freq >= f1_min) & (fft_freq <= f1_max)
   
    freq_f1 = fft_freq[f1_range]
    env_f1 = envelope_log_spectrum[f1_range]
    
    # ピーク検出
    peaks, properties = find_peaks(
        env_f1,
        prominence=1.0,   # ピークらしさ。必要に応じて調整
        distance=20       # 近すぎるピークを避ける
    )
    # print("The number of peaks is",len(peaks))
    
    if len(peaks) == 0:
        print("F1 candidate was not found")
        F1 = np.nan
    else:
        # 低周波側から最初のピークをF1とする
        first_peak_idx = peaks[0]
        second_peak_idx = peaks[1]
        F1 = freq_f1[first_peak_idx]
        F2 = freq_f1[second_peak_idx]
        # print(f"F1 = {F1:.1f} Hz")
        # print(f"F2 = {F2:.1f} Hz")
    
        # 確認プロット
        plt.figure(figsize=(8, 1))
        plt.plot(fft_freq, Spectrum_db, label='Original Spectrum')
        plt.plot(fft_freq, envelope_log_spectrum, label='Spectral Envelope', linewidth=2)
        plt.scatter(F1, env_f1[first_peak_idx], color='red', zorder=5, label=f'F1 = {F1:.1f} Hz')
        plt.scatter(F2, env_f1[second_peak_idx], color='red', zorder=5, label=f'F2 = {F2:.1f} Hz')
        plt.xscale("log")
        plt.xlabel("Frequency, Hz")
        plt.ylabel("Spectrum, dB")
        plt.grid(True)
        plt.legend(fontsize=7,loc="upper right")
        plt.show()

フォルマント周波数の結果をまとめると次のようになる.

	F1	F2
/a/	1045	2760
/i/	327.5	1150
/u/	320	1307.5
/e/	407.5	2067.5
/o/	460	1405

冒頭で述べたフォルマント周波数と母音の関係は次のようなものだった.

フォルマント周波数と母音の関係

この結果から, 第一フォルマントは上記のような傾向が確認できた.
しかし第二フォルマントにその傾向はみられなかった.
考えられる原因は
適切なケフレンシーがよくわからなかったことがあげられる.
低ケフレンシっていったいどれくらいのことなのか.

また実はこのデータは２回目の収録である.１回目に収録したときは第一フォルマントの傾向さえ確認できなかったため, もっとはっきりと母音を出して録音してみたのが今回のデータである.
そのため, やっていて再現性あるのかは不安になった.
もっとこうしたほうがいいという案があれば是非いただきたいです.

ちなみにiphoneで収録したためMP3ファイルによって高周波帯域がカットされているためこのような結果になったのではないかと思ったが、よく調べてみるとMP3がカットする周波数は16khz付近とかなり高い. フォルマントはせいぜい数千Hzと予想されるため, その説はなさそうである.

参考文献

[1]ビジュアル音声学,川原繁人, 三省堂, 12018

母音のフォルマント解析