特許 6152690 音響解析装置

(19)【発行国】日本国特許庁(JP)

(12)【公報種別】特許公報(B2)

(11)【特許番号】6152690

(24)【登録日】2017年6月9日

(45)【発行日】2017年6月28日

(54)【発明の名称】音響解析装置

(51)【国際特許分類】

   G10L 25/90 20130101AFI20170619BHJP

   G10L 25/27 20130101ALI20170619BHJP

   G10L 25/45 20130101ALI20170619BHJP

【ＦＩ】

   G10L25/90

   G10L25/27

   G10L25/45

【請求項の数】5

【全頁数】17

(21)【出願番号】特願2013-97088(P2013-97088)

(22)【出願日】2013年5月2日

(65)【公開番号】特開2014-219481(P2014-219481A)

(43)【公開日】2014年11月20日

【審査請求日】2016年3月22日

(73)【特許権者】

【識別番号】000004075

【氏名又は名称】ヤマハ株式会社

(74)【代理人】

【識別番号】100125689

【弁理士】

【氏名又は名称】大林 章

(74)【代理人】

【識別番号】100121108

【弁理士】

【氏名又は名称】高橋 太朗

(72)【発明者】

【氏名】フェルナンド ヴィラヴィセンシオ

【審査官】 菊池 智紀

(56)【参考文献】

【文献】 特開平０６−１６１４９４（ＪＰ，Ａ）

【文献】 特開２００２−１７５０９９（ＪＰ，Ａ）

【文献】 特開２００４−１０９７４２（ＪＰ，Ａ）

【文献】 米国特許第６１３８０９２（ＵＳ，Ａ）

(58)【調査した分野】（Int.Cl.，ＤＢ名）

Ｇ１０Ｌ １９／００−２５／９０

(57)【特許請求の範囲】

【請求項1】

音響信号の基本周波数に対応する周期の整数倍に相当する基準長を含む探索範囲内の複数の候補長の各々について、前記音響信号の特定点を中心として当該候補長にわたる候補区間内における調波成分の優勢度を示す調波強度指標を算定する指標算定手段と、
前記複数の候補長のうち前記指標算定手段が候補長毎に算定した調波強度指標に応じて選択した候補長から解析長を特定する解析長特定手段と、
前記音響信号のうち前記解析長特定手段が特定した解析長に対応する解析区間を解析することで、前記音響信号の調波成分が優勢な調波帯域と調波帯域の高域側の非調波帯域との境界である調波境界周波数を特定する境界特定手段と
を具備する音響解析装置。

【請求項2】

前記解析長特定手段は、前記複数の候補長のうち前記調波強度指標が示す調波成分の優勢度が最大となる候補長を前記解析長として選択する
請求項１の音響解析装置。

【請求項3】

前記境界特定手段は、
前記音響信号のうち前記解析長特定手段が特定した解析長に対応する解析区間を解析することで、前記調波帯域と前記非調波帯域との境界に対応する調波成分の次数である調波次数を算定する調波次数解析手段と、
前記調波次数解析手段が算定した調波次数と前記音響信号の基本周波数とに応じて前記調波境界周波数を算定する周波数算定手段とを含む
請求項１または請求項２の音響解析装置。

【請求項4】

前記解析長特定手段が特定した解析長に応じて前記音響信号の基本周波数を補正する周波数補正手段を具備し、
前記周波数算定手段は、前記調波次数解析手段が算定した調波次数と前記周波数補正手段による補正後の基本周波数とに応じて前記調波境界周波数を算定する
請求項３の音響解析装置。

【請求項5】

前記調波次数解析手段は、前記特定点を包含するように時間軸上の位置を相違させた前記解析長に対応する複数の解析区間の各々について、前記調波帯域と前記非調波帯域との境界に対応する調波成分の次数である暫定調波次数を算定する第１演算手段と、
前記複数の解析区間について前記第１演算手段が算定した複数の暫定調波次数から前記調波次数を算定する第２演算手段とを含む
請求項３または請求項４の音響解析装置。

【発明の詳細な説明】

【技術分野】

【0001】

本発明は、音響信号を解析する技術に関する。

【背景技術】

【0002】

音響信号のスペクトルのうち調波成分が優勢な低域側の調波帯域と非調波成分が優勢な高域側の非調波帯域との境界となる周波数（以下「調波境界周波数」という）を音響信号から推定する技術が従来から提案されている。例えば非特許文献１には、音響信号のうち基本周波数に対応した周期（以下「ピッチ周期」という）に相当する時間長の解析区間を対象として、調波帯域内の各調波成分の合計強度（累積エネルギー）と非調波帯域内の非調波成分の合計強度との加重和が最大となるように調波境界周波数（VCO：Voicing Cut-Off frequency）を推定する技術が開示されている。

【先行技術文献】

【非特許文献】

【0003】

【非特許文献1】Kris Hermus, et. al., "ESTIMATION OF THE VOICING CUT-OFF FREQUENCY CONTOUR OF NATURAL SPEECH BASED ON HARMONIC AND APERIODIC ENERGIES", ICASSP 2008, p.4473-p.4476, 2008

【発明の概要】

【発明が解決しようとする課題】

【0004】

しかし、特許文献１の技術では、音響信号から推定された基本周波数に対応する唯１個の解析区間を対象として解析処理が実行されるから、基本周波数の推定精度が低い場合に適切な調波境界周波数を高精度に推定することが困難であるという問題がある。以上の事情を考慮して、本発明は、音響信号の基本周波数の推定誤差に対して頑健に調波境界周波数を特定することを目的とする。

【課題を解決するための手段】

【0005】

以上の課題を解決するために、本発明の音響解析装置は、音響信号の基本周波数に対応する周期の整数倍に相当する基準長（例えば基準長Ｎ0(k)）を含む探索範囲内の複数の候補長（例えばＧ個の候補長ＮC(1)〜ＮC(G)）の各々について、音響信号の特定点（例えば特定点Ｔ(k)）を中心として当該候補長にわたる候補区間内における調波成分の優勢度を示す調波強度指標（例えば調波強度指標Ｈ(g)）を算定する指標算定手段（例えば指標算定部３２）と、複数の候補長のうち指標算定手段が候補長毎に算定した調波強度指標に応じて選択した候補長から解析長（例えば解析長Ｎ(k)）を特定する解析長特定手段（例えば解析長特定部３４）と、音響信号のうち解析長特定手段が特定した解析長に対応する解析区間（例えば解析長の半分の区間）を解析することで、音響信号の調波成分が優勢な調波帯域と調波帯域の高域側の非調波帯域との境界である調波境界周波数（例えば調波境界周波数ＦC(k)）を推定する境界特定手段（例えば境界特定部２８）とを具備する。以上の構成では、複数の候補長のうち調波強度指標に応じて選択された候補長から解析長が特定されるから、音響信号の基本周波数の推定誤差に対して頑健に調波境界周波数を特定することが可能である。なお、音響信号の特定点は、時間軸上の特定の時点を意味する。例えば音響信号の単位波形（基本周波数に対応する周期の波形）毎に設定されたピッチマークが特定点の好適例である。

