音声符号化について考えてみる」実践編その4です。
いよいよボコーダーの基礎まで来ました。
ピッチや基本周波数の抽出、フーリエ変換やケプストラム分析、音声再合成を行います。
…今回、呼び出している各種ライブラリが何をしているのか、自分自身あまり理解できていません。
ご指摘・ご指導、是非お待ちしております。

目的

ボコーダーの基本操作をマスターする。

環境

  • python 3.6.3
  • numpy 1.14.2
  • scipy 1.0.0
  • pysptk 0.1.10
  • matplotlib 2.1.1

コード解説

こことかこことか参考にしまくってます。ありがとうございます。

synthesis.py
from scipy.io import wavfile as io
from pysptk.sptk import swipe, rapt
from pysptk import mgc2sp, mgcep, excite, mgc2b
from pysptk.synthesis import  MGLSADF, Synthesizer
from scipy.fftpack import fft, fftfreq
import numpy as np
import matplotlib.pyplot as plt


def main():
    """
    メイン処理
    ----------
    """

    use_wav_path = "file/say.wav"
    n_frame = 2**16
    sr, wavdata = io.read(use_wav_path)
    wavdata = wavdata.astype(np.float64)
    wavdata = wavdata[:n_frame]

    # ピッチ抽出
    hopsize = int(sr / 100)     # 0.01秒間隔で分析
    min = 60.0                  # 最小周波数
    max = 4000.0                # 最大周波数
    pitch = process_pitch(wavdata, sr, hopsize, min, max)

    # 基本周波数抽出
    process_f0(wavdata, sr, hopsize, min, max)

    # 高速フーリエ変換(振幅スペクトル抽出)
    hammingWindow = np.hamming(n_frame)     # ハミング窓
    window_data = hammingWindow * wavdata
    process_fft(window_data, sr, n_frame)

    # ケプストラム分析
    ub_cep, lpc, mel_cep, mel_lpc = process_cepstrum(window_data, n_frame)

    plt.show()

    # 再合成
    process_synthesis(ub_cep, lpc, mel_cep, mel_lpc, hopsize, pitch, sr)



def process_pitch(wavdata, sr, hopsize, min, max):
    """
    ピッチ処理
    ----------
    :param wavdata: 音声データ
    :param sr:      サンプリング周波数
    :param hopsize: 抽出間隔
    :param min:     最小周波数
    :param max:     最大周波数
    """

    swipe_pitch = swipe(wavdata, fs=sr, hopsize=hopsize, min=min, max=max, otype="pitch")
    rapt_pitch = rapt(wavdata.astype(np.float32), fs=sr, hopsize=hopsize, min=min, max=max, otype="pitch")

    plt.subplot(2,2,1)
    plt.title("pitch")
    plt.plot(swipe_pitch, label="swipe")
    plt.plot(rapt_pitch, label="rapt")
    plt.legend(loc="lower left")
    plt.tight_layout()

    return swipe_pitch
    # return rapt_pitch


def process_f0(wavdata, sr, hopsize, min, max):
    """
    基本周波数処理
    -------------
    :param wavdata: 音声データ
    :param sr:      サンプリング周波数
    :param hopsize: 抽出間隔
    :param min:     最小周波数
    :param max:     最大周波数
    """

    swipe_f0 = swipe(wavdata, fs=sr, hopsize=hopsize, min=min, max=max, otype="f0")
    rapt_f0 = rapt(wavdata.astype(np.float32), fs=sr, hopsize=hopsize, min=min, max=max, otype="f0")

    plt.subplot(2,2,2)
    plt.title("f0")
    plt.plot(swipe_f0, label="swipe")
    plt.plot(rapt_f0, label="rapt")
    plt.legend(loc="lower left")
    plt.tight_layout()


def process_fft(window_data, sr, n_frame):
    """
    高速フーリエ変換処理
    -------------------
    :param window_data: 窓関数適用済み音声データ
    :param sr:          サンプリング周波数
    :param n_frame:     フレーム数
    :return:
    """

    fft_data = fft(window_data)
    amplitude_spectrum = [np.sqrt(c.real ** 2 + c.imag ** 2) for c in fft_data]
    freq = fftfreq(n_frame, 1/sr)

    plt.subplot(2,2,3)
    plt.title("amplitude spectrum")
    plt.plot(freq[1:int(n_frame/2)], np.abs(amplitude_spectrum[1:int(n_frame/2)]),
             linewidth=1.0)
    plt.tight_layout()


