Pythonで正弦波モデルを実装

Pythonを使って正弦波モデルを実装しました。ソースコードがネット上にあまりないため、頑張って自力で書きました。 今回は、M&Qアルゴリズムとよばれる方法で実装を行っています。

参考文献は以下の2つです。

[1] 小坂 直敏、「サウンドエフェクトのプログラミング―Cによる音の加工と音源合成」、オーム社、2012.

[2] R.J.McAulay and T.F.Quatieri , “Speech Analysis/Synthesis Based on a Sinusoidal Representation”, IEEE Trans. on Acoust., Speech, and Singnal Processing, vol. ASSP-34, No. 4, Aug. 1986.

正弦波モデルとは

正弦波モデルとは、正弦波に似ている信号を足し合わせて、音声や楽音などを表現する方法です。

音声のスペクトログラムの例を図1に示します。

図1:音声のスペクトログラムの例
図1:音声のスペクトログラムの例

音声のスペクトログラムでは、周波数方向に間隔をあけてエネルギーが多い箇所があり、こういう特徴を調波構造と呼んだりします。正弦波モデルでは、こういうエネルギーが多い箇所を振幅と周波数が滑らかに時間変化する正弦波で近似します。

正弦波モデルの用途としては、元信号よりも少ないデータで信号を近似するので、音声圧縮で使われたりします(MP3とかはこの方法を使いませんが...)。また、正弦波モデルによって調波成分ごとに操作が可能となるため、信号処理の前段階としても使われている印象です。

プログラム

記述したプログラムについて説明します。正直いって、論文の内容で見逃している点やバグがあるかもしれないので、参考程度にしていただければ幸いです。

プログラムは一つのファイルに記述しています。また、モノラル音源のみに対応しており、ステレオには対応していません。

メイン処理

メインの処理プログラム箇所は以下です。

import soundfile as sf
import numpy as np
import scipy.signal as sg
import cmath
import datetime
from librosa import resample, stft
import matplotlib.pyplot as plt

PEAK        = 1
CONNECT     = 2
DEATH       = 3
BIRTH       = 4
DEATH_BIRTH = 5

# フレームごとの極大値を検知
def DetectPeak(X):
        ~ 略 ~

# フレーム間ピークマッチング
def PeakMatching(peak_plane):
        ~ 略 ~

# 音声合成(位相補間と振幅補間)
def SynthesizeSpeech(X, peak_plane, track_plane):
        ~ 略 ~

# パラメータ
path = "voice.wav"
sr_o = 10000
n_fft = 512
hop_length = 256
delta_freq = 50
window = sg.windows.get_window('hamming', n_fft)
window = window / np.sum(window)
    
# WAVファイル読み込みとリサンプリング
x, sr = sf.read(path)
x = resample(x, sr, sr_o)
sf.write("sample.wav", x, sr_o, subtype='PCM_16')

# 必要な値をあらかじめ計算
delta_bin  = int(delta_freq / (sr_o/n_fft))
T = hop_length/sr_o
omega_b = 2.0 * np.pi *(sr_o/n_fft)
x_shape = x.shape

# 計測スタート
start_time = datetime.datetime.now()

# 短時間フーリエ変換
X = stft(x, n_fft, hop_length, None, window)

# フレームごとの極大値を検知
peak_plane = DetectPeak(X)

# フレーム間ピークマッチング
track_plane = PeakMatching(peak_plane)

# 音声合成(位相補間と振幅補間)
y = SynthesizeSpeech(X, peak_plane, track_plane)

# 処理時間の出力
elapsed_time = datetime.datetime.now() - start_time
print("time: %s" % (elapsed_time))

# WAVファイル出力
sf.write("soundout.wav", y, sr_o, subtype='PCM_16')

メインの処理では、以下の手順で音声を再合成しています。

  1. 音声を読み込む
  2. リサンプリングする
  3. 短時間フーリエ変換をする
  4. フレームごとに極大値(ピーク)を検知
  5. フレームどうしのピークをつなげる
  6. フレーム間の位相補間と振幅補間
  7. 音声を出力する

注意:窓関数  w[n] については以下の式を満たすように補正しています。

\displaystyle{
\hspace{3em} \sum_{n=-N/2}^{N/2} w[n] = 1 
}

ここで、 N は窓関数の幅です。

なぜ、このようなことをするかというと、補正していない窓関数を使うと短時間フーリエ変換後の各要素のパワーと元信号のパワーを単純比較できなくなるからです。

補正に関しては以下の方のブログがわかりやすいです。

フレームごとに極大値の検知

フレームごとに極大値の検知をする DetectPeak の関数定義は以下になります。

