MUSIC法による音源定位

本記事ではMUSIC(MUltiple SIgnal Classification)法Pythonで実装しました。MUSIC法はSchmidtにより1986年に提案されており、音源の位置を推定する技術(音源定位)では最も代表的な手法ではないかと思います。今回は、シミュレーション環境で録音した音を用いて、MUSIC法の有効性を確かめました。

本記事は前回の記事(音声・音響処理で使う主成分分析)を読んでいることを前提に書かれています。

主に参考にした本は以下です。

MUSIC法

まず、MUSIC法を理解するのに必要な部分空間の直交性について説明します。それから、狭帯域の信号と広帯域の信号に対するMUSIC法について説明したいと思います。

部分空間の直交性

前回の記事(音声・音響処理で使う主成分分析)でも説明したように、アレイマニフォールドベクトルによって張られる空間は信号部分空間と呼び、空間相関行列の固有ベクトル \boldsymbol{w_1, w_2} については信号部分空間の正規直交基底となっています(図1)。

図1:アレイマニフォールドベクトルと固有ベクトルの幾何学的関係
図1:アレイマニフォールドベクトルと固有ベクトル幾何学的関係

一方、音源の数より大きい番号の固有ベクトル \boldsymbol{w}_3 については以下の式が成り立ちます。

\displaystyle{
\boldsymbol{A}^H \boldsymbol{w}_3 = \boldsymbol{0}
}

音源の数を J 、マイクの数を Mとすると、一般的には以下の式が成り立ちます。

\displaystyle{
\boldsymbol{A}^H \boldsymbol{w}_i = \boldsymbol{0} \hspace{1em} i=J+1, \cdots, M
}

\boldsymbol{A}=[\boldsymbol{A}_1, \cdots, \boldsymbol{A}_J]\boldsymbol{A}_j=[A_{1j}, \cdots, A_{Mj} ]^T とすると、

\displaystyle{
\boldsymbol{A}_j^H \boldsymbol{w}_i = \boldsymbol{0} \hspace{1em} i=J+1, \cdots, M\hspace{1em} j=1,\cdots,J
}

つまり、アレイマニフォールドベクトルと音源の数より大きい番号の固有ベクトルは直交するという性質があります。

空間スペクトル

MUSIC法では次式の空間スペクトル P_k(\theta) \theta 方向から到来する信号のパワー)によって音源の到来方向を求めます。

\displaystyle{
P_{k}(\theta) = \frac{||\boldsymbol{a}_k(\theta)||^2}{\sum_{i=J+1}^{M}|\boldsymbol{a}_k^H(\theta)\boldsymbol{w}_i|^2} = \frac{\boldsymbol{a}_k^H(\theta)\boldsymbol{a}_k(\theta)}{\boldsymbol{a}_k^H(\theta)\boldsymbol{W}\boldsymbol{W}^H \boldsymbol{a}_k(\theta)}
}

ここで、

\displaystyle{
\begin{equation}
\boldsymbol{a}_k(\theta)
=
\begin{bmatrix}
\exp{\left(-j 2\pi \frac{k}{N}f_s \tau_1(\theta)\right)} \\
\vdots \\
\exp{\left(-j 2\pi \frac{k}{N}f_s \tau_M(\theta)\right)} 
\end{bmatrix} 
\end{equation}
}

は周波数ビン番号 k、\theta 方向の仮のアレイマニフォールドベクトルで、線状アレイの場合、遅延時間は以下のようになります。

\displaystyle{
\begin{equation}
\tau_m(\theta)=\left((m-1)-\frac{M-1}{2}\right) \frac{d_x}{c} \sin{\theta}  
\end{equation}
}

 d_x はマイクの間隔、 c は音速です。遅延時間の求め方を知りたい方はビームフォーミングで特定方向の音源を強調という記事をご覧ください。

また、

\displaystyle{
\boldsymbol{W}=[\boldsymbol{w}_{J+1}, \cdots, \boldsymbol{w}_{M}]
}

です。

 \theta = 0°\sim 180° の空間スペクトルを計算して、空間スペクトルのピークの位置が音源の到来方向となります。ピークが音源の数よりも多くある場合は、空間スペクトルが大きいピークの位置を音源の到来方向とします。

なぜMUSIC法ではこのような空間スペクトルを使うかというと、部分空間の直交性から、 \boldsymbol{a}_k(\theta)正しいアレイマニフォールドベクトルの場合、分母は 0 となり、空間スペクトルはピークを持つからです。分子の ||\boldsymbol{a}(\theta)||^2 はアレイマニフォールドベクトルの正規化のための項です。あってもなくてもピークの位置は変わらないため、計算しなくてもしてもいいです。

広帯域信号の場合

さきほどの方法は狭帯域信号に対する方法でした。 音楽信号などの広帯域信号の場合、周波数ビンごとの情報を統合して、最終的な推定結果を得る必要があります。

周波数ビンごとの情報を統合する手法はいくつかあるようですが、この記事では空間スペクトルを平均します。

周波数ビンごとの空間スペクトルを次式のように重み付き平均して、最終的な空間スペクトルを推定します。

\displaystyle{
\bar{P}(\theta) = \frac{1}{K} \sum_{k=k_l}^{k_h} \beta_k P_k(\theta)
}

\displaystyle{
\beta_k = \left[\sum_{i=1}^J \lambda_i[k] \right]^\alpha
}

ここで、K は周波数ビンの数、 \lambda_i は空間相関行列の固有値\alpha は定数で、\alpha=1\alpha=1/2 とします。

空間相関行列の固有値は信号の平均パワーに相当しますので、パワーが大きい周波数ビンほど重みづけするということです

プログラム

MUSIC法によって空間スペクトルと音源の方向を求めるプログラムは以下のmusic.pyです。

import soundfile as sf
import numpy as np
import scipy.signal as sg
import matplotlib.pyplot as plt

# パラメータ
dir_name = "mic/"  # 録音データがあるディレクトリ
png_name = "music.png"  # 出力するグラフの画像名
n_mic = 8          # マイクの数
n_src = 3          # 音源の数
N     = 512        # 窓の大きさ
window = "hann"    # 窓の種類
d     = 0.01       # マイクの間隔[m]
c     = 340        # 音速[m/s]
freq_l = 800       # 空間スペクトルを計算する周波数の下限
freq_h = 3000      # 空間スペクトルを計算する周波数の下限

# WAVファイル名のリスト作成
wav_list = []
for i in range(n_mic):
    wav_list.append("mic"+str(i)+".wav")

# WAVファイルを読み込む
for i, wav_name in enumerate(wav_list):
    x, fs = sf.read(dir_name+wav_name)
    if i==0:
        audio=np.zeros((n_mic,len(x)))
        audio[0,:] = x 
    else:
        audio[i,:] = x

# 周波数ビンに変換
k_l = int(freq_l/(fs/N))
k_h = int(freq_h/(fs/N))

# 短時間フーリエ変換(STFT)を行う X.shape=(n_mic, n_bin, n_frame)
f, t, X = sg.stft(audio, fs, window=window, nperseg=N)
n_bin = X.shape[1]

# 空間相関行列を求める
XH = np.conjugate(X)
XXH = np.einsum("mki,nki->mnki", X, XH)
R = np.mean(XXH, axis=3)

# 固有値を求めて大きい順に並べる
eig_val = np.zeros((n_mic, n_bin), dtype=np.complex64)
eig_vec = np.zeros((n_mic, n_mic, n_bin), dtype=np.complex64)
beta = np.zeros(n_bin)
for k in range(n_bin):
    eig_val_k, eig_vec_k = np.linalg.eig(R[:,:,k]) # 固有値、固有ベクトルを求める
    sort = np.argsort(-1.0*np.abs(eig_val_k))      # 大きい順にargsort
    eig_val[:,k] = eig_val_k[sort]                 # 固有値を大きい順に並べる
    eig_vec[:,:,k] = eig_vec_k[:,sort]             # 固有ベクトルも並べ替える
    
W = eig_vec[:,n_src:,:]   # 音源の数より大きい番号の固有ベクトル    
beta = np.sum(np.abs(eig_val[:n_src,:]), axis=0)  # 重みβを求める

# 空間スペクトル作成
theta = np.linspace(-90.0, 90.0, 181)
P_MU  = np.zeros_like(theta)
a = np.zeros((n_mic, 1, n_bin), dtype=np.complex64)
for i, th in enumerate(theta):
    th = np.radians(th) # deg -> rad 変換
    # 直線状アレイのアレイ・マニフォールド・ベクトルの作成
    for k in range(n_bin):
        fk = fs*k/N  # 周波数
        for m in range(n_mic):
            delay = (m-(n_mic-1)/2)*d*np.sin(th)/c # 遅延時間
            a[m,0,k] = np.exp(-1j*2*np.pi*fk*delay)
    aH = np.conjugate(a)
    aHa = np.einsum("mik,mjk->ijk", aH, a)
    P_MU_th = np.einsum("mik,mek->iek", aH, W)
    P_MU_th = np.einsum("iek,jek->ijk", P_MU_th, np.conjugate(P_MU_th))
    P_MU_th = np.abs(aHa)/np.abs(P_MU_th)
    P_MU[i] = np.mean(beta[k_l:k_h+1]*P_MU_th[0,0,k_l:k_h+1])

# グラフを出力
fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(theta, 10*np.log10(P_MU), c="red")
ax.set_xlim([-90,90])
ax.set_xlabel("Direction [deg]") 
ax.set_ylabel("Spatial spectrum [dB]")
fig.savefig(png_name)