def process_cepstrum(window_data, n_frame):
    """
    ケプストラム処理
    ---------------
    :param window_data: 窓関数適用済み音声データ
    :param n_frame:     フレーム数
    :return:
    """

    # 対数スペクトル
    sp = 20 * np.log10(np.abs(np.fft.rfft(window_data)))

    # unbiased ケプストラム
    alpha = 0
    gamma = 0
    ub_cep = mgcep(window_data, alpha=alpha, gamma=gamma)
    env_ub_cep = mgc2sp(ub_cep, alpha, gamma, n_frame).real

    # LPC
    alpha = 0
    gamma = -1
    lpc = mgcep(window_data, alpha=alpha, gamma=gamma)
    env_lpc = mgc2sp(lpc, alpha, gamma, n_frame).real

    # mel ケプストラム
    alpha = 0.55
    gamma = 0
    mel_cep = mgcep(window_data, alpha=alpha, gamma=gamma)
    env_mel_cep = mgc2sp(mel_cep, alpha, gamma, n_frame).real

    # mel LPC
    alpha = 0.55
    gamma = -1
    mel_lpc = mgcep(window_data, alpha=alpha, gamma=gamma)
    env_mel_lpc = mgc2sp(mel_lpc, alpha, gamma, n_frame).real

    plt.subplot(2,2,4)
    plt.plot(sp, linewidth=0.1, label="Original log spectrum 20log|X(w)|")
    plt.plot(20.0/np.log(10)*env_ub_cep, linewidth=1.0, label="unbiased cepstrum")
    plt.plot(20.0/np.log(10)*env_lpc, linewidth=1.0, label="LPC")
    plt.plot(20.0/np.log(10)*env_mel_cep, linewidth=1.0, label="mel cepstrum")
    plt.plot(20.0/np.log(10)*env_mel_lpc, linewidth=1.0, label="mel LPC")
    plt.xlim(0, len(sp))
    plt.xlabel("frequency bin")
    plt.ylabel("log amplitude")
    plt.legend()

    return ub_cep, lpc, mel_cep, mel_lpc


def process_synthesis(ub_cep, lpc, mel_cep, mel_lpc, hopsize, pitch, sr):
    """
    音声再合成処理
    -------------
    :param ub_cep:  unbiased メルケプストラム
    :param lpc:     LPC ケプストラム
    :param mel_cep: メルケプストラム
    :param mel_lpc: メルLPCケプストラム
    :param hopsize: 分析間隔
    :param pitch:   ピッチ
    :param sr:      サンプリング周波数
    """

    # シンセサイザー作成
    ub_cep_synthesizer = Synthesizer(MGLSADF(alpha=0, stage=100), hopsize)
    lpc_synthesizer = Synthesizer(MGLSADF(alpha=0, stage=1), hopsize)
    mel_cep_synthesizer = Synthesizer(MGLSADF(alpha=0.55, stage=100), hopsize)
    mel_lpc_synthesizer = Synthesizer(MGLSADF(alpha=0.55, stage=1), hopsize)

    # フィルタ係数作成
    b_ub_cep = mgc2b(ub_cep, alpha=0, gamma=0)
    b_lpc = mgc2b(lpc, alpha=0, gamma=-1)
    b_mel_cep = mgc2b(mel_cep, alpha=0.55, gamma=0)
    b_mel_lpc = mgc2b(mel_lpc, alpha=0.55, gamma=-1)
    source_excitation = excite(pitch, hopsize)

    # 合成実行
    data_from_ub_cep = ub_cep_synthesizer.synthesis(source_excitation, np.array([b_ub_cep]))
    data_from_lpc = lpc_synthesizer.synthesis(source_excitation, np.array([b_lpc]))
    data_from_mel_cep = mel_cep_synthesizer.synthesis(source_excitation, np.array([b_mel_cep]))
    data_from_mel_lpc = mel_lpc_synthesizer.synthesis(source_excitation, np.array([b_mel_lpc]))

    # 音声書き出し&表示
    data_from_ub_cep = data_from_ub_cep.astype(np.int16)
    io.write("file/ubcep.wav", sr, data_from_ub_cep)
    plt.subplot(2, 2, 1)
    plt.plot(data_from_ub_cep)
    data_from_lpc = data_from_lpc.astype(np.int16)
    io.write("file/lpc.wav", sr, data_from_lpc)
    plt.subplot(2, 2, 2)
    plt.plot(data_from_lpc)
    data_from_mel_cep = data_from_mel_cep.astype(np.int16)
    io.write("file/melcep.wav", sr, data_from_mel_cep)
    plt.subplot(2, 2, 3)
    plt.plot(data_from_mel_cep)
    data_from_mel_lpc = data_from_mel_lpc.astype(np.int16)
    io.write("file/mellpc.wav", sr, data_from_mel_lpc)
    plt.subplot(2, 2, 4)
    plt.plot(data_from_mel_lpc)

    plt.tight_layout()
    plt.show()