【0006】

本発明の好適な態様において、解析長特定手段は、複数の候補長のうち調波強度指標が示す調波成分の優勢度が最大となる候補長を解析長として選択する。以上の態様では、複数の候補長のうち調波強度指標が示す調波成分の優勢度が最大となる候補長が解析長として選択されるから、解析長を簡便に特定できるという利点がある。

【0007】

本発明の好適な態様において、境界特定手段は、音響信号のうち解析長特定手段が特定した解析長の解析区間を解析することで、調波帯域と非調波帯域との境界に対応する調波成分の次数である調波次数（例えば調波次数ｈC(k)）を算定する調波次数解析手段（例えば調波次数解析部４４）と、調波次数解析手段が算定した調波次数と音響信号の基本周波数とに応じて調波境界周波数を算定する周波数算定手段（例えば周波数算定部４６）とを含む。更に好適な態様に係る音響解析装置は、解析長特定手段が特定した解析長に応じて音響信号の基本周波数を補正する周波数補正手段（例えば周波数補正部３６）を具備し、周波数算定手段は、調波次数解析手段が算定した調波次数と周波数補正手段による補正後の基本周波数とに応じて調波境界周波数を算定する。以上の態様では、解析長特定手段が特定した解析長に応じて基本周波数が補正され、補正後の基本周波数に応じて調波境界周波数が算定されるから、調波境界周波数を高精度に特定できるという利点がある。

【0008】

本発明の好適な態様において、調波次数解析手段は、特定点を包含するように時間軸上の位置を相違させた解析長に対応する複数の解析区間（例えばＭ個の解析区間Ｓ(1)〜Ｓ(M)）の各々について、調波帯域と非調波帯域との境界に対応する調波成分の次数である暫定調波次数（例えば暫定調波次数ｈA(m)）を算定する第１演算手段（例えば第１演算部５２）と、複数の解析区間について第１演算手段が算定した複数の暫定調波次数から調波次数を算定する第２演算手段（例えば第２演算部５４）とを含む。具体的には、第１演算手段は、暫定調波次数と、複数の候補次数の各々が当該暫定調波次数に該当する確度（例えば確度ｐ(m,j)）とを、複数の解析区間の各々について時間軸上の特定点毎に算定し、第２演算手段は、１個の特定点に対応した複数の解析区間にわたる暫定調波次数の代表値である代表調波次数（例えば代表調波次数ｈB(k)）と、当該複数の解析区間にわたる各候補次数の確度の代表値である候補次数毎の代表確度（例えば代表確度λ(k,j)）とを、時間軸上の特定点毎に順次に算定する代表値算定手段（例えば代表値算定部５４２）と、代表値算定手段が時間軸上の特定点毎に算定した代表調波次数と各候補次数の代表確度とを利用して、代表調波次数の時系列と比較して変動が平滑化された調波次数の時系列を算定する平滑処理手段（例えば平滑処理部５４４）とを含む。以上の構成の具体例において、代表値算定手段は、複数の解析区間にわたる暫定調波次数の平均値（例えば数式(14)）を代表調波次数として特定点毎に算定し、複数の解析区間にわたる各候補次数の確度の中央値（例えば数式(15)）を候補次数毎の代表確度として算定する。

【0009】

以上の各態様に係る音響解析装置は、音響信号の解析に専用されるＤＳＰ（Digital Signal Processor）などのハードウェア（電子回路）によって実現されるほか、ＣＰＵ（Central Processing Unit）等の汎用の演算処理装置とプログラムとの協働によっても実現される。本発明のプログラムは、コンピュータを、音響信号の基本周波数に対応する周期の整数倍に相当する基準長を含む探索範囲内の複数の候補長の各々について、音響信号の特定点を中心として当該候補長にわたる候補区間内における調波成分の優勢度を示す調波強度指標を算定する指標算定手段と、複数の候補長のうち指標算定手段が候補長毎に算定した調波強度指標に応じて選択した候補長から解析長を特定する解析長特定手段と、音響信号のうち解析長特定手段が特定した解析長に対応する解析区間を解析することで、音響信号の調波成分が優勢な調波帯域と調波帯域の高域側の非調波帯域との境界である調波境界周波数を特定する境界特定手段として機能させる。以上に例示したプログラムによれば、本発明の音響解析装置と同様の効果が実現される。なお、本発明のプログラムは、コンピュータが読取可能な記録媒体に格納された形態で提供されてコンピュータにインストールされ得る。記録媒体は、例えば非一過性（non-transitory）の記録媒体であり、ＣＤ-ＲＯＭ等の光学式記録媒体（光ディスク）が好例であるが、半導体記録媒体や磁気記録媒体等の公知の任意の形式の記録媒体を包含し得る。また、例えば、本発明のプログラムは、通信網を介した配信の形態で提供されてコンピュータにインストールされ得る。

【図面の簡単な説明】

【0010】

【図1】本発明の実施形態に係る音響解析装置のブロック図である。

【図2】音響信号の特定点および候補長の説明図である。

【図3】解析区間設定部のブロック図である。

【図4】指標算定部が実行する指標算定処理のフローチャートである。

【図5】複数の候補長を選択する処理の説明図である。

【図6】候補長と調波強度指標との関係を示すグラフである。

【図7】境界特定部のブロック図である。

【図8】複数の解析区間を設定する処理の説明図である。

【図9】調波次数解析部のブロック図である。

【図10】第１演算部が実行する調波解析処理のフローチャートである。

【図11】実施形態の効果を説明するためのグラフである。

【図12】実施形態の効果を説明するためのグラフである。

【発明を実施するための形態】

【0011】

