卒業論文

機械学習を用いた音声の伝達特性からの声道断面積の自動抽出

北海道科学大学 工学部
情報工学科

2-15-3-033
指導教員

松﨑 博季

2019 年 (平成 31 年)2 月

目次

第 1 章 緒言1

第 2 章 音声生成機構

2

2.1

声道断面積関数および声道伝達特性、音源フィルタ理論について . . . . . . .

2

2.1.1

音源フィルタ理論 . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

2

2.1.2

声道断面積関数 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

2

2.1.3

声道伝達特性 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

3

第 3 章 声道形状を考慮した多量の声道断面積関数の自動作成について

4

3.1

まえがき . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

4

3.2

声道断面積関数のセクション数について . . . . . . . . . . . . . . . . . . . .

4

声道形状の決定要因 . . . . . . . . . . . . . . . . . . . . . . . . . . .

4

声道断面積関数の導出方法 . . . . . . . . . . . . . . . . . . . . . . . . . . . .

4

3.3.1

調音位置の導出方法 . . . . . . . . . . . . . . . . . . . . . . . . . . .

5

3.3.2

調音位置における声道断面積. . . . . . . . . . . . . . . . . . . . . .

5

3.3.3

多量の声道断面積関数の求めるプログラム . . . . . . . . . . . . . . .

5

第 4 章 電気等価回路モデルを用いた声道断面積関数からの声道伝達特性の導出方法

7

4.1

まえがき . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

7

4.2

声道伝達特性の導出方法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

7

電気的等価回路モデル . . . . . . . . . . . . . . . . . . . . . . . . . .

7

3.2.1
3.3

4.2.1

第 5 章 機械学習による声道断面積の出力方法および結果と考察

9

5.1

まえがき . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

95.2

ニューラルネットワークの設定 . . . . . . . . . . . . . . . . . . . . . . . . .

9

5.2.1

ニューラルネットワークの概略図 . . . . . . . . . . . . . . . . . . . .

9

5.2.2

教師あり学習 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10

5.2.3

ノード数、学習係数、学習回数およびバッチサイズ . . . . . . . . . . 10

5.2.4

活性化関数 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10i

5.2.5重みの初期値 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10

5.2.6

パラメータ更新 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11

5.2.7

Batch Normalization . . . . . . . . . . . . . . . . . . . . . . . . . . . 11

5.3

実験結果 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11

5.4

考察 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11

5.4.1

学習精度が低い原因 . . . . . . . . . . . . . . . . . . . . . . . . . . . 125.4.2出力された声道断面積関数の誤差が大きい原因 . . . . . . . . . . . . . 13

5.4.3出力された声道断面積関数に負数が含まれる原因 . . . . . . . . . . . 13

第 6 章 結論

14

参考文献

15

付 録 A 声道断面積関数プログラム

17

付 録 B 声道伝達特性のプログラム

21

付 録 C ニューラルネットワークのプログラム

46

謝辞

54

ii

1

第1章

緒言

音声の音韻性（声道伝達特性）は音源フィルタ理論 [1] より発話時の声道形状により一意的に特徴付けられる。一方、音声から声道形状を求めることは同じ音声波形であっても声道形状が異なる場合があり得るため、一意的に決まらず容易ではない。近年、機械学習により多量のデータを学習させることで、これまで得ることが困難であったデータを得ることが可能になった。そこで、この技術を用いることで音声から声道形状を求めることができれば、声
道形状の制御を行う仕組みが求められる。この仕組みを用いて人間と同じ発話機構を有す
る発話ロボットなどの開発が可能となる。これらのことから本研究では、機械学習を用い、音声情報としての声道伝達特性から声道形状を得ることを試みた結果について報告する。3次元の声道断面積関数を声道形状として用いると情報量が大きく、機械学習による計算コス
トが高くなり、学習結果を得るまでに多大な時間がかかる。概ね 3Hz の母音であれば、1 次元の近似で声道伝達特性を求められることから、最初の試みとして本報告では、情報量が小さく機械学習による計算コストの低い、1 次元の声道断面積関数を声道形状として扱う。また、学習に必要な多量の声道断面積関数を生成する方法についても述べる。

2

第2章

音声生成機構声道断面積関数および声道伝達特性、音源フィルタ理論について

2.1

音声の生成過程は、声門側の片端を閉鎖し、口唇側となるもう一方の端が開いた音響管でモデル化される。声門で発生した音（音源）は声道内で声道形状によって決定される共鳴特性（伝達特性）により変化し、口唇から空気中に放射される。以上より本章では、音声の生成過程を線形モデルで表現した音源フィルタ理論、および声道断面積関数と声道伝達特性について述べる。

2.1.1

音源フィルタ理論

人間は、喉頭における発声（音源）と声道における調音（フィルタ）とが独立であると仮定できる。音源が声道フィルタに入力されると共鳴特性が反映される。その結果出力されるものが音声であると考えることができる。このような線型システムによるモデル化を音源フィルタ理論 [1] と いう。もし、声道形状が一定ならば、声道フィルタは線形時不変となり、入力信号 x(t)、インパルス応答 s(t)、出力信号 h(t) としたとき以下の式で表せる。ここで、tは時間を表す。

h(t) = s(t) ∗ x(t)(2.1)

ここで*は、たたみ込み演算を表す。上の式を周波数領域において表現すると、音源のス
ペクトル X(ω)、声道フィルタのスペクトル H(ω)、声道の伝達関数 T (ω) と口唇からの放射
特性 R(ω) としたとき以下の式で表せる。ここで ω は周波数を表す。

H(ω) = [T (ω)R(ω)]X(ω)

2.1.2

(2.2)

声道断面積関数

声門から口唇に向かって伝搬する軸（声道中心線）と中心線に垂直な面と声道の輪郭との交線により声道の断面積形状を近似できる。声道断面積関数とは、この断面積を声道から口第 2 章 音声生成機構

3

唇に至る中心線上の距離を関数として表したものである [1]。また、声道断面積関数は声道の伝達特性を求める重要なパラメータである。

2.1.3

声道伝達特性

声道伝達特性とは、音源が声道フィルタに入力されると反映される共鳴特性である [1]。
その結果出力されるものが音声であると考えることができるため、音声を構成する要素とも
言える。

4

第3章

声道形状を考慮した多量の声道断面積
関数の自動作成について

3.1

まえがき

機械学習には、大量のデータが必要である。しかし、声道断面積関数を求める際に、MRI等を用いて多数の人から得るのは困難である。そこで、声道形状を考慮してランダムに声道
断面積関数を求めることができれば、容易に多量のデータを得ることができる。そこで本章
では、声道形状を考慮した制約条件よりランダムに求めた多量の声道断面積関数の自動作成について示す。

3.2

声道断面積関数のセクション数について

本手法では、成人男性の声道長 17cm [2] を 0.5cm ごとに区切った 34 セクションを求める。
セクション数に関して、神山の論文 [3] より、母音を区別するために必要な３つのホルマントと呼ばれる伝達特性内の特徴を求めるのに影響しないセクション数は 20 セクション以上あればよいとされている。

3.2.1

声道形状の決定要因

声道形状は、唇、口蓋、咽頭壁、舌、喉頭蓋、披裂軟骨、声帯などの発声器官の位置によって決まる [1]。このことから、本研究では、調音位置として声門、歯列、喉頭、咽頭、軽口
蓋、口腔、口唇の 7 つの位置を求め、声道の形状を求める。

3.3

声道断面積関数の導出方法

導出手順として、まず声道断面積関数の計算のため、文献 [6] の図 69 などの断面積の模式図を参考に 7 つの調音位置と声道断面積関数を求める。調音位置以外のセクションの声道断
面積関数はスプライン関数 [4, 5] を用いて、先に求めた 7 つの調音位置が滑らかにつながるよう求める。

第 3 章 声道形状を考慮した多量の声道断面積関数の自動作成について

5

調音位置の導出方法

3.3.1

声門および口唇の位置を固定し、喉頭、咽頭、軟口蓋、口腔、歯列の位置を乱数により、

1〜2 セクション程前後させ 7 セクションの調音位置を求める。また、歯列は出口から 1cmの位置で固定している。以下にそれぞれの位置の導出方法を示す。乱数は 0.0〜1.0 の間の値である。

glottis(声門) = 0
larynx（喉頭） =velum（軟口蓋） =pharynx（咽頭） =teeth（歯列） =oral（口腔） =lips（口唇） =11+seclng seclng ∗ 乱数
セクション数−1
2(velum − larynx − 1) ∗ 乱数
larynx +2
1
セクション数 − 1 −seclng(teeth − velum − 2) ∗ 乱数
velum +
2セクション数 − 1

(3.1)
(3.2)
(3.3)
(3.4)
(3.5)
(3.6)
(3.7)

上式より、喉頭 (larynx) を例に seclng=0.5 とした場合、larynx = 2 + 2*乱数 となる。
よって喉頭のセクションは 2〜3 という結果が求められる。その他 6 つの調音位置についても同様に求められる。

3.3.2

調音位置における声道断面積

ここでは、先に求めた 7 つのセクションの声道断面積を求める。声道断面積の出力範囲は
神山の文献 [3] より 0.45029cm2 〜12.67086cm2 とする。また、声門の最大値は同文献 [3] より 1.23364cm2 とする。

3.3.3

多量の声道断面積関数の求めるプログラム

本章 2.1〜2.4 を元に作成したプログラムを用いて多量の声道断面積関数を求める。付録

A にプログラムの詳細を記述する。図 3.1 に本プログラムより求めた声道断面積関数の例を示す。
以上より、声門より 14cm の位置での声道断面積はいずれも絞られていることが分かる。

第 3 章 声道形状を考慮した多量の声道断面積関数の自動作成について

図 3.1: 声道断面積関数の例

6

7

第4章

電気等価回路モデルを用いた声道断面
積関数からの声道伝達特性の導出方法

4.1

まえがき

本章では、文献 [3] の電気的等価回路モデルを参考にして声道断面積関数から声道伝達特性を求める。

4.2

声道伝達特性の導出方法

第 3 章より求めた声道断面積関数から電気的等価回路モデル [3] を用いて声道伝達特性を導出する方法を以下に示す。

4.2.1

電気的等価回路モデル

文献 [1] の電気的等価回路より、声道伝達特性を求めるには、声道断面積関数は階段関数により近似し、等しい長さ l で分割（本研究では 34 分割）する。その時、1 区間を伝送線路
による等価回路として扱った場合、全区間を縦続接続音響管として扱うことができる。このような処理により、声道断面積関数から声道伝達特性を求めることができる。
図 4.1 に上記の手順で求めた声道伝達特性の例を示す。

第 4 章 電気等価回路モデルを用いた声道断面積関数からの声道伝達特性の導出方法