if __name__ == '__main__':
    main()

最近、if __name__ == '__main__':の中は変数がグローバルになると知って、main()関数をファイルの先頭に作成してそこでメイン処理を書くようにしています。

n_frame = 2**16

フレーム数が2の累乗じゃないと、ケプストラム処理でエラーが出るので注意。

min = 60.0                  # 最小周波数
max = 4000.0                # 最大周波数

ピッチ抽出・基本周波数抽出に用いられる最小・最大周波数ですが、具体的にどう効いているのか理解できていません。
最大周波数を使用する音声データのナイキスト周波数未満にしておかないとエラーになるのですが…

hammingWindow = np.hamming(n_frame)     # ハミング窓

キレイにフーリエ変換するための窓関数はハミング窓を使用しています。
他にもいくつか種類がありますが、使い分けについてはよくわかりません!sorry!!

def process_pitch(wavdata, sr, hopsize, min, max):
    """
    ピッチ処理
    ----------
    :param wavdata: 音声データ
    :param sr:      サンプリング周波数
    :param hopsize: 抽出間隔
    :param min:     最小周波数
    :param max:     最大周波数
    """

    swipe_pitch = swipe(wavdata, fs=sr, hopsize=hopsize, min=min, max=max, otype="pitch")
    rapt_pitch = rapt(wavdata.astype(np.float32), fs=sr, hopsize=hopsize, min=min, max=max, otype="pitch")

ピッチ抽出処理です。
swipeとraptの2種類がありますが、こちらも使い分けについてはよくわかっていません。
とりあえずswipeの方が上手く抽出できる?的なことをどこかで見たような気がしたのですが…

def process_f0(wavdata, sr, hopsize, min, max):
    """
    基本周波数処理
    -------------
    :param wavdata: 音声データ
    :param sr:      サンプリング周波数
    :param hopsize: 抽出間隔
    :param min:     最小周波数
    :param max:     最大周波数
    """

    swipe_f0 = swipe(wavdata, fs=sr, hopsize=hopsize, min=min, max=max, otype="f0")
    rapt_f0 = rapt(wavdata.astype(np.float32), fs=sr, hopsize=hopsize, min=min, max=max, otype="f0")

基本周波数抽出処理です。
ピッチと基本周波数似た概念ですが、ピッチが心理量、基本周波数が物理量なようで。

def process_fft(window_data, sr, n_frame):
    """
    高速フーリエ変換処理
    -------------------
    :param window_data: 窓関数適用済み音声データ
    :param sr:          サンプリング周波数
    :param n_frame:     フレーム数
    :return:
    """

    fft_data = fft(window_data)
    amplitude_spectrum = [np.sqrt(c.real ** 2 + c.imag ** 2) for c in fft_data]
    freq = fftfreq(n_frame, 1/sr)

高速フーリエ変換処理です。
フーリエ変換を行うことで、時間軸の音声データを周波数軸の音声データに変換できます。
つまり、その音声データに含まれる周波数とその強さがわかる形になるんですね。