図１は、本発明の好適な態様に係る音響解析装置１００のブロック図である。音響解析装置１００は、音響信号ｓ(n)のスペクトルのうち調波成分が優勢な低域側の調波帯域と非調波成分が優勢な高域側の非調波帯域との境界となる調波境界周波数（VCO：Voicing Cut-Off frequency）を特定する信号処理装置であり、図１に示すように、演算処理装置１２と記憶装置１４とを具備するコンピュータシステムで実現される。

【0012】

記憶装置１４は、演算処理装置１２が実行するプログラムや演算処理装置１２が使用する各種のデータを記憶する。例えば記憶装置１４は、音響解析装置１００による解析対象となる音響信号ｓ(n)を記憶する。音響信号ｓ(n)は、調波音を含有する音響の波形を示す時間領域のサンプル系列である。記号ｎは、時間軸上の各サンプルの番号を意味する。調波音は、基本周波数（ピッチ）に相当する間隔で複数の調波成分（すなわち基音成分および複数の倍音成分）を周波数軸上に配列した調波構造が観測される音響成分である。人間が発声した音声を収録した音響信号（音声信号）ｓ(n)を以下の説明では想定する。半導体記録媒体または磁気記録媒体等の公知の記録媒体や複数種の記録媒体の組合せが記憶装置１４として任意に採用される。

【0013】

演算処理装置１２は、記憶装置１４に記憶されたプログラムを実行することで、音響信号ｓ(n)の調波境界周波数を特定するための複数の機能（基本周波数推定部２２，特定点設定部２４，解析区間設定部２６，境界特定部２８）を実現する。なお、演算処理装置１２の各機能を複数の装置に分散した構成や、専用の電子回路（例えばＤＳＰ）が一部の機能を実現する構成も採用され得る。演算処理装置１２が特定した調波境界周波数は、音響信号ｓ(n)に対する各種の音響処理に適用される。例えば、音響信号ｓ(n)のうち調波境界周波数で規定される調波帯域内の成分に対して選択的に音響処理を実行する構成や、音響信号ｓ(n)のうち調波境界周波数で規定される調波帯域内の成分を解析する（例えば特徴量を抽出する）構成が採用され得る。

【0014】

基本周波数推定部２２は、音響信号ｓ(n)の基本周波数ＦAを所定の周期で順次に推定する。基本周波数ＦAの推定には公知の技術が任意に採用される。基本周波数推定部２２が順次に推定する基本周波数ＦAの時系列は、時間軸上で折線として表現される。

【0015】

特定点設定部２４は、音響信号ｓ(n)について時間軸上の特定点Ｔ(k)を順次に設定する（ｋ＝１,２,３,……）。記号ｋは各特定点Ｔ(k)の番号を意味する。各特定点Ｔ(k)は、図２に示すように、基本周波数ＦAに対応するピッチ周期毎（すなわち音響信号ｓ(n)の単位波形毎）に設定された時間軸上の時点（ピッチマーク）である。例えば、音響信号ｓ(n)が示す音声の発音時に発声者の声門が閉塞する時点（GCI: Glottal Closure Instant）または声門が開口する時点（GOI: Glottal Opening Instant）が特定点Ｔ(k)として設定される。声門が閉塞または開口する時点の推定には、例えば、Thomas Drugman, et. al., "Glottal Closure and Opening Instant Detection from Speech Signals", INTERSPEECH 2009, p.2891-p.2894, 2009に開示された技術が好適に利用される。特定点設定部２４が設定した特定点Ｔ(k)毎（すなわち音響信号ｓ(n)の単位波形毎）に調波境界周波数ＦC(k)が順次に推定される。

【0016】

図１の解析区間設定部２６は、音響信号ｓ(n)のうち調波境界周波数ＦC(k)の解析に適用される解析区間を規定する時間長（以下「解析長」という）Ｎ(k)を設定する。解析長（フレームサイズ）Ｎ(k)は、特定点設定部２４が設定した特定点Ｔ(k)毎に個別に設定され、音響信号ｓ(n)のサンプルの個数で表現される。図３は、解析区間設定部２６のブロック図である。図３に示すように、解析区間設定部２６は、指標算定部３２と解析長特定部３４と周波数補正部３６とを含んで構成される。

【0017】

指標算定部３２は、解析長Ｎ(k)の候補であるＧ個（Ｇは２以上の自然数）の候補長ＮC(1)〜ＮC(G)の各々について、当該候補長ＮC(g)（ｇ＝１〜Ｇ）を解析長Ｎ(k)として採択する妥当性の指標となる調波強度指標Ｈ(g)（Ｈ(1)〜Ｈ(G)）を算定する。特定点設定部２４が設定した特定点Ｔ(k)毎にＧ個の調波強度指標Ｈ(1)〜Ｈ(G)が順次に算定される。

【0018】

図４は、特定点設定部２４が設定した第ｋ番目の特定点Ｔ(k)について指標算定部３２が各候補長ＮC(g)の調波強度指標Ｈ(g)（Ｈ(1)〜Ｈ(G)）を算定する指標算定処理のフローチャートである。指標算定処理は、例えば利用者からの指示（調波境界周波数ＦC(k)の算定指示）を契機として特定点Ｔ(k)毎に順次に実行される。

【0019】

指標算定処理を開始すると、指標算定部３２は、音響信号ｓ(n)のうち特定点設定部２４が設定した特定点Ｔ(k)での基本周波数Ｆ0(k)を算定する（ＳA1）。具体的には、指標算定部３２は、以下の数式(1)で表現される通り、基本周波数推定部２２が推定した複数の基本周波数ＦAの時系列から特定点Ｔ(k)での基本周波数Ｆ0(k)を算定する。例えば、基本周波数ＦAが算定された各時点の間隔内に位置する特定点Ｔ(k)での基本周波数Ｆ0(k)が、特定点Ｔ(k)の前後の各時点の基本周波数ＦAを利用して補間される。基本周波数Ｆ(k)の算定（補間）には公知の補間技術が任意に採用される。

【数1】

【0020】

指標算定部３２は、以下の数式(2)の演算を実行することで、調波強度指標Ｈ(g)の算定の基礎（各候補長ＮC(g)の基準）となる基準長Ｎ0(k)を算定する（ＳA2）。

【数2】

数式(2)の記号ＦSは、音響信号ｓ(n)のサンプリング周波数である。数式(2)から理解される通り、基準長Ｎ0(k)は、音響信号ｓ(n)のうち時間軸上の特定点Ｔ(k)での基本周波数Ｆ0(k)に対応するピッチ周期（１/Ｆ0(k)）の整数倍に対応した時間長として定義され、当該時間長の区間内に存在するサンプルの個数として表現される。本実施形態では、数式(2)から理解される通り、基準長Ｎ0(k)を音響信号ｓ(n)のピッチ周期の２倍（偶数倍）とした場合を例示する。