# フレームごとの極大値を検知
def DetectPeak(X):

    # 全てゼロのnumpy配列の作成
    peak_plane = np.zeros(X.shape, dtype=int)
    # 極大値の箇所を検知
    index = sg.argrelextrema(X, np.greater, axis=0, order=1)
    # 極大値の箇所に1を代入
    peak_plane[index] = PEAK

    return peak_plane

極大値がある箇所にPEAK(1)を代入して、極大値ではない箇所は0とする2次元配列 peak_plane を作成します。

DetectPeak の出力は図2のようになります。

図2:DetectPeak の出力
図2:DetectPeak の出力

補足Pythonでfor文を使うと、基本的に処理が遅いので、scipy.signal にある argrelextrema を使用して、極大値を検知しました。

フレーム間ピークマッチング

フレーム間ピークマッチングする PeakMatching の関数定義は以下になります。

# フレーム間ピークマッチング
def PeakMatching(peak_plane):

    n_bin = peak_plane.shape[0]    # ビン数
    n_frame = peak_plane.shape[1]  # フレーム数
    link_flag = False              # 接続先があるかのフラグ 
    candidate_match = n_bin        # 接続先の候補1
    next_candidate_match = n_bin   # 接続先の候補2

    # 全てゼロのnumpy配列の作成
    track_plane = np.ones(peak_plane.shape, dtype=int) * (-1)

    for i in range(n_frame-1):
        for k in range(n_bin):

            # 極大値の箇所を検知
            if peak_plane[k,i]==PEAK or peak_plane[k,i]==CONNECT:

                # 前の極大値で接続先があったとき
                if link_flag == True:
                    # step2(c)
                    if abs(pre_peak-candidate_match) < abs(k-candidate_match):
                        track_plane[pre_peak,i] = candidate_match
                        peak_plane[candidate_match, i+1] = CONNECT
                    else:
                        # step2(d)
                        if next_candidate_match > candidate_match:
                            track_plane[pre_peak,i] = pre_peak     
                            peak_plane[pre_peak,i+1] = DEATH
                        # step2(e)
                        else:
                            track_plane[pre_peak,i] = next_candidate_match
                            peak_plane[next_candidate_match, i+1] = CONNECT

                # 初期化
                link_flag = False
                candidate_match = n_bin
                next_candidate_match = n_bin

                # 接続先候補を探す
                if k < delta_bin:
                    for j in range(-k,delta_bin+1):
                        if peak_plane[k+j,i+1]==PEAK:
                            link_flag = True
                            if abs(j) < abs(candidate_match):
                                next_candidate_match = candidate_match 
                                candidate_match = k+j  #step1(b)
                elif k >= n_bin-delta_bin: 
                    for j in range(-delta_bin,n_bin-k):
                        if peak_plane[k+j,i+1]==PEAK:
                            link_flag = True
                            if abs(j) < abs(candidate_match):
                                next_candidate_match = candidate_match 
                                candidate_match = k+j  #step1(b)
                else:
                    for j in range(-delta_bin,delta_bin+1):
                        if peak_plane[k+j,i+1]==PEAK:
                            link_flag = True  
                            if abs(j) < abs(candidate_match):
                                next_candidate_match = candidate_match 
                                candidate_match = k+j  #step1(b)
                
                # 接続先がなかった場合
                if link_flag==False:
                    track_plane[k,i]  = k 
                    peak_plane[k,i+1] = DEATH

                # 極大値の位置の保存
                pre_peak = k
    
            # trackがDEATHのとき
            elif peak_plane[k,i]==DEATH:
                track_plane[k,i] = (-2)  # step1(a) DEATHとする

        # 一番高い周波数ビンのピークを接続させる
        if link_flag == True:
            track_plane[pre_peak,i] = candidate_match
            peak_plane[candidate_match, i+1] = CONNECT

        # i+1フレームに接続されてない極大値がないか確認
        for k in range(n_bin):
            # step3(f) BIRTHとする
            if peak_plane[k,i+1]==PEAK:
                track_plane[k,i] = k
                peak_plane[k,i+1] = CONNECT      
                if peak_plane[k,i]==DEATH:
                    peak_plane[k,i] = DEATH_BIRTH
                else:
                    peak_plane[k,i] = BIRTH

    return track_plane

PeakMatching は図3のようにピークの接続先のビン番号を記録した2次元配列 track_planeを出力します。ピークではない箇所は (-1) となっており、接続先がないときは (-2) とします。

図3:PeakMatching の出力の模式図
図3:PeakMatching の出力の模式図

ピークの接続は以下の手順で決めます。

STEP1

図4:STEP1のピークマッチングの模式図
図4:STEP1のピークマッチングの模式図