# 候補を出力
index=sg.argrelmax(P_MU) # ピークのあるインデックスを求める
sort=np.argsort(-1.0*P_MU[index])   # ピークの中でargsort
for i in range(n_src):
    print(theta[index[0][sort[i]]]) # 大きい順に方向を表示

6~16行目:録音データがあるディレクトリ、出力するグラフの画像名、マイクの数、窓関数の種類などを指定する。

18~30行目:指定したディレクトリの中にあるmicX.wavという名前のWAVファイルをマイクの数だけ読み込む。

32~34行目:空間スペクトルを平均する周波数の下限と上限を周波数ビン番号に変換する。

36~38行目短時間フーリエ変換を行う。

40~43行目:周波数ビンごとに空間相関行列を求める。

45~53行目固有値を大きい順に並べかえ、対応する固有ベクトルも同じように並べかえる。

55~56行目\boldsymbol{W}, \beta_k を求める。

58~75行目 -90^\circ \leq \theta \leq 90^\circ の空間スペクトルP_k を求める。

77~84行目:空間スペクトルのグラフを出力する。

86~90行目:空間スペクトルのピークの位置を到来方向の候補としてターミナルに出力する。

シミュレーション実験

室内音響シミュレーションで作成した音に対してMUSIC法を使用して、音源の方向を推定しました。

方法

シミュレーション環境の録音方法は前回の記事:音声・音響処理で使う主成分分析での方法と同様です。

室内音響シミュレーションにはPythonのライブラリであるPyRoomAcousticsを用いました。

シミュレーション上で図2のように音源とマイクを配置して、8個のマイクで音楽データを録音しました。

図2:シミュレーションにおける音源とマイクの配置
図2:シミュレーションにおける音源とマイクの配置

残響時間は0.2秒、収録時のサンプリング周波数は16kHz 、SN比は90dB(ほとんど雑音なし)に設定しました。

以下をクリックすると、シミュレーションに用いた Python のコード room_sim.py が展開されます。

結果

プログラムmusic.pyで処理した結果は図3のようになりました。

図3:空間スペクトルのグラフ1
図3:空間スペクトルのグラフ1

また、ターミナルに出力された到来方向の推定結果は表1のようになりました。

表1:到来方向の推定結果1
推定 正解
シンセ 53.0 60.0
ボーカル -13.0 -15.0
ドラム -47.0 -45.0

シンセの方向については7度ずれていましたが、だいたい推定できていました。

補足:シンセの方向を30度、ボーカルの方向を-15度、ドラムの方向を-30度としたとき、図4のような結果となりました。

図4:空間スペクトルのグラフ2
図4:空間スペクトルのグラフ2

また、ターミナルに出力された到来方向の推定結果は表2のようになりました。

表2:到来方向の推定結果2
推定 正解
シンセ -51.0 30.0
ボーカル -13.0 -15.0
ドラム -33.0 -30.0

音源どうしが近すぎるためか、空間スペクトルにいくつもピークが生じてしまい、シンセの方向については推定できませんでした。一応、30度付近にピークはありますが、他のピークのほうが大きいため、推定結果がずれてしまいました。

原因としては、反射・残響により部分空間の直交性が成り立たなくなったことでビーム幅が広がり、近すぎる音源どうしを分離できなかったのではないかと思っています。

おわりに

本記事ではMUSIC法をPythonで実装して、MUSIC法の有効性を確認しました。他の音源定位の手法も実装できたら、手法間で音源定位の精度を比較しようかなと思います。

参考文献
[1] 浅野太、”音のアレイ信号処理 音源の定位・追跡と分離”、コロナ社、2011.

音声・音響処理で使う主成分分析

本記事では、音声・音響処理で使われている主成分分析 (PCA: Principal Component Analysis) を紹介します。自分が知る限り、主成分分析はアレイ信号処理によく使われています。

音声・音響処理する人向けに記述しましたので、データ分析する人向けではないことをご承知おきください。

本記事で主に参考にした本は以下です。

主成分分析とは

主成分分析では図1のようにx1とx2の散布図があるときw1のような新しい軸(変数)を作ります。

図1:x1 とx2 の第1主成分
図1:x1 とx2 の第1主成分

w1の軸は、データを軸に射影したとき、その成分の分散が最大となるように決定します。このw1の軸は第1主成分と呼びます。

第2主成分、第3主成分、・・・については、これまでに決定した主成分と直交して、かつ射影した成分の分散が最大となるように決定します。x1 と x2 の場合、直交する軸は一つしかないので、自動的に第2主成分 w2 が図2のように決定します。

図2:x1 と x2 の第2主成分
図2:x1 と x2 の第2主成分

主成分分析はデータの次元を減らしてデータを近似したり、データ分析に使われたりします。音のアレイ信号処理では、音源の位置を推定する音源定位に使われたりします。

主成分の求め方

以下の空間相関行列 \boldsymbol{R}固有ベクトルが主成分の方向となります。

\displaystyle{
\boldsymbol{R}=E[\boldsymbol{x[n]x^H[n]}] 
}

 E[\cdot] は期待値演算、\boldsymbol{x}[n]=[x_1[n], \cdots, x_m[n], \cdots, x_M[n]]^T \{\cdot\}^H はエルミート転置です。

固有値  \lambda_i固有ベクトル  \boldsymbol{w_i} は次式を満たすものとなっており、M個(マイクの個数)あります。

\displaystyle{
\boldsymbol{R w_i} = \lambda_i \boldsymbol{w_i}
}

相関行列の最大の固有値第1主成分の分散となり、最大の固有値に対応する固有ベクトル第1主成分の方向となります。

それから、2番目に大きい固有値が第2主成分の分散になり、その固有値に対応する固有ベクトルが第2主成分の方向、3番目に大きい固有値が第3主成分の分散になり、その固有値に対応する・・・となっていきます。

期待値演算については、複数の信号の平均(集合平均)を求める演算ですが、普通は1回分の観測信号しか得られないことがほとんどのため、以下のエルゴード性を仮定して解析を行います。

\displaystyle{
E[\boldsymbol{x[n]x^H[n]}] = \frac{1}{N} \sum_{n=0}^{N-1}\boldsymbol{x[n]x^H[n]}
}

 N はサンプル数です。

つまり、集合平均と時間平均が同じになると仮定して計算を行います。

固有ベクトルが主成分の方向となる理由

相関行列の固有ベクトルが主成分の方向となることはいろいろな文献で紹介されていますが、なぜ固有ベクトルが主成分の方向となるかの説明はあまりないと思います。今回、文献をいろいろ探して、参考文献 [1] にその説明がありましたので紹介します。

観測値を \boldsymbol{x} (複素ベクトル)とすると、単位ベクトル  \boldsymbol{w_1} 上に射影された  \boldsymbol{x} の成分である y_1 は以下のようになります。

\displaystyle{
\begin{align}
y_1 &= \frac{\boldsymbol{w_1^H x}}{\boldsymbol{w_1^H w_1}} \\
  &= \boldsymbol{w_1^H x}
\end{align}
}

 E[\boldsymbol{x}]=\boldsymbol{0} であるとすると、 y_1 の分散は以下のように計算されます。

\displaystyle{
\begin{align}
E[|y_1|^2] -  (E[y])^2 &= E[y_1y_1^*] -  0 \\
  &= \boldsymbol{w_1^H} E[\boldsymbol{xx^H}] \boldsymbol{w_1} \\
  &= \boldsymbol{w_1^H R w_1}
\end{align}
}

 ||\boldsymbol{w_1}||=1 (単位ベクトル)という拘束条件を満たしながら、 y_1 の分散(平均パワー)が最大となるように \boldsymbol{w_1} を決定します。

この問題はラグランジュの未定乗数法によって、次式を最大化する問題に置き換えられます。

\displaystyle{
 J = \boldsymbol{w_1^H R w_1} + \lambda (1-\boldsymbol{w_1^H w_1})
}

 \boldsymbol{w_1^*}偏微分して、 \boldsymbol{0} とおくことにより次式のようになります。

\displaystyle{
\begin{align}
 & \frac{\partial J}{\partial \boldsymbol{w_1^*}} = \boldsymbol{R w_1} - \lambda \boldsymbol{w_1} = \boldsymbol{0} \\
 & \boldsymbol{R w_1} = \lambda \boldsymbol{w_1}
\end{align}
}

ここでは、以下の複素ベクトルの偏微分の公式を用いています。

\displaystyle{
\begin{align}
& \frac{\partial}{\partial \boldsymbol{w^*}} \boldsymbol{w^H R w} = \boldsymbol{R w} \\
& \frac{\partial}{\partial \boldsymbol{w^*}} \boldsymbol{w^H w} = \boldsymbol{w}
\end{align}
}

つまり、 \boldsymbol{R w_1} = \lambda \boldsymbol{w_1}固有値問題となるので、相関行列の固有ベクトルが主成分の方向となります。

また、さきほどの式に左から  \boldsymbol{w_1}^H をかけることで、

\displaystyle{
\boldsymbol{w_1^H R w_1} = \lambda 
}

左辺が  y_1 の分散となるので、主成分の分散は相関行列の固有値となります。

さらに、固有ベクトルは複数ありますが、第1主成分は分散が最大となるように決定しますので、最大の固有値に対応する固有ベクトルが第1主成分の方向となります。

第2主成分、第3主成分については、これまでに求めた主成分に直交するという条件の下で、J が最大となるように決定します。固有ベクトルは直交するので、2番目、3番目に大きい固有値に対応した固有ベクトルが第2主成分、第3主成分となります。