【0021】

指標算定部３２は、基準長Ｎ0(k)を含む所定の数値範囲（以下「探索範囲」という）からＧ個の候補長ＮC(1)〜ＮC(G)を特定する（ＳA3）。具体的には、図５から理解される通り、基準長Ｎ0(k)を中心とする探索範囲Ｒ内から間隔δでＧ個の候補長ＮC(1)〜ＮC(G)が抽出される。間隔δは偶数（例えばδ＝２）に設定される。前述の通り基準長Ｎ0(k)は奇数であるから、Ｇ個の候補長ＮC(1)〜ＮC(G)の各々は奇数である。図５に例示される通り、探索範囲Ｒは、基準長Ｎ0(k)から所定値ｗを減算した数値｛Ｎ0(k)−ｗ｝を下限値とし、基準長Ｎ0(k)に所定値ｗを加算した数値｛Ｎ0(k)＋ｗ｝を上限値とする数値範囲（幅２ｗ）である。所定値ｗは、例えば基準長Ｎ0(k)に所定の正数（例えば０.１５）を乗算した数値に設定される。以上の説明から理解される通り、図２に示すように、相異なる候補長ＮC(g)に対応するＧ個の区間（以下「候補区間」という）Ｃ(1)〜Ｃ(G)が特定点Ｔ(k)を中心として設定される。

【0022】

指標算定部３２は、探索範囲Ｒから選択したＧ個の候補長ＮC(1)〜ＮC(G)の各々について調波強度指標Ｈ(g)（Ｈ(1)〜Ｈ(G)）を算定する（ＳA4）。具体的には、指標算定部３２は、以下の数式(3)の演算で各候補長ＮC(g)の調波強度指標Ｈ(g)を算定する。

【数3】

【0023】

数式(3)の記号Ｑ(g,i)は、音響信号ｓ(n)のうち特定点Ｔ(k)を中心とした候補長ＮC(g)の候補区間（フレーム）Ｃ(g)内のスペクトルにおける各周波数成分の強度（パワー）を意味する。記号ｉは、周波数軸上に離散的に設定された複数の周波数（周波数ビン）のうち任意の１個の周波数を指示する変数である。具体的には、指標算定部３２は、特定点Ｔ(k)を中心とする候補長ＮC(g)の窓幅の窓関数（典型的には矩形窓）を時間領域で音響信号ｓ(n)に乗算したうえで周波数領域のスペクトルに変換することで各周波数成分の強度Ｑ(g,i)を算定する。時間領域から周波数領域への変換には、例えば高速フーリエ変換が好適に利用される。なお、各周波数成分の振幅を強度Ｑ(g,i)として算定することも可能である。

【0024】

数式(3)から理解される通り、調波強度指標Ｈ(g)は、候補長ＮC(g)の候補区間Ｃ(g)内の音響信号ｓ(n)を構成する複数の周波数成分のうち周波数軸上の奇数番目の周波数成分の強度Ｑ(g,i)（ｉ＝３,５,７,……）の合計値と全部の周波数成分にわたる強度Ｑ(g,i)の合計値との強度比に相当する。周波数軸上の奇数番目の各周波数成分（ｉ＝３,５,７,……）は調波成分を優勢に包含し、偶数番目の各周波数成分（ｉ＝２,４,６,……）は非調波成分を優勢に包含する。したがって、数式(3)の調波強度指標Ｈ(g)は、候補長ＮC(g)の候補区間Ｃ(g)内の音響信号ｓ(n)の合計強度（数式(3)の右辺の分母）に対する各調波成分の合計強度（数式(3)の右辺の分子）の相対比である。以上の説明から理解される通り、調波強度指標Ｈ(g)は、音響信号ｓ(n)のうち候補長ＮC(g)の候補区間Ｃ(g)内における調波成分の優勢度（非調波成分に対する調波成分の強度の優勢度）の指標に相当する。すなわち、音響信号ｓ(n)のうち候補長ＮC(g)の候補区間Ｃ(g)内で調波成分が非調波成分と比較して優勢であるほど調波強度指標Ｈ(g)は大きい数値に設定される。以上に説明した調波強度指標Ｈ(g)がＧ個の候補長ＮC(1)〜ＮC(G)の各々について算定される。図６に示すように、調波強度指標Ｈ(g)は候補長ＮC(g)の数値（横軸）に応じて変動する。

【0025】

図３の解析長特定部３４は、Ｇ個の候補長ＮC(1)〜ＮC(G)に応じて解析区間（フレーム）の解析長Ｎ(k)を特定する。具体的には、Ｇ個の候補長ＮC(1)〜ＮC(G)のうち各候補長ＮC(g)の調波強度指標Ｈ(g)に応じて選択された候補長ＮC(g)から解析長Ｎ(k)が特定される。実施形態の解析長特定部３４は、図６に例示される通り、Ｇ個の候補長ＮC(1)〜ＮC(G)のうち調波強度指標Ｈ(g)が最大となる候補長（すなわち、調波成分の優勢度が最大となる候補長）ＮC(g)を解析長Ｎ(k)として選択する。すなわち、解析長Ｎ(k)は、以下の数式(4)で表現される。

【数4】

以上の説明から理解される通り、特定点Ｔ(k)を中心とするＧ個の候補区間Ｃ(1)〜Ｃ(G)のうち音響信号ｓ(n)の調波成分が非調波成分と比較して優勢である１個の候補区間Ｃ(g)が解析長Ｎ(k)の解析区間として特定される。図６に矢印で図示されるように、解析長特定部３４による解析長Ｎ(k)の特定は、基本周波数推定部２２が推定した基本周波数ＦAに対応する基準長Ｎ0(k)を調波成分が優勢な解析長Ｎ(k)に補正する処理とも換言され得る。

【0026】

図３の周波数補正部３６は、解析長特定部３４が特定した解析長Ｎ(k)に応じて特定点Ｔ(k)での音響信号ｓ(n)の基本周波数Ｆ0(k)を補正する。具体的には、周波数補正部３６は、特定点Ｔ(k)での音響信号ｓ(n)の基本周波数Ｆ0(k)を以下の数式(5)の演算で算定する。基本周波数Ｆ0(k)の補正は特定点Ｔ(k)毎（解析長Ｎ(k)毎）に実行される。

