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

音声分析合成

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

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

音声音響信号処理の基礎と実践 - フィルタ,ノイズ除去,音響エフェクトの原理 - (次世代信号情報処理シリーズ 2)

作者:川村新
コロナ社

Amazon

ヘリウムボイスの原理
ヘリウムボイスの作成方法
- スペクトル包絡の伸縮
- ヘリウムボイス変換の流れ
プログラム
処理結果
おわりに

ヘリウムボイスの原理

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

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

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

$\displaystyle{ c = f \lambda }$

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

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

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

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

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

スペクトル包絡の伸縮

ヘリウムボイスの原理からフォルマント周波数をｒ倍、つまり対数スペクトル包絡 $H[k]$ を図２のように周波数方向に 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 としたときが実際のヘリウムボイスに近いみたいです。

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

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

フレームごとにケプストラム分析で対数スペクトル包絡を求めて、ｒ倍に伸縮しています。抽出していなかった高次のケプストラムについても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/

2021-10-27

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

音声分析合成

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

ケプストラム分析
プログラム
ケプストラム分析の結果
おわりに

ケプストラム分析

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

スペクトル包絡とは

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

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

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

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

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

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

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

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

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

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

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

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

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

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

一方、声道フィルタの対数振幅スペクトルは緩やかな形状を持ちますので、さらにフーリエ変換（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したものをケプストラムとします。

リフタリング

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

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

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

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

プログラム

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

解析するWAVファイルをソースコードと同じディレクトリに入れて、wav_nameとstartを書き換えて、「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] に記載されていた男性33人のフォルマント周波数 F1～F3 の平均です。

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

おわりに

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

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

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

2021-10-24

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

音声・音響処理

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

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

頭部伝達関数（HRTF）
頭の周りをまわる音源の作成
処理結果
おわりに

頭部伝達関数（HRTF）

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

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

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

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

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

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

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

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

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

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

音源の位置の計算

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

音源が図２のように角速度 $\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$ の２つの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 とする。

プログラム

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と４に対応する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行目：ノーマライズして、波形データをファイルに書き込んでいます。

処理結果

頭部伝達関数を使って作成した８D立体音響は以下のようになります。仰角は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

2021-10-11

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

音源分離・音源定位

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

音のアレイ信号処理 (音響テクノロジーシリーズ)

作者:浅野　太
コロナ社

Amazon

遅延和ビームフォーマ
プログラム
シミュレーション実験
- 方法
- 結果
おわりに
参考文献
使用した楽曲について

遅延和ビームフォーマ

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

概要

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

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

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

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

直線状アレイの遅延

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

図２の $\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} {\bf 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} }$

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

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

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

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

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