まず、パラメータとして Δ(delta_bin) を決めておきます。フレーム i、ビン番号 kにピークがある場合、次のフレームの [k-Δ, k+Δ] の範囲にピークがなければ、接続先はないものとします。接続先がないときは、図4(a) のように track_plane[k,i]=ktrack_plane[k,i+1]=(-2) を代入します。

一方、 [k-Δ, k+Δ] の範囲にピークがあれば、接続先の第一候補のビン番号(candidate_match)と第二候補のビン番号(next_candidate_match)を保存しておいて、STEP2 に移行します。

ピークの接続先がなくなった箇所(track_plane[k,i+1]=(-2)を代入した箇所)については、peak_plane[k,i+1]=DEATH(3)を代入します。

STEP2

図5:STEP2のピークマッチングの模式図
図5:STEP2のピークマッチングの模式図

接続先候補のピークがある場合、次のピークとkビンのピークのどちらが接続先第一候補(candidate_match)と近いか比較します。

kビンの方が近い場合(σ1>σ2)、track_plane[k,i]=candidate_match を代入します(図5(c))。

次のビンの方が近くて(σ1<σ2)、candidate_match<next_candidate_matchの場合、図5(d)のように track_plane[k,i]=ktrack_plane[k,i+1]=(-2) を代入します。

次のビンの方が近くて(σ1<σ2)、candidate_match>next_candidate_matchの場合、図5(e)のように track_plane[k,i]=next_candidate_match を代入します。

前のフレームのピークと接続されたピークの場所については、
peak_plane[k,i]=CONNECT(2) を代入します。
また、STEP1と同様にtrack_plane[k,i+1]=(-2)を代入した箇所については、
peak_plane[k,i+1]=DEATH(3)を代入します。

STEP3

図6:STEP3のピークマッチングの模式図
図6:STEP3のピークマッチングの模式図

フレーム i の全ピークの接続先が決まったら、フレーム(i+1) で接続されていないピークがあれば、track_plane[k,i]=kpeak_plane[k,i]=BIRTH(4) を代入します(図6(f))。

また、既にpeak_plane[k,i]=DEATH(2)となっていた場合、peak_plane[k,i]=DEATH_BIRTH(5)とします。

位相補間と振幅補間

位相補間と振幅補間をする SynthesizeSpeech の関数定義は以下になります。

# 音声合成(位相補間と振幅補間)
def SynthesizeSpeech(X, peak_plane, track_plane):

    n_bin = X.shape[0]                # ビン数
    n_frame = X.shape[1]              # フレーム数
    y = np.zeros(x_shape)             # 出力信号
    n = np.arange(hop_length)         # サンプル番号
    t = np.arange(hop_length)/sr_o    # 時間
    r_arr = np.zeros(hop_length)      # 振幅の一時保存
    theta_arr = np.zeros(hop_length)  # 位相の一時保存

    for i in range(n_frame-1):
        for k in range(n_bin):
            
            # 極大値、DEATH、BIRTHを判定
            if peak_plane[k,i] == CONNECT or \
               peak_plane[k,i] == BIRTH or \
               peak_plane[k,i] == DEATH_BIRTH:
                
                # 振幅と位相と角周波数
                if peak_plane[k,i] == BIRTH or peak_plane[k,i] == DEATH_BIRTH:
                    r1 = 0
                else:
                    if k==0 or k==n_bin-1:
                        r1 = abs(X[k,i])
                    else:
                        r1 = 2.0 * abs(X[k,i])
                theta1 = cmath.phase(X[k,i])
                omega1 = k * omega_b
                if peak_plane[k,i+1] == DEATH or peak_plane[k,i+1] == DEATH_BIRTH:
                    r2 = 0
                else:
                    if track_plane[k,i]==0 or track_plane[k,i]==n_bin-1:
                        r2 = abs(X[track_plane[k,i],i+1])
                    else:
                        r2 = 2.0 * abs(X[track_plane[k,i],i+1])
                theta2 = cmath.phase(X[track_plane[k,i],i+1])
                omega2 = track_plane[k,i] * omega_b
                
                # Mの計算
                M = ((theta1+omega1*T-theta2)+(omega2-omega1)*T/2.0)/(2.0*np.pi)
                M = round(M)

                # alphaとbeta
                matrixup = theta2-theta1-omega1*T+M*2.0*np.pi
                matrixdw = omega2-omega1
                alpha = matrixup*(3.0/(T**2))-matrixdw/T
                beta  = matrixup*(-2.0/(T**3))+matrixdw/(T**2)

                # r, theta計算とyに足し合わせる
                r_arr = r1+(r2-r1)*n/(hop_length)
                theta_arr = theta1+omega1*t+alpha*(t**2)+beta*(t**3)
                y[i*hop_length:(i+1)*hop_length] += r_arr*np.cos(theta_arr)

フレーム間の位相補間と振幅補間を行い、音声を再合成します。