【数5】

【0027】

以上が解析区間設定部２６の具体例である。図１の境界特定部２８は、音響信号ｓ(n)のうち解析区間設定部２６（解析長特定部３４）が特定した解析長Ｎ(k)に対応する時間長の解析区間（フレーム）を解析することで調波境界周波数ＦC(k)を特定点Ｔ(k)毎（解析長Ｎ(k)毎）に特定する。第１実施形態の解析区間は、解析長Ｎ(k)の半分の時間長（Ｎ(k)/２）に設定される。数式(2)を参照して前述した通り、基準長Ｎ0(k)は音響信号ｓ(n)のピッチ周期の２倍に相当するから、解析区間の時間長は、音響信号ｓ(n)のピッチ周期の１個分に概略的には近似する。図７は、境界特定部２８のブロック図である。図７に示すように、境界特定部２８は、区間シフト部４２と調波次数解析部４４と周波数算定部４６とを含んで構成される。

【0028】

区間シフト部４２は、時間軸上の位置が相違するＭ個の解析区間Ｓ(1)〜Ｓ(M)を特定点Ｔ(k)毎に設定する。時間軸上の１個の特定点Ｔ(k)に対応するＭ個の解析区間Ｓ(1)〜Ｓ(M)の各々は、図８に示すように、解析長特定部３４が特定した共通の解析長Ｎ(k)にわたる区間であり、当該特定点Ｔ(k)を区間内に包含するように時間軸上の相異なる位置に設定される。具体的には、音響信号ｓ(n)のサンプルの１個分を単位として解析長Ｎ(k)に応じた解析区間を時間軸上で順次に移動（シフト）させることでＭ個の解析区間Ｓ(1)〜Ｓ(M)が設定される。図８から理解される通り、Ｍ個の解析区間Ｓ(1)〜Ｓ(M)の各々の始端（時間軸上の左端）は範囲ＡL内に分布し、Ｍ個の解析区間Ｓ(1)〜Ｓ(M)の各々の終端（時間軸上の右端）は範囲ＡR内に分布する。

【0029】

各解析区間Ｓ(m)（ｍ＝１〜Ｍ）の始端の範囲ＡLは、以下の数式(6)で表現される時点ａL1から時点ａL2までの時間範囲であり、各解析区間Ｓ(m)の終端の範囲ＡRは、以下の数式(7)で表現される時点ａR1から時点ａR2までの時間範囲である。数式(6)および数式(7)の記号τは所定の正数（例えばτ＝０.２５）に設定される。

【数6】

【数7】

【0030】

図７の調波次数解析部４４は、音響信号ｓ(n)の各解析区間Ｓ(m)を解析することで調波次数ｈC(k)を特定点Ｔ(k)毎に順次に算定する。調波次数ｈC(k)は、音響信号ｓ(n)の複数の調波成分のうち調波帯域と非調波帯域との境界に対応（近似または合致）する調波成分の次数（周波数軸上の低域側から複数の調波成分を計数した場合の各調波成分の番号）である。

【0031】

図７の周波数算定部４６は、調波次数解析部４４が算定した調波次数ｈC(k)と音響信号ｓ(n)の基本周波数Ｆ0(k)とに応じて調波境界周波数ＦC(k)を算定する。具体的には、周波数算定部４６は、調波次数解析部４４が算定した調波次数ｈC(k)と周波数補正部３６による補正後の基本周波数Ｆ0(k)とに応じて特定点Ｔ(k)毎に調波境界周波数ＦC(k)を算定する。例えば、以下の数式(8)で表現される通り、調波次数ｈC(k)と補正後の基本周波数Ｆ0(k)との乗算値が調波境界周波数ＦC(k)として算定される。

【数8】

【0032】

図９を参照して調波次数解析部４４の具体例を説明する。図９に示すように、調波次数解析部４４は、第１演算部５２と第２演算部５４とを含んで構成される。第１演算部５２は、区間シフト部４２が設定したＭ個の解析区間Ｓ(1)〜Ｓ(M)の各々について暫定調波次数ｈA(m)（ｈA(1)〜ｈA(M)）と確度系列Ｐ(m)（Ｐ(1)〜Ｐ(M)）とを算定する。暫定調波次数ｈA(m)は、音響信号ｓ(n)の解析区間Ｓ(m)内のスペクトルに観測されるＪ個（Ｊは２以上の自然数）の調波成分のうち調波帯域と非調波帯域との境界に対応する１個の調波成分の次数に相当する。確度系列Ｐ(m)は、周波数軸上の相異なる次数（以下「候補次数」という）に対応するＪ個の確度ｐ(m,1)〜ｐ(m,J)の系列（すなわちＪ次元ベクトル）である。第ｊ番目の候補次数の確度ｐ(m,j)は、周波数軸上のＪ個の候補次数のうち第ｊ番目の候補次数が暫定調波次数ｈA(m)（または調波次数ｈC(k)）に該当する可能性の指標（確率または尤度）である。すなわち、概略的には、確度系列Ｐ(m)内の相異なる候補次数に対応するＪ個の確度ｐ(m,1)〜ｐ(m,J)の最大値に対応する候補次数が暫定調波次数ｈA(m)に相当する。

【0033】

図１０は、第１演算部５２が暫定調波次数ｈA(m)および確度系列Ｐ(m)を算定する調波解析処理のフローチャートである。区間シフト部４２が設定したＭ個の解析区間Ｓ(1)〜Ｓ(M)の各々について図９の調波解析処理が実行される。調波解析処理を開始すると、第１演算部５２は、暫定調波次数ｈA(m)の候補として検討すべき候補次数の最大値（以下「最大次数」という）Ｊを設定する（ＳB1）。最大次数Ｊは、暫定調波次数ｈA(m)（または調波次数ｈC(k)）に該当する可能性がある調波成分の範囲（帯域）を規定する変数である。

【0034】