${\bf 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を用いました。

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

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

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

room_sim.py（クリックで展開）

import numpy as np
import pyroomacoustics as pra
import soundfile as sf
from librosa import resample

fs = 16000  # リサンプリング後の周波数 
rt60_tgt = 0.6    # 残響時間
room_size = [10, 10, 3.5]  # 部屋の大きさ
SNR = 90          # SN比
dir_src = "src/"  # 音源のフォルダ
dir_mic = "mic/"  # 録音したデータのフォルダ
src_list = ["Vox.wav", "Synth.wav", "drum.wav"]

# 音源とマイクの位置を決定
center = [5.0,5.0]  # 中心の位置 
s_num  = 3     # 音源の数
phi0   = 90    # 一個目の音源の方向 [度]
radius = 3     # 音源の中心から距離 [m]
m_num  = 64    # マイクの数
d      = 0.01  # マイクの間隔 [m]
# 音源を円状アレイで配置
phi0 = np.radians(phi0) # deg -> rad 変換
s_locs = pra.circular_2D_array(center, s_num, phi0, radius)
s_locs = np.insert(s_locs, 2, 1.25, axis=0)
# マイクを線状アレイで配置
m_locs = pra.linear_2D_array(center, m_num, 0, d)
m_locs = np.insert(m_locs, 2, 1.25, axis=0)

# WAVファイルを読み込む
audio = []
for i, src_name in enumerate(src_list):
    x, fs_tmp = sf.read(dir_src+src_name)
    x = resample(x, fs_tmp, fs)
    audio.append(x)
audio = np.array(audio)

# 残響時間と部屋の大きさから吸音率と反射回数を決定
e_absorption, max_order = pra.inverse_sabine(rt60_tgt, room_size)

# 部屋の作成
room = pra.ShoeBox(
    room_size, fs=fs, materials=pra.Material(e_absorption), max_order=max_order
)

# 部屋に音源を配置
for i in range(s_num):
    room.add_source(s_locs[:,i], signal=audio[i,:], delay=0)

# 部屋にマイクを配置
room.add_microphone_array(
    pra.MicrophoneArray(m_locs, fs=fs)
)

# シミュレーションする
room.simulate(snr=SNR)

# マイクで録音した音を書き出す
mic_data = room.mic_array.signals
mic_data = mic_data/np.max(np.abs(mic_data)) # ノーマライズ
for i in range(m_num):
    name = dir_mic+"mic"+str(i)+".wav"
    sf.write(name, mic_data[i,:], fs, subtype='PCM_16')

結果

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

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

処理結果

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

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

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

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

おわりに

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

参考文献

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

使用した楽曲について

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

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

2021-09-27

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

音声分析合成

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

サウンドエフェクトのプログラミング　Cによる音の加工と音源合成

作者:東京電機大学小坂直敏
オーム社

Amazon

ソースフィルタモデル
音声合成の手順
- 音源信号を作成
- 声道フィルタを作成
プログラム
結果
おわりに
参考文献

ソースフィルタモデル

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

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

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

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

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

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

音声合成の手順

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

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

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

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

音源信号を作成

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

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

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

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

声道フィルタを作成

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

声道フィルタはフォルマント周波数F1～F3を通す３つのバンドパスフィルタ（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）の周波数特性を図４に示します。

プログラム

ソースフィルタモデルによって ”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" を音声合成しました。思っている以上の母音が合成できて正直驚いています。

参考文献

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

2021-09-21

双一次変換によるIIRフィルタ設計

音声・音響処理

カットオフ周波数を決めれば、単純な計算だけでIIRフィルタの係数が求まる方法はないかと思っていたら、双一次変換と呼ばれる方法があったので紹介します。

今回の記事は以下の本を主に参考にしています。

ディジタルフィルタ原理と設計法 (設計技術シリーズ)

作者:陶山健仁
科学情報出版株式会社

Amazon

双一次変換によるIIRフィルタ設計手順
プログラム
周波数特性について
おわりに
参考文献

双一次変換によるIIRフィルタ設計手順

双一次変換によるIIRフィルタの設計手順は以下です。

アナログ基準ローバスフィルタを設計する。
周波数変換で所望のカットオフ周波数 $\omega_c$ にする。
アナログフィルタに双一次変換を行い、伝達関数 $H(z)$ を求める。

今回は、手始めに以下のような2次のIIRフィルタの伝達関数を作ります。

$\displaystyle{ H(z) = \frac{b_0+b_1z^{-1}+b_2z^{-2}}{a_0+a_1z^{-1}+a_2z^{-2}} }$

アナログ基準ローバスフィルタの設計

カットオフ角周波数が $\omega_c=1$ であるローバスフィルタ（LPF）を基準LPFと呼び、最初にアナログ基準LPFの伝達関数を設計します。

ちなみに、デジタル信号処理の世界では暗黙のうちに、周波数 $f'$ をサンプリング周波数 $f_s$ で割った正規化周波数 $f'/f_s$ を周波数 $f$ と呼びます。同様に、角周波数についても正規化角周波数 $2\pi f'/f_s$ を角周波数 $\omega$ と呼びます。

アナログフィルタの設計では、 $s=j\omega$ を代入すると以下の理想LPFの周波数特性を近似するような有理関数 $H(s)$ を求めます。

$\displaystyle{ H_{id}(\omega) = \left\{ \begin{array}{clc} 1 & \omega \lt 1 \\ 0 & \omega \geq 1 \end{array} \right. }$

有理関数というのは以下のように表せる関数です。

$\displaystyle{ H(s) = \frac{b_0s^m+b_1s^{m-1}+\cdots+b_m}{s^{n}+a_1s^{n-1}+\cdots + a_n} }$

直接周波数特性を近似するのは難しい（周波数特性は複素数になる）ので、まず理想の2乗振幅特性を有理関数で近似します。アナログフィルタの代表であるバタワースフィルタでは以下の $|H_n(\omega)|^2$ で近似させています。

$\displaystyle{ |H_n(\omega)|^2 = \frac{1}{1+\omega^{2n}} }$

図１をみると、確かに理想LPFに近い特性となっています。

2乗振幅特性は有理関数で表せたので、つぎに伝達関数 $H(s)$ を求めていきます。

$|H_n(\omega)|^2$ の平方根を求めて $\omega=s/j$ と置き換えればいいと思いますが、有理関数にならないので工夫が必要です。

$|H_n(\omega)|^2$ が以下のように表せることを利用します。

$\displaystyle{ \begin{align} |H_n(\omega)|^2 &= H_n(\omega) H_n^{*}(\omega) \\ &= H_n(\omega) H_n(-\omega) \end{align} }$

$H_n^{*}(\omega)=H_n(-\omega)$ は実数をフーリエ変換したときの性質となっています。この性質がないと周波数特性を逆フーリエ変換した結果が複素数になってしまいます。

$|H_n(\omega)|^2$ に $\omega=s/j$ とおくと、

$\displaystyle{ \begin{align} |H(s)|^2 &= \frac{1}{1+(s/j)^{2n}} \\ &= \frac{1}{1+(-1)^n s^{2n}} \end{align} }$

$|H(\omega)|^2$ の極を求めてみると、

n が奇数のときは

$\displaystyle{ \begin{align} 1 - s^{2n} &= 0 \\ s &= e^{j2\pi k/2n} \end{align} }$

n が偶数のときは

$\displaystyle{ \begin{align} 1 + s^{2n} &= 0 \\ s &= e^{j(2k+1)\pi/2n} \end{align} }$

よって、n=2、n=3のときの $|H(s)|^2$ の極の位置は図２のようになります。

$|H(s)|^2=H(s)H(-s)$ という条件があるので、この条件を満たすように $H(s)$ の極を選ぶ必要があります。加えて、安定なフィルタを設計するには s 平面の左半平面に極が位置する必要があるので、 $H(s)$ の極は左半平面にある極全部となります。

したがって、n=2のときの伝達関数 $H_2(s)$ は以下のようになります。

$\displaystyle{ \begin{align} H_2(s) &= \frac{1}{(s-e^{j3\pi/4})(s-e^{j5\pi/4})} \\ &= \frac{1}{s^2+\sqrt{2}s + 1} \end{align} }$

n=1 のときの伝達関数 $H_1(s)$ も求めてみると以下のようになります。

$\displaystyle{ H_1(s) = \frac{1}{s+1} }$

これで、アナログ基準 LPF の設計ができました。

周波数変換

フィルタのタイプは図３のように４つあります。

カットオフ角周波数が $\omega_c$ である各タイプのフィルタを設計したい場合は表のように s を変換すればよいです。

タイプ	変換
LPF	$s\rightarrow \frac{s}{\omega_c}$
HPF	$s\rightarrow \frac{\omega_c}{s}$
BPF	$s\rightarrow \frac{s^2+\omega_o^2}{s\omega_b}$
BRF	$s\rightarrow \frac{s\omega_b}{s^2+\omega_o^2}$

ここで、 $\omega_b=\omega_{c2}-\omega_{c1}$ です。

$H_2 (s)$ を以上のLPF変換、HPF変換すると以下のようになります。

$\displaystyle{ \begin{align} H_{LPF}(s) &= \frac{\omega_c^2}{s^2+\sqrt{2}\omega_c s + \omega_c^2} \\ H_{HPF}(s) &= \frac{s^2}{s^2+\sqrt{2}\omega_c s + \omega_c^2} \end{align} }$

また、 $H_1 (s)$ を以上のBPF変換、BRF変換すると以下のようになります。

$\displaystyle{ \begin{align} H_{BPF}(s) &= \frac{\omega_b s}{s^2+\omega_b s + \omega_o^2} \\ H_{BRF}(s) &= \frac{s^2+\omega_o^2}{s^2+\omega_b s + \omega_o^2} \end{align} }$

双一次変換

双一次変換は s 変数を以下の式で z 変数にする変換をいいます。

$\displaystyle{ s=2\frac{1-z^{-1}}{1+z^{-1}} }$

このとき、図４のようにアナログ角周波数 $\omega_a : [-\infty, \infty]$ とデジタル角周波数 $\omega_d : [-\pi, \pi]$ が対応するようになっています。

双一次変換の式に、 $s=j\omega_a$ 、 $z=e^{j\omega_d}$ を代入すると関係式が以下のようにわかります。

$\displaystyle{ \omega_a=2 \tan{\frac{\omega_d}{2}} }$

バタワースフィルタの伝達関数を双一次変換することで、 $H(z)$ を得ることができます。ただし、双一次変換によりカットオフ周波数の位置もずれてしまうため、アナログフィルタ設計の時点で予め考慮する必要があります。

$\omega_a$ と $\omega_d$ の関係式より、デジタルフィルタのカットオフ角周波数を $\omega_c$ としたいときは、アナログフィルタの設計時ではカットオフ角周波数を以下の $\omega'_c$ にすればよいです。

$\displaystyle{ \omega'_c=2 \tan{\frac{\omega_c}{2}} }$

このことはプリワーピングと呼ばれています。

フィルタ係数の決定

双一次変換をバタワースフィルタの伝達関数に行った結果、各フィルタ係数は以下のようになります。 $\omega_c$ はプリワーピングしたカットオフ角周波数です。

表：LPFのフィルタ係数
$b_0$	$\omega_c^2$
$b_1$	$2\omega_c^2$
$b_2$	$\omega_c^2$
$a_0$	$\omega_c^2+2\sqrt{2}\omega_c+4$
$a_1$	$2\omega_c^2-8$
$a_2$	$\omega_c^2-2\sqrt{2}\omega_c+4$

表：HPFのフィルタ係数
$b_0$	$4$
$b_1$	$-8$
$b_2$	$4$
$a_0$	$\omega_c^2+2\sqrt{2}\omega_c+4$
$a_1$	$2\omega_c^2-8$
$a_2$	$\omega_c^2-2\sqrt{2}\omega_c+4$

表：BPFのフィルタ係数
$b_0$	$2\omega_b$
$b_1$	$0$
$b_2$	$-2\omega_b$
$a_0$	$\hspace{15px} \omega_o^2+2\omega_b+4 \hspace{15px}$
$a_1$	$2\omega_o^2-8$
$a_2$	$\omega_o^2-2\omega_b+4$

表：BRFのフィルタ係数
$b_0$	$\omega_o^2+4$
$b_1$	$2\omega_o^2-8$
$b_2$	$\omega_o^2+4$
$a_0$	$\hspace{15px} \omega_o^2+2\omega_b+4 \hspace{15px}$
$a_1$	$2\omega_o^2-8$
$a_2$	$\omega_o^2-2\omega_b+4$

プログラム

Python で記述したIIRフィルタの係数と周波数特性を求めるプログラムは以下です。

import numpy as np
import scipy.signal as sg
import matplotlib.pyplot as plt
import sys

# コマンドラインの引数を取得
args = sys.argv
if len(args)!=4:
    print("\nusage: python main.py  fs  fc  type")
    print("   fs        : sampling frequency")
    print("   fc        : cutoff frequency")
    print("   type      : filter type (lpf, hpf, bpf, brf)")
    raise Exception("Argument error ")

fs = int(args[1])    # サンプリング周波数
fc = float(args[2])  # カットオフ周波数
ftype = args[3]      # フィルタの種類　LPF, HPF, BPF, BRF
fb = 5000            # バンドパスの幅

omega_b  = 2*np.pi*(fb/fs)  # バンドパスの幅
omega_c  = 2*np.pi*(fc/fs)    # カットオフ正規化角周波数
omega_c  = 2*np.tan(0.5*omega_c)  # プリワーピング
omega_c1 = omega_c-0.5*omega_b
omega_c2 = omega_c+0.5*omega_b
omega_o  = np.sqrt(omega_c1*omega_c2) 
a = [0,0,0]
b = [0,0,0] 

# フィルタ係数決定
if(ftype=="lpf"):
    b[0] = omega_c**2
    b[1] = 2*omega_c**2
    b[2] = omega_c**2
    a[0] = omega_c**2+2*np.sqrt(2)*omega_c+4
    a[1] = 2*omega_c**2-8
    a[2] = omega_c**2-2*np.sqrt(2)*omega_c+4
elif(ftype=="hpf"):
    b[0] = 4
    b[1] = -8
    b[2] = 4
    a[0] = omega_c**2+2*np.sqrt(2)*omega_c+4
    a[1] = 2*omega_c**2-8
    a[2] = omega_c**2-2*np.sqrt(2)*omega_c+4
elif(ftype=="bpf"):
    b[0] = 2*omega_b
    b[1] = 0
    b[2] = -2*omega_b
    a[0] = omega_o**2+2*omega_b+4
    a[1] = 2*omega_o**2-8
    a[2] = omega_o**2-2*omega_b+4
else: # brf
    b[0] = omega_o**2+4
    b[1] = 2*omega_o**2-8
    b[2] = omega_o**2+4
    a[0] = omega_o**2+2*omega_b+4
    a[1] = 2*omega_o**2-8
    a[2] = omega_o**2-2*omega_b+4

# 周波数特性の出力
w, h = sg.freqz(b, a)
fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(w, np.abs(h), c="red")
ax.set_ylim([0,1.2])
ax.set_xlabel("omega") 
ax.set_ylabel("Amplitude [dB]")
fig.savefig(ftype+".png")

以下のようにコマンドを実行することで、サンプリング周波数 fs=48kHz、カットオフ周波数 fc=12kHz の LPF のフィルタ係数と振幅特性が求まります。BPF や BRF を指定したときはバンド幅は 5kHz で、第二引数の値がバンドの中心となります。

  python main.py 48000 12000 lpf

周波数特性について

サンプリング周波数 fs=48kHz、カットオフ周波数 fc=12kHz を指定してさきほどのプログラムを実行した結果、図５のような特性になりました。

あまり急峻な特性ではないですが、２次の IIR フィルタなので仕方ないですかね。

白色雑音に対して各フィルタを畳み込んだときのスペクトログラムは図６のようになります。

やはり急峻ではないですが、上手く機能しているのがわかると思います。

おわりに

長くなりましたが、双一次変換によるIIRフィルタ設計についてまとめました。FPGAでフィルタを作成するときなどに利用しようと思います。

参考文献

陶山健仁、”ディジタルフィルタ原理と設計法”、科学情報出版株式会社、2018．
尾知博、”ディジタル・フィルタ設計入門各種フィルタの原理とDSPによる実現”、CQ出版社、1990．

2021-09-12

RX651用のマウス基板の設計とRX631との違い

マイクロマウス

いままでRX631を載せたマウス基板を用いて、デバイスの動作確認をおこなっていましたが、AD変換が上手くいかなかったので、マウス基板をまた作りなおしました。

RX631-48ピンが手に入りづらくなっていたことに加えて、RAMが多い高性能のマイコンを使いたいと思ったので、今回はRX651-64ピン用のマウス基板を設計しました。

RX651用のマウス基板の回路設計
- 回路図の変更点
- 配線図の変更点
届いたマウス基板
RX651とRX631の違い
センサーホルダーでの失敗
おわりに
参考文献

RX651用のマウス基板の回路設計

赤外線センサーのAD変換結果がおかしな値となっていたため、RX651のユーザーズマニュアルハードウェア編を参考にして回路を再設計をしました。

回路図の変更点

再設計した回路図は図１です。

回路図PDF版

今回の回路図で前回の回路図から変更した主な箇所は以下の３つです。

マイコンをRX631からRX651に変更
1つのFETで1つの赤外線LEDをスイッチング
マウス基板とAS5047P基板をつなげるコネクタを4x2から3x2に変更

冒頭で述べたようにRX631からRX651に変更しています。

加えて、いままで1つのFETで2つの赤外線LEDをスイッチングしていましたが、プログラムが書きにくかったので、1つのFETで1つの赤外線LEDをスイッチングするようにしました。

あとは、AS5047PのA端子とB端子はマイコンと接続しなくてもよいことがわかりましたので、コネクタを4x2から3x2に変更しています。

DRV8836やMPU6000も手に入りづらくなっているので、変更しようかなと思いましたが、使い方を再び調査しないといけないので辞めときました。

配線図の変更点

再設計した配線図は以下です。

配線図PDF版（縮尺：2.5倍）

今回の配線図で前回の配線図から変更した主な箇所は以下の３つです。

マイコンのAVSSや赤外線センサにつながっているGNDはレギュレータ付近のコンデンサのGNDと一点接続する
赤外線センサーの出力信号とデジタル信号線を交差させない
モータドライバーのGNDはレギュレータ付近のコンデンサのGNDと一点接続する

変更点１と２についてはRX651のマニュアルp.2548にある「53.6.12 ボード設計上の注意」に書いてあることを参考にしています。いままでは、プリント基板の2層ともベタGNDにしていましたが、1層だけベタGNDにして、GNDを１点接続できるようにしました。

また、モータドライバーからのノイズで問題が生じるというブログ記事をいくつか拝見しましたので、モータドライバーのGNDについてもレギュレータ付近のコンデンサのGNDと一点接続させました。

届いたマウス基板

Elecrowに発注して届いたマウス基板が図３です。

文字は少し乱れていましたが、フットプリントがとても綺麗にできていて良かったです。

ハンダ付けして組み立てたマウスが図４となります。

前回のマウス基板では、センサーの出力はオシロコープで正常だったにもかかわらず、AD変換の結果がずっと4095になっていたり、電圧が上がると逆にAD変換の結果が下がったりして、謎の動作をしていました。

今回のマウス基板では、AD変換の結果が正常な値でしたので、とても良かったです。

RX651とRX631の違い

RXマイコン（RX631）のソフト開発で記述したプログラムを書き換えました。その過程で気づいたRX631と違った点を述べていきます。

性能

前回購入したRX631と今回購入したRX651の性能の違いは以下の表のようになっています。システムクロック（ICLK）というのはRXマイコンのCPU動作クロックの呼び名です。

	RX631	RX651
型名	R5F5631PDDFL	R5F5651EHDFM
ピンの数	48	64
RAM容量	64KB	640KB
ROM容量	512KB	2MB
ICLK周波数	100MHz	120MHz

全体的にメモリの容量が多くなっているので、プログラムや計算結果などを多く保存できます。

また、 ICLK周波数も120MHzに上がっているので、計算速度が速くなってます。

書き込みにメインクロックがいらない

以下の記事で説明しましたが、RX631では書き込みにメインクロックが必要でした。

しかし、RX651ではSCI (UART) で書き込む場合、メインクロックがいらなくなりました。試しにセラロックをつけずに書き込んだところ、無事書き込めました。

また、PLL回路の入力クロックソースに高速オンチップオシレータ（HOCO）を選択できるようになりましたので、メインクロックがなくてもICLKを最大周波数で動作できます。

これでRXマイコンの書き込みのハードルも下がったかなと思います。

プログラム的な変更点

RX631からRX651でプログラムを書き換えたときの主な変更点を説明します。

ROMWTレジスタの設定

RX651からICLKを50MHz以上で動かす場合、ROM ウェイトサイクル設定レジスタ (ROMWT)の設定が必要となりました（マニュアルp.312「9.2.2 ROM ウェイトサイクル設定レジスタ (ROMWT)」を参照）。

そのため、メインクロックの設定とICLKの設定の間にROMWAITレジスタの設定を以下のように加えました。

/*---- MainCLK, SUBCLK and RTC Initialization ----*/
void init_clock(void){
     ～略～
    // Main CLK setting
    RTC.RCR4.BIT.RCKSEL     = 1;    // Select main CLK as count source
    SYSTEM.MOSCCR.BIT.MOSTP = 0;    // Main CLK ON
    for(i=0;i<1000;i++);            // Wait about 50ms

    // ROMWAIT setting
    SYSTEM.ROMWT.BIT.ROMWT  = 2;
    while(SYSTEM.ROMWT.BIT.ROMWT!=2);  // ROM wait cycle = 2cycle

    // SYSTEM CLK setting
    SYSTEM.PLLCR.BIT.PLIDIV = 0x0;  // Main CLK divided by 1
    SYSTEM.PLLCR.BIT.STC    = 0x17; // Main CLK multiplied by 12
     ～略～
}

MTU3とRSPIがPCLKAを使用

RX631ではMTU3とRSPIはクロックにPCLKBを使用していましたが、RX651ではPCLKAを使用しています。それなので、MTU3とRSPIがPCLKBで動いていると勘違いしないように注意する必要があります。

スピーカーを鳴らしたときに音が高くなっていたので、気づくことができました。

センサーホルダーでの失敗

図５のセンサーホルダーを3Dプリンタで作成したのですが、失敗したので報告します。

赤外線LEDを光らせて、センサーホルダーを付けたときのAD変換値をみてください。

表：センサーホルダーありのＡＤ変換値
障害物	左前	左横	右前	右横
無	3607	3634	3587	3564
有	3689	3747	3661	3673

障害物がセンサーの前に無いときと有るときでAD変換値の差がほとんどありません。

おそらく、原因としてはセンサーホルダーによる散乱光だと考えられます。それにしても、散乱光が強すぎると思いますが....。

センサーホルダーを使う人で、散乱光の対策はどうしているのかと気になったところ以下の記事が見つかりました。

見つけた記事：うむ夫の歩み Nucleo-32boardを使ったクラシックマウスの開発その8 ～センサーホルダー～

どうやら、熱収縮チューブを付けてホルダーに挿入すれば良いみたいです。ただ、作成したホルダーの穴の大きさが熱収縮チューブの分を考慮していないので、作り直す必要があります。

試しにセンサーホルダーを外して、赤外線LEDに熱収縮チューブを付けた場合のAD変換値は以下です。

表：センサーホルダーなし、熱収縮チューブありのＡＤ変換値
障害物	左前	左横	右前	右横
無	350	343	278	574
有	1916	2078	2316	2951

AD変換の値に差があっていい感じです。

センサーホルダーを付けなくてもいいかなと少し思いましたが、基板にホルダーを差し込む穴もあるので、もう一度作成しなおそうと思います。

おわりに

今回、RX651用のマウス基板の設計とRX651とRX631の違いについて話しました。これで全てのデバイスの動作確認ができましたので、次回から本格的にマイクロマウスのプログラムを作成していきたいと思います。