部分空間法

部分空間法というのは、参考文献 [1] から引用すると以下のようなことです。

部分空間法は、観測信号を空間相関行列の固有ベクトルが張る固有空間に変換し、解析・処理する方法である。

つまり、主成分分析で得られたものを使った解析や処理のことです。

マイクの数が音源の数より多いときは、信号を低次の部分空間で表すことができます。 さらに、部分空間法を利用するMUSIC法などは精度の良い音源定位を実現できます。

ここでは、信号を低次の部分空間で表せることを説明します。

観測信号を短時間フーリエ変換したとき、それは以下のように表せます。

\displaystyle{
\boldsymbol{X}[i,k] = \boldsymbol{A}[k] \boldsymbol{S}[i,k] \hspace{3em} (A)
}

ここで、 \boldsymbol{X}[i,k]=[X_1[i,k] \cdots X_m[i,k] \cdots X_M[i,k]]^T は観測信号を短時間フーリエ変換したもの、 \boldsymbol{S}[i,k]=[S_1[i,k] \cdots S_j[i,k] \cdots S_J[i,k]]^T は音源信号を短時間フーリエ変換したもの、 i はフレーム番号、k は周波数ビン番号です。\boldsymbol{A}[k]アレイマニフォールド行列で以下のように表せます。

\displaystyle{
\boldsymbol{A}[k] = 
\begin{bmatrix}
A_{11}[k] & \cdots & A_{1j}[k] & \cdots & A_{1J}[k] \\
\vdots & \ddots & & & \vdots \\
A_{m1}[k] & & A_{mj}[k] & & A_{mJ}[k] \\
 \vdots & & & \ddots & \vdots \\
A_{M1}[k] & \cdots & A_{Mj}[k] & \cdots & A_{MJ}[k] 
\end{bmatrix}
}

アレイマニフォールド行列の要素は、各音源から各センサまでの経路の伝達関数となっています。

ビームフォーミングで特定方向の音源を強調 では音源が1つでしたので、アレイマニフォールドベクトルでしたが、今回は音源が複数ある場合を考えるので、アレイマニフォールド行列となっています。

音源の数を  J=3、マイクの数を M=2 とします。このとき、観測信号は以下のように表せます。

\displaystyle{
\begin{bmatrix}
X_{1}  \\
X_{2}  \\
X_{3} 
\end{bmatrix}
= 
\begin{bmatrix}
A_{11} & A_{12}  \\
A_{21} & A_{22}  \\
 A_{31}& A_{32} 
\end{bmatrix}
\begin{bmatrix}
S_{1}  \\
S_{2}  \\
\end{bmatrix}
}

この式を整理して、 S_1, S_2 を消去すると、以下のような式となります。

\displaystyle{
(A_{31}A_{22}-A_{21}A_{32})X_{1}+(A_{11}A_{32}-A_{31}A_{12})X_{2}+(A_{21}A_{12}-A_{11}A_{22})X_{3} = 0
}

これは平面の式になりますので、観測信号は3次元空間において平面上に存在します。したがって、主成分分析で第2主成分の方向まで求めれば、それを使って観測信号を2次元(音源の数の次元)で表現できるということです。

また、アレイマニフォールド行列の列ベクトルであるアレイマニフォールドベクトル \boldsymbol{A_1, A_2}固有ベクトル幾何学的関係は図3のようになります。

図:アレイマニフォールドベクトルと固有ベクトルの幾何学的関係
図:アレイマニフォールドベクトルと固有ベクトル幾何学的関係

アレイマニフォールドベクトルによって張られる空間は信号部分空間と呼びます。固有ベクトル \boldsymbol{w_1, w_2}については信号部分空間の正規直交基底となっています。

実験

信号を低次の部分空間で表せることを実験で確認します。

方法

室内音響シミュレーションにはPythonのライブラリであるPyRoomAcousticsを用いました。

シミュレーション上で図4のように音源とマイクを配置して、8個のマイクで音楽データを録音しました。

図4:シミュレーションにおける音源とマイクの配置
図4:シミュレーションにおける音源とマイクの配置

残響時間は0.2秒、収録時のサンプリング周波数は16kHz 、SN比は90dB(ほとんど雑音なし)に設定しました。また、反射・残響あり、反射・残響なしの2つの場合について調べました。

以下をクリックすると、シミュレーションに用いた Python のコード room_sim.py が展開されます。


信号を低次元(音源の数の次元)で表せる場合、1~3番目の固有値は大きくなり、4番目以降の固有値は急激に小さくなるはずです。

プログラム

各周波数ビンで主成分分析を行うプログラムが以下の pca.py です。

import soundfile as sf
import numpy as np
import scipy.signal as sg
import matplotlib.pyplot as plt

# パラメータ
dir_name = "mic/"  # 録音データがあるディレクトリ
png_name = "eigen.png" # 出力するグラフのファイル名
n_mic = 8          # マイクの数
N     = 512        # 窓の大きさ
window = "hann"    # 窓の種類
f_l   = 800        # 固有値を合計する周波数の下限
f_h   = 3000       # 固有値を合計する周波数の上限

# WAVファイル名のリスト作成
wav_list = []
for i in range(n_mic):
    wav_list.append("mic"+str(i)+".wav")

# WAVファイルを読み込む
for i, wav_name in enumerate(wav_list):
    x, fs = sf.read(dir_name+wav_name)
    if i==0:
        audio=np.zeros((n_mic,len(x)))
        audio[0,:] = x 
    else:
        audio[i,:] = x

# 周波数ビンを求める
k_l = int(f_l/(fs/N))
k_h = int(f_h/(fs/N))

# 短時間フーリエ変換(STFT)を行う X.shape=(n_mic, n_bin, n_frame)
f, t, X = sg.stft(audio, fs, window=window, nperseg=N)
n_bin = X.shape[1]

# 空間相関行列を求める
XH = np.conjugate(X)
XXH = np.einsum("mki,nki->mnki", X, XH)
R = np.mean(XXH, axis=3)

# 固有値を求めて昇順に並べる
eig_val = np.zeros((n_mic, n_bin), dtype=np.complex64)
eig_vec = np.zeros((n_mic, n_mic, n_bin), dtype=np.complex64)
for k in range(n_bin):
    # 固有値と固有ベクトルを求める
    eig_val_k, eig_vec_k = np.linalg.eig(R[:,:,k])
    # 固有値の大きい順に固有値と固有ベクトルを並べる
    sort = np.argsort(np.abs(eig_val_k))[::-1]
    eig_val[:,k] = eig_val_k[sort]
    eig_vec[:,:,k] = eig_vec_k[:,sort]

# 固有値の絶対値の合計を求める
eig_sum = np.sum(np.abs(eig_val[:,k_l:k_h]), axis=1)
eig_val_a = eig_sum/eig_sum[0]  #一番目の固有値の合計で正規化
i = np.arange(1,n_mic+1)

# グラフを出力
fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(i, 20*np.log10(eig_val_a), c="red", marker='o')
ax.set_xlabel("Eigenvalue number") 
ax.set_ylabel("Eigenvalue [dB]")
fig.savefig(png_name)

# ターミナルに出力
print(20*np.log10(eig_val_a))

このプログラムで 800Hzから3000Hz の固有値の絶対値の合計を以下のように求めています。

\displaystyle{
\Lambda_i = \sum_{k=k_l}^{k_h-1} |\lambda_i[k]|
}

k_l は800Hzに対応する周波数ビン番号、k_h は3000Hz に対応する周波数ビン番号です。

結果

番号ごとの固有値の合計は以下のようになりました。

図:固有値分布
図:固有値分布

反射・残響なしの場合は予想したように、4番目の固有値から急激に小さくなっています。しかし、反射・残響ありの場合は、急激に小さくなってはいません。

これは、反射・残響ありの場合、観測信号が式 (A) のように表せないからです。残響時間が長い場合、前のフレームの音源信 \boldsymbol{S}[i-1,k] によっても観測信号は影響されるので、低次元で表せなくなるということです。

おわりに

長くなりましたが本記事では主成分分析を紹介しました。次回は、主成分分析を利用した音源定位であるMUSIC法について紹介しようと思います。

参考文献
[1] 浅野太、”音のアレイ信号処理 音源の定位・追跡と分離”、コロナ社、2011.

信号処理でヘリウムボイスに音声変換

f:id:Setoti:20211029223035p:plain:w250

ヘリウムガスを吸い込んだときの特徴的な声はヘリウムボイスと呼ばれています。本記事では、信号処理でヘリウムボイスを擬似的に実現してみました。前回の記事で紹介したケプストラム分析を用いてスペクトル包絡の加工をします。

本記事で主に参考にした本は以下です。

ヘリウムボイスの原理

参考文献 [1] によると、ヘリウムボイスは、ヘリウムガスを吸い込むことで音速が増大し、声道の共鳴周波数が変化することによって生じる現象のようです。

片方が閉じた管で人間の声道をモデル化して考えます。図1のように管の長さが L [m] のとき、波長が  \lambda=4L となる周波数、およびその奇数倍の周波数で共鳴が起こります。

図1:片側が閉じた管で共鳴する波長
図1:片側が閉じた管で共鳴する波長

音速 c [m/s]、周波数 f [Hz]、波長 λ [m] には次式のような関係がありますので、

\displaystyle{
c = f \lambda
}

音速が r 倍になれば、共鳴する周波数は以下のようになります。

\displaystyle{
f = \frac{rc}{4L} n \hspace{1em} (n=1,3,5,\cdots)
}