具体的には、第１演算部５２は、周波数補正部３６による補正後の基本周波数Ｆ0(k)に応じて最大次数Ｊを特定点Ｔ(k)毎に可変に設定する。例えば、第１演算部５２は、補正後の基本周波数Ｆ0(k)で所定値（以下「境界上限係数」という）ＦCmaxを除算した数値（ＦCmax/Ｆ0(k)）を整数化することで最大次数Ｊを算定する。境界上限係数ＦCmaxは、調波境界周波数ＦC(k)に想定される最大値を上回るように選定される。具体的には、調波境界周波数ＦC(k)の用途に応じて境界上限係数ＦCmaxは実験的または統計的に選定される。例えば、音響信号ｓ(n)のスペクトルのうち充分に明確な調波構造が観測される調波帯域を推定する必要がある用途では境界上限係数ＦCmaxが比較的に小さい数値に設定され、僅かでも調波構造が観測される調波帯域を推定する必要がある用途では境界上限係数ＦCmaxが比較的に大きい数値に設定される。

【0035】

第１演算部５２は、音響信号ｓ(n)のうち解析区間Ｓ(m)内の各周波数成分について正規化強度ＱN(m,i)を算定する（ＳB2）。正規化強度ＱN(m,i)は、周波数軸上の第１番目（ｉ＝１）の周波数の強度Ｑ(m,1)が所定値（以下の例示では１）となるように、音響信号ｓ(n)の解析区間Ｓ(m)のスペクトルにおける各周波数成分の強度Ｑ(m,i)を正規化した強度である。具体的には、第１演算部５２は、解析区間Ｓ(m)に対応する窓関数（典型的には解析長Ｎ(k)の窓幅の矩形窓）を音響信号ｓ(n)に乗算したうえで周波数領域のスペクトルに変換することで各周波数成分の強度Ｑ(m,i)を算定し、各強度Ｑ(m,i)を適用した以下の数式(9)の演算で正規化強度ＱN(m,i)を算定する。時間領域から周波数領域への変換には、例えば高速フーリエ変換が好適に利用される。なお、各周波数成分の振幅を強度Ｑ(m,i)として算定することも可能である。

【数9】

【0036】

第１演算部５２は、Ｊ個の候補次数の各々について累算調波強度Ｅh(m,j)（Ｅh(m,1)〜Ｅh(m,J)）と累算非調波強度Ｅa(m,j)（Ｅa(m,1)〜Ｅa(m,J)）とを算定する（ＳB3）。累算調波強度Ｅh(m,j)は、第ｊ番目の調波成分の低域側（すなわち周波数ｊ・Ｆ0(k)以下の帯域内）に位置するｊ個の調波成分の合計強度（累算エネルギー）の指標である。具体的には、以下の数式(10)で表現される通り、調波成分に相当する第ｌ番目（ｌ＝３,５,……）の周波数成分の正規化強度ＱN(m,l)と非調波成分に相当する第(ｌ−１)番目の周波数成分の正規化強度ＱN(m,l-1)との差分をｊ個の調波成分にわたり累算することで累算調波強度Ｅh(m,j)が算定される。数式(10)の記号max｛0,a｝は、数値ａを非負範囲に制限する演算を意味する。

【数10】

【0037】

他方、累算非調波強度Ｅa(m,j)は、第ｊ番目の調波成分の高域側（すなわち周波数ｊ・Ｆ0(k)以上の帯域内）に位置する非調波成分の合計強度（累算エネルギー）の指標である。具体的には、以下の数式(11)で表現される通り、非調波成分に相当する第ｌ番目（ｌ＝２j，２j＋２，……）の周波数成分の正規化強度ＱN(m,l)を累算することで累算非調波強度Ｅa(m,j)が算定される。

【数11】

【0038】

以上の演算でＪ個の候補次数の各々について累算調波強度Ｅh(m,j)と累算非調波強度Ｅa(m,j)とを算定すると、第１演算部５２は、相異なる候補次数に対応するＪ個の確度ｐ(m,j)（ｐ(m,1)〜ｐ(m,J)）を累算調波強度Ｅh(m,j)と累算非調波強度Ｅa(m,j)とに応じて算定する（ＳB4）。すなわち、Ｊ個の確度ｐ(m,1)〜ｐ(m,J)を含む確度系列Ｐ(m)が生成される。第ｊ番目の候補次数に対応する確度ｐ(m,j)は、例えば以下の数式(12)で表現される通り、累算調波強度Ｅh(m,j)と累算非調波強度Ｅa(m,j)との加重和として算定される。

【数12】

数式(12)の加重値ｂは所定の整数（例えばｂ＝２）に設定される。前述の通り、累算調波強度Ｅh(m,j)は第ｊ番目の候補次数の調波成分の低域側に位置する各調波成分の合計強度の指標であり、累算非調波強度Ｅa(m,j)は第ｊ番目の候補次数の調波成分の高域側に位置する非調波成分の合計強度の指標である。したがって、調波帯域と非調波帯域との境界（調波境界周波数ＦC(k)）に近い調波成分の候補次数に対応する確度ｐ(m,j)ほど大きい数値となる。

【0039】

以上の傾向を考慮して、第１演算部５２は、解析区間Ｓ(m)の暫定調波次数ｈA(m)を各確度ｐ(m,j)に応じて特定する（ＳB5）。具体的には、第１演算部５２は、以下の数式(13)で表現される通り、調波成分のＪ個の候補次数のうち確度ｐ(m,j)が最大となる候補次数（すなわち、調波帯域と非調波帯域との境界に最も近い調波成分の次数）を解析区間Ｓ(m)の暫定調波次数ｈA(m)として選択する。

【数13】

【0040】

以上に説明した調波解析処理を第１演算部５２が実行することで、１個の特定点Ｔ(k)に対応するＭ個の解析区間Ｓ(1)〜Ｓ(M)の各々について暫定調波次数ｈA(m)（ｈA(1)〜ｈA(M)）と確度系列Ｐ(m)（Ｐ(1)〜Ｐ(M)）とが算定される。図９の第２演算部５４は、相異なる解析区間Ｓ(m)について第１演算部５２が算定したＭ個の暫定調波次数ｈA(1)〜ｈA(M)とＭ個の確度系列Ｐ(1)〜Ｐ(M)とに応じて調波次数ｈC(k)を特定点Ｔ(k)毎に算定する。図９に示すように、第２演算部５４は、代表値算定部５４２と平滑処理部５４４とを含んで構成される。

【0041】