8

図 4.1: 声道断面積関数から求めた声道伝達特性の例（左図：声道断面積関数、右図：声道
伝達特性

9

第5章

機械学習による声道断面積の出力方法
および結果と考察

5.1

まえがき

本研究では、4 章を元に作成した声道伝達特性データを入力とし、3 層ニュートラルネッ
トワークを用いて声道断面積関数の出力を試みる。

5.2

ニューラルネットワークの設定

本研究で用いるニューラルネットワークの設定を本項に示す。

5.2.1

ニューラルネットワークの概略図

ニューラルネットワークの構造を図 5.1 に示す。

図 5.1: 学習に用いるニューラルネットワークの構造
上図は 3 層ニューラルネットワークを示し、1 層目を Aﬃne、Batch Normalization とし、2 層目を Sigmoid、3 層目を Aﬃne、Cross entropy とした図である。Aﬃne は入力値をそのまま出力する処理であり、Cross entropy は損失関数と呼ばれるものであり、ニューラルネッ
トワークの精度を測定するものである。その他の処理についての詳細は次項以降に示す。

第 5 章 機械学習による声道断面積の出力方法および結果と考察

5.2.2

10

教師あり学習ニューラルネットによる出力を試みるにあたり、ニューラルネットワークの学習手法の一種である教師あり学習を用いた手法を提案する。ニューラルネットワークとは、機械学習の一種であり、適切な重みパラメータをデータから自動で学習するものである。教師あり学習とは、入力データと教師データ（入力データに対しラベル付けしたデータ）を対としたデー
タセットを学習させる手法である。本手法では、入力データに声道伝達特性、教師データに
声道断面積関数を用いて学習を試みる。

5.2.3

ノード数、学習係数、学習回数およびバッチサイズ

ノード数は入力層、中間層、出力層の順に 1024、50、34 とする。学習係数は 0.001 とし、学習回数は 10000 回とする。バッチサイズを 100 とする。学習係数とは、1 回の学習でどれだけパラメータを更新するかを定めるものである。バッチとは、あるデータセットをひとまとめにして処理する単位である。バッチサイズを入力データ数より小さくすることで学習速度を高めることができる [7]。

5.2.4

活性化関数

活性化関数として式 (5.1) で表される Sigmoid 関数を用い、出力層では恒等関数を用いる。
活性化関数とは、入力信号の総和を出力信号に変換する関数のことを表す。ニューラルネッ
トワークでは活性化関数に非線形関数を用いる必要がある [7] ため、非線形関数の一種であり古くからニューラルネットワークの学習に用いられてきた Sigmoid 関数を使用する。恒等関数とは、入力信号をそのまま出力信号とする関数である。

h(x) =

5.2.5

1
1 + exp(−x)(5.1)

重みの初期値

重みの初期値として Xavier らの論文 [8] で推奨される Xvaier の初期値を用いる。重みの初期値とは、ニューラルネットワークの学習の成否が分かれるものとされており [7]、活性化関数として Sigomid 関数を用いる場合は、Xvaier の初期値が推奨されている [7]。Xvaier の初期値は、前層のノード数を n とした場合、√1n の標準誤差をもつ分布を使うことで、各層のアクティベーション分布（活性化関数の出力）の広がりを同程度にすることができる。この
結果、Sigmoid 関数の表現力を制限することなく学習を効率的に進めることが期待できる。

第 5 章 機械学習による声道断面積の出力方法および結果と考察

5.2.6

11

パラメータ更新

重みパラメータの更新には Adam [9] という手法を用いる。この手法は、他のパラメータ更新手法である Momenturm の最適化の更新経路 [7] と AdamGrad [10] の学習係数の衰退という利点を併せ持つ手法であり、ハイパーパラメータの偏りを補正するという特徴を持つ。

5.2.7

Batch Normalization

学習を効率的に進めることを目的に、以下の利点が期待できる Batch Normaliztion [11]を用いる。

• 学習を速く進めることができる（学習係数を大きくできる）• 初期値にそれほど依存しない• 過学習を抑制する
過学習とは、あるデータセットにだけ過度に対応した汎化能力が失われている状態のことを言う。Batch Normalization は重みの初期値を強制的にアクティベーション分布を同程度にするよう調節する手法である。具体的には、学習を行う際のミニバッチを単位としてミニ
バッチごとに、データの分布の平均が０、分散が１になるように正規化する。

5.3

実験結果

下記に、学習後のニュートラルネットワークの認識精度と自動抽出より得た声道断面積を
図 5.2、図 5.3 に示す。
図 5.2 はニュートラルネットワークの学習精度の推移を示しており、学習後の精度は約 35％前後となっている。この学習により得た重みパラメータを用いて出力した声道断面積関数が図 5.3 である。出力された声道断面積関数とテストデータとして使用した声道断面積関数の平均誤差は 140cm2 であり、最小誤差は 32cm2 、最大誤差は 266cm2 であった。また、出力された声道断面積関数には負数が含まれていた。

5.4

考察

実験より以下、3 点の問題があることがわかった。ニュートラルネットワークの学習精度が低い点と出力した声道断面積関数の誤差が大きい点および負数が含まれる点。本章では、これらの問題の考察を示す。

第 5 章 機械学習による声道断面積の出力方法および結果と考察

12

図 5.2: ニュートラルネットワークの精度の推移（実線：テストデータの認識精度、破線：
訓練データの認識精度）

図 5.3: テストデータの声道断面積データと出力された声道断面積データ（実線：テストデー
タとして用いた声道断面積関数、破線：学習より得た声道断面積関数）

5.4.1

学習精度が低い原因

学習精度が低い原因の 1 つとして以下の要素が考えられる• 訓練データに過度に対応している（オーバーフィッティング）• 訓練データに過度に対応していない（アンダーフィッティング）
上記が起こる原因として以下のようなことがある。
オーバーフィッティング

アンダーフィッティング

訓練データ数が少ない

訓練データ数が多い

特徴量の種類が多い

特徴量の種類が少ない

第 5 章 機械学習による声道断面積の出力方法および結果と考察

13以上より、オーバーフィッティングによる訓練データ数の過小が原因と仮定し、追加実験を行った。追加実験は訓練データおよびテストデータ数のみ変更する。実験内容は訓練デー
タ数を 35000、テストデータを 15000 に変更しそれ以外は同条件にて学習を行った。図 5.4に追加実験の学習結果を示す。

図 5.4: 追加実験結果：ニュートラルネットワークの精度の推移（実線：テストデータの認
識精度、破線：訓練データの認識精度）以上の結果から、データ数を増やした結果、学習精度が低下することが分かる。これらから、訓練データ数およびテストデータ数は当初の設定で適当であり、特徴量の種類が学習精度低下の要因ではないかと考える。

5.4.2

出力された声道断面積関数の誤差が大きい原因

出力結果にニューラルネットワークの学習後に得た重みパラメータを用いていることから出力された声道断面積関数とテストデータの声道断面積関数との誤差が大きい原因としてニューラルネットワークの学習精度が低いことが関与していると考えられる。

5.4.3

出力された声道断面積関数に負数が含まれる原因

出力された声道断面積関数に負数が含まれる原因として、出力結果の誤差が大きい原因と同様、ニューラルネットワークの学習精度に関係があると考える。

14

第6章

結論

本研究では、3 層ニューラルネットワークを用いて伝達特性データから声道断面積関数データの出力を行った。ニュートラルネットワークの学習精度は出力に十分なほど上がらず、適切な声道断面積関数データの出力には至らなかった。原因としては、ニューラルネットワークの認識精度が低いため、出力した結果の精度も低いものなったと考える。本研究ではニュー
トラルネットワークを用いたが、機械学習では新たな技術の開発が進んでいる。今後はそれらの技術が出力に利用できないか実験と考察を行いたい。

15

参考文献

[1]

正木信夫, 元木邦俊, 松﨑博季, 北村達也, 音声生成の計算モデルと可視化, コロナ社,

2010[2] Goldstein, U.G. An articulatory model for the vocal tracts of growing children. Ph.D.
dissertation, Massachusetts Institute of Technology, Cambridge, MA, pp.184–197,
June. 1980.
[3] 神山直久, 音声生成過程における音響特徴量抽出と実体的声道モデルに関する研究, 北海道大学博士論文, March. 1993.

[4] PolynomialSplineFunction(Commons Math 2.2 API), https://commons.apache.
org/proper/commons-math/javadocs/api-2.2/org/apache/commons/math/analysis/polynomials/PolynomialSplineFunction.html[5] SplineInterpolator(Commons

Math

3.0

API),

https://commons.apache.

org/proper/commons-math/javadocs/api-3.0/org/apache/commons/math3/
analysis/interpolation/SplineInterpolator.html[6] 千葉勉, 梶山正登, 母音-その性質と構造-, 岩波書店, 2003
[7] 斎藤泰毅, ゼロから作る Deep Learning − Python で学ぶディープラーニングの理論と実装, O’REILLY, 2016

[8] Xvaier Glorot, Yoshua Bengio, Understanding the diﬃculty of training deep feedforward neural nerworks, AISTATS2010, 2010[9] Diederik Kingma, Jimmy Ba, Adam:A Method for Stochastic Optimization,
arXiv:1414.6980, December. 2014[10] Jhon Duchi, Elad Hazan, Yoram Singer, Adaptive Subgradient Methods for Online Learning and Stochastic Optimization, Journal of Machine Learning Reasearch,

文献

16

pp.2121–2159. July. 2011[11] Sergey Ioﬀe, Christian Szegedy, Batch Normalization: Accelerating Deep Network
Training by Reducing Internal Covariate Shift, arXiv:1502.03167, February. 2015

17

付 録A

声道断面積関数プログラム

声道断面積関数を求めるプログラムを以下に示す。
1

package jp.ac.hus.mhlab;222

import java.util.Random;
import static java.lang.Math.*;import org.apache.commons.math3.analysis.interpolation.
SplineInterpolator;
import org.apache.commons.math3.exception.DimensionMismatchException;
import org.apache.commons.math3.exception.NonMonotonicSequenceException;
import org.apache.commons.math3.exception.NumberIsTooSmallException;
import org.apache.commons.math3.exception.OutOfRangeException;;
import org.apache.commons.math3.analysis.polynomials.
PolynomialSplineFunction;
import java.io.File;
import java.io.FileWriter;
import java.io.PrintWriter;
import java.io.BufferedWriter;
import java.io.IOException;/* 各部位の断面積 [cm2]（千葉＆梶山）
gl la P v o t lip/a/ 1.18 1.18 7.45 1.15 3.70 11.30 5.65 8.40
/i/ 1.18 1.55 7.94 6.81 7.10 0.80 1.76 4.55
/u/ 1.18 1.55 6.59 6.59 2.00 4.75 0.75 1.65
/e/ 1.18 1.18 6.59 6.61 3.10 1.90 3.20 5.00
/o/ 1.18 1.18 7.30 1.35 2.95 10.65 2.00 1.40

23

*/

