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

ビームフォーミングで特定方向の音源を強調しました。ビームフォーミングとは所定の方向に波(電波、音波など)の指向性を高める技術のことです(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

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

© 2021 Setoti All rights reserved.