代表値算定部５４２は、１個の特定点Ｔ(k)に対応するＭ個の解析区間Ｓ(1)〜Ｓ(M)にわたる暫定調波次数ｈA(m)の代表値（以下「代表調波次数」という）ｈB(k)と、Ｍ個の解析区間Ｓ(1)〜Ｓ(M)にわたる確度系列Ｐ(m)の各確度ｐ(m,j)の候補次数毎の代表値（以下「代表確度」という）λ(k,j)（λ(k,1)〜λ(k,J)）とを算定する。例えば、以下の数式(14)で表現される通り、相異なる解析区間Ｓ(m)に対応するＭ個の暫定調波次数ｈA(1)〜ｈA(M)の平均値（mean）が代表調波次数ｈB(k)として特定点Ｔ(k)毎に算定される。なお、Ｍ個の暫定調波次数ｈA(1)〜ｈA(M)の平均値以外の代表値（例えば中央値）を代表調波次数ｈB(k)として算定することも可能である。

【数14】

【0042】

また、以下の数式(15)で表現される通り、相異なる解析区間Ｓ(m)に対応するＭ個の確度ｐ(1,j)〜ｐ(M,j)の中央値（median）が代表確度λ(k,j)として特定点Ｔ(k)毎に算定される。なお、Ｍ個の確度ｐ(1,j)〜ｐ(M,j)の中央値以外の代表値（例えば平均値）を代表確度λ(k,j)として算定することも可能である。

【数15】

代表確度λ(k,j)は、Ｊ個の候補次数にわたる合計値（λ(k,1)＋λ(k,2)＋……＋λ(k,J)）が１となるように正規化される。したがって、代表確度λ(k,j)は、音響信号ｓ(n)の解析長Ｎ(k)にわたる解析区間内のＪ個の調波成分のうち第ｊ番目の調波成分が、調波境界周波数ＦC(k)に近似または合致する調波成分に該当する確度の指標（確率または尤度）として利用される。

【0043】

以上に説明した通り、特定点Ｔ(k)毎に代表調波次数ｈB(k)と各候補次数の代表確度λ(k,j)（λ(k,1)〜λ(k,J)）とが算定される。平滑処理部５４４は、代表調波次数ｈB(k)の時系列と比較して変動が平滑化された調波次数ｈC(k)の時系列を算定する。具体的には、平滑処理部５４４は、代表値算定部５４２が時間軸上の特定点Ｔ(k)毎に算定した代表調波次数ｈB(k)と各候補次数の代表確度λ(k,j)とを利用して特定点Ｔ(k)毎の調波次数ｈC(k)の時系列を算定する。

【0044】

各調波次数ｈC(k)の算定（代表調波次数ｈB(k)の時系列の平滑化）には公知の動的計画法（DP：Dynamic Programming）が好適に利用される。例えば、数式(16)で表現される通り、評価関数（スコア関数）Ｚが最大となる代表調波次数ｈB(k)を時間軸上の特定点Ｔ(k)毎（解析区間毎）に選択した経路を探索するビタビアルゴリズム（constrained Viterbi algorithm）で調波次数ｈC(k)の時系列を算定することが可能である。

【0045】

【数16】

【0046】

数式(16)の評価関数Ｚは、例えば以下の数式(17)で表現される。

【数17】

数式(17)の記号ｔ_{hB(j-1),hB(j)}は、相前後する各特定点Ｔ(k)で代表調波次数ｈB(j-1)から代表調波次数ｈB(j)に変化する確度（すなわち、調波成分の個数がｈB(j-1)個からｈB(j)個に変化する遷移確率）を意味する。ただし、相前後する各特定点Ｔ(k)の間の調波成分の個数差は所定個（例えば５個）に制限される。また、数式(17)の代表確度λ(k,hB(j))は、代表調波次数ｈB(j)が調波帯域と非調波帯域との境界に該当する確度（状態確率）を意味する。以上に説明した処理で特定点Ｔ(k)毎に算定された調波次数ｈC(k)が、図７の周波数算定部４６による調波境界周波数ＦC(k)の算定（数式(8)）に適用される。

【0047】

以上に説明した形態では、数式(4)を参照して前述した通り、Ｇ個の候補長ＮC(1)〜ＮC(G)のうち調波強度指標Ｈ(g)に応じて選択された候補長ＮC(g)から解析長Ｎ(k)が特定される。したがって、例えば解析区間の時間長を基準長Ｎ0(k)（すなわち基本周波数ＦAのみに依存する時間長）に固定する構成（以下「構成例１」という）と比較して、音響信号ｓ(n)の基本周波数ＦAの推定誤差に対して頑健に調波境界周波数ＦC(k)を特定できるという利点がある。実際の実験結果を参照しながら以上の効果を説明する。

【0048】

図１１の部分(A)は、構成例１で算定される調波次数ｈC(k)の時間変化であり、図１１の部分(B)は、本実施形態で算定される調波次数ｈC(k)の時間変化である。図１１の部分(A)および部分(B)では、調波境界周波数ＦC(k)（調波次数ｈC(k)）に対する基本周波数ＦAの推定誤差の影響を評価するために、音響信号ｓ(n)から推定された基本周波数ＦAに対して相異なる程度（０.５％，１％，２％，５％，１０％）の誤差を人為的に付与した複数の場合について調波次数ｈC(k)の時間変化が併記されている。

【0049】

図１１の部分(A)から把握される通り、構成例１では、基本周波数ＦAの誤差の程度に応じて調波次数ｈC(k)が大幅に変動するという傾向がある。他方、図１１の部分(B)から把握される通り、本実施形態では、基本周波数ＦAに相異なる誤差が付与された各場合の調波次数ｈC(k)の時系列は実質的に合致する。すなわち、基本周波数ＦAの推定誤差に対して頑健に調波境界周波数ＦC(k)を特定できるという前述の効果が図１１から確認できる。

【0050】

また、本実施形態では、１個の特定点Ｔ(k)を包含するように時間軸上の位置を相違させたＭ個の解析区間Ｓ(1)〜Ｓ(M)の各々について算定された暫定調波次数ｈA(m)（ｈA(1)〜ｈA(M)）を利用して調波次数ｈC(k)（調波境界周波数ＦC(k)）が算定される。したがって、解析長Ｎ(k)の１個の解析区間のみについて調波次数を算定する構成（以下「構成例２」という）と比較して、音響信号ｓ(n)に対する解析区間の位置の変動（特定点Ｔ(k)の位置の変動）に対して頑健に調波境界周波数ＦC(k)を特定できるという利点がある。実際の実験結果を参照しながら以上の効果を説明する。

【0051】