3
4
5

6
7
8
9
10

11
12
13
14
15
16
17
18
19
20
21

24
25
26
27

28

29

30
31
32
33

public class AreaFunctionMaker {private final static double minArea = 0.45029; // 最小面積 [cm2]。
神山さんのD 論から
private final static double maxArea = 12.67086; // 最大面積 [cm2]。
神
山さんのD 論から
private final static double maxInitArea = 1.23364; // 入力面の最大面
積 [cm2]。神山さんのD 論から
private int nsec; //セクション数
private double length; // 軸長
private double[] area; // 面積
private double[] dist; // 入力からの距離

付 録A

34
35
36
37
38
39
40
41
42

声道断面積関数プログラム

18

private double seclng;
private int iglottis;
private int ilarynx;
private int ipharynx;
private int ivelum;
private int ioral;
private int iteeth;private int ilips;
private int dataset = 10000;

4345

public AreaFunctionMaker() {this(34, 17.0); // 1セクション 0.5cm とすると 17cm で 34 セクション

46

}

44

47
48
49
50
51
52

public AreaFunctionMaker(int nsec, double length) {this.nsec = nsec;
this.length = length;area = new double[nsec];
dist = new double[nsec];53

seclng = length / nsec;

54
55

58

for(int i=0; i<nsec; i++) {dist[i] = seclng * i;//System.out.println("deit[" + i + "]=" + dist[i]);

59

}

56
57

60

setSection();

61
62

}

63

76

private void setSection(){iglottis = 0;
Random rand = new Random();
ilarynx = (int) (1.0 / seclng) + (int) (1.0 / seclng *random());//ipharynx = (int) (length / 3.33333 / seclng) - 1;
ivelum = nsec / 2 - 1;//ipharynx = ilarynx + (int)((ivelum - ilarynx - 1)*random()) + 1;
ipharynx = ilarynx + (int)((ivelum - ilarynx - 1)*random())/2;//ioral = (int) (length / 1.41421356 / seclng) - 1;
iteeth = nsec - 1 - (int) (1.0 / seclng); // 歯は出口から 1
cm の位置で固定//ioral = ivelum + (int)((iteeth - ivelum - 2)*random()) + 1;
ioral = ivelum + (int)((iteeth - ivelum - 2)*random())/2;
ilips = nsec - 1;

77

}

64
65
66
67
68
69
70
71
72
73

74
75

78
79
80

public void calculate(){boolean flag = true;

付 録A

81
82
83
84
85
86
87
88
89
90
91

声道断面積関数プログラム

19

int count=0;
while(count<dataset){
double[] x = new double[7];double[] y = new double[7];x[0]=dist[iglottis];x[1]=dist[ilarynx];x[2]=dist[ipharynx];x[3]=dist[ivelum];x[4]=dist[ioral];x[5]=dist[iteeth];x[6]=dist[ilips];

92

102

y[0]=1.0 + 0.2*random();
y[1]=1.2 + (1.55 - 1.2) * random();y[2]=1.15 + (7.94 - 1.15) * random();y[3]=2.0 + (7.1 - 2.0) * random();
y[4]=0.8 + (11.3 - 0.8) * random();y[5]=0.75 + (5.65 - 0.75) * random();y[6]=1.4 + (8.4 - 1.4) * random();//for(int i=0; i< x.length; i++){// System.out.print("x[" + i + "]=" + x[i]);// System.out.println(", y[" + i + "]=" + y[i]);

103

//}

93
94
95
96
97
98
99
100
101

104

112

SplineInterpolator si = new SplineInterpolator();try {PolynomialSplineFunction psf = si.interpolate(x, y);for (int i=0; i<area.length; i++){area[i] = psf.value(dist[i]);if(area[i]<0.0) {flag = false;
break;

113

}

105
106
107
108
109
110
111

114
115
116
117
118
119
120
121
122
123

//System.out.println(dist[i] + " " + area[i]);//System.out.println(i + " " + area[i]);}} catch (OutOfRangeException ore) {//ore.printStackTrace();flag = false;} catch (Exception ex) {flag = false;//ex.printStackTrace();
}

124
125
126
127
128

if(flag) {String dataname = "" + count;
try{String dir = System.getProperty("dir","C:\\Users\\Masataka

付 録A

声道断面積関数プログラム\\Documents\\NetBeansProjects\\Peaks\\area");
File file = new File(dir,dataname);
FileWriter fw = new FileWriter(file);
BufferedWriter bw = new BufferedWriter(fw);
PrintWriter pw = new PrintWriter(bw);

129
130
131
132
133

141

pw.format("%f",17.0);
pw.format("\n");pw.format("%d", area.length);pw.format("\n");for (int i=0; i<area.length; i++){
System.out.println(area[i]);pw.format(" %5f ",area[i]);pw.format("\n");

142

}

134
135
136
137
138
139
140

bw.close();

143

}

144

catch(IOException ex){
System.out.println(ex);

145
146

}

147

System.out.println(count);
count++;

148
149

}

150
151

flag=true;
setSection();

152
153

}

154

}

155
156

164

public void show(){
System.out.println("iglottis: " + iglottis);
System.out.println("ilarynx: " + ilarynx);
System.out.println("ipharynx: " + ipharynx);
System.out.println("ivelum: " + ivelum);
System.out.println("ioral: " + ioral);
System.out.println("iteeth: " + iteeth);
System.out.println("ilips: " + ilips);

165

}

157
158
159
160
161
162
163

166

public static void main(String args[]) {AreaFunctionMaker af = new AreaFunctionMaker();//af.show();
af.calculate();

167
168
169
170

}

171
172

}

20

21

付 録B

声道伝達特性のプログラム

声道伝達特性を求めるプログラムを以下に示す。
1/******************************************************

2
3

peaks : 声道伝達関数計算プログラム

4
5

CC peaks.C calf.C radz.C news.C -lcomplex -o peaks

6
7

pea.com でコンパイルできる。

89

ヘッダーファイルの例は ex.hed。

10
11
12

(C)opyright NAOHISA KAMIYAMA DEC. 1991

13
14******************************************************/

15

package jp.ac.hus.mhlab;
import java.io.IOException;
import java.io.File;
import java.io.FileReader;import java.io.PrintWriter;
import java.io.BufferedReader;import java.io.BufferedWriter;
import java.io.FileWriter;
import static java.lang.Math.*;import org.apache.commons.math3.transform.*;
import org.apache.commons.math3.complex.*;import org.apache.commons.math3.exception.*;

16
17
18
19
20
21
22
23
24
25
26
27
28

import java.util.Scanner;

29
30
31
32
33
34
35
36
37
38

public class Peaks {private WallImpedance ZW, ZWN;
private double freq, frest, finalf, dit, Na1, Ndit;
private int cal_tag, Ef_tag, nasal_tag;
private int nareav, icntlr, Nport, Nsec;
private String areav, arean, filekp, fileV, fileN;
private double amp_1, amp_2, amp_3, fin_1, fin_2, fin_3;
private double vleng, vlg, Vvleng, Nvlg;
private double[] radin, radout, sfact, Vradin, Vradout, Vsfact;

付 録B

39
40
41
42
43
44
45

声道伝達特性のプログラム

private double[] fdat = new double[4000]; //*double → int
private double[] abdat = new double[4000]; //*double → int
private double[] Nabdat = new double[4000]; //*double → int
private double[] Sabdat = new double[4000]; //*double → int//private Complex gam, fimp, Vgam, Vfimp;
private int ntimes;
private Complex hvtra, Nhvtra;

46
47

private double temper = 25.0;

48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86

public Peaks(File file, int a) throws IOException{//try {
FileReader fr = new FileReader(file);
BufferedReader br = new BufferedReader(fr);String data;
while((data = br.readLine()) != null) {String[] tmp = data.split("[ ]+");
for(int i=0; i<tmp.length; i++){//System.out.println(tmp[i]);switch(tmp[i]) {case "FinalFreq":
finalf = Double.parseDouble(tmp[++i]);
break;
case "CalSpeed":
cal_tag = Integer.parseInt(tmp[++i]);
break;
case "AreaName":areav = tmp[++i];//bufferedWriter.write("1");//bufferedWriter.close();break;
case "Ef_tag":
Ef_tag = Integer.parseInt(tmp[++i]);
break;
case "LossMode":
icntlr = Integer.parseInt(tmp[++i]);
break;
case "ZwallParam":
ZW = new WallImpedance(Double.parseDouble(tmp[++i]),
Double.parseDouble(tmp[++i]),
Double.parseDouble(tmp[++i]),
Double.parseDouble(tmp[++i]),
Double.parseDouble(tmp[++i]));break;
case "Thickness":
dit = Double.parseDouble(tmp[++i]);
break;
case "NasalFile":

22

付 録B

声道伝達特性のプログラム

107

arean = tmp[++i];
break;
case "ConnecSec":
Nport = Integer.parseInt(tmp[++i]);
break;
case "ConnecArea":
Na1 = Double.parseDouble(tmp[++i]);
break;
case "ZwnPara":
ZWN = new WallImpedance(Double.parseDouble(tmp[++i]),
Double.parseDouble(tmp[++i]),
Double.parseDouble(tmp[++i]),
Double.parseDouble(tmp[++i]),
Double.parseDouble(tmp[++i]));break;
case "ThickN":
Ndit = Double.parseDouble(tmp[++i]);
break;
case "Ofile":
filekp = tmp[++i];
break;

108

}

87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106

}

109
110

}

111

fr.close();//} catch (IOException e) {// e.printStackTrace();

112
113

//}
/*String[] tmp = new String[3];finalf = 4096; //FinalFreq
cal_tag = 4; //CalSpeed
areav = "0"; //AreaName
tmp[0] = areav;
Ef_tag = 0; //Ef_tag
icntlr = 0; //LossMode
ZW = new WallImpedance(3700.0, 1800.0, 2900.0, 2.2, 410.0); //

114
115
116
117
118
119
120
121
122

ZwallParamdit = 0.0; //Thickness
arean = "no"; //NasalFile
tmp[1] = arean;
Nport = 18; //ConncSec
Na1 = 0.4; //ConnecArea
ZWN = new WallImpedance(3700.0, 1800.0, 2900.0, 2.2, 410.0); //ZwnParaNdit = 2.0; //ThickN
filekp = "ddd"; //Ofiletmp[2] = filekp;

123
124
125
126
127
128

129
130
131
132

*/

23

付 録B

133
134
135
136
137
138
139
140
141
142

声道伝達特性のプログラム

if(cal_tag==4){ // added by Matsuzaki
frest = 4.0; // added by Matsuzaki
finalf+=1.0;} else { // added by Matsuzaki
frest= 20.0;
if(finalf<=0.0) finalf=8000.0;
if(cal_tag<=0) cal_tag=2;
finalf+=40.0;
CALTAG(cal_tag);} // added by Matsuzaki

143
144
145

if((dit<0.0)||(icntlr !=2)) dit=0.0;
ZW.setGw(ZW.getGw()*1000.0);

146
147
148
149
150
151
152

153
154
155

156
157
158
159
160
161
162
163
164
165
166
167
168

////////////////////////////////////////
if(icntlr==3){ // 裏技
System.out.println(">Wall Impedance Parameter");
System.out.println(" [R],[L],[K], [800, 2.1, 84500] ===> ");Scanner scan = new Scanner(System.in);//cin>> ZW.R1; if(ZW.R1<=0.0){ ZW.R1=800.;ZW.Lo=2.1 ;ZW.Gw=84500.;}// else { cin >> ZW.Lo; cin >> ZW.Gw;}
ZW.setR1(scan.nextDouble());
if(ZW.getR1()<=0.0){ ZW.setR1(800.0); ZW.setLo(2.1); ZW.setGw(84500.0);}
else {
ZW.setLo(scan.nextDouble());
ZW.setGw(scan.nextDouble());}}////////////////////////////////////////
readAreaFunctionData();if(!arean.equals("no")) {readNasalAreaFunctionData();
nasal_tag = 1;}
else
nasal_tag = 0;

169
170
171
172

/***** 出力ファイル *****/
fileV = filekp + ".v";
fileN = filekp + ".n";

173
174
175
176
177
178

FileWriter fw = new FileWriter(file);
BufferedWriter bw = new BufferedWriter(fw);
PrintWriter pw =new PrintWriter(bw);
pw.format("FinalFreq");
pw.format(" 4096");

24

付 録B

声道伝達特性のプログラム

25

pw.format("\n");pw.format("CalSpeed");pw.format(" %s", cal_tag);pw.format("\n");pw.format("AreaName");pw.format(" %s", ++a);
pw.format("\n");pw.format("Ef_tag");
pw.format(" %s", Ef_tag);pw.format("\n");pw.format("LossMode");pw.format(" %s", icntlr);
pw.format("\n");pw.format("ZwallParam");pw.format(" 3700.0 1800.0 2900.0 2.2 410.0");
pw.format("\n");pw.format("Thickness");
pw.format(" %s", dit);pw.format("\n");pw.format("NasalFile");pw.format(" no");pw.format("\n");pw.format("ConnecSec");
pw.format(" %s", Nport);
pw.format("\n");pw.format("ConnecArea");
pw.format(" %s", Na1);
pw.format("\n");pw.format("ZwnPara");pw.format(" 3700.0 1800.0 2900.0 2.2 410.0");
pw.format("\n");pw.format("ThickN");pw.format(" %s", Ndit);pw.format("\n");pw.format("Ofile");pw.format(" %s", a);pw.format("\n");

179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216

bw.close();

217
218

}

219
220
221

public void showHelp() {System.out.println(" o..... P E A K S (C) N.KAMIYAMA 1991 .....o%n

");
222
223
224
225

System.out.println("******* Hedder Format *******");
System.out.println("FianlFreq [freq]");
System.out.println("Calculate Speed [1~3]");

付 録B

声道伝達特性のプログラム

System.out.println("FileName [Area-curcum format file]");System.out.println("Ef Mode [1/0]");
System.out.println("LossMode [0:Lossless 1:Air 2:Zwall]");System.out.println("Zwall Parameter [R1][R2][R3][Lo][Gw(xe3)]");
System.out.println("WallThickness [dit]");System.out.println("NasalName [name]");System.out.println("ConnectSec [int]");System.out.println("ConnectArea [Area]");System.out.println("Zwn Param (if Loss mode = 2) [R1][R2][R3][Lo][Gw(xe3)]");
System.out.println("WallThickness [dit]");System.out.println("Outputfile [name]");System.out.println("*****************************");

226
227
228
229
230
231
232
233
234

235
236
237
238

26

}

239
240
241
242
243
244
245
246

247

public void showStatus(){
System.out.println(" * ************** STATUS *************** *");
System.out.printf("\tInput File : %s(%dsec.)%n", areav, nareav);
System.out.printf("\tOutput Fie : %s%n", filekp);
System.out.printf("\tLoss Mode : %d%n", icntlr);
if(icntlr == 2) {System.out.printf("\t Zw R1:%5.1f R2:%5.1f R3:%5.1f%n", ZW.getR1(), ZW.getR2(), ZW.getR3());
System.out.printf("\t Lo:%5.1f Gw:%f %n",ZW.getLo(), ZW.getGw());
System.out.printf("\t Thickness : %f%n",dit);

248
249

}

250

System.out.printf("\tElongation : %d%n",Ef_tag);if(!arean.equals("no")){System.out.printf("\tNasal File : %s(%dsec.)%n", arean, Nsec);
System.out.printf("\t Connect Sec : %d%n", Nport);
System.out.printf("\t Connect Area: %5.2f%n", Na1);
if(icntlr == 2){
System.out.printf("\t Zwn R1:%5.1f R2:%5.1f R3:%5.1f%n", ZWN.
getR1(), ZWN.getR2(), ZWN.getR3());
System.out.printf("\t Lo:%5.1f Gw:%f %n", ZWN.getLo(), ZWN.
getGw());
System.out.printf("\t Thickness : %f%n", Ndit);

251
252
253
254
255
256

257

258

}

259
260

}

261

else System.out.println("\tNasal File : <None>");System.out.println(" * **************************************

262

*");
263

}

264
265
266
267

private void CALTAG(int n) {switch (n) {case 1:

付 録B

声道伝達特性のプログラム

amp_1= 8.0; fin_1=78.0;
amp_2=14.0; fin_2=48.0;
amp_3=18.0; fin_3=16.0;
break;
case 2:
amp_1= 5.0; fin_1=30.0;
amp_2=10.0; fin_2=20.0;
amp_3=16.0; fin_3=10.0;
break;
case 3:
amp_1= -10.0; fin_1=16.0;
amp_2= 3.0; fin_2=8.0;
amp_3=10.0; fin_3=4.0;
break;

268
269
270
271
272
273
274
275
276
277
278
279
280
281

}

282
283

}

284
285

286
287
288
289
290

291
292
293
294
295

//* ---------------- READ AREA FUNCTION DATA------------------------*/
public void readAreaFunctionData(){readAreaFunctionData(areav);}public void readAreaFunctionData(String filename){
String dir = System.getProperty("dir","C:\\Users\\Masataka\\
Documents\\NetBeansProjects\\Peaks\\area");
File file = new File(dir,filename);if(!file.exists()) {System.out.println("ファイルが存在しません。");
System.exit(1);
}

296
297
298
299
300
301
302

try {
FileReader fileReader = new FileReader(file);
BufferedReader bufferedReader = new BufferedReader(fileReader);String data;//while((data = bufferedReader.readLine()) != null) {// System.out.println(data);

303

//}

304

data = bufferedReader.readLine();vleng = Double.parseDouble(data);
data = bufferedReader.readLine();nareav = Integer.parseInt(data);

305
306
307
308
309
310
311

radin = new double[nareav];radout = new double[nareav];sfact = new double[nareav];

312
313

double[] area = new double[nareav];

27

付 録B

声道伝達特性のプログラム

28

double[] circ = new double[nareav];for(int i=0; i<nareav; i++){data = bufferedReader.readLine();String[] tmp = data.split("[ ]+");//for(int j=0; j<tmp.length; j++){// System.out.println("i=" + i + "tmp[" + j + "]=" + tmp[j]);

314
315
316
317
318
319
320

//}

321

area[i] = Double.parseDouble(tmp[1]);//circ[i] = Double.parseDouble(tmp[2]);//System.out.println(i + " " + area[i] + " " + circ[i]);

322
323
324

}

325

fileReader.close();

326

329

for(int i=0; i<nareav; i++) {radin[i] = sqrt(area[i]/PI);
radout[i] = radin[i] + dit;

330

}

327
328

331

if(Ef_tag == 1)
for(int i=0; i<nareav; i++) sfact[i]= circ[i]/(2.0*sqrt(PI*area[i]));else
for(int i=0; i<nareav; i++) sfact[i]=1.0;

332
333

334
335
336

/*===========================================*/
revar(radin); revar(radout);/*===========================================*/
vlg=vleng/nareav;
System.out.println("vlg=" + vlg);} catch (IOException e){e.printStackTrace();
System.exit(1);
}

337
338
339
340
341
342
343
344
345
346

}

347
348
349
350
351
352
353
354
355
356
357

/*---------------------- Nasal Part ---------------------------*/
public void readNasalAreaFunctionData(){
readNasalAreaFunctionData(arean);}public void readNasalAreaFunctionData(String filename){
File file = new File(filename);if(!file.exists()) {System.out.println("ファイルが存在しません。");
System.exit(1);
}

358
359
360

try {
FileReader fileReader = new FileReader(file);

付 録B

声道伝達特性のプログラム

364

BufferedReader bufferedReader = new BufferedReader(fileReader);String data;//while((data = bufferedReader.readLine()) != null) {// System.out.println(data);

365

//}

366

data = bufferedReader.readLine();Vvleng = Double.parseDouble(data);
data = bufferedReader.readLine();Nsec = Integer.parseInt(data);

361
362
363

367
368
369

29

370
371
372
373

Vradin = new double[Nsec];
Vradout = new double[Nsec];
Vsfact = new double[Nsec];

374
375
376
377
378
379
380
381

double[] area = new double[Nsec];
double[] circ = new double[Nsec];
for(int i=0; i<Nsec; i++){data = bufferedReader.readLine();String[] tmp = data.split("[ ]+");//for(int j=0; j<tmp.length; j++){// System.out.println("i=" + i + "tmp[" + j + "]=" + tmp[j]);

382

//}

383

area[i] = Double.parseDouble(tmp[1]);//circ[i] = Double.parseDouble(tmp[2]);//System.out.println(i + " " + area[i] + " " + circ[i]);

384
385
386

}

387

fileReader.close();

388
389
390

area[0]=Na1; /* まだArea[cm2]. */Nport = nareav - Nport;

391

394

for(int i=0; i<Nsec; i++){
Vradin[i]=sqrt(area[i]/PI);
Vradout[i]= Vradin[i]+ Ndit;

395

}

396

if(Ef_tag == 1)
for(int i=0; i<Nsec; i++) Vsfact[i]= circ[i]/(2.0*sqrt(PI*area[i]));else
for(int i=0; i<Nsec; i++) Vsfact[i]=1.0;

392
393

397

398
399
400
401
402
403

/*-------------------------------------*/
revar(Vradin); revar(Vradout);/*-------------------------------------*/

404
405

System.out.printf("Vrad %f%n", Vradin[Nsec-1]);

406
407

Nvlg=Vvleng/Nsec;

付 録B

声道伝達特性のプログラム} catch (IOException e){e.printStackTrace();
System.exit(1);
}

408
409
410
411
412

}

413
414
415
416
417
418

/************* revar.c ****************/
private void revar(double[] data){
int N = data.length;
double[] tmp = new double[N];419

for(int i=0 ; i<N ; i++) tmp[i]=data[i];for(int i=0 ; i<N ; i++) data[i]=tmp[N-i-1];

420
421
422

}

423
424
425

public void calculate(){System.out.println("\t!!!!! NOW CALCULATING !!!!!");

426
427
428
429

freq=frest;
ntimes=0;
hvtra = null; Nhvtra = null;

430
431
432
433
434
435
436
437
438
439
440
441
442
443
444

while(freq<finalf){// ------------------------------------------if(nasal_tag ==1) {VocalAndNose(hvtra, Nhvtra); /* NASAL VERSION */Nabdat[ntimes]= (20.0*log10(Nhvtra.abs()));
abdat[ntimes] = (20.0*log10(hvtra.abs()));
Sabdat[ntimes]= (20.0*log10(Nhvtra.abs() + hvtra.abs()));}
else {//Vocal(hvtra); /* VOCAL VERSION */Vocal(); /* VOCAL VERSION */Sabdat[ntimes]= (20.0*log10(hvtra.abs()));}// --------------------------------------------

445
446

fdat[ntimes]= freq; //*double → int

447
448
449
450
451
452
453
454
455

// if(((ntimes)%40)==0) printf("%4.0f %n",freq);if(cal_tag==4){ // added by Matsuzaki
freq += 4.0; // added by Matsuzaki} else { // added by Matsuzakiif(Sabdat[ntimes]< amp_1) freq+=fin_1;
else{ if(Sabdat[ntimes]< amp_2) freq+=fin_2;
else{ if(Sabdat[ntimes]< amp_3) freq+=fin_3;
else freq+=2.0;}}

30

付 録B

声道伝達特性のプログラム

456

}

457

ntimes+=1;

}

458
459

31

}

460
461

public void saveFile(){/* ファイルへの書き込み */try{if(nasal_tag ==1) {String dir = System.getProperty("dir","C:\\Users\\Masataka\\
Documents\\NetBeansProjects\\Peaks\\peaks");
File file_V = new File(dir,fileV);
PrintWriter pwV = new PrintWriter(file_V);
pwV.format("%d%n", ntimes);for(int i=0; i<ntimes ; i++) pwV.format("%f%n", abdat[i]);pwV.format("%n");
pwV.close();

462
463
464
465

466
467
468
469
470
471
472

File file_N = new File(dir,fileN);
PrintWriter pwN = new PrintWriter(file_N);
pwN.format("%d%n", ntimes);for(int i=0; i<ntimes ; i++) pwN.format("%f%n", Nabdat[i]);pwN.format("%n");
pwN.close();
System.out.println("File: " + fileV);
System.out.println("File: " + fileN);

473
474
475
476
477
478
479
480

}

481
482

String dir = System.getProperty("dir","C:\\Users\\Masataka\\
Documents\\NetBeansProjects\\Peaks\\peaks");
File file = new File(dir,fileV);
PrintWriter pwkp = new PrintWriter(file);//pwkp.format("%d%n", ntimes);for(int i=0; i<ntimes; i++){//fdat[i]=fdat[i] * 0.001;
fdat[i]=fdat[i];pwkp.format("%f%n", Sabdat[i]);

483

484
485
486
487
488
489
490

}
//pwkp.format("%n");
pwkp.close();
System.out.println("File: "+ filekp + " is made !");} catch (IOException ex){ex.printStackTrace();
}

491
492
493
494
495
496
497
498

}

499
500

///*-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-*/

付 録B

501
502

503
504
505
506
507
508
509
510
511
512

513
514
515
516
517
518
519
520
521
522
523
524

525
526
527
528
529
530
531
532
533

声道伝達特性のプログラム///* KEEP IN THE LOG FILE "Log.Peaks"///*-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-*/////fp=fopen("Log.Peaks","a");//// fprintf(fp," * ************** STATUS *************** *%n");// fprintf(fp,"\tInput File : %s%n",areav);// fprintf(fp,"\tOutput Fie : %s%n",filekp);// fprintf(fp,"\tLoss Mode : %d%n",icntlr);// if(icntlr == 2)// {// fprintf(fp,"\t Zw R1:%5.1f R2:%5.1f R3:%5.1f%n",ZW.R1, ZW.R2, ZW.
R3);// fprintf(fp,"\t Lo:%5.1f Gw:%f %n",ZW.Lo, ZW.Gw);// fprintf(fp,"\t Thickness : %f%n",dit);// }// fprintf(fp,"\tElongation : %d%n",Ef_tag);// if(strcmp(arean,"no") != 0)// {// fprintf(fp,"\tNasal File : %s%n",arean);// fprintf(fp,"\t Connect Sec : %d%n",Nport);// fprintf(fp,"\t Connect Area: %5.2f%n",Na1);// if(icntlr == 2)// {// fprintf(fp,"\t Zw R1:%5.1f R2:%5.1f R3:%5.1f%n",ZWN.R1, ZWN.R2,
ZWN.R3);// fprintf(fp,"\t Lo:%5.1f Gw:%f %n",ZWN.Lo, ZWN.Gw);// fprintf(fp,"\t Thickness : %f%n",Ndit);
// }
// }
// else fprintf(fp,"\tNasal File : <None> %n%n");// fclose(fp);//////--------------- End of main. --------------------------

534
535
536
537
538
539
540
541
542
543
544
545

32

private void VocalAndNose(Complex hvtra, Complex Nhvtra){double area1=radin[0]*radin[0]*PI;
double Varea1=Vradin[0]*Vradin[0]*PI;
Complex[] gam = new Complex[nareav];Complex[] fimp = new Complex[nareav];Complex[] Vgam = new Complex[Nsec];
Complex[] Vfimp = new Complex[Nsec];
Complex Vzl, zl;
Complex Al, Bl, Cl, Dl;
Complex Av, Bv, Cv, Dv;
Complex An, Bn, Cn, Dn;

付 録B

546

声道伝達特性のプログラム

33

Complex KKnose, KKlips, FC, FD;

547

550

if(icntlr <= 1){
newsec01(radin, nareav, gam, fimp, sfact);
newsec01(Vradin, Nsec, Vgam, Vfimp, Vsfact);

551

}

552

else{
newsec23(radin, radout, nareav, gam, fimp, sfact);
newsec23(Vradin, Vradout, Nsec, Vgam, Vfimp, Vsfact);

548
549

553
554
555
556
557
558
559
560
561
562
563

}
// ----------------------------------------------RadiationZ radimpV = new RadiationZ(Varea1, freq, temper);radimpV.calculate();Vzl = Complex.valueOf(radimpV.getRadr(), radimpV.getRadi());//radimp.setArea(area1);
RadiationZ radimp = new RadiationZ(area1, freq, temper);radimp.calculate();zl = Complex.valueOf(radimp.getRadr(), radimp.getRadi());

564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582

CalcF calfl = new CalcF(gam, fimp, 0, Nport-1, vlg);
calfl.calculate();Al = calfl.getA();
Bl = calfl.getB();
Cl = calfl.getC();
Dl = calfl.getD();
CalcF calfv = new CalcF(gam, fimp, Nport, nareav-1, vlg);
calfv.calculate();Av = calfv.getA();
Bv = calfv.getB();
Cv = calfv.getC();
Dv = calfv.getD();
CalcF calfn = new CalcF(Vgam, Vfimp, 0, Nsec-1, Nvlg);
calfn.calculate();An = calfn.getA();
Bn = calfn.getB();
Cn = calfn.getC();
Dn = calfn.getD();

583
584
585

586

KKnose = ((Cn.multiply(Vzl)).add(Dn)).divide((An.multiply(Vzl)).add(Bn));
KKlips = ((Cl.multiply(zl)).add(Dl)).divide((Al.multiply(zl)).add(
Bl));

587
588
589
590
591

/* ----------------------------------H(w)=Hv(w)+sum(n=0,to nareav-1) Hn(w)----------------------------------- */付 録B

声道伝達特性のプログラム

34

FC = (Al.multiply(Cv)).add(Dv.multiply((KKnose.multiply(Al)).add(Cl

592

)));
FD = (Bl.multiply(Cv)).add(Dv.multiply((KKnose.multiply(Bl)).add(Dl)));
hvtra = Complex.valueOf(1.0).divide((FC.multiply(zl)).add(FD));

593

594
595
596

FC = (An.multiply(Cv)).add(Dv.multiply((KKlips.multiply(An)).add(Cn597

FD = (Bn.multiply(Cv)).add(Dv.multiply((KKlips.multiply(Bn)).add(Dn598

Nhvtra = Complex.valueOf(1.0).divide((FC.multiply(zl)).add(FD));

)));
)));
599

}

600
601
602
603
604
605
606
607

//private void Vocal(Complex hvtra) {private void Vocal() {double area1=radin[0]*radin[0]*PI;
Complex[] gam = new Complex[nareav];Complex[] fimp = new Complex[nareav];Complex Vzl, zl;
Complex Cv, Dv;

608

if(icntlr <=1)
newsec01(radin, nareav, gam, fimp, sfact);
else
newsec23(radin, radout, nareav, gam, fimp, sfact);

609
610
611
612
613

RadiationZ radimp = new RadiationZ(area1, freq, temper);radimp.calculate();zl = Complex.valueOf(radimp.getRadr(), radimp.getRadi());

614
615
616
617

CalcF calf = new CalcF(gam, fimp, 0, nareav-1, vlg);
calf.calculate();Cv = calf.getC();
Dv = calf.getD();

618
619
620
621

/* ----------------------------------H(w)=Hv(w)+sum(n=0,to nareav-1) Hn(w)----------------------------------- */hvtra = Complex.valueOf(1.0).divide((Cv.multiply(zl)).add(Dv));

622
623
624
625
626

}

627
628
629

630
631
632
633
634

private void newsec01(double[] radin, int nsec, Complex[] gamma,
Complex[] z, double[] sfact){
double Ai,Li,Gi,Ri,Ci,ome,ome2,S;Complex Z2, gam;double mu=1.86e-4;
double lambda=5.5e-5;
double cp=0.24;

付 録B

635
636
637

声道伝達特性のプログラム

35

double eta=1.4;// const double c0=34359.0;// const double rho=1.2046e-3;

638
639
640

double c0=3.3145e4+60.7*temper;
double rho=1.293e-3/(1.0+0.00367*temper);

641
642
643

ome=PI*2.0*freq;
ome2=ome*ome;

644
645
646
647

for(int i=0; i<nsec ;i++){Ai = radin[i] * radin[i] * PI;S = 2.0 * PI * radin[i] * sfact[i];

648

652

if(Ai < 1.0e-4){
System.out.println("Too small ! -> A[" + i + "]" + Ai);
gamma[i] = Complex.valueOf(1.0e-10);z[i] = Complex.valueOf(1.0e10);

653

}

654

else{
switch (icntlr){
case 0:z[i] = Complex.valueOf(rho*c0/Ai,0.0);
gamma[i] = Complex.valueOf(0.0,ome/c0);
break;

649
650
651

655
656
657
658
659
660
661
662
663
664
665

666

667
668

669
670
671
672
673
674
675
676
677

case 1:
Li = rho/Ai;
Ci = Ai/(rho*c0*c0);
Ri = (S/(Ai*Ai))*sqrt(ome*rho*mu/2.0);
Gi = (S*(eta-1.0)/(rho*c0*c0))*sqrt(lambda*ome/(2.0*cp*rho));
gam = Complex.valueOf(Ri,ome*Li).multiply(Complex.valueOf(Gi,ome*Ci));gamma[i] = gam.sqrt();
Z2 = Complex.valueOf(Ri,ome*Li).divide(Complex.valueOf(Gi,ome*Ci));z[i] = Z2.sqrt();break;
default:
System.out.println("icntlr is not 0 nor 1...");
break;} //<-switch} //<-else} //<-for} // ----------------------- End of newsec01 -------

678
679

///************************************

付 録B

680
681
682
683
684

685
686
687
688
689
690

声道伝達特性のプログラム

36// newsec23 calculates Zi, ri in the vocal tract// from various parameters (R1, G, etc...)// ********************************************///
private void newsec23(double[] radin, double[] radout, int nsec,
Complex[] gamma, Complex[] z, double[] sfact){
double Ai,Li,Gi,Ri,Ci,ome,ome2,Lw,S,dit,xx,yy1,yy,bunbo,bunb2;
Complex Yw, Z2, gam;double mu=1.86e-4;
double lambda=5.5e-5;
double cp=0.24;
double eta=1.4;

691
692
693

double c0 = 3.3145e4 + 60.7*temper;
double rho = 1.293e-3 / (1.0+0.00367*temper);

694
695
696

ome=PI*2.0*freq;
ome2=ome*ome;

697
698
699
700
701
702

for(int i=0; i<nsec ;i++){dit = (radout[i]) - (radin[i]);Ai = radin[i] * (radin[i]) * PI;S = 2.0 * PI * (radin[i]) * (sfact[i]);Lw = ZW.getLo();

703

707

if(Ai < 1.0e-4){
System.out.println("Too small ! -> A[" + i + "]" + Ai);
gamma[i] = Complex.valueOf(1.0e-10);z[i] = Complex.valueOf(1.0e10);

708

}

709

else{
switch (icntlr){
case 2:/******* Wall impedance *******/
if(dit==0.0)
System.out.println("Thickness =0 ! Something is wrong

704
705
706

710
711
712
713
714

715

716
717

718
719

...");
bunbo = (ZW.getGw() * ZW.getGw()) + ome2 * ZW.getR1() * ZW.
getR1();
bunb2 = ZW.getR3() * ZW.getR3() + ome2 * Lw * Lw;
xx=((ZW.getGw() * ZW.getGw() * ZW.getR1() / bunbo) + ZW.
getR2()) / dit + ome2 * Lw * Lw * ZW.getR3() / bunb2;
yy1 = (-1.0) * ome * ZW.getGw() * ZW.getR1() * ZW.getR1();
yy = yy1 / bunbo / dit + ome * Lw * ZW.getR3() * ZW.getR3()/ bunb2;

720
721
722

Yw = Complex.valueOf(S,0.0).divide(Complex.valueOf(xx,yy));
Li= rho / Ai;

付 録B

声道伝達特性のプログラム

37

Ci = Ai / (rho * c0 * c0);
Ri=(S / (Ai * Ai)) * sqrt(ome * rho * mu / 2.0);
Gi=(S * (eta - 1.0)/(rho * c0 * c0)) * sqrt(lambda * ome /(2.0 * cp * rho));gam = Complex.valueOf(Ri, ome * Li).multiply(Complex.valueOf(
Gi, ome * Ci).add(Yw));
gamma[i] = gam.sqrt();
Z2 = Complex.valueOf(Ri, ome * Li).divide(Complex.valueOf(Gi,ome * Ci).add(Yw));z[i] = Z2.sqrt();break;

723
724
725

726

727
728

729
730
731
732
733
734
735
736
737
738

739

740
741

742
743

case 3:xx=ZW.getR1(); yy = ome * ZW.getLo() - ZW.getGw() / ome;
Yw=Complex.valueOf(S,0.0).divide(Complex.valueOf(xx, yy));
Li = rho / Ai;
Ci = Ai / (rho * c0 * c0);
Ri = (S / (Ai * Ai)) * sqrt(ome * rho * mu / 2.0);
Gi = (S * (eta - 1.0) / (rho * c0 * c0)) * sqrt(lambda * ome/ (2.0 * cp * rho));gam = Complex.valueOf(Ri, ome * Li).multiply(Complex.valueOf(
Gi, ome * Ci).add(Yw));
gamma[i] = gam.sqrt();
Z2 = Complex.valueOf(Ri, ome * Li).divide(Complex.valueOf(Gi,ome * Ci).add(Yw));z[i] = Z2.sqrt();break;

744
745
746
747
748
749
750
751

default:
System.out.println(" icntlr is not 2 nor 3 ...");
break;} //<-switch} //<-else} //<-for} // ----------------------- End of newsec23 -------

752
753
754
755
756
757
758
759
760
761

public static void main(String[] args){
int dataset = 10000;
for(int i=0; i<dataset; i++){
System.out.println(args);
File file = new File(args[0]);if(!file.exists()) {System.out.println("ファイルが存在しません。");
return;

}

762
763
764

try {

付 録B

声道伝達特性のプログラム

772

Peaks peaks = new Peaks(file, i);if (args.length == 2)
if(args[1].equals("-h")) peaks.showHelp();
peaks.showStatus();peaks.calculate();peaks.saveFile();} catch (IOException e) {e.printStackTrace();

773

}

765
766
767
768
769
770
771

}

774

}

775
776

}

以下は声道伝達特性を求めるプログラムの補足プログラム（CalcF.java)
1
2

package jp.ac.hus.mhlab;

3
4
5
6
7
8

import java.io.*;import static java.lang.Math.*;import org.apache.commons.math3.transform.*;
import org.apache.commons.math3.complex.*;import org.apache.commons.math3.exception.*;

9
10
11
12
13
14

public class CalcF {private Complex A, B, C, D;
private Complex[] gamma, fimp;
private int nbgn, nend;
private double vlg;

15

21

public CalcF(Complex[] gamma, Complex[] fimp, int nbgn, int nend,
double vlg) {this.gamma = gamma;
this.fimp = fimp;
this.nbgn = nbgn;
this.nend = nend;
this.vlg = vlg;

22

}

16

17
18
19
20

23

25

public Complex getA(){
return A;

26

}

27

24

28

public Complex getB(){
return B;

29

}

30
31

public Complex getC(){
return C;

32

}

38

付 録B

声道伝達特性のプログラム

33
34

public Complex getD(){
return D;

35

}

39

36
37
38
39
40
41
42
43
44
45
46
47
48

public void calculate(){Complex tanal,pk2;
Complex[] a = new Complex[nend+1];
Complex[] b = new Complex[nend+1];
Complex[] c = new Complex[nend+1];
Complex[] d = new Complex[nend+1];
Complex[] fa = new Complex[nend+1];
Complex[] fb = new Complex[nend+1];
Complex[] fc = new Complex[nend+1];
Complex[] fd = new Complex[nend+1];
double ei,er;
int ise,iid;49

for(int i=nbgn; i<=nend; i++){er = (gamma[i].multiply(vlg)).getReal();ei = (gamma[i].multiply(vlg)).getImaginary();

50
51
52
53

tanal = Complex.valueOf(tanh(er)*cos(ei),sin(ei)).divide(Complex.
valueOf(cos(ei),tanh(er)*sin(ei)));//tanal=(std::complex<double>(tanh(er)*cos(ei),sin(ei))/std::
complex<double>(cos(ei),tanh(er)*sin(ei)));pk2=Complex.valueOf(1.0).divide(tanal.sqrt1z());//pk2=1.0/sqrt(1.0-tanal*tanal);

54

55

56
57
58

a[i]=pk2;
b[i]=pk2.multiply(fimp[i].multiply(tanal));c[i]=pk2.multiply(tanal.divide(fimp[i]));d[i]=pk2; //Initialize...
//a[i]=pk2 ; b[i]=pk2*fimp[i]*tanal;//c[i]=pk2*tanal/fimp[i] ; d[i]=pk2; //Initialize...

59
60
61
62
63
64
65
66

}

67
68
69
70

/*****-------------------------Calculate Whole F_Matrix--------------------------********/

71
72
73
74
75

fa[nbgn]=a[nbgn];
fb[nbgn]=b[nbgn];
fc[nbgn]=c[nbgn];
fd[nbgn]=d[nbgn];

76
77
78

for(ise=nbgn, iid = nbgn+1 ; iid<= nend ; ise++,iid++){fa[iid]=(a[iid].multiply(fa[ise])).add(b[iid].multiply(fc[ise]));

付 録B

声道伝達特性のプログラム

40

81

fb[iid]=(a[iid].multiply(fb[ise])).add(b[iid].multiply(fd[ise]));fc[iid]=(c[iid].multiply(fa[ise])).add(d[iid].multiply(fc[ise]));fd[iid]=(c[iid].multiply(fb[ise])).add(d[iid].multiply(fd[ise]));

82

}

79
80

83

A = fa[nend];
B = fb[nend];
C = fc[nend];
D = fd[nend];

84
85
86
87

}

88
89

}

以下は声道伝達特性を求めるプログラムの補足プログラム（RadiationZ.java)
1
2
3

package jp.ac.hus.mhlab;
import static java.lang.Math.*;

4
5
6
7

public class RadiationZ {private static final double SMALL = -1.0e36;
private double area, freq, temper, radr, radi;

8

14public RadiationZ(double area, double freq, double temper){this.area = area;
this.freq = freq;
this.temper = temper;
this.radr = 0.0;
this.radi = 0.0;

15

}

9
10
11
12
13

16

18

public double getRadr(){return radr;

19

}

17

20

22

public double getRadi(){return radi;

23

}

21

24

26

public void setArea(double area){this.area = area;

27

}

25

28
29
30
31
32
33
34

public void calculate(){double[] p = new double[10000];
double x, conv, d1, d2, df, d=0.0;
double xmyu, xreal, ximag, p1=0.0, p2=0.0, order=0.0, pp;
int m, norder;
final double C0 = 33145.0 + 60.7*temper; // cm

付 録B

35
36

声道伝達特性のプログラム

41

final double RHO = 1.293e-3 /(1.0 + 0.00367*temper); // cm
final int NPX = 41;

37
38
39
40

double rhead=8.0; // cm
x=sqrt(rhead*rhead-area/PI)/rhead; /* x: cos[theta0] */conv = 0.5/(1.0-x); /*const=1/{4*sin(theta0/2)*sin(theta0/2)}*/

41
42
43
44
45
46
47
48
49

/* Generate p[ ] */p[1] =1.0 ; p[2] =1.0 ; p[3] =x;for(int i=3 ; i <= NPX; i++){d1 = x*p[i] ;d2 = d1-p[i-1];
df = (double)i - 1.0;
p[i+1] = d1+d2 - d2/df ;
}

50
51

xmyu=2.0*PI*freq*rhead/C0;

52
53
54
55
56
57
58
59
60
61
62
63

64
65
66
67
68
69

/* myu = 2*pai*a/lamda= 2*pai*a*freq./c
a: = rhead*/for(m=0, xreal=0.0, ximag=0.0; m<=26 ; m++ ) {if(m == 0) {p1= - snbes(1,xmyu);
p2= sjbes(1,xmyu);}
else {
p1=(double)m*snbes(m-1,xmyu)-((double)m + 1.0)*snbes(m+1,xmyu);
p2=((double)m +1.0)*sjbes(m+1,xmyu)-(double)m*sjbes(m-1,xmyu);}if(abs(p1)>1.0e16 || abs(p2)>1.0e16) break;norder=(int)(log10(max(abs(p1),abs(p2))));
order=pow(10.0,(-norder));p1=p1*order; p2=p2*order;70

pp=p1*p1+p2*p2;
d=(p[m+1]-p[m+3])*(p[m+1]-p[m+3]);

71
72
73

/* P[m-1] --> p(m+1)
P[m+1] --> p(m+3) */

74
75
76

xreal =xreal+(2.0*(double)m +1.0)*d/pp*order*order;
ximag =ximag+d*(sjbes(m,xmyu)*p2-snbes(m,xmyu)*p1)/pp*order;

77
78
79

}

80
81

//System.out.printf("d p1 p2= %f, %f, %f %f %n", d,p1,p2,order) ;付 録B

声道伝達特性のプログラム

42

82
83
84
85
86
87
88

radr =conv*xreal/(xmyu*xmyu);radi =conv*ximag;
radr *= (conv*RHO*C0/(PI*rhead*rhead)); /* Radiation impedance */radi *= (conv*RHO*C0/(PI*rhead*rhead)); /* Radiation impedance */}/* =-=-=-=-=-=- End of radimp -=-=-=-=-=-=-=-=-=-=- */

89
90
91
92

private double sjbes(int n, double x) {double t1=0.0,t2,t3=0.0,dx,y,w,z,sj=0.0,dy,d55,xdum;
int lme=0,i,nm,k;93

97

xdum=x - 3.0e4;
if( ((n-30000) > 0) ||(xdum >= 0.0)||(x < 0.0) || (n < 0)){
System.out.println("Bad Arg.1!");
return 0.0;

98

}

94
95
96

99
100
101
102

d55=1.0e-36;dx=x;
xdum=x-7.0e-4;

103
104
105
106
107
108
109

if(xdum < 0.0){
w=1.0;
if(n == 0) return(w);if(dx <= 1.0e-37) return(0.0);
dy =1.0e-37/dx;
t1=3.0 ; t2 =1.0;

110

115

for(i=1; i<=n ; i++) {t3=t2*dx/t1;
t1=t1+2.0;
t2=t3;
if( abs(t3) <= t1*dy ) return(0.0);

116

}

111
112
113
114

return(t3);

117
118

}

119
120
121
122
123

if( (x-0.2) < 0.0 ) { y = dx*dx;
w = 1.0-y*(1.0-0.05*y)/6.0;}
else w = sin(dx)/dx;if( n == 0) return(w);

124
125
126
127
128
129

if( (x-100.0) >= 0.0) lme=(int)(0.5*x+5.0);else{ if((x-10.0) >= 0.0) lme=(int)(0.1*lme+10.0);else{ if( (x-1.0) >= 0.0) lme=(int)(0.02*x+18.0);
else lme=5; }}付 録B

声道伝達特性のプログラム

nm=max(n,(int)x)+ lme;
z =1.0/dx;
t3=0.0;t2=1.0e-36;

130
131
132
133
134

for(i=1 ; i<=nm ; i++){k=nm-i;t1=(2.0*(double)k+3.0)*z*t2-t3;
if( n-k == 0) sj=t1;

135
136
137
138
139
140

if((abs(t1)-1.0e36) >= 0.0) {t1=t1*d55;
t2=t2*d55;
if((n-k) >= 0) sj=sj*d55;}

141
142
143
144
145

t3=t2; t2=t1;}/* End of for*/

146
147
148

151

if(abs(w) <= 0.01){
w=w/dx-cos(dx)/dx;
return (sj/t3*w);

152

}

153

return (w/t1*sj);

149
150

154
155
156
157

}
/* ===== End of sjbes ===================### snbes ###
*/

158
159
160
161

private double snbes(int n,double x) {double dx,z,qn0,qn1,qn2=0.0,w,fdum;
int ndum,i,m;

162

166

ndum=n-30000;
if((ndum >= 0)||(x <= 0.0)||(n < 0)){
System.out.println("Bad Arg.2!");
return(0.0);

167

}

163
164
165

168
169
170

dx=x;
z=1.0/dx;

171
172
173
174
175

if( (x-7.0e-4) <= 0.0) {if(n == 0) return (-z);
m=(int)(30.0/log10(z));
if((n-m+1) > 0) return(SMALL);

176
177

for(i=1,qn0=z,qn1=1.0; i<=n; i++){

43

付 録B

声道伝達特性のプログラム

qn2=qn0*qn1*z;
qn1+=2.0;
qn0=qn2;

178
179
180

}

181

return (-qn2);

182

}

183
184

qn0= -z*cos(dx);if(n == 0) return(qn0);
qn1=z*(qn0-sin(dx));

185
186
187
188

ndum=n-1;
if( ndum < 0) return (qn0);
if( ndum == 0) return(qn1);

189
190
191
192

for(i=2; i <=n ; i++){fdum=abs(qn1)-1.0e36;
if(fdum >= 0.0) {w=qn1*1.0e-10;qn2=(2.0*(double)i-1.0)*z*w-qn0*1.0e-10;
if(abs(qn2) >= 1.0e36){
System.out.println("Bad ..3");
return(SMALL);

193
194
195
196
197
198
199
200

}

201

qn2*=1.0e10;

202
203

}

204

else qn2=(2.0*(double)i-1.0)*z*qn1-qn0;

205

qn0=qn1; qn1=qn2;

206

}

207

return(qn2);

208

}

209
210

}

以下は声道伝達特性を求めるプログラムの補足プログラム（WallImpedance.java)
1
2

package jp.ac.hus.mhlab;

345public class WallImpedance {private double R1, R2, R3, Lo, Gw;

6
7

8
9
10
11

public WallImpedance(double R1, double R2, double R3, double Lo,
double Gw) {this.R1 = R1;
this.R2 = R2;
this.R3 = R3;
this.Lo = Lo;

44

付 録B

声道伝達特性のプログラム

this.Gw = Gw;

12

}

13
14

16

public double getR1(){
return R1;

17

}

18

15

19

public double getR2(){
return R2;

20

}

21
22

public double getR3(){
return R3;

23

}

24
25

public double getLo(){return Lo;

26

}

27
28

public double getGw(){
return Gw;

29

}

30
31

public void setR1(double R1){
this.R1 = R1;

32

}

33
34

public void setLo(double Lo){this.Lo = Lo;

35

}

36
37

public void setGw(double Gw){
this.Gw = Gw;

38

}

39

}

45

46

付 録C

ニューラルネットワークのプログ
ラム

ニューラルネットワークのプログラムを以下に示す。
1
2
3
4
5
6
7
8
9

# coding: utf-8
import sys, os
sys.path.append(os.pardir)
import matplotlib.pyplot as plt
import numpy as np
import pickle
from dataset import load_data #要変更
from two_layer_net import TwoLayerNetfrom common.optimizer import Adam

データ取得

10
11
12

# データの読み込み(x_train, t_train), (x_test, t_test) = load_data()

13
14
15

network = TwoLayerNet(input_size=1024, hidden_size=50, output_size=34)
optimizer = Adam(lr=0.001)

16
17
18
19
20

iters_num = 10000
train_size = x_train.shape[0]test_size = x_test.shape[0]batch_size = 100

21
22
23
24
25

train_loss_list = []test_loss_list = []train_acc_list = []test_acc_list = []26
27

iter_per_epoch = max(train_size / batch_size, 1)

28
29
30
31
32

for i in range(iters_num):
batch_mask = np.random.choice(train_size, batch_size)
x_batch = x_train[batch_mask]t_batch = t_train[batch_mask]

33
34
35
36

batch_mask_for_test = np.random.choice(test_size, batch_size)
x_t_batch = x_test[batch_mask_for_test]
t_t_batch = t_test[batch_mask_for_test]

付 録C

ニューラルネットワークのプログラム

47

37
38
39
40

# 勾配
grad = network.gradient(x_batch, t_batch)
optimizer.update(network.params, grad)

41
42
43

loss = network.loss(x_batch, t_batch)
loss_test = network.loss(x_t_batch, t_t_batch)

44
45
46

train_loss_list.append(loss)
test_loss_list.append(loss_test)

47
48
49
50
51
52

if i % iter_per_epoch == 0:
train_acc = network.accuracy(x_train, t_train)
test_acc = network.accuracy(x_test, t_test)
train_acc_list.append(train_acc)
test_acc_list.append(test_acc)

53
54

55

with open(’C:￥￥Users￥￥Masataka￥￥Desktop￥￥akemi￥￥grad_sample.pkl’,mode=’wb’) as f:
pickle.dump(grad, f)

56
57
58
59
60
61
62

plt.figure()
plt.xlabel(’iterations’)
plt.ylabel(’loss’)
plt.plot(train_loss_list, label = ’train’, linestyle = ’dashed’)
plt.plot(test_loss_list, label = ’test’, linestyle = ’solid’)
plt.legend()

63
64
65
66
67
68
69
70

plt.figure()
plt.xlabel(’epochs’, fontsize = 18)
plt.ylabel(’accuracy’, fontsize = 18)
plt.plot(train_acc_list, label = ’train’, linestyle = ’dashed’)
plt.plot(test_acc_list, label = ’test’, linestyle = ’solid’)
plt.legend()
plt.show()

以下はニューラルネットワークの補足プログラム（dataset.py）
1234

import glob
import os
import numpy as np
import random

5
6
7
8

def load_data():train_peaks = []for filepath in glob.glob(’C:￥￥
Users￥￥Masataka￥￥Documents￥￥NetBeansProjects￥￥Peaks￥￥peaks￥￥*.v’):

付 録C

9
10
11
12

ニューラルネットワークのプログラム

48

with open(filepath, ’r’) as f:l = [s.strip() for s in f.readlines()]l_f = [float(s) for s in l]train_peaks.append(l_f)

13
14
15

random.seed(10)t_p = np.array(random.sample(train_peaks, 7000))

16
17
18

train_area = []for filepath in glob.glob(’C:￥￥
Users￥￥Masataka￥￥Documents￥￥NetBeansProjects￥￥Peaks￥￥area￥￥

*’):
19
20
21
22
23

with open(filepath, ’r’) as f:l = [s.strip() for s in f.readlines()]l_f = [float(s) for s in l]
del l_f[0:2]
train_area.append(l_f)

24
25
26

random.seed(10)t_a = np.array(random.sample(train_area, 7000))

27
28
29

30
31
32
33

test_peaks = []for filepath in glob.glob(’C:￥￥
Users￥￥Masataka￥￥Documents￥￥NetBeansProjects￥￥Peaks￥￥peaks￥￥*.v’):
with open(filepath, ’r’) as f:l = [s.strip() for s in f.readlines()]l_f = [float(s) for s in l]
test_peaks.append(l_f)

34
35
36

random.seed(1)ts_p = np.array(random.sample(test_peaks, 3000))

37
38
39

test_area = []for filepath in glob.glob(’C:￥￥
Users￥￥Masataka￥￥Documents￥￥NetBeansProjects￥￥Peaks￥￥area￥￥

*’):
40
41
42
43
44

with open(filepath, ’r’) as f:l = [s.strip() for s in f.readlines()]l_f = [float(s) for s in l]
del l_f[0:2]
test_area.append(l_f)

45
46
47

random.seed(1)ts_a = np.array(random.sample(test_area, 3000))

48
49

return (t_p, t_a),(ts_p, ts_a)

付 録C

ニューラルネットワークのプログラム

49

以下はニューラルネットワークの補足プログラム（two layer net.py）
1
2
3

4
5
6
7

# coding: utf-8
import sys, os
sys.path.append(os.pardir) # 親ディレクトリのファイルをインポートするための設定
import numpy as np
from common.layers import *from common.gradient import numerical_gradientfrom collections import OrderedDict89class TwoLayerNet:

10
11

def __init__(self, input_size, hidden_size, output_size,
weight_init_std = 0.01):

12# 重みの初期化

13

self.params = {}self.params[’W1’] = np.random.randn(input_size, hidden_size) /np.sqrt(input_size)
self.params[’b1’] = np.zeros(hidden_size)
self.params[’W2’] = np.random.randn(hidden_size, output_size) /np.sqrt(hidden_size)
self.params[’b2’] = np.zeros(output_size)

14

15
16

17
18
19
20
21

22
23
24

#レイヤの生成self.layers = OrderedDict()
self.layers[’Affine1’] = Affine(self.params[’W1’], self.params[’b1’])
self.layers[’BatchNormalization’] = BatchNormalization(1, 0)
self.layers[’Sigmoid1’] = Sigmoid()
self.layers[’Affine2’] = Affine(self.params[’W2’], self.params[’b2’])

25
26

self.lastLayer = IdentityWithLoss()

27
28
29
30

def predict(self, x):
for layer in self.layers.values():x = layer.forward(x)3132

return x

33
34
35
36
37

# x:入力データ, t:教師データ
def loss(self, x, t):
y = self.predict(x)return self.lastLayer.forward(y, t)

38
39
40
41

def accuracy(self, x, t):
y = self.predict(x)y = np.argmax(y, axis=1)

付 録C

42

ニューラルネットワークのプログラム

if t.ndim != 1 : t = np.argmax(t, axis=1)

43
44
45

accuracy = np.sum(y == t) / float(x.shape[0])return accuracy

46
47
48
49

# x:入力データ, t:教師データ
def numerical_gradient(self, x, t):
loss_W = lambda W: self.loss(x, t)

50
51
52
53
54
55

grads = {}grads[’W1’] = numerical_gradient(loss_W, self.params[’W1’])
grads[’b1’] = numerical_gradient(loss_W, self.params[’b1’])
grads[’W2’] = numerical_gradient(loss_W, self.params[’W2’])
grads[’b2’] = numerical_gradient(loss_W, self.params[’b2’])

56
57

return grads

58
59
60
61

def gradient(self, x, t):# forward
self.loss(x, t)

62
63
64
65

# backward
dout = 1
dout = self.lastLayer.backward(dout)

66
67
68
69
70

layers = list(self.layers.values())
layers.reverse()
for layer in layers:
dout = layer.backward(dout)

71
72
73
74

75

# 設定
grads = {}grads[’W1’], grads[’b1’] = self.layers[’Affine1’].dW, self.
layers[’Affine1’].db
grads[’W2’], grads[’b2’] = self.layers[’Affine2’].dW, self.
layers[’Affine2’].db

76
77

return grads

以下はニューラルネットワークの補足プログラム（optimizer.py）
1# coding: utf-8

2

import numpy as np

3
4class SGD:

567"""確率的勾配降下法（Stochastic Gradient Descent）"""

50

付 録C

ニューラルネットワークのプログラム

8

def __init__(self, lr=0.01):
self.lr = lr

9

51

10
11
12
13

def update(self, params, grads):
for key in params.keys():params[key] -= self.lr * grads[key]

14
15
16

class Momentum:

17
18"""Momentum SGD"""

19
20
21
22
23

def __init__(self, lr=0.01, momentum=0.9):
self.lr = lr
self.momentum = momentum
self.v = None

24
25
26
27
28
29

def update(self, params, grads):
if self.v is None:self.v = {}for key, val in params.items():
self.v[key] = np.zeros_like(val)

30
31
32

for key in params.keys():self.v[key] = self.momentum*self.v[key] - self.lr*grads[key]33params[key] += self.v[key]

34
35
36

class Nesterov:

37
38"""Nesterov’s Accelerated Gradient (http://arxiv.org/abs/1212.0901)"""

39
40
41
42
43

def __init__(self, lr=0.01, momentum=0.9):
self.lr = lr
self.momentum = momentum
self.v = None

44
45
46
47
48
49

def update(self, params, grads):
if self.v is None:self.v = {}for key, val in params.items():
self.v[key] = np.zeros_like(val)

50
51
52
53

for key in params.keys():self.v[key] *= self.momentum
self.v[key] -= self.lr * grads[key]付 録C

ニューラルネットワークのプログラム

52

params[key] += self.momentum * self.momentum * self.v[key]params[key] -= (1 + self.momentum) * self.lr * grads[key]

54
55
56
57
58

class AdaGrad:

59
60"""AdaGrad"""

61
62
63
64

def __init__(self, lr=0.01):
self.lr = lr
self.h = None

65
66
67
68
69
70

def update(self, params, grads):
if self.h is None:self.h = {}for key, val in params.items():
self.h[key] = np.zeros_like(val)

71
72
73
74

for key in params.keys():self.h[key] += grads[key] * grads[key]params[key] -= self.lr * grads[key] / (np.sqrt(self.h[key])+ 1e-7)75
76
77

class RMSprop:

78
79"""RMSprop"""

80
81
82
83
84

def __init__(self, lr=0.01, decay_rate = 0.99):
self.lr = lr
self.decay_rate = decay_rate
self.h = None

85
86
87
88
89
90

def update(self, params, grads):
if self.h is None:self.h = {}for key, val in params.items():
self.h[key] = np.zeros_like(val)

91
92
93
94

95

for key in params.keys():self.h[key] *= self.decay_rate
self.h[key] += (1 - self.decay_rate) * grads[key] * grads[
key]
params[key] -= self.lr * grads[key] / (np.sqrt(self.h[key])+ 1e-7)96
97
98

class Adam:

付 録C

ニューラルネットワークのプログラム

53

99
100

"""Adam (http://arxiv.org/abs/1412.6980v8)"""

101
102
103
104
105
106
107
108

def __init__(self, lr=0.001, beta1=0.9, beta2=0.999):
self.lr = lr
self.beta1 = beta1
self.beta2 = beta2
self.iter = 0
self.m = Noneself.v = None

109
110
111
112
113
114
115

def update(self, params, grads):
if self.m is None:self.m, self.v = {}, {}for key, val in params.items():
self.m[key] = np.zeros_like(val)self.v[key] = np.zeros_like(val)

116
117
118

self.iter += 1
lr_t = self.lr * np.sqrt(1.0 - self.beta2**self.iter) / (1.0 self.beta1**self.iter)

119
120
121

122

123

for key in params.keys():#self.m[key] = self.beta1*self.m[key] + (1-self.beta1)*
grads[key]#self.v[key] = self.beta2*self.v[key] + (1-self.beta2)*(
grads[key]**2)
self.m[key] += (1 - self.beta1) * (grads[key] - self.m[key])124

self.v[key] += (1 - self.beta2) * (grads[key]**2 - self.v[
key])

125
126

params[key] -= lr_t * self.m[key] / (np.sqrt(self.v[key]) +
1e-7)

127
128

129

130

#unbias_m += (1 - self.beta1) * (grads[key] - self.m[key])# correct bias
#unbisa_b += (1 - self.beta2) * (grads[key]*grads[key] self.v[key]) # correct bias
#params[key] += self.lr * unbias_m / (np.sqrt(unbisa_b) + 1
e-7)

54

謝辞

本研究を進めるにあたり，終始，御指導，御鞭撻をいただいた北海道科学大学工学部 情報工学科 音声処理分野, 松﨑教授に心より感謝致します．
本研究の遂行にあたり様々な御助言，御協力をいただいた北海道科学大学工学部 情報工学科

松﨑ゼミの皆さんに深く感謝致します．

最後に本研究の遂行にあたり著者を常に支援し応援してくれた両親，並びに親族のみなさ
まに心より感謝申しあげます．