# [論文]時間領域差分法に基づく1次元の声道音響モデルを読んでまとめる

## リンク

[[論文]時間領域差分法に基づく1次元の声道音響モデルで使う定数](/hM0BXlh8QlG4MWHN6VleLg)

## 論文の理解




### はじめに
声帯音源波が声道を通過して音声が生成される過程を時間領域でシミュレーションする手法はあるが、アルゴリズムが複雑でプログラミングが難しい

よって、声道形状からフィルタを生成し，声帯音源波に畳み込んで音声を生成することにする

声門から口唇までの断面積の変化を表す声道断面積関数に基づき，声道を縦続音響管でモデル化する。
次に，各管の入口と出口の音圧と体積速度を周波数領域で順次計算し，最終的に声道の伝達関数を計算する(Transmission Line Model: 以降TLM法)。
そして，これを逆離散フーリエ変換してインパルス応答に変換し，フィルタとして声帯音源波に畳み込む。
この方法は，周波数依存性のある壁振動，熱交換，粘性抵抗などによる損失を考慮できるため，声道の伝達関数を精度よく計算できる。

↑これは本来やる計算法（今回関係ない）

だが、これも難しいので、母音の生成過程に限れば，==1次元の時間領域差分法==(Finite-Difference Time-Domain method:FDTD法)で簡単に実現することにする。

### モデル

#### 支配方程式とは
連続の式と運動方程式で立てられる。
音圧・粒子速度・音速・密度・k=1/pc^2^・減衰定数

#### 支配方程式の差分化
1次元(x方向のみ)の音響 FDTD 法の支配方程式は,

- 音圧をp
- 粒子速度をu
- 音速をc
- 密度をρとして、κ = 1/ρc^2^
- 時間離散化間隔 ∆t
- 空間離散化間隔 ∆x
- 減衰定数 α,α^∗^
- 例）𝛼 = 2.0 × 10^−3^ ，𝛼^∗^ = 0

とおく。

δ = デルタ , κ = カッパ , ρ = ロー
κ = 1/ρc^2^

-κ＊δ音圧/時間離散化間隔-減衰定数＊音圧＝δ粒子速度/空間離散化間隔 (1)連続の式
-密度δ＊δ粒子速度/時間離散化間隔-減衰定数^＊^ ＊ 粒子速度＝δ音圧/空間離散化間隔 (2)運動方程式

式(1), (2) を差分化すると式(3), (4)が得られる。
カエル跳び差分アルゴリズム

#### 入力体積速度


音圧[n]=

式(3)が示すように,声道の入力端である第1セクションの音圧p[1]の計算には**定義されていないU[0]が必要**なので,これを入力値として外部から与える。
第m計算ステップで入力する体積速度をU~in~[m]としてU[0]に代入する。
ここで,シミュレーションのステップ数をMとすると, ==m = 1, 2, ... , M== である。

声道の伝達関数を計算する場合にはインパルスを入力すればよい。
mは 1以上で定義されているので, ==U~in~ [m]はm = 1で 1,それ以外では 0== とする。これによりインパルス応答が計算できる。


そして,これを離散フーリエ変換することにより伝達関数が得られる。
なお,母音を生成する場合にはU~in~として Rosenberg波などの声帯の音源波形を与えればよい。

#### 放射インピーダンス
式(4)で示すように,声道の出力端である第Nセクションの体積速度U[N]の計算には,**定義されていないp[N + 1]が必要**となるため,このままでは計算できない。
そこで,第Nセクションの体積速度U[N]と音圧p[N]を関係付ける放射インピーダンスを導入する。

- 出力端の断面積をA~out~
- 出力端の音圧をp~out~
- 出力端の体積速度をU~out~


第m計算ステップでは,
- A~out~ = A[N]
- U~out~[m] = U[N]
- p~out~[m] = p[N]

なお,声道は時不変なので==A~out~はmによらない定数==である。

よって, U~out~[m],p~out~[m]は最も単純な近似計算として以下の式で表される。

### 計算手順
以上を踏まえて,計算手順の概略を述べる。
プログラムには,

- 声道断面積関数A,
- 空間離散化間隔(断面積の計測間隔)∆x
- 時間離散化間隔∆t
- 入力体積速度U~in~

を入力するものとする。

なお,これらを入力する前に,シミュレーションの安定条件として ==∆t ≤ c⁄∆xを満たすように∆tを設定==する。
また,ミュレーションの==継続時間をd s== とすると時間離散化間隔∆t,==ステップ数M = d⁄∆t==でU~in~を作成しておく必要がある。

#### 1. プログラムの冒頭で音速などの定数を設定する。
#### 2. 断面積関数を読み込み,必要があれば単位の変換を行う。
#### 3. 計算に必要なpとU(セクション数=N), p~out~(ステップ数=M)のメモリを確保する。→zeros()
#### 4. 以下の計算を実行する。

for m=1~M
    式(3)によるUin[m]から p[1]の計算
    
    for n=2~ N
        式(3)によるp[n]の計算
    end
    
    for n=1~ N − 1
        式(4)によるU[n]の計算
    end
    
    式(6)によるUout[m] = U[N]の計算
end

これを MATLAB でコーディングしたところ,空行,コメント行を含めても 65 行,含めなければ 34 行であった。

## わかったこと
- ==1次元のFDTD法を用いて縦続音響管でモデル化した声道内における音波の伝搬をシミュレーションする。==
- 手法は短いコードでインパルス応答を計算し,声帯の音源波から音声を生成することができる。
- ある特定の周波数で励振したときの瞬時体積速度・音圧の分布を可視化できる
- 声帯振動モデルと結合させて声帯と声道の相互作用も検討できる。

## わからないこと
- インパルス応答(音響工学でやった気がするけど覚えてない)

## 足りていない情報・材料
- MRIデータから抽出した声道断面積関数
- U~in~(体積速度) がない
- Rosenberg波などの声帯の音源波形を与えれば良いらしい
- 音速
- 毎秒340メートル？