フレーム i、ビン番号 k にピークがあるとき、接続先のビン番号を k' とすると、フレームごとの振幅と位相と周波数を以下のように計算します。

\displaystyle{
 r_1 = 2 |X[k,i]| \\
 r_2 = 2 |X[k',i]| \\
 \theta_1 = {\rm arg}(X[k,i])\\
 \theta_2 =  {\rm arg}(X[k',i])\\
 \omega_1 = 2 \pi k \frac{f_s}{N_F}\\
 \omega_2 = 2 \pi k' \frac{f_s}{N_F}
}

ここで、 f_s はサンプリング周波数、 N_FFFT点数です。

ただし、k=0 または k=n_bin-1(一番大きい周波数ビン番号) のときは、

\displaystyle{
 r_1 = |X[k,i]|
}

k'=0 または k'=n_bin-1 のときは、

\displaystyle{
 r_2 = |X[k',i]|
}

また、peak_plane[k,i] == BIRTH または peak_plane[k,i] == DEATH_BIRTH のときは、

\displaystyle{
 r_1 = 0
}

peak_plane[k',i] == DEATH または peak_plane[k',i] == DEATH_BIRTH のときは、

\displaystyle{
 r_2 = 0
}

とします。

その後、ピークごとにフレーム間の信号を再合成します。

ピークのインデックスを p とすると、プレーム間の各信号の振幅  A_p [n] は以下のように線形補間で求めます。

\displaystyle{
 A_p[n] = r_1 + \frac{r_2 - r_1}{S} n \hspace{1em}(n=0,1,\cdots, S-1)  
}

ここで、 S はフレームシフト点数です。

フレーム間の各信号の位相  \theta_p [n] は以下のように3次補間で求めます。

\displaystyle{
 \theta_p[n] = \theta_1 + \omega_1 \frac{n}{f_s} + \alpha (M) \left(\frac{n}{f_s}\right)^2 + \beta (M) \left(\frac{n}{f_s}\right)^3 \hspace{1em}(n=0,1,\cdots, S-1) 
}

 \alpha (M)、\beta (M) は以下のように求めます。

\displaystyle{
 \begin{bmatrix}
\alpha(M)\\
\beta(M)
\end{bmatrix}
=
 \begin{bmatrix}
\frac{3}{T^2}&\frac{-1}{T}\\
\frac{-2}{T^3}&\frac{1}{T^2}
\end{bmatrix}
 \begin{bmatrix}
\theta_2 -\theta_1-\omega_1 T + 2\pi M\\
\omega_2-\omega_1
\end{bmatrix}
}

 M については以下のように求めます。

\displaystyle{
 M = {\rm round} \left( \frac{1}{2\pi} \left[(\theta_1+\omega_1 T - \theta_2) + (\omega_2 - \omega_1)\frac{T}{2}\right] \right) 
}

 A_p[n] と \theta_p[n] を求めたら、再合成信号  y [n] を以下の式で求めます。

\displaystyle{
 y[n] = \sum_{p} A_p[n] \cos{(\theta_p[n])}
}

プログラムの実行

プログラム中のパラメータを表のように設定して、プログラムを実行しました。

パラメータ名 記号 変数名 設定値
サンプリング周波数 f_s sr_o 10 kHz
FFT点数 N_F n_fft 512
フレームシフト点数  S hop_length 256
接続先の周波数範囲 \Delta f delta_freq 50 Hz
窓関数 w window ハミング窓

プログラムの実行は記述したスクリプトmain.pyと音声ファイルvoice.wavを同じフォルダに置き、以下のようにして実行できます。

python main.py

歌声を再合成した結果は以下のようになります。

10kHz にリサンプリング後の歌声

正弦波モデルで再合成した歌声

リサンプリング後の歌声と再合成後の歌声の波形とスペクトログラムは図7です。

図7:歌声の波形とスペクトログラム(上:リサンプリング後、下:再合成後)
図7:歌声の波形とスペクトログラム(上:リサンプリング後、下:再合成後)

少しノイズがありますが、歌声も聴きとれるので多分上手くできてるのかな?

おわりに

今回は、正弦波モデルを実装しました。論文では、窓関数の幅を基本周波数に応じて変えたり、1フレームの最大ピーク数を決めていたりしているので、それらを実装することでさらに良質な歌声を再合成できるかもしれません。

使用した楽曲について

この記事で信号処理した楽曲は Cambridge Music Technology で提供されている Actions の Devil's Words を使用させていただきました。

【楽曲を提供している Cambridge Music Technology のページ】
https://www.cambridge-mt.com/ms/mtk

【ActionsのFacebook
https://www.facebook.com/actionsuk

お問い合わせフォーム プライバシーポリシー

© 2021 Setoti All rights reserved.