図１２の部分(A)は、構成例２で算定される調波次数ｈC(k)の時間変化であり、図１１の部分(B)は、本実施形態で算定される調波次数ｈC(k)の時間変化である。図１１の部分(A)および部分(B)では、調波境界周波数ＦC(k)（調波次数ｈC(k)）に対する解析区間の位置の誤差の影響を評価する観点から、特定点Ｔ(k)を基準位置（誤差０％）として相異なる程度（０.５％，１％，２％，５％，１０％）の誤差を解析区間の時間軸上の位置に対して人為的に付与した複数の場合における調波次数ｈC(k)の時間変化が併記されている。

【0052】

図１２の部分(A)から把握される通り、構成例２では、解析区間（特定点Ｔ(k)）の時間軸上の位置に応じて調波次数ｈC(k)が大幅に変動するという傾向がある。すなわち、調波次数ｈC(k)の推定精度に対する特定点Ｔ(k)（例えば発声者の声門が開口または閉塞する時点）の推定誤差の影響が顕著である。他方、図１２の部分(B)から把握される通り、本実施形態では、解析区間の位置に相異なる誤差が付与された各場合の調波次数ｈC(k)の時系列が実質的に合致する。すなわち、音響信号ｓ(n)の特定点Ｔ(k)の推定誤差に対して頑健に調波境界周波数ＦC(k)を特定できるという利点がある。

【0053】

また、本実施形態では、解析長特定部３４が特定した解析長Ｎ(k)に応じて音響信号ｓ(n)の基本周波数Ｆ0(k)が補正される。したがって、数式(1)で算定された補正前の基本周波数Ｆ0(k)を数式(8)に適用する構成と比較して、調波境界周波数ＦC(k)を高精度に特定できるという利点がある。ただし、数式(1)で算定された基本周波数Ｆ0(k)を数式(8)に適用して調波境界周波数ＦC(k)を算定することも可能である。したがって、周波数補正部３６は省略され得る。

【0054】

＜変形例＞
以上の形態は多様に変形され得る。具体的な変形の態様を以下に例示する。以下の例示から任意に選択された２以上の態様は適宜に併合され得る。

【0055】

（１）境界特定部２８が調波境界周波数ＦC(k)（調波次数ｈC(k)）を算定する方法は適宜に変更される。例えば、前述の実施形態では、時間軸上の位置が相違するＭ個の解析区間Ｓ(1)〜Ｓ(M)の各々について算定された暫定調波次数ｈA(m)を利用して調波次数ｈC(k)（調波境界周波数ＦC(k)）を算定したが、解析長Ｎ(k)の１個の解析区間について図１０の調波解析処理を実行することで調波次数ｈC(k)（前述の形態における暫定調波次数ｈA(m)に相当する）を算定することも可能である。以上の説明から理解される通り、図１２の説明で想定した構成例２は本発明の範囲に包含される。

【0056】

また、平滑処理部５４４による平滑処理の具体的な方法は前述の例示（ビタビアルゴリズム）に限定されない。例えば、時間軸上の特定点Ｔ(k)毎の代表調波次数ｈB(k)の時系列を平滑化する（例えば移動平均を算定する）ことで平滑化後の調波次数ｈC(k)の時系列を算定することも可能である。

【0057】

（２）前述の形態では、Ｇ個の候補長ＮC(1)〜ＮC(G)のうち調波強度指標Ｈ(g)に応じて選択された候補長ＮC(g)から解析長Ｎ(k)が特定する構成（以下「特徴Ａ」という）と、時間軸上の位置が相違するＭ個の解析区間Ｓ(1)〜Ｓ(M)の各々の暫定調波次数ｈA(m)（ｈA(1)〜ｈA(M)）を利用して調波次数ｈC(k)（調波境界周波数ＦC(k)）を算定する構成（以下「特徴Ｂ」という）との双方を具備する音響解析装置１００を例示したが、特徴Ａおよび特徴Ｂの一方にとって他方は必須の要件ではなく各々が独立に成立し得る。例えば特徴Ｂを具備する構成において解析長Ｎ(k)を特定する方法は任意であり、Ｇ個の候補長ＮC(1)〜ＮC(G)から解析長Ｎ(k)を選択する処理は省略され得る。

【0058】

（３）前述の形態では、Ｇ個の候補長ＮC(1)〜ＮC(G)のうち調波強度指標Ｈ(g)が最大となる１個の候補長ＮC(g)を解析区間を規定する解析長Ｎ(k)として選択したが、各候補長ＮC(g)の調波強度指標Ｈ(g)を利用して解析長Ｎ(k)を特定する方法は以上の例示に限定されない。例えば、調波強度指標Ｈ(g)の降順で選択された所定個の候補長ＮC(g)の代表値（例えば平均値）や加重和を解析長Ｎ(k)として算定する構成や、調波強度指標Ｈ(g)の降順で上位に位置する所定個の候補長ＮC(g)を除外したうえで調波強度指標Ｈ(g)が最大となる候補長ＮC(g)を解析長Ｎ(k)として選択する構成も採用され得る。以上の説明から理解される通り、解析長特定部３４は、Ｇ個の候補長ＮC(1)〜ＮC(G)のうち各候補長ＮC(g)の調波強度指標Ｈ(g)に応じて選択した候補長ＮC(g)から解析長Ｎ(k)を特定する要素として包括的に表現される。

【0059】

（４）携帯電話機等の端末装置と通信するサーバ装置で音響解析装置１００を実現することも可能である。例えば、音響解析装置１００は、端末装置から受信した音響信号ｓ(n)から調波境界周波数ＦC(k)を算定して端末装置に送信する。なお、音響信号ｓ(n)の基本周波数ＦAを音響信号ｓ(n)とともに端末装置から受信する構成（例えば端末装置が基本周波数推定部２２を具備する構成）では音響解析装置１００から基本周波数推定部２２が省略され、音響信号ｓ(n)の特定点Ｔ(k)を端末装置から受信する構成（例えば端末装置が特定点設定部２４を具備する構成）では音響解析装置１００から特定点設定部２４が省略される。

【符号の説明】

【0060】

１００……音響解析装置、１２……演算処理装置、１４……記憶装置、２２……基本周波数推定部、２４……特定点設定部、２６……解析区間設定部、２８……境界特定部、３２……指標算定部、３４……解析長特定部、３６……周波数補正部、４２……区間シフト部、４４……調波次数解析部、４６……周波数算定部、５２……第１演算部、５４……第２演算部、５４２……代表値算定部、５４４……平滑処理部。

【図1】

【図2】

【図3】

【図4】

【図5】

【図6】

【図7】

【図8】

【図9】

【図10】

【図11】

【図12】