したがって、音速が増大すれば共鳴周波数(フォルマント周波数)も大きくなるので、ヘリウムボイスはあのような特徴的な音声になります。

ヘリウムボイスの作成方法

ヘリウムボイスを信号処理で擬似的に実現する方法について説明したいと思います。

スペクトル包絡の伸縮

ヘリウムボイスの原理からフォルマント周波数をr倍、つまり対数スペクトル包絡 H[k] を図2のように周波数方向に r 倍伸縮すればヘリウムボイスを実現できます。

図2:対数スペクトル包絡のr倍の伸縮
図2:対数スペクトル包絡のr倍の伸縮

数式で表現すると以下のようになります。

\displaystyle{
\tilde{H}[k]=H\left[\frac{k}{r}\right] 
}

しかし、k/r が整数になることはあまりないため、以下のような線形補間で \tilde{H}[k] を求めます。

\displaystyle{
\tilde{H}[k]=(1-\alpha_k)H[k'] +\alpha_k H[k'+1] 
}

ただし、

\displaystyle{
\begin{align}
\hspace{2em} k' &= \left\lfloor \frac{k}{r} \right\rfloor \\
\hspace{2em} \alpha_k &= \frac{k}{r} - k'
\end{align}
}

となっています。

参考文献 [1] によると、ヘリウムは空気の約3倍の音速ですが、市販のヘリウムガスは安全のため空気が混入しているそうです。そのため、1.5<r<2.0 としたときが実際のヘリウムボイスに近いみたいです。

ヘリウムボイス変換の流れ

音声をヘリウムボイスに変換するためのブロック図を図3に示します。

図3:ヘリウムボイスに変換するブロック図
図3:ヘリウムボイスに変換するブロック図

フレームごとにケプストラム分析で対数スペクトル包絡を求めて、r倍に伸縮しています。抽出していなかった高次のケプストラムについてもFFTして、伸縮した対数スペクトルに加算しています。それから、log の逆関数であるexpによって振幅を再合成します。位相については観測信号の位相をそのまま用いています。

プログラム

音声をヘリウムボイスに変換するソースコードは以下の helium_vc.py となっています。

import soundfile as sf
import numpy as np
from scipy.fft import rfft, irfft
import scipy.signal as sg

# パラメータ
wav_name = "ATR_PM00.wav"  # 読み込むWAVデータの名前
out_name = "soundout.wav"  # 出力するWAVデータの名前
window = "hann"     # 窓関数の種類
N = 1024            # FFT点数
r = 1.5             # スペクトル包絡の伸縮率

# WAVファイルを読み込む
x, fs = sf.read(wav_name)

# 短時間フーリエ変換(STFT)を行う X.shape=(n_bin, n_frame)
_, _, X = sg.stft(x, fs, window=window, nperseg=N)
X_phase = np.angle(X)   # 観測信号の位相
n_bin   = X.shape[0]    # ビン数
n_frame = X.shape[1]    # フレーム数

# 各numpy配列を準備
ceps_l = np.zeros(N)       # 低次のケプストラム用の配列
ceps_h = np.zeros(N)       # 高次のケプストラム用の配列
H_tilde = np.zeros(n_bin)  # 伸縮後のスペクトル包絡用の配列
Y_abs = np.zeros(X.shape, dtype=np.float64) # 出力信号の振幅用の配列
eps = np.finfo(np.float64).eps  # マシンイプシロン

# フレームごとにr培に伸縮したスペクトル包絡を求める
for i in range(n_frame):
    spec_log = np.log(np.abs(X[:,i])+eps)  # 対数変換
    ceps = irfft(spec_log)     # IFFTしてケプストラムを求める
    lifter = 72                # 低次のケプストラムを72点まで抽出
    ceps_l[0:lifter] = ceps[0:lifter]        # 低次の抽出(前半)
    ceps_l[N-lifter+1:] = ceps[N-lifter+1:]  # 低次の抽出(後半)
    ceps_h[lifter:N-lifter+1] = ceps[lifter:N-lifter+1]  # 高次の抽出
    H = np.real(rfft(ceps_l))  # FFTして実部だけ取り出す
    G = np.real(rfft(ceps_h))  # FFTして実部だけ取り出す
    # 対数スペクトル包絡をr倍に伸縮
    for k in range(n_bin):
        k2 = int(k/r)
        alpha = k/r - k2
        if k2 < n_bin-1:
            H_tilde[k] = (1-alpha)*H[k2] + alpha*H[k2+1] # 線形補間
        else:  # k2がn_binを超えた場合 
            H_tilde[k] = np.log(eps)  # -∞ に近いものを代入
    Y_abs[:,i] = np.exp(H_tilde+G) # 振幅スペクトルを求める

# 位相と振幅でスペクトログラムを合成
Y = Y_abs * np.exp(X_phase)

# 逆短時間フーリエ変換(ISTFT)を行う
_, y = sg.istft(Y, fs=fs, window=window, nperseg=N)

# ファイルに書き込む
y = y/np.max(np.abs(y)) # ノーマライズ
sf.write(out_name, y, fs, subtype="PCM_16")

6~11行目:読みこむWAVデータ名、出力するデータ名、窓関数の種類、FFTの点数、スペクトル包絡の伸縮率などを指定しています。

13~14行目:WAVデータを読み込んでいます。

16~20行目短時間フーリエ変換をしています。また、観測信号の位相を取得しています。

22~26行目:あらかじめ各numpy配列を用意しています。

31~38行目:フレームごとにケプストラム分析をして、対数スペクトル包絡 H を求めています。また、高次抽出をして微細構造 G も求めています。子音を発するときなどは、ケプストラムに基本周波数のピークがないので、低次抽出する点数は72点で固定しています。