def process_cepstrum(window_data, n_frame):
    """
    ケプストラム処理
    ---------------
    :param window_data: 窓関数適用済み音声データ
    :param n_frame:     フレーム数
    :return:
    """

    # 対数スペクトル
    sp = 20 * np.log10(np.abs(np.fft.rfft(window_data)))

    # unbiased ケプストラム
    alpha = 0
    gamma = 0
    ub_cep = mgcep(window_data, alpha=alpha, gamma=gamma)
    env_ub_cep = mgc2sp(ub_cep, alpha, gamma, n_frame).real

    # LPC
    alpha = 0
    gamma = -1
    lpc = mgcep(window_data, alpha=alpha, gamma=gamma)
    env_lpc = mgc2sp(lpc, alpha, gamma, n_frame).real

    # mel ケプストラム
    alpha = 0.55
    gamma = 0
    mel_cep = mgcep(window_data, alpha=alpha, gamma=gamma)
    env_mel_cep = mgc2sp(mel_cep, alpha, gamma, n_frame).real

    # mel LPC
    alpha = 0.55
    gamma = -1
    mel_lpc = mgcep(window_data, alpha=alpha, gamma=gamma)
    env_mel_lpc = mgc2sp(mel_lpc, alpha, gamma, n_frame).real

ケプストラム処理です。
こ れ が 難 し い
alphaは-1~1の値を、gammaは-1~0の値を取ります。
alphaの値によってmel度合いが決定されるそうですが、一体いくらに設定すればよいのか…
gammaの値は-1に近いほどLPCという分析手法に近づくのですかね…?
そもそもLPCとそうでない手法の違いや、この実装で本当にLPC等の実装になっているのかなど、分からないことだらけ…

def process_synthesis(ub_cep, lpc, mel_cep, mel_lpc, hopsize, pitch, sr):
    """
    音声再合成処理
    -------------
    :param ub_cep:  unbiased メルケプストラム
    :param lpc:     LPC ケプストラム
    :param mel_cep: メルケプストラム
    :param mel_lpc: メルLPCケプストラム
    :param hopsize: 分析間隔
    :param pitch:   ピッチ
    :param sr:      サンプリング周波数
    """

    # シンセサイザー作成
    ub_cep_synthesizer = Synthesizer(MGLSADF(alpha=0, stage=100), hopsize)
    lpc_synthesizer = Synthesizer(MGLSADF(alpha=0, stage=1), hopsize)
    mel_cep_synthesizer = Synthesizer(MGLSADF(alpha=0.55, stage=100), hopsize)
    mel_lpc_synthesizer = Synthesizer(MGLSADF(alpha=0.55, stage=1), hopsize)

    # フィルタ係数作成
    b_ub_cep = mgc2b(ub_cep, alpha=0, gamma=0)
    b_lpc = mgc2b(lpc, alpha=0, gamma=-1)
    b_mel_cep = mgc2b(mel_cep, alpha=0.55, gamma=0)
    b_mel_lpc = mgc2b(mel_lpc, alpha=0.55, gamma=-1)
    source_excitation = excite(pitch, hopsize)

    # 合成実行
    data_from_ub_cep = ub_cep_synthesizer.synthesis(source_excitation, np.array([b_ub_cep]))
    data_from_lpc = lpc_synthesizer.synthesis(source_excitation, np.array([b_lpc]))
    data_from_mel_cep = mel_cep_synthesizer.synthesis(source_excitation, np.array([b_mel_cep]))
    data_from_mel_lpc = mel_lpc_synthesizer.synthesis(source_excitation, np.array([b_mel_lpc]))

大 問 題 が あ り ま す
実行してみると、正しく合成されていないみたいなんですよね…
どう正しくないのかは下記のグラフを見ていただくとして、なぜこうなるのか、不明です。
そのそもこの実装自体ほとんどここのを丸パクりしただけですので、どこで間違えたのか…
正しいシンセサイザーを用いて、正しい励振信号?と正しいフィルタ係数を合成すればいいはずなんですけどね…

というわけで、最後に各種値をプロットしてみましょう。
Figure_1.png
上図はケプストラム処理までの図。見た感じ正しく処理できていそうですね。

Figure_2.png
そしてこちらが問題の再合成後の音声波形…
4種法どの音声波形も、先頭付近のフレームだけ信号が立っていて、以降は信号なし。
Help Meィィ!

おわりに

いや~…音声信号処理って難しいですね!
やる夫で学ぶディジタル信号処理を読んで出直してきます…
そもそも文系出身の筆者にはここは難しかったので一度ならず何度か読んでは諦めているんですけども。

Sign up for free and join this conversation.
Sign Up
If you already have a Qiita account log in.