39~46行目:対数スペクトル包絡を r 倍に伸縮しています。r<1.0 のときは k2(k') にスペクトル包絡がない場合があるので、そのときはとても大きいマイナスの数を代入しています。

47行目:対数スペクトル包絡 H と微細構造 G を足し合わせて、指数関数 exp で振幅スペクトルを求めています。

49~50行目:再合成した振幅と観測信号の位相を用いて、出力信号のスぺクトログラムを求めています。

52~53行目:逆短時間フーリエ変換で波形データに戻しています。

55~57行目:オーバーフローしたときのためにノーマライズして、WAVデータで出力しています。

処理結果

ヘリウムボイスに変換した男性と女性の音声が以下です。窓関数はハン窓、FFT点数は1024点、オーバーラップは1/2、スペクトル包絡の伸縮率は r=1.5 倍としました。一応、r=0.5 としたときの処理結果についても載せました。音速が遅い気体を吸い込んだらおそらくr=0.5のときのような声になると思います。

男声の音声

男声の音声の処理結果(r=1.5)

男声の音声の処理結果(r=0.5)

女声の音声

女声の音声の処理結果(r=1.5)

女声の音声の処理結果(r=0.5)

信号処理によるノイズ(ミュージカルノイズ)はありますが、ヘリウムガスを吸い込んだときのような声になりました。r=0.5のときは、野太い人の声になってますね。

おわりに

本記事ではケプストラム分析を用いて音声をヘリウムボイスに変換しました。よくあるボイスチェンジャーがどのようなアルゴリズムになっているか理解できてよかったです。

参考文献
[1] 川村新、”音声音響信号処理の基礎と実践”、コロナ社、2021.

使用したデータについて
この記事で信号処理した音声は髙橋弘太研究室の話速バリエーション型音声データベース(SRV-DB)で提供されているATR25文の読み上げ音声を使用させていただきました。

【音声データベースを提供している髙橋弘太研修室のページ】
http://www.it.cei.uec.ac.jp/SRV-DB/

ケプストラム分析によるスペクトル包絡推定

ケプストラム分析を用いてスペクトル包絡推定を行いました。参考文献 [1] によると、ケプストラム分析は1964年にNoll らにより提案された方法です。ちなみにケプストラム(cepstrum)という名称は、spectrumのアナグラム(specだけを逆から書くとcepstrum)となっています。

ケプストラム分析

ケプストラム分析について説明します。スペクトル包絡について説明してから、ケプストラム分析の考え方、リフタリングについて説明しようと思います。

スペクトル包絡とは

ソースフィルタモデルで母音を音声合成 で説明しましたように、ソースフィルタモデルではz変換された音声 Y は以下のように表せます。

\displaystyle{
Y(z) = H(z) X(z)
}

ここで、H(z) は声道の伝達関数、X(z) は声帯振動で得られた音源をz変換したものです。

つまり、音声の周波数特性 Y(ω) は以下のように表せるということです。

\displaystyle{
Y(\omega) = H(\omega) X(\omega)
}

H(ω) は声道フィルタの周波数特性、X(ω) は音源の周波数特性です。

声道フィルタの振幅特性 |H(ω)| についてはスペクトル包絡とも呼ばれています。

今回はケプストラム分析を用いて、実際の音声からこのスペクトル包絡を推定します。

ケプストラム分析の考え方

ケプストラム分析の考え方は意外に単純です。「あ」を発音した音声データの対数振幅スペクトルを図1に示します。ちなみに0dBに対応するのはフルスケールの正弦波のパワーです。

図1:「あ」の対数振幅スペクトル
図1:「あ」の対数振幅スペクトル

ソースフィルタモデルから音声の対数振幅スペクトルは以下のように表せます。

\displaystyle{
\begin{align}
\log{(|Y(\omega)|)} &= \log{(|H(\omega) X(\omega)|)}\\
&=  \log{(|H(\omega)|)} +  \log{(|X(\omega)|)}
\end{align}
}

つまり、音声の対数振幅スペクトルは声道フィルタの対数振幅スペクトルと音源の対数振幅スペクトルの加算で表せるということです。

声帯振動波形はパルス列のようになりますので、音源の対数振幅スペクトルも基本周波数の間隔でパルスが立つパルス列となります。つまり、図1が周期的にギザギザになっているのは音源の特性が原因と考えられます。

一方、声道フィルタの対数振幅スペクトルは緩やかな形状を持ちますので、さらにフーリエ変換(DFT)することで声道特性と音源特性が低域と高域に分離できるのではないかというのがケプストラム分析の考え方です。

音声の対数振幅スペクトルをさらにDFTしたものをケプストラムと呼びます。また、ケプストラムの横軸についてはケフレンシー(quefrency はfrequecyのアナグラムです)軸と呼びます。

補足:参考文献[1][2] では、対数振幅スペクトルをDFTではなくIDFTしたものをケプストラムとしています。理由としては、対数振幅スペクトルが偶関数であるため、以下のようにIDFTとDFTが係数を除いて、同じになるからです。

\displaystyle{
\begin{align}
{\rm IDFT}[\log{|Y[k]|}]&= \frac{1}{N} \sum_{k=0}^{N-1} \log{|Y[k]|} e^{j 2\pi kn/N} \\
&=  \frac{1}{N} \sum_{k=0}^{N-1} \log{|Y[N-k]|} e^{-j 2\pi kn/N} \\
&=  \frac{1}{N} \sum_{k=0}^{N-1} \log{|Y[-k]|} e^{-j 2\pi kn/N} \\
&= \frac{1}{N} \sum_{k=0}^{N-1} \log{|Y[k]|} e^{-j 2\pi kn/N}  \\
&= \frac{1}{N} {\rm DFT}[\log{|Y[k]|}]
\end{align}
}

また、スペクトルを逆フーリエ変換したと考えると横軸は時間 [s] とみなせるので、音声分野では逆フーリエ変換したものをケプストラムと呼ぶ説明が多いようです。これ以降の説明ではIDFTしたものをケプストラムとします。

リフタリング

「あ」のケプストラムは図2のようになります。

図2:「あ」のケプストラム
図2:「あ」のケプストラム

音源特性は図2のように基本周期  T_o にピークがあり、声道特性は低ケフレンシーに集中するため、ケプストラムの低ケフレンシー部だけを抽出し、これをDFTすることでスペクトル包絡が推定できます。

ケプストラムの低ケフレンシー部だけを抽出することをリフタリング(liftering は filtering のアナグラムです)と呼びます。

ケプストラムを抽出する範囲については、参考文献やWEBに明確な記載はなかったので、私は基本周期の半分まで抽出しました。

プログラム

ケプストラムと推定したスペクトル包絡を出力するプログラムは以下のcepstrum.py です。

解析するWAVファイルをソースコードと同じディレクトリに入れて、wav_namestartを書き換えて、「python cepstrum.py」と実行すれば、ceps.png(ケプストラムのグラフ) と sp_env.png(スペクトル包絡のグラフ)が出力されます。

import soundfile as sf
import numpy as np
from scipy.fft import rfft, irfft
import scipy.signal as sg
import matplotlib.pyplot as plt

# パラメータ
wav_name = "aiueo.wav"   # 読み込むWAVデータの名前
ceps_name = "ceps.png"   # 出力するPNGデータの名前
spec_name = "sp_env.png" # 出力するPNGデータの名前
window = "hann"          # 窓関数の種類
start = 10000            # DFTする開始点
N = 2**(14)              # DFT点数

# WAVファイルを読み込む
x, fs = sf.read(wav_name)
w = sg.get_window(window, N) # 窓関数の作成
w = w / np.sum(w)            # 窓関数の補正

# 対数振幅スペクトルを求める
eps = np.finfo(np.float64).eps         # マシンイプシロン
spec = rfft(x[start:start+N]*w)        # DFTする
spec_log = np.log10(np.abs(spec)+eps)  # 対数変換
freq  = np.arange(N//2+1) * (fs / N)   # 周波数軸

# ケプストラムを求める
ceps = irfft(spec_log)       # IDFTする
t = np.arange(N)*1000 / fs   # ケフレンシー軸

# ケプストラムのグラフ出力
To_l = int(fs/800) # 基本周波数の推定範囲の上限を800Hzとする
To_h = int(fs/40)  # 基本周波数の推定範囲の下限を40Hzとする
ylim = np.max(ceps[To_l:To_h]) # 基本周期のピーク
fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(t, ceps, c="red")
ax.set_ylim([-0.02,ylim+0.02]) # 基本周期のピーク+0.2まで表示
ax.set_xlim([0,30])            # 30msまで表示
ax.set_xlabel("Quefrency [ms]") 
ax.set_ylabel("log amplitude")
fig.savefig(ceps_name)  # ケプストラムを出力

# リフタリング
To = np.argmax(ceps[To_l:To_h])+To_l  # 基本周期の点数を求める
lifter = To//2                        # 基本周期の半分まで抽出
ceps[lifter:N-lifter+1] = 0           # リフタリング
sp_env = np.real(rfft(ceps))          # DFTして実部だけ取り出す

# スペクトル包絡のグラフ出力
fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(freq, 10*spec_log, c="red")
ax.plot(freq, 10*sp_env, c="blue")
ax.set_xlim([0,4000])   # 4kHzまで表示
ax.set_ylim([-50,-10]) 
ax.set_xlabel("Frequency [Hz]") 
ax.set_ylabel("Amplitude [dB]")
fig.savefig(spec_name)  # スペクトル包絡を出力

7~13行目:読み込むWAVファイルや出力ファイルの名前、窓関数の種類、データの解析する位置、FFT点数などのパラメータを決めています。

15~18行目:WAVファイルの読み込みと窓関数の作成を行っています。窓関数の補正も行っています。

20~24行目:窓関数をかけてDFTしています。さらに対数振幅スペクトルを求めています。対数変換する際は、計算エラーが発生しないように微小量(マシンイプシロン)を加算してから対数変換しています。

26~28行目:対数振幅スペクトルをIDFTして、ケプストラムを求めています。

30~41行目:ケプストラムのグラフを出力しています。いい感じにグラフが表示されるようにy軸の上限を基本周期のピーク+0.2にしています。

43~47行目:基本周期の点数を求めて、基本周期の半分までケプストラムを抽出しています。基本周期は1.25ms (800Hz) ~25ms (40Hz) の間の最大値を探すことで求めています。また、リフタリングする際は、ケプストラムはN/2で偶対称となっているので、(N-lifter+1)~(N-1) のケプストラムについても残しておきます。

49~58行目:対数振幅スペクトルとスペクトル包絡のグラフを出力しています。

ケプストラム分析の結果

以下の「あいうえお」を発音した音声をケプストラム分析しました。

あいうえおの音声

変な音声になっているのは、記事末尾に記載しているデータベースの音声から「あ」「い」「う」「え」「お」を抜き出して、つなぎ合わせてるからです。

「あいうえお」をケプストラム分析してスペクトル包絡を推定した結果は図3のようになりました。

図3:「あいうえお」のスペクトル包絡
図3:「あいうえお」のスペクトル包絡

緑の点線は、参考文献 [3] に記載されていた男性33人のフォルマント周波数 F1~F3 の平均です。

緑の点線がある箇所にスペクトル包絡の山がだいたいあるので、上手く推定できているかな? 正直いって、これが上手く推定できていると断定はできませんが、多分上手くできているでしょう。

おわりに

ケプストラム分析を用いてスペクトル包絡の推定を行いました。実は今回ケプストラム分析を紹介したのは、次回の記事でスペクトル包絡を用いて、音声変換がしたかったためです。次回、音声変換について書きたいと思います。

参考文献
[1] 森勢将雅、”音声分析合成”、コロナ社、2018.
[2] 川村新、”音声音響信号処理の基礎と実践”、コロナ社、2021.
[3] 電子情報通信学会、”聴覚と音声”、コロナ社、1973.

使用したデータについて
この記事で信号処理した音声は髙橋弘太研究室の話速バリエーション型音声データベース(SRV-DB)で提供されているATR25文の読み上げ音声を使用させていただきました。

【音声データベースを提供している髙橋弘太研修室のページ】
http://www.it.cei.uec.ac.jp/SRV-DB/

頭部伝達関数を使って8D立体音響を実装

頭部伝達関数(HRTF: Head-Related Transfer Function)を使って8D立体音響をPythonで実装してみました。8D立体音響というのは正直言って私も厳密な定義はわかりませんが、この記事で言っているのはYoutubeによくある頭の周りをクルクルまわる音源のことです。「左右のチャネルの音量を変えれば作れるのでは?」という声を聞きますが、それではできず、頭部伝達関数というものを使う必要があります。

注意点として、8D立体音響を作るソフトがどのように処理しているかはわからないため、細かいところは異なることをご承知おきください。

頭部伝達関数(HRTF)

頭部伝達関数とはヒトが音の方向を知覚するために重要な役割を果たすものです。

参考文献 [1] では、頭部伝達関数を以下のように定義しています。

音波は鼓膜に届く直前に頭や耳介,あるいは胴体の影響を受ける。このような,頭部周辺による入射音波の物理特性の変化を周波数領域で表現したものを頭部伝達関数という。

つまり、頭や耳、鼻などが与える音への影響のことです。

例えば、図1のように左側に音源があるとき、音源と左耳の間には何もないため、音はあまり変化なく左耳に伝わります。一方、音源と右耳の間には頭や鼻という障害物があるため、右耳に伝わる音は減衰したり、左耳と比べると遅れて音が届いたりします。

図1:音が伝わる模式図
図1:音が伝わる模式図

このように左右の耳で異なる音の伝わり方をすることで、ヒトは音の位置や方向を知覚できるそうです。冷静に考えると驚くハナシで、信号処理では音が到来する2次元方向を把握するのも結構難しいので、音源の3次元の方向や位置がわかる人間はすごいなと思います。

8D立体音響を作成する場合は、測定した頭部伝達関数を音データに作用させることで左右前後に音源があると錯覚させることができます。

頭の形などは個人差があるので、他人の頭部伝達関数を使うと前後方向については上手く錯覚させることはできませんが、左右方向については上手く錯覚させることができます。

頭の周りをまわる音源の作成

Youtube によくある頭の周りをまわる音源を作成したいと思います。

音源の位置の計算

音源が頭の周りをまわっている様子を図2に示します。

図2:音源の位置の計算
図2:音源の位置の計算

音源が図2のように角速度  \omega [deg/s] でまわっているとき、サンプル番号 n のときの音源方向角度  \theta [deg] は以下のように求まります。

\displaystyle{
\theta = \omega \frac{n}{f_s}
}

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

音源方向角度が360度以上の場合、プログラム上扱いにくいので、0度以上360度未満になるように360度の倍数で引きます。

頭部伝達関数の補間

今回使用するHRTFのデータベースは5度間隔でHRTFが測定されています。そのため、音源方向角度が7度や8度のような中途半端な場合、どのようなHRTFを使えばいいかという問題があります。

今回は、参考文献 [2] に記述されている線形2点補間を用いて、HRTFの補間を行います。計算は簡単で以下のようにHRTFを補間します。

\displaystyle{
H[k] = rH_1[k] + (1-r)H_2[k]
}

ここで、r 0\leq r \leq 1 の2つのHRTFの内分比を表します。例えば、音源方向角度が8度の場合、 H_1は5度のHRTF、 H_2 は10度のHRTFを用いて、 r は0.4とします。

補足:今回行う線形2点補間は参考文献 [2] のやり方とは少し異なっています。参考文献 [2] ではHRTFの振幅応答を補間して求めて、位相応答はHRTFが最小位相フィルタとなるように求めています。私の補間の仕方では位相が打ち消しあって、変な音になる可能性がありますが、処理した音はそこまで問題がなかったのでこのままにしています。

オーバーラップ加算法

いままでデータベースの中にはHRTFが入っているような言い方をしていましたが、厳密には頭部インパルス応答(HRIR:Head-Related Impulse Response)が入っています。HRIRをFFTすることで、HRTFになります。

HRIRを波形データに畳み込めば、クルクルまわる音源ができますが、かなり時間がかかってしまいます。また、for文の処理が遅いPythonではなおさら計算時間がかかります。

そこで、参考文献 [1] に記述されているオーバーラップ加算法 を行います。オーバーラップ加算法の手順は以下です。

(1) 音源信号 x からchg_len ずつずらしながらHRIRの長さ512点を取り出していく。

(2) 取り出した音源信号とHRIRに512点の0埋めを行う

(3) FFT、複素乗算、逆FFTを施してデータ長 1024点の結果を得る。

(4) 処理した結果を加算していって、出力信号 y とする。

図3:オーバーラップ加算法
図3:オーバーラップ加算法

プログラム

8D立体音響を作成するプログラム 3Dsound.py は以下です。

import soundfile as sf
import numpy as np
from scipy.fft import rfft, irfft
import scipy.signal as sg

# パラメータ
wav_name = "sound_mono.wav"  # 読み込むWAVデータの名前
out_name = "soundout.wav"    # 出力するWAVデータの名前
elev = 0         # 仰角
hrtf_dir = "hrtfs/elev"+str(elev)+"/"  # HRTFがあるディレクトリ
N       = 512    # HRTFの点数
chg_len = 128    # HRTFを変える間隔
omega   = 30     # 角速度 [deg/s]

# WAVファイルを読み込む
x, fs = sf.read(wav_name)
if len(x.shape) == 2:  # ステレオのときは左チャネルだけ読む
    x = x[:,0]
x_len = len(x)
x = np.pad(x,[0,N*2], "constant")

# 方位角0°~355°のHRIRを読み込む
hrir_l = np.zeros((72, N))
hrir_r = np.zeros((72, N))
for i, angle in enumerate(range(0,360,5)):
    angle_str = str(angle)
    angle_str = angle_str.zfill(3)
    # 左耳のHRIRを読み込む
    path = hrtf_dir + "L"+str(elev)+"e" + angle_str + "a.dat"
    with open(path) as f:
        for j, s_line in enumerate(f):
            hrir_l[i,j] = float(s_line)
    # 右耳のHRIRを読み込む
    path = hrtf_dir + "R"+str(elev)+"e" + angle_str + "a.dat"
    with open(path) as f:
        for j, s_line in enumerate(f):
            hrir_r[i,j] = float(s_line)

# FFT をして HRTF 作成
HRTF_L = np.zeros((72, N+1), dtype=np.complex128)
HRTF_R = np.zeros((72, N+1), dtype=np.complex128)
for m in range(72):
    h = np.pad(hrir_l[m,:], [0,N], 'constant')  # 0埋め
    HRTF_L[m,:] = rfft(h)
    h = np.pad(hrir_r[m,:], [0,N], 'constant')  # 0埋め
    HRTF_R[m,:] = rfft(h)

# フレームごとに異なるHRTFを掛ける
y = np.zeros((len(x),2))
n_frame = x_len // chg_len + 1  # フレーム数
for i in range(n_frame):
    # 移動音源がどの角度にあるか計算
    theta = omega * i * chg_len / fs
    while int(theta) > 359:   # 0<theta<360にする
        theta = theta - 360
    # HRTFを線形2点補間するためのパラメータを求める
    m = theta / 5
    m1 = int(m)
    m2 = m1 + 1
    if m2 == 72:
        m2 = 0
    r2 = m - int(m)
    r1 = 1.0 - r2
    # 取り出した x を FFT する
    x_N = np.pad(x[i*chg_len:i*chg_len+N], [0,N], 'constant') # 0埋め
    X = rfft(x_N)
    # 補間したHRTF と X を掛ける
    YL = X * (r1*HRTF_L[m1,:]+r2*HRTF_L[m2,:])
    YR = X * (r1*HRTF_R[m1,:]+r2*HRTF_R[m2,:])
    # 逆FFT をして足し合わせる
    y[i*chg_len:i*chg_len+2*N, 0] += irfft(YL)
    y[i*chg_len:i*chg_len+2*N, 1] += irfft(YR)

# ファイルに書き込む
y = y/np.max(y)  # ノーマライズ
sf.write(out_name, y[:x_len,:], fs, subtype="PCM_16")

6~13行目:使用するHRTFの仰角やchg_len、角速度などのパラメータを入力しています。

15~20行目:WAVファイルから波形データを読み込んでいます。WAVファイルがステレオデータの場合は、左チャネルだけ取り出します。

22~37行目:0度から355度のHRIRをnumpy配列に読み込んでいます。

39~46行目:読み込んだHRIRを0埋めして、FFTすることでHRTFを作成しています。

52~55行目:フレームごとに移動音源の角度を求めています。

56~63行目:HRTFを線形2点補間するためのパラメータを求めています。例えば、音源方向角度が17度の場合、17/5=3.4 でm1=3, m2=4 となります。インデックス3と4に対応するHRIRは15度と20度です。また、r2=3.4-3=0.4 、r1=1-r2=0.6 となり、線形補間のパラメータが求まります。

64~66行目:x から512点取り出して、0埋めをして、FFTをしています。

67~72行目:補間して求めたHRTFとXを乗算して、逆FFTをして、加算しています。

74~76行目ノーマライズして、波形データをファイルに書き込んでいます。

処理結果

頭部伝達関数を使って作成した8D立体音響は以下のようになります。仰角は0度、chg_lenは128点で作成しました。ヘッドホンやイヤホンでお聴きください。

入力音源

処理結果(角速度ω:30 [deg/s])

処理結果(角速度ω:120 [deg/s])

そんなに雑音もなく作成できていると思います。角速度ω:120 [deg/s] はさすがに少し気持ち悪くなりますね。

おわりに

HRTFを使って8D立体音響を実装してみました。Youtubeの8D立体音響ではボーカルだけこのような処理をして、他の伴奏とかは普通にミックスしているのかもしれません。あとはコンサートホールなどのインパルス応答も畳み込んだりして、いい感じにしているのだろうと思います。

参考文献
[1] 飯田 一博、”頭部伝達関数の基礎と3次元音響システムへの応用”、コロナ社、2017.
[2] 西野隆典, 梶田将司, 武田一哉, 板倉文忠, "水平方向及び仰角方向に関 する頭部伝達関数の補間," 日本音響学会誌, 57巻, 11号, pp.685-692, 2001.

使用したデータについて
この記事で信号処理した楽曲は Cambridge Music Technology で提供されている The Balazs Daniel Boogie Woogie Trio の Own Way To Boogie を使用させていただきました。

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

頭部伝達関数につきましては西野隆典研究室のホームページにある頭部伝達関数データベースからHRTF data (2)を使用させていただきました。

【頭部伝達関数データベースのページ】
http://www.sp.m.is.nagoya-u.ac.jp/HRTF

ビームフォーミングで特定方向の音源を強調

ビームフォーミングで特定方向の音源を強調しました。ビームフォーミングとは所定の方向に波(電波、音波など)の指向性を高める技術のことです(Wikipediaより)。今回は、一番基本的なビームフォーマである 遅延和ビームフォーマ(DSBF:Delay and Sum Beam Former)Python で実装して、室内音響シミュレーションで作った音に作用させました。

遅延和ビームフォーマ

はじめに、DSビームフォーマの概要について説明していきます。それから、周波数領域のハナシや行列形式のハナシをしていこうと思います。

概要

DSビームフォーマでは複数のマイクを並べて音を録音します。マイクの代表的な並べ方には直線状円状がありますが、今回は直線状に並べて考えていきます。

1つの平面波だけ存在すると考えると、図1のようにマイクで録音される信号(観測信号)に遅延が発生します。

図1:DSビームフォーマのシステム
図1:DSビームフォーマのシステム

そこで、DSビームフォーマでは  \tau_m だけ時間を進ませて時間遅れをなくします。そうすると、すべてのチャネルで信号の位相がそろいますので、信号を足し合わせることで平面波が到来する方向の音を強調することができます。一方、他の方向から到来する信号は位相がずれて足し合わされるので、減衰します。

以上がDSビームフォーマの概要となっています。

直線状アレイの遅延

図2のように直線状にマイクを配置した場合の中心からの遅延時間について考えていきます。ビームフォーマでは絶対的な遅延時間ではなく、相対的な遅延時間差が重要となりますので、音源ではなく別の場所からの遅延時間を考えても問題ないです。

図2:直線状アレイの遅延時間
図2:直線状アレイの遅延時間

図2の  \theta 方向から平面波が到来した場合、赤い点線の距離  D_m は以下のように表せます。

\displaystyle{
\begin{equation}
D_m=\left((m-1)-\frac{M-1}{2}\right) d_x \sin{\theta}  
\end{equation}
}

ここで、 d_x はマイクの間隔 [m]、 M はマイクの個数です。

そのため、 c を音速 [m/s] とすると、m番目のマイクにおける遅延時間は次式のように計算されます。

\displaystyle{
\begin{equation}
\tau_m=\left((m-1)-\frac{M-1}{2}\right) \frac{d_x}{c} \sin{\theta}  
\end{equation}
}

周波数領域

概要で説明したことを周波数領域で考えていきます。また、デジタル信号で考えていきます。時間領域でもDSビームフォーマはできますが、周波数領域ではアナログな遅延時間を扱えるというメリットがあります。

時間領域のデジタル信号で時間進み  \tau_m のある信号は左辺、それをSTFTしたものは右辺のように表せます。ここで、fsはサンプリング周波数、i はフレーム番号、k は周波数ビン番号、NはFFT点数です。

\displaystyle{
\begin{equation}
x_m[n+\tau_m f_s] \xrightarrow{STFT} \exp{\left(j 2\pi \frac{k}{N}f_s \tau_m\right)}X_m[i, k] 
\end{equation}
}

時間領域の左辺では  \tau_m f_s は整数ではないので、小数点以下を四捨五入などして整数にする必要があり、きちんと他の観測信号と位相を揃えることができません。

一方、短時間フーリエ変換(STFT)した右辺の周波数領域では、 \tau_m f_s が整数でなくても扱えますので、きちんと他の観測信号と位相を揃えることができます。

このようなメリットがあるので、STFTした信号を扱っていきます。

アレイマニフォールドベクトル

源信号と観測信号の周波数領域における関係を行列で表現します。

源信号の周波数領域の表現を S[i,k]とすると、観測信号の周波数領域の表現 X_m[i,k]は以下のように表せます。

\displaystyle{
\begin{align}
x_m[n] &= s[n-f_s\tau_m] \\ 
& ↓ {\small STFT} \\
X_m[i, k] &= \exp{\left(-j 2\pi \frac{k}{N}f_s \tau_m\right)}S[i, k] 
\end{align}
}

これを行列で表現すると以下のように表せます。

\displaystyle{
\begin{equation}
\begin{bmatrix}
X_1[i,k] \\
\vdots \\
X_M[i,k]
\end{bmatrix}
=
\begin{bmatrix}
\exp{\left(-j 2\pi \frac{k}{N}f_s \tau_1\right)} \\
\vdots \\
\exp{\left(-j 2\pi \frac{k}{N}f_s \tau_M\right)} 
\end{bmatrix}
S[i, k] 
\end{equation}
}

ここで、

\displaystyle{
\begin{equation}
\boldsymbol{a}[k]
=
\begin{bmatrix}
\exp{\left(-j 2\pi \frac{k}{N}f_s \tau_1\right)} \\
\vdots \\
\exp{\left(-j 2\pi \frac{k}{N}f_s \tau_M\right)} 
\end{bmatrix} 
\end{equation}
}

アレイマニフォールドベクトルと呼ばれ、音源定位や音源分離で重要なものとなります。

このアレイマニフォールドを打ち消したいので、次のようなベクトル \boldsymbol{w}^Hを考えます。

\displaystyle{
\begin{equation}
Y[k,i] = \boldsymbol{w}^H[k] X_m[i, k] = S[i, k]
\end{equation}
}

この  \boldsymbol{w}^H は以下のように計算できます。

\displaystyle{
\begin{equation}
\boldsymbol{w}^H[k] =\frac{\boldsymbol{a}^H[k]}{M}
\end{equation}
}

 \boldsymbol{w}^Hステアリングベクトルと呼ばれ、ビーム方向(強調する方向)を決定するものとなっていますので、取り出したい音源方向のステアリングベクトルを計算して、観測信号との行列積を求めれば、その方向の音源だけ強調することができます。

プログラム

作成したDSビームフォーマのプログラム dsbf.py は以下です。

import soundfile as sf
import numpy as np
import scipy.signal as sg

# パラメータ
dir_name = "mic/"  # 録音データがあるディレクトリ
n_mic = 64         # マイクの数
theta = 120        # 強調する方向 [度]
d     = 0.01       # マイク間の距離 [m]
N     = 256        # 窓の大きさ
c     = 340        # 音速 [m]
window = "hann"    # 窓の種類

# WAVファイル名のリスト作成
wav_list = []
for i in range(n_mic):
    wav_list.append("mic"+str(i)+".wav")

# WAVファイルを読み込む
for i, wav_name in enumerate(wav_list):
    x, fs = sf.read(dir_name+wav_name)
    if i==0:
        audio=np.zeros((n_mic,len(x)))
        audio[0,:] = x 
    else:
        audio[i,:] = x

# 短時間フーリエ変換(STFT)を行う X.shape=(M, K, I)
f, t, X = sg.stft(audio, fs, window=window, nperseg=N)

# 直線状アレイのアレイ・マニフォールド・ベクトル
n_bin = N//2+1  # ビンの数
a = np.zeros((n_bin, n_mic, 1), dtype=np.complex64)
theta = np.radians(theta) # deg -> rad 変換
for k in range(n_bin):
    fk = fs*k/N  # 周波数
    for m in range(n_mic):
        delay = (m-(n_mic-1)/2)*d*np.sin(theta)/c # 遅延時間
        a[k,m,0] = np.exp(-1j*2*np.pi*fk*delay)

# ステアリングベクトル
w = np.zeros((n_bin, 1, n_mic), dtype=np.complex64)
for k in range(n_bin):
    w[k,:,:] = np.conj(a[k,:,:].T)/n_mic

# ステアリングベクトルをかける
Y  = np.einsum("ksm,mki->ski", w, X)

# 逆短時間フーリエ変換(ISTFT)を行う
t, y = sg.istft(Y, fs=fs, window=window, nperseg=N)

# ファイルに書き込む
sf.write("soundout.wav", y[0], fs, subtype="PCM_16")

5~12行目:マイクの個数や音を強調する方向、短時間フーリエ変換(STFT)のパラメータなどを設定しています。

19~26行目:マイクの個数分のWAVファイルを読み込んで、numpy配列を作成しています。

28~29行目:波形データを短時間フーリエ変換しています。scipy signal の stft で マイクの個数分の波形データを一気にスぺクトログラムに変換しています。フレームシフト点数は指定しない場合、窓関数の大きさの半分となっています。

31~39行目:周波数ビンごとに遅延時間を計算して、直線状アレイのアレイマニフォールドベクトルを作成しています。

41~44行目:アレイマニフォールドベクトルからステアリングベクトルを作成しています。

46~47行目:周波数ビンごとに観測信号をSTFTしたものとステアリングベクトルの行列積を計算しています。計算にはアインシュタインの縮約記法を使えるnumpy einsum を用いています。

49~50行目:逆短時間フーリエ変換で時間波形に戻しています。

52~53行目:波形データをWAVファイルに書き出しています。

シミュレーション実験

室内音響シミュレーションで作成した音にDSビームフォーマを作用させ、ボーカルの音だけ強調してみました。

方法

室内音響シミュレーションにはPythonのライブラリであるPyRoomAcousticsを用いました。

シミュレーション上で図3のように音源とマイクを配置して、64個のマイクで音楽データを録音しました。

図3:シミュレーションにおける音源とマイクの配置
図3:シミュレーションにおける音源とマイクの配置

残響時間は0.4秒、収録時のサンプリング周波数は16kHz に設定しました。

一応、以下をクリックすると、シミュレーションに用いた Python のコード room_sim.py が展開されます。

結果

DSビームフォーマを用いてボーカルを強調した結果は以下のようになりました。参考としてm=32のマイクで録音されたデータも載せておきます。

m=32のマイクの録音データ

処理結果

ボーカルの音だけ取り出されて、ドラムの音が減衰しているのがわかると思います。

補足:はじめ、マイクの個数 M=8 で実験を行ったのですが、効果は以下のようになんとなくわかる程度でした。

マイクの個数 M=8 の処理結果

人間の感覚の大きさは、受ける刺激の強さの対数に比例するというウェーバー・フェヒナーの法則があるので、音源を8倍大きくしても感じる音の大きさは8倍にはならないということかなと思います。

おわりに

DSビームフォーマを用いて特定方向の音源を強調してみました。DSビームフォーマの有効性を実感するのに64個のマイクを使う必要があったことには驚いています。DSビームフォーマは64個のマイクを1㎝間隔で配置する必要があるので、実際の環境で用いるのはかなり難しそうですね。

参考文献

  1. 浅野太、”音のアレイ信号処理 音源の定位・追跡と分離”、コロナ社、2011.
  2. 戸上真人、”機械学習実践シリーズ Pythonで学ぶ音源分離”、インプレス、2020.

使用した楽曲について

この記事で信号処理した楽曲は Cambridge Music Technology で提供されている AM Contra の Heart Peripheral を使用させていただきました。

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

ソースフィルタモデルで母音を音声合成

ソースフィルタモデルによる古典的な方法で母音 ”ieaou” を音声合成しました。今回の記事は主に以下の本を参考にしています。

ソースフィルタモデル

ソースフィルタモデルは参考文献 [1] の引用では、以下のようなことです。

ソースフィルタモデルは音声を、声帯振動あるいは呼気流などによる音源を入力として、声道フィルタによる系を経た出力信号とみなす。

図1は簡単な人の口の模式図ですが、声帯の振動を音源、声道による影響をフィルタとして考えるということです。このフィルタは声道フィルタと呼ばれています。

図1:口の模式図
図1:口の模式図

数式で説明すると、z変換された音声 Y は以下のように表せるということです。

\displaystyle{
Y(z) = H(z) X(z)
}

ここで、Xは声帯振動で得られた音源をz変換したもの、Hは声道フィルタの伝達関数です。

音声合成の手順

ソースフィルタモデルによる音声合成手順は以下です。

  1. 声帯振動を模擬した音源信号を作成
  2. 声道フィルタを作成
  3. 源信号に声道フィルタを畳み込む

今回は、参考文献 [1] を基に図2のブロック図で音声を合成します。

図2:母音合成のブロック図
図2:母音合成のブロック図

声帯振動は周波数F0の周期波形でモデル化します。また、3つのバンドパスフィルタ(BPF)を並列に置いたものを声道フィルタとします。

源信号を作成

今回、声帯振動波形は図3のような周期T0のパルス列でモデル化します。

図3:パルス列の声帯振動波形
図3:パルス列の声帯振動波形

声帯振動波形の周期の逆数は基本周波数と呼ばれ、参考文献 [3] によると、人間が知覚する声の高さにだいたい対応しているそうです。そのため、イントネーションの解析などに利用される重要なパラメータとなります。

今回使用した各母音の基本周波数は参考文献 [2] に記載されている女性28人の平均値となっています。その値は以下の表のようになります。

表:使用した基本周波数
/i/ /e/ /a/ /o/ /u/
235 Hz 223 Hz 212 Hz 216 Hz 231 Hz

声道フィルタを作成

声道で共振がおこる周波数はフォルマント周波数と呼ばれており、音声の音韻らしさを形作る重要なものとなっています。

声道フィルタはフォルマント周波数F1~F3を通す3つのバンドパスフィルタ(BPF)を並列に置いて作成します。

1つのバンドバスフィルタの伝達関数  H_n は以下のようになります。

\displaystyle{
H_n(z) = A_n \frac{\frac{\alpha^2+\omega_n^{2}}{\omega_n}z^{-1}}{1-2e^{-\alpha}\cos(\omega_n)z^{-1}+e^{-2\alpha}z^{-2}}
}

\displaystyle{
\begin{align}
\omega_n &=2\pi F_n/f_s \\
\alpha &=  \pi f_{bn}/f_s \\
A_n &= 10^{L_n/20}
\end{align}
}

ここで、 L_n は第 n フォルマントの振幅 [dB]、 f_{bn} は第 n フォルマントのバンド幅 [Hz]、 F_n は第 n フォルマント周波数 [Hz]、 f_s はサンプリング周波数です。

今回使用した各母音のフォルマント周波数、フォルマントの振幅についても参考文献 [2] に記載されている女性28人の平均値となっています。その値は以下の表のようになります。

表:使用したフォルマント周波数
/i/ /e/ /a/ /o/ /u/
F1 310 Hz 610 Hz 850 Hz 590 Hz 370 Hz
F2 2790 Hz 2330 Hz 1220 Hz 920 Hz 950 Hz
F3 3310 Hz 2990 Hz 2810 Hz 2710 Hz 2670 Hz
表:使用したフォルマント振幅
/i/ /e/ /a/ /o/ /u/
F1 -4 dB -2 dB -1 dB 0 dB -3 dB
F2 -24 dB -17 dB -5 dB -7 dB -19 dB
F3 -28 dB -27 dB -28 dB -34 dB -43 dB

また、使用したフォルマントのバンド幅については参考文献 [2] に記載されていた各母音の全平均値となっています。その値は以下の表のようになります。

表:用いたフォルマントバンド幅
F1 F2 F3
49.7 Hz 64.0 Hz 115.2 Hz

音韻/a/のときの声道フィルタ(並列に3つ並んだBPF)の周波数特性を図4に示します。

図4:/a/ の声道フィルタの周波数特性
図4:/a/ の声道フィルタの周波数特性

プログラム

ソースフィルタモデルによって ”ieaou” を合成するプログラムは以下です。

参考文献 [1] に掲載されているプログラムを基に作成しています。

from math import sin, cos, exp
import soundfile as sf
import numpy as np
import scipy.signal as sg

fs = 48000        # サンプリング周波数
mora_per_sec = 1  # 話す速度(モーラ/秒)
gain = 0.5        # 音の大きさ  

sec = (1/mora_per_sec)*5  # 音源の長さ
T   = 1.0/fs           # サンプリング周期
nsample = int(fs*sec)  # データ数

# フォルマント周波数 i,e,a,o,u
F = [[235,   223,  212,  216,  231], # F0
     [310,   610,  850,  590,  370], # F1
     [2790, 2330, 1220,  920,  950], # F2
     [3310, 2990, 2810, 2710, 2670]] # F3

# フォルマント振幅 [dB] i,e,a,o,u
L = [[  0,   0,   0,   0,   0],
     [ -4,  -2,  -1,   0,  -3],  # F1
     [-24, -17,  -5,  -7, -19],  # F2
     [-28, -24, -28, -34, -43]]  # F3

F = np.array(F, dtype=np.float64)
L = np.array(L, dtype=np.float64)
L = L + 12  # 振幅特性が小さすぎたので、大きくする
L = 10**(L/20.0)  # log->linear

band  = np.array([0, 50.0, 64.0, 115.0]) # 帯域幅
alpha = band*np.pi*T

sample_per_mora=int(fs/mora_per_sec) # 1モーラのサンプル数
x = np.zeros(sample_per_mora)
y = np.zeros(nsample, dtype=np.float64)

for i in range(5):

    # 声帯振動を作成
    T0 = int(fs/F[0,i])  # 基本周期(サンプル数)
    x[:] = 0
    x[0::T0] = gain

    # 音韻の始まりと終わり
    start = sample_per_mora*i
    end   = sample_per_mora*(i+1)

    # フォルマント周波数の共振フィルタをかける
    for k in range(1,4):
        w  = 2.0*np.pi*F[k,i]*T
        b = np.array([0.0, 0.0, 0.0])
        a = np.array([1.0, 0.0, 0.0])
        b[1] = ((alpha[k]**2+w**2)/w)*sin(w)*exp(-alpha[k])
        b = b * L[k,i]
        a[1] = -2*exp(-alpha[k])*cos(w)
        a[2] = exp(-2.0*alpha[k])           
        y[start:end] = y[start:end] + sg.lfilter(b, a, x)
    
sf.write("ieaou.wav", y, fs, subtype='PCM_16')

6~8行目:サンプリング周波数、1秒あたりに話すモーラ数、声帯振動波形のパルスの大きさを設定しています。

28行目:合成信号が小さくなりすぎたので、フィルタの振幅特性を12dB(4倍)大きくしています。

29行目:記載している振幅値が log 表記なので、linear に変換しています。

40~43行目:パルス列の声帯振動波形を作成します。

49~58行目:3つのBPFを作成して、58行目でBPFを畳み込んだものを足し合わせています。

60行目:合成した音声をwavで出力しています。

結果

パルス列でモデル化した声帯振動信号と声道フィルタを畳み込んで母音合成した信号は以下です。

声帯振動信号

合成音声(声道フィルタを畳み込んだ信号)

たしかに、"ieaou" に聞こえるようになりました!!!

おわりに

ソースフィルタモデルを使って、母音 "ieaou" を音声合成しました。思っている以上の母音が合成できて正直驚いています。

参考文献

  1. 小坂直敏,”サウンドエフェクトのプログラミング Cによる音の加工と音源合成”,株式会社オーム社,2012.
  2. 電子情報通信学会,”聴覚と音声”,pp.350-371,コロナ社,1973.
  3. 森勢将雅,”音声分析合成”,コロナ社,2018.

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

© 2021 Setoti All rights reserved.