特開 2024-37379 周波数帯毎に境界を計算する波動計算プログラム、装置及び方法

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

(12)【公報種別】公開特許公報(A)

(11)【公開番号】P2024037379

(43)【公開日】2024-03-19

(54)【発明の名称】周波数帯毎に境界を計算する波動計算プログラム、装置及び方法

(51)【国際特許分類】

   G01H 17/00 20060101AFI20240312BHJP

   G10K 15/00 20060101ALN20240312BHJP

   G06F 30/23 20200101ALN20240312BHJP

【ＦＩ】

G01H17/00 Z

G10K15/00 M

G06F30/23

【審査請求】未請求

【請求項の数】11

【出願形態】ＯＬ

(21)【出願番号】P 2022142194

(22)【出願日】2022-09-07

(71)【出願人】

【識別番号】000208891

【氏名又は名称】ＫＤＤＩ株式会社

(74)【代理人】

【識別番号】100135068

【弁理士】

【氏名又は名称】早原 茂樹

(74)【代理人】

【識別番号】100141313

【弁理士】

【氏名又は名称】辰巳 富彦

(72)【発明者】

【氏名】大久保 翔太

(72)【発明者】

【氏名】堀内 俊治

【テーマコード（参考）】

2G064

5B146

【Ｆターム（参考）】

2G064AB01

2G064AB02

2G064AB13

2G064CC43

2G064DD02

5B146DJ01

5B146DJ07

5B146DJ11

(57)【要約】

【課題】周波数帯によって生じる、計算値の実測値からの乖離を抑制可能な波動計算プログラムを提供する。
【解決手段】本プログラムは、所定の境界を有する空間における波動に係る量である波動関連量を、有限差分時間領域法を用いて計算する波動計算プログラムであり、波動に係る周波数範囲に含まれる複数の周波数帯における周波数帯毎に、当該周波数帯について設定された境界計算式及びこの波動に係る離散化式を用いて、この空間における当該周波数帯での波動の計算に係る計算対象量を算出する周波数帯毎計算手段と、周波数帯毎に算出されたこの空間における当該周波数帯での計算対象量から、この空間における波動関連量を算出する波動関連量算出手段としてコンピュータを機能させる。ここで、波動は音波であって、計算対象量は音圧に係る量であり、波動関連量は、設定された音源位置と受音位置との間におけるインパルス応答とすることもできる。
【選択図】図１

【特許請求の範囲】

【請求項1】

所定の境界を有する空間における波動に係る量である波動関連量を、有限差分時間領域法を用いて計算するプログラムであって、
当該波動に係る周波数範囲に含まれる複数の周波数帯における周波数帯毎に、当該周波数帯について設定された境界計算式と、当該波動に係る離散化式とを用いて、当該空間における当該周波数帯での当該波動の計算に係る計算対象量を算出する周波数帯毎計算手段と、
周波数帯毎に算出された当該空間における当該周波数帯での当該計算対象量から、当該空間における当該波動関連量を算出する波動関連量算出手段と
してコンピュータを機能させることを特徴とする波動計算プログラム。

【請求項2】

前記周波数帯毎計算手段は、周波数帯毎に、当該境界の吸収又はインピーダンスに係る量である境界量をパラメータとして含む当該境界計算式に対して、当該周波数帯について設定された当該境界量の値を適用し、当該境界量の値を適用した当該周波数帯での当該境界計算式と、当該離散化式とを用いて、当該境界に隣接する当該境界外の位置における当該周波数帯での当該計算対象量と、当該境界外の位置での当該計算対象量を用いて計算される当該境界での当該計算対象量を含む、当該空間における当該周波数帯での当該計算対象量と、を算出することを特徴とする請求項１に記載の波動計算プログラム。

【請求項3】

前記周波数帯毎計算手段は、周波数帯毎に、当該周波数帯に合わせて設定されたサンプリング周波数に対応した当該境界計算式及び当該離散化式を用いて、当該空間における当該周波数帯での当該計算対象量を算出することを特徴とする請求項１又は２に記載の波動計算プログラム。

【請求項4】

前記波動関連量算出手段は、周波数帯毎に、当該周波数帯での当該計算対象量から当該周波数帯での当該波動関連量を算出し、算出した各周波数帯での当該波動関連量を合算して、又は所定の周波数帯に係る当該波動関連量に対しより大きい若しくは小さい重みを付与した上で合算して、当該空間における当該波動関連量を算出することを特徴とする請求項１又は２に記載の波動計算プログラム。

【請求項5】

前記波動関連量算出手段は、周波数帯毎に、当該周波数帯での当該計算対象量から当該周波数帯での当該波動関連量を算出し、算出した各周波数帯での当該波動関連量に対し、当該周波数帯に合わせて通過帯域が設定されたバンドパスフィルタ処理を施し、当該バンドパスフィルタ処理を施した各周波数帯での当該波動関連量を合算して、又は所定の周波数帯に係る当該波動関連量に対しより大きい若しくは小さい重みを付与した上で合算して、当該空間における当該波動関連量を算出することを特徴とする請求項１又は２に記載の波動計算プログラム。

【請求項6】

前記波動関連量算出手段は、周波数帯毎に、当該周波数帯での当該計算対象量から当該周波数帯での当該波動関連量を算出し、算出した各周波数帯での当該波動関連量に対しフーリエ変換処理を施し、当該フーリエ変換処理を施した各波数帯での当該波動関連量から、当該空間におけるフーリエ変換された当該波動関連量を決定し、決定した当該フーリエ変換された当該波動関連量に対し逆フーリエ変換処理を施して、当該空間における当該波動関連量を算出することを特徴とする請求項１又は２に記載の波動計算プログラム。

【請求項7】

当該有限差分時間領域法はＣＥ-ＦＤＴＤ（Compact-Explicit Finite-Difference Time-Domain）法であって、当該境界計算式は、複数の参照方向における参照方向毎に設定されており、当該離散化式は、当該複数の参照方向に対応した式となっていることを特徴とする請求項１又は２に記載の波動計算プログラム。

【請求項8】

当該有限差分時間領域法は、支配方程式における目的の変数にベクトル量を含む方式、又は支配方程式における目的の変数がスカラ量である方式であることを特徴とする請求項１又は２に記載の波動計算プログラム。

【請求項9】

当該波動は音波であって、当該計算対象量は音圧に係る量であり、当該波動関連量は、設定された音源位置と受音位置との間におけるインパルス応答であることを特徴とする請求項１又は２に記載の波動計算プログラム。

【請求項10】

所定の境界を有する空間における波動に係る量である波動関連量を、有限差分時間領域法を用いて計算する装置であって、
当該波動に係る周波数範囲に含まれる複数の周波数帯における周波数帯毎に、当該周波数帯について設定された境界計算式と、当該波動に係る離散化式とを用いて、当該空間における当該周波数帯での当該波動の計算に係る計算対象量を算出する周波数帯毎計算手段と、
周波数帯毎に算出された当該空間における当該周波数帯での当該計算対象量から、当該空間における当該波動関連量を算出する波動関連量算出手段と
を有することを特徴とする波動計算装置。

【請求項11】

所定の境界を有する空間における波動に係る量である波動関連量を、有限差分時間領域法を用いて計算する方法であって、
当該波動に係る周波数範囲に含まれる複数の周波数帯における周波数帯毎に、当該周波数帯について設定された境界計算式と、当該波動に係る離散化式とを用いて、当該空間における当該周波数帯での当該波動の計算に係る計算対象量を算出するステップと、
周波数帯毎に算出された当該空間における当該周波数帯での当該計算対象量から、当該空間における当該波動関連量を算出するステップと
を有することを特徴とする、コンピュータによって実施される波動計算方法。

【発明の詳細な説明】

【技術分野】

【0001】

本発明は、空間を伝播する波動についての数値シミュレーション技術に関する。

【背景技術】

【0002】

コンピュータによる数値解析としての数値シミュレーションによって、空間を伝播する波動に係る量である波動関連量を計算する技術が、数多く開発されている。例えば、音波の伝播する音場のインパルス応答は、空間音響特性を表す重要な指標であるが、このインパルス応答を数値シミュレーションによって計算する技術の開発が盛んに進んでいる。

【0003】

このような音場の数値シミュレーション技術の例として、特許文献１や非特許文献１には、音線法や虚像法を用いてインパルス応答を計算する手法が開示されている。これらの手法では、音波の波動性が考慮されないため、特に低周波数帯において、計算値と実測値との乖離が顕著となってしまう。

【0004】

これに対し、非特許文献２や非特許文献３には、有限差分時間領域（ＦＤＴＤ, Finite Difference Time Domain）法を用いてインパルス応答を計算する技術が開示されている。ここでＦＤＴＤ法は、例えば特許文献２に開示されているように、電磁場解析で盛んに用いられてきた数値シミュレーション手法である。非特許文献２や非特許文献３では、通常のＦＤＴＤ法と比較してより多方向を参照した離散化式を適用するＣＥ（Compact explicit）-ＦＤＴＤ法が用いられている。これにより、計算結果の精度が向上し、例えばサンプリング周波数までのインパルス応答の計算も可能になるのである。

【先行技術文献】

【特許文献】

【0005】

【特許文献1】特開平９－２９２８８５号公報

【特許文献2】特開２００６－１３９７２３号公報

【非特許文献】

【0006】

【非特許文献1】石田康二, 「虚音源によるホール音響の可視化」, 騒音制御, 22巻, 1号, pp.31-33, １９９８年

【非特許文献2】土屋隆生, 「音場シミュレーションと音空間レンダリング」, 電子情報通信学会 基礎・境界ソサイエティ Fundamentals Review, 10巻, 3号, pp.206-218, ２０１７年

【非特許文献3】Konrad Kowalczyk and Maarten van Walstijn, “Room Acoustics Simulation Using 3-D Compact Explicit FDTD Schemes”, IEEE Transactions on Audio, Speech, and Language Processing, Vol.19, Issue.1, pp.34-36, ２０１１年

【発明の概要】

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

【0007】

しかしながら、上記の非特許文献２や非特許文献３に開示されたような従来技術では、ＣＥ-ＦＤＴＤ法を用いて精度向上を図っているにもかかわらず、依然、周波数帯によっては再現精度（実測値との一致度）が低下してしまう問題が生じていた。

【0008】

実際、本願発明者等は、ある防音室についてＣＥ-ＦＤＴＤ法を用いてインパルス応答の残響時間を計算し、この計算値が、境界条件の設定に用いた周波数帯よりも低い周波数帯及び高い周波数帯において実測値と乖離し、特に低い周波数帯でより大きく乖離してしまうことを確認している。

【0009】

ここで、このような実測値との乖離は、通常のＦＤＴＤ法を用いた音波関連量の計算結果や、さらには電磁波等、音波以外の波動に係るＦＤＴＤ法の計算結果においても、しばしば認められるものとなっている。

【0010】

そこで、本発明は、周波数帯によって生じる、計算値の実測値からの乖離を抑制可能な波動計算プログラム、装置及び方法を提供することを目的とする。

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

【0011】

本発明によれば、所定の境界を有する空間における波動に係る量である波動関連量を、有限差分時間領域法を用いて計算するプログラムであって、
当該波動に係る周波数範囲に含まれる複数の周波数帯における周波数帯毎に、当該周波数帯について設定された境界計算式と、当該波動に係る離散化式とを用いて、当該空間における当該周波数帯での当該波動の計算に係る計算対象量を算出する周波数帯毎計算手段と、
周波数帯毎に算出された当該空間における当該周波数帯での当該計算対象量から、当該空間における当該波動関連量を算出する波動関連量算出手段と
してコンピュータを機能させる波動計算プログラムが提供される。

【0012】

この本発明による波動計算プログラムにおいて、周波数帯毎計算手段は具体的に、周波数帯毎に、当該境界の吸収又はインピーダンスに係る量である境界量をパラメータとして含む当該境界計算式に対して、当該周波数帯について設定された当該境界量の値を適用し、当該境界量の値を適用した当該周波数帯での当該境界計算式と、当該離散化式とを用いて、当該境界に隣接する当該境界外の位置における当該周波数帯での当該計算対象量と、当該境界外の位置での当該計算対象量を用いて計算される当該境界での当該計算対象量を含む、当該空間における当該周波数帯での当該計算対象量と、を算出することも好ましい。

【0013】

また、本発明による波動計算プログラムの一実施形態として、周波数帯毎計算手段は、周波数帯毎に、当該周波数帯に合わせて設定されたサンプリング周波数に対応した当該境界計算式及び当該離散化式を用いて、当該空間における当該周波数帯での当該計算対象量を算出することも好ましい。

【0014】

さらに、本発明による波動計算プログラムの他の実施形態として、波動関連量算出手段は、周波数帯毎に、当該周波数帯での当該計算対象量から当該周波数帯での当該波動関連量を算出し、算出した各周波数帯での当該波動関連量を合算して、又は所定の周波数帯に係る当該波動関連量に対しより大きい若しくは小さい重みを付与した上で合算して、当該空間における当該波動関連量を算出することも好ましい。

【0015】

また、本発明による波動計算プログラムの更なる他の実施形態として、波動関連量算出手段は、周波数帯毎に、当該周波数帯での当該計算対象量から当該周波数帯での当該波動関連量を算出し、算出した各周波数帯での当該波動関連量に対し、当該周波数帯に合わせて通過帯域が設定されたバンドパスフィルタ処理を施し、当該バンドパスフィルタ処理を施した各周波数帯での当該波動関連量を合算して、又は所定の周波数帯に係る当該波動関連量に対しより大きい若しくは小さい重みを付与した上で合算して、当該空間における当該波動関連量を算出することも好ましい。

【0016】

さらに、本発明による波動計算プログラムの更なる他の実施形態として、波動関連量算出手段は、周波数帯毎に、当該周波数帯での当該計算対象量から当該周波数帯での当該波動関連量を算出し、算出した各周波数帯での当該波動関連量に対しフーリエ変換処理を施し、当該フーリエ変換処理を施した各波数帯での当該波動関連量から、当該空間におけるフーリエ変換された当該波動関連量を決定し、決定した当該フーリエ変換された当該波動関連量に対し逆フーリエ変換処理を施して、当該空間における当該波動関連量を算出することも好ましい。

【0017】

また、本発明による波動計算プログラムにおける当該有限差分時間領域法は、ＣＥ-ＦＤＴＤ（Compact-Explicit Finite-Difference Time-Domain）法であって、当該境界計算式は、複数の参照方向における参照方向毎に設定されており、当該離散化式は、当該複数の参照方向に対応した式となっていることも好ましい。

【0018】

さらに、本発明による波動計算プログラムにおける当該有限差分時間領域法は、支配方程式における目的の変数にベクトル量を含む方式、又は支配方程式における目的の変数がスカラ量である方式であることも好ましい。

【0019】

また、本発明による波動計算プログラムにおいて、当該波動は音波であって、当該計算対象量は音圧に係る量であり、当該波動関連量は、設定された音源位置と受音位置との間におけるインパルス応答であることも好ましい。

【0020】

本発明によれば、また、所定の境界を有する空間における波動に係る量である波動関連量を、有限差分時間領域法を用いて計算する装置であって、
当該波動に係る周波数範囲に含まれる複数の周波数帯における周波数帯毎に、当該周波数帯について設定された境界計算式と、当該波動に係る離散化式とを用いて、当該空間における当該周波数帯での当該波動の計算に係る計算対象量を算出する周波数帯毎計算手段と、
周波数帯毎に算出された当該空間における当該周波数帯での当該計算対象量から、当該空間における当該波動関連量を算出する波動関連量算出手段と
を有する波動計算装置が提供される。

【0021】

本発明によれば、さらに、所定の境界を有する空間における波動に係る量である波動関連量を、有限差分時間領域法を用いて計算する方法であって、
当該波動に係る周波数範囲に含まれる複数の周波数帯における周波数帯毎に、当該周波数帯について設定された境界計算式と、当該波動に係る離散化式とを用いて、当該空間における当該周波数帯での当該波動の計算に係る計算対象量を算出するステップと、
周波数帯毎に算出された当該空間における当該周波数帯での当該計算対象量から、当該空間における当該波動関連量を算出するステップと
を有することを特徴とする、コンピュータによって実施される波動計算方法が提供される。

【発明の効果】

【0022】

本発明の波動計算プログラム、装置及び方法によれば、周波数帯によって生じる、計算値の実測値からの乖離を抑制することができる。

【図面の簡単な説明】

【0023】

【図1】本発明による波動計算装置の一実施形態における機能構成を示す機能ブロック図である。

【図2】ＣＥ-ＦＤＴＤ法における参照方向を説明するための模式図である。

【図3】ＣＥ-ＦＤＴＤ法における境界計算式を説明するための模式図である。

【図4】ベクトル変数型ＦＤＴＤ法におけるスタッガードグリッド及びコロケードグリッドを示した模式図である。

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

【0024】

以下、本発明の実施形態について、図面を用いて詳細に説明する。

【0025】

［波動計算装置］
図１は、本発明による波動計算装置の一実施形態における機能構成を示す機能ブロック図である。

【0026】

図１に示した、本発明の一実施形態としての波動計算装置１は、所定の境界を有する空間における波動に係る量である「波動関連量」を、有限差分時間領域（ＦＤＴＤ）法を用いて計算する装置である。ここで本実施形態の波動計算装置１は、波動として音波を採用して、波動計算に係る「計算対象量」としての音圧を計算し、さらにこの音圧の計算結果を用いて、この音波が伝播する音場における「波動関連量」としてのインパルス応答を算出することが可能となっている。

【0027】

このような波動計算処理を実施すべく、波動計算装置１は具体的に、
（Ａ）波動（本実施形態では音波）に係る周波数範囲に含まれる複数の周波数帯における周波数帯毎に、当該周波数帯について設定された「境界計算式」と、この波動（音波）に係る「離散化式」とを用いて、所定の境界を有する空間における当該周波数帯での「計算対象量」（本実施形態では音圧）を算出する周波数帯毎（計算対象量）計算部１１１と、
（Ｂ）周波数帯毎に算出されたこの空間における当該周波数帯での「計算対象量」（音圧）から、この空間における「波動関連量」（本実施形態ではインパルス応答）を算出する波動関連量算出部１１２と
を有している。

【0028】

ここで従来、ＦＤＴＤ法を用いて導出された波動関連量の計算値は、境界条件の設定に用いた周波数帯よりも低い周波数帯及び高い周波数帯において実測値と乖離してしまう傾向にあった。これに対し波動計算装置１は、上記の構成（Ａ）及び（Ｂ）から分かるように、周波数帯毎に設定された「境界計算式」を用いて、各周波数帯での「計算対象量」を算出した上で、目的である所定の境界を有する空間における「波動関連量」を算出している。その結果、周波数帯によって生じる、計算値の実測値からの乖離を抑制し、波動計算処理の精度向上を図ることが可能となるのである。

【0029】

なお本実施形態において、上述したように波動は音波であって、計算対象量及び波動関連量はそれぞれ、音圧（又は音圧に係る量）及びインパルス応答となっているが勿論、本発明に係る波動、計算対象量及び波動関連量はこれらに限定されるものではない。例えば、音圧と粒子速度を計算対象量として、波動関連量としての音響インテンシティを算出することもできる。さらに波動として電磁波を採用し、電磁場強度（電場強度や磁場強度等）を計算対象量として、波動関連量としての電磁場勾配や電磁場エネルギー（ポインティングベクトル）を算出することも可能である。

【0030】

いずれにしても本発明による波動計算処理は、その他種々様々な波動、計算対象量及び波動関連量を取り扱うことができるが、特に、境界条件が周波数によって相当に変化する空間を伝播する波動に係る波動関連量を、好適に計算することが可能となっている。

【0031】

また、波動の伝播する空間は勿論、３次元に限定されるものではなく、例えば２次元空間（平面）を伝播する波動のＦＤＴＤ法も使用可能である。さらに、本実施形態において使用されるＦＤＴＤ法は後述するように、通常のＦＤＴＤ法と比較してより多方向を参照した離散化式を適用するＣＥ（Compact explicit）-ＦＤＴＤ法となっている。しかしながら勿論、本発明による波動計算処理は、通常の（多方向参照ではない）ＦＤＴＤ法を用いて実施することもできる。

【0032】

さらにまた、本発明による波動計算処理においては、
（ａ）「支配方程式における目的の変数にベクトル量を含む方式」（以下、ベクトル変数型と略称）、及び
（ｂ）「支配方程式における目的の変数がスカラ量である方式」（以下、スカラ変数型と略称）
のいずれのＦＤＴＤ法を採用することも可能である。ちなみにＣＥ-ＦＤＴＤ法は、スカラ変数型ＦＤＴＤ法の多方向参照による改良版となっている。

【0033】

ここで以下、上記（Ａ）の「境界計算式」及び「離散化式」を含め、これらのＦＤＴＤ法について、図面を用いて説明を行う。ちなみに、ＣＥ-ＦＤＴＤ法を含めＦＤＴＤ法については、非特許文献２や、非特許文献：Osamu Yamashita, Takao Tsuchiya, Yukio Iwaya, Makoto Otani, and Yasushi Inoguchi, “Reflective boundary condition with arbitrary boundary shape for compact-explicit finite-difference time-domain method”, Japanese Journal of Applied Physics, Volume 54, Number 7S1, ２０１５年において詳細に解説されている。

【0034】

＜ＦＤＴＤ法の説明＞
最初に、波動を音波、計算対象量を音圧とした場合における、ＣＥ-ＦＤＴＤ法（スカラ変数型）を用いた、所定の境界を有するある媒質（例えば空気）の３次元空間における音圧（の分布）のシミュレーションを説明する。

【0035】

この場合、支配方程式として、次式
（１） c^-2×∂²p/∂t²＝∇²p
を用いる。ここで、pは音圧であり、cは（その媒質での）音速である。ちなみに上式（１）は、状態方程式p＝c²×ρ（ρは媒質密度）を用いて、後述する音波（音場）の波動方程式及び連続方程式から導出される、音圧の波動方程式となっている。

【0036】

次いで、上式（１）を離散化することによって、音圧についてのＣＥ-ＦＤＴＤ法の離散化式
（２） p_i,j,k ⁿ⁺¹＝
d₁×(p_i+1,j,k ⁿ＋p_i-1,j,k ⁿ＋p_i,j+1,k ⁿ＋p_i,j-1,k ⁿ＋p_i,j,k+1 ⁿ＋p_i,j,k-1 ⁿ)
＋d₂×(p_i+1,j+1,k ⁿ＋p_i+1,j-1,k ⁿ＋p_i+1,j,k+1 ⁿ＋p_i+1,j,k-1 ⁿ＋p_i,j+1,k+1 ⁿ
＋p_i,j+1,k-1 ⁿ＋p_i,j-1,k+1 ⁿ＋p_i,j-1,k-1 ⁿ＋p_i-1,j+1,k ⁿ＋p_i-1,j-1,k ⁿ
＋p_i-1,j,k+1 ⁿ＋p_i-1,j,k-1 ⁿ)
＋d₃×(p_i+1,j+1,k+1 ⁿ＋p_i+1,j-1,k+1 ⁿ＋p_i+1,j+1,k-1 ⁿ＋p_i+1,j-1,k-1 ⁿ
＋p_i-1,j+1,k+1 ⁿ＋p_i-1,j-1,k+1 ⁿ＋p_i-1,j+1,k-1 ⁿ＋p_i-1,j-1,k-1 ⁿ)
＋d₄×p_i,j,k ⁿ－p_i,j,k ^n-1
ここで、d₁＝χ²×(1－4a＋4b)
d₂＝χ²×(a－2b)
d₃＝χ²×b
d₄＝2×(1－3×χ²＋6a×χ²－4b×χ²)
が得られる。

【0037】

ここで上式（２）において、p_i,j,k ⁿは、グリッド間隔をΔ、時間ステップをΔtとして、空間内の３次元グリッド位置(i×Δ,j×Δ,k×Δ)における時点n×Δtでの音圧となっている。また、χは、次式
（３） χ＝c×Δt／Δ
によって定義される、安定条件を表すＣＦＬ（Courant-Friedrichs-Lewy）数である。さらに、a及びbは精度を制御するためのパラメータであり、手法に応じて適切な値に設定される。例えばa＝1／4、b＝1／16とすると、カットオフ周波数がナイキスト周波数に一致し精度が向上する。

【0038】

また上式（２）において、
・d₁の項は、グリッド空間における座標軸方向の参照（以下、ＳＬＦ(Standard Leapflog)と略称）に対応しており、
・d₂の項は、グリッド空間における（面内）対角線方向の参照（以下、ＣＣＰ(Cubic Close packed)と略称）に対応しており、
・d₃の項は、グリッド空間における立体対角線方向の参照（以下、ＯＣＴＡ(OCTAhedral)と略称）に対応している。

【0039】

ここで図２に、ＳＬＦ（図２（Ａ））、ＣＣＰ（図２（Ｂ））及びＯＣＴＡ（図２（Ｃ））における、グリッド空間内での参照方向を示している。このように、ＣＥ-ＦＤＴＤ法の離散化式（上式（２））は、設定された複数の参照方向に対応した式となっていることが分かる。

【0040】

ちなみにａ及びbをゼロ（a＝0，b＝0）とした場合、d₂もd₃もゼロとなり、その結果、上式（２）は、ＳＬＦに対応したd₁の項を含む、通常の（参照方向が１種類の）スカラ変数型ＦＤＴＤ法の離散化式となる。

【0041】

次に、この３次元空間の境界における音圧の計算方法について説明する。図３は、ＣＥ-ＦＤＴＤ法における境界計算式を説明するための模式図である。

【0042】

図３に示したように、音圧を計算すべきグリッド空間（図３では実線で示した空間）における境界位置（例えば図３の黒色印の位置）での音圧値を算出するには、この境界に隣接する境界外の仮想グリッド位置（例えば図３の白色印の位置）での（推定値としての）音圧値も必要となる。そこで、これらの仮想グリッド位置での音圧も含めた境界計算式を、インピーダンス境界条件を適用して認定し、この境界計算式を用いて境界での音圧を計算するのである。

【0043】

具体的に、ＳＬＦ（図３では白丸印の参照）に対応した境界計算式の１つとして、次式
（４） p_i+1,j,k ⁿ⁺¹＝p_i,j,k ⁿ＋A×(p_i,j,k ⁿ⁺¹－p_i+1,j,k ⁿ)
ここで、A＝(Z^*×χ－1)／(Z^*×χ＋1)
Z^*＝Z／ρc
を用いることができる。ここで上式（４）は（図３の白丸印のうちの１つの位置での）音圧p_i+1,j,k ⁿ⁺¹についての計算式となっているが、他のＳＬＦ相当位置での音圧についても同様の境界計算式を規定することができる。ちなみに、通常の（参照方向が１種類の）スカラ変数型ＦＤＴＤ法の境界計算式は、このようなＳＬＦに対応したもののみで構成されることになる。

【0044】

なお上式（４）において、Z^*は境界における規格化インピーダンスであって、Zは境界における音響インピーダンスとなっている。この音響インピーダンスZと、境界の垂直入射吸音率αや同反射率rとの間では次式
（５） r＝((Z－Z₀)／(Z＋Z₀))²
（６） α＝1－r
の関係が成立する。ここでZ₀は、媒質（例えば空気）の音響インピーダンスである。したがって、境界条件として境界の音響インピーダンスZ、又は吸音率α若しくは反射率r、を指定してやれば上記のパラメータAが決まり、これにより境界計算式（４）が使用可能となるのである。

【0045】

さらに、ＣＣＰ（図３では白四角印の参照）に対応した境界計算式の１つとして、次式
（７） p_i+1,j+1,k ⁿ⁺¹＝p_i,j,k ⁿ
＋(2^0.5×Z^*χ－1)^-1×(p_i+1,j,k ⁿ＋p_i,j+1,k ⁿ－p_i+1,j,k ⁿ⁺¹－p_i,j+1,k ⁿ⁺¹)
＋(2^0.5×Z^*χ－1)×(2^0.5×Z^*χ－1)^-1×(p_i,j,k ⁿ⁺¹－p_i+1,j+1,k ⁿ)
を用いることができる。ここで上式（７）は（図３の白四角印のうちの１つの位置での）音圧p_i+1,j+1,k ⁿ⁺¹についての計算式となっているが、他のＣＣＰ相当位置での音圧についても同様の境界計算式を規定することができる。

【0046】

さらにまた、ＣＯＴＡ（図３では白二重丸印の参照）に対応した境界計算式の１つとして、次式
（８） p_i+1,j+1,k+1 ⁿ⁺¹＝p_i,j,k ⁿ
＋(3^0.5－Z^*χ)×(3^0.5＋3Z^*χ)^-1×(p_i,j+1,k+1 ⁿ＋p_i+1,j+1,k ⁿ
＋p_i+1,j,k+1 ⁿ－p_i+1,j,k ⁿ⁺¹－p_i,j+1,k ⁿ⁺¹－p_i,j,k+1 ⁿ⁺¹)
＋(3^0.5－Z^*χ)×(3^0.5＋3Z^*χ)^-1×(p_i+1,j,k ⁿ＋p_i,j+1,k ⁿ
＋p_i,j,k+1 ⁿ－p_i+1,j+1,k ⁿ⁺¹－p_i+1,j,k+1 ⁿ⁺¹
－p_i,j+1,k+1 ⁿ⁺¹)
＋(3^0.5×Z^*χ－1)×(3^0.5×Z^*χ＋1)^-1×
×(p_i,j,k ⁿ⁺¹－p_i+1,j+1,k+1 ⁿ)
を用いることができる。ここで上式（８）は（図３の白二重丸印の位置での）音圧p_i+1,j+1,k+1 ⁿ⁺¹についての計算式となっているが、他のＣＯＴＡ相当位置での音圧についても同様の境界計算式を規定することができる。

【0047】

以上、ＣＥ-ＦＤＴＤ法の境界計算式（例えば上式（４）,（７）,（８））は、複数の参照方向における参照方向毎に設定されており、また、所定の媒質（すなわち所定のρ値やc値）の空間において、境界量、すなわち音響インピーダンスZ、又は境界の垂直入射吸音率α若しくは同反射率r、を所定値に設定してやれば、一意に決定されることが理解される。なお、音響インピーダンスZの設定によって境界計算式が一意に決まることは、当然ながら、通常の（参照方向が１種類の）ＦＤＴＤ法にも当てはまる事項となっている。

【0048】

次に、上述したベクトル変数型のＦＤＴＤ法によるシミュレーションを説明する。以下具体的に、波動を音波、波動関連量を音圧（スカラ量）及び粒子速度（ベクトル量）としている。また音波の伝播する空間（音場）を、式表記の簡略化のため、所定の境界を有する２次元空間としている。この場合、支配方程式として、次式
（９） ∂p/∂t＋ρc²×∇u＝0
（１０） ∂u/∂t＋ρ^-1×∇ρ＝0
が用いられる。ちなみに上式（９）及び（１０）はそれぞれ、音波（音場）の連続方程式及び波動方程式であり、また式中のuは、粒子速度ベクトルとなっている（u＝(u_x, u_y)）。

【0049】

次いで支配方程式（９）及び（１０）を離散化することによって、音圧p及び粒子速度uについてのベクトル変数型ＦＤＴＤ法の離散化式
（１１） p_i,j ⁿ⁺¹＝p_i,j ⁿ－(ρc²Δt／Δ)×(u_x_i+0.5,j ^n+0.5
－u_x_i-0.5,j ^n+0.5＋u_y_i,j+0.5 ^n+0.5－u_y_i,j-0.5 ^n+0.5)
（１２） u_x_i+0.5,j ^n+0.5＝u_x_i+0.5,j ^n-0.5－(Δt／ρΔ)×(p_i+1,j ⁿ－p_i,j ⁿ)
（１３） u_y_i,j+0.5 ^n+0.5＝u_y_i,j+0.5 ^n-0.5－(Δt／ρΔ)×(p_i,j+1 ⁿ－p_i,j ⁿ)
が得られる。 さらに、境界計算式は、例えば次式
（１４） u_x_i+0.5,j ^n+0.5＝p_i,j ⁿ／Z
（１５） u_y_i,j+0.5 ^n+0.5＝p_i,j ⁿ／Z
のように求められる。

【0050】

ちなみに、以上に述べた離散化式は、図４（Ａ）に示したスタッガードグリッド上で離散化した結果となっている。このスタッガードグリッドでは、音圧pと粒子速度uとが、空間方向において順次、Δ×0.5だけずれて交互に定義され、時間方向においても順次、Δt×0.5だけずれて交互に定義される。ベクトル変数型ＦＤＴＤ法では通常、精度の低下を防止するべく、このスタッガードグリッドが用いられる。しかしながら勿論、図４（Ｂ）に示したコロケードグリッドを用いて、ベクトル変数型ＦＤＴＤ法を実施することも可能である。また、このコロケードグリッドを用いたベクトル変数型ＦＤＴＤ法の境界計算式も、媒質密度ρ、媒質での音速c、及び境界の音響インピーダンスZを含む式となるのである。

【0051】

このように、ベクトル変数型ＦＤＴＤ法の境界計算式（例えば上式（１４）,（１５））も、所定の媒質（すなわち所定のρ値やc値）の空間において、境界量、すなわち境界の音響インピーダンスZ、又は境界の垂直入射吸音率α若しくは同反射率r、を所定値に設定してやれば、一意に決定されることが理解される。

【0052】

以上、種々のタイプのＦＤＴＤ法について、その境界条件を含めて詳細に説明したが、本発明による波動計算処理においては、そのいずれのタイプのＦＤＴＤ法も採用可能となっている。すなわち、いずれのタイプのＦＤＴＤ法を採用するにしても、周波数帯毎に、当該周波数帯における音響インピーダンスZ、又は境界の垂直入射吸音率α若しくは同反射率r、を用いて境界計算式を設定し、この周波数帯毎に設定した境界計算式と、離散化式とを用いて計算対象量、ひいては波動関連量を算出するのである。

【0053】

［装置構成，波動計算プログラム・方法］
以下、本発明の一実施形態としての波動計算装置１の機能構成について、より詳細に説明を行う。同じく図１の機能ブロック図において、波動計算装置１は、通信インタフェース部１０１と、入力データ保存部１０２と、計算結果保存部１０３と、入出力インタフェース（ＩＦ）部１０４と、プロセッサ・メモリ（メモリ機能を備えた演算処理系）とを有する。ここで、プロセッサ・メモリは、本発明による波動計算プログラムを保存しており、また、コンピュータ機能を有していて、この波動計算プログラムを実行することによって波動計算処理を実施する。

【0054】

またこのことから波動計算装置１は、波動計算処理専用の装置であってもよいが、本発明による波動計算プログラムを搭載した、汎用のクラウドサーバや非クラウド型サーバであってもよく、さらにはパーソナルコンピュータ（ＰＣ）、ノート型若しくはタブレット型コンピュータや、スマートフォン、さらにはヘッドマウントディスプレイ（ＨＭＤ）といったようなウェアラブルデバイス等とすることも可能である。

【0055】

また、プロセッサ・メモリは、機能構成部として、
（ａ）サンプリング周波数設定部１１１ａ、境界計算式設定部１１１ｂ、離散化式設定部１１１ｃ、及び境界・空間内計算部１１１ｄを含む周波数帯毎（計算対象量）計算部１１１と、
（ｂ）音源・受音位置設定部１１２ａ及びバンドパスフィルタ処理部１１２ｂを含む波動関連量算出部１１２と、
（ｃ）入出力制御部１２１と
を有する。なお、これらの機能構成部は、プロセッサ・メモリに保存された、本発明による波動計算プログラムの実行によって具現する機能と捉えることができる。また、図１の機能ブロック図における波動計算装置１の機能構成部間を矢印で接続して示した処理の流れは、本発明による波動計算方法の一実施形態としても理解される。

【0056】

ちなみに、以下説明する実施形態では、所定の室壁（境界）で取り囲まれた対象室（空間，音場）を伝播する音波の音圧（計算対象量）、ひいてはインパルス応答（波動関連量）を、ＣＥ-ＦＤＴＤ法を用いて計算している。しかしながら上述したように、本発明の波動計算処理は勿論、これらの条件・対象・手法に限定されるものではない。

【0057】

＜周波数帯毎計算手段＞
同じく図１の機能ブロック図において、周波数帯毎計算部１１１は本実施形態において、３次元のＣＥ-ＦＤＴＤ法に基づき、音波に係る複数の周波数帯における周波数帯毎に、当該周波数帯について設定された境界計算式と、離散化式とを用いて、対象室における当該周波数帯での音圧を算出する。

【0058】

ここで本実施形態において、上記の複数の周波数帯は、当該分野で常用されている１２５Ｈｚ帯、２５０Ｈｚ帯、５００Ｈｚ帯、１０００Ｈｚ帯、２０００Ｈｚ帯、４０００Ｈｚ帯、８０００Ｈｚ帯、及び１６０００Ｈｚ帯に設定されている。さらに本実施形態において、周波数帯毎計算部１１１は、例えば外部ＰＣ２による指示の下、外部ＤＢ３から、通信インタフェース部１０１及び入出力制御部１２１を介して、各周波数帯での室壁の（垂直入射）吸音率αの値を取得する。またこれらの値から、次式
（１６） Z＝ρc×(r_p＋1)／(r_p－1)
r_p＝(1－α)^0.5
を用いて、各周波数帯での室壁の音響インピーダンスZを算出し、例えば入力データ保存部１０２に保存するのである。ここで、r_pは（垂直入射）音圧反射率であり、上式（１６）は、（垂直入射）吸音率αから音響インピーダンスを求めるための公式となっている。

【0059】

ちなみに、上式（１６）を用いて各周波数帯での音響インビーダンスZを算出する場合、媒質（例えば空気）の音響インピーダンスZ₀についても、各周波数帯での値を用いることが好ましい。また、吸音率αのＤＢとしては、例えば非特許文献：「吸音率データ」，[online]，［令和４年８月３０日検索］，インターネット＜URL: http://www.ymec.com/products/ymcad/doc2.htm＞が挙げられる。

【0060】

さらに、周波数帯毎計算部１１１は本実施形態において、対象室・室壁（音場）の形状・サイズに係るデータや、その他ＣＥ-ＦＤＴＤ法を実施するのに必要となるパラメータ（例えばΔ, Δt, ρ, c等）値を、例えば入出力ＩＦ部１０４や外部ＰＣ２から（通信インタフェース部１０１や入出力制御部１２１を介し）取得し、例えば入力データ保存部１０２に保存した上で、実施する波動計算処理に適用する。なお、適用されるパラメータとしての時間ステップΔtは、使用されるサンプリング周波数fsの逆数として、周波数帯毎計算部１１１のサンプリング周波数設定部１１１ａで設定されてもよい。

【0061】

同じく図１の機能ブロック図において、周波数帯毎計算部１１１の境界計算式設定部１１１ｂは、本実施形態において、使用すべき境界計算式（例えば上式（４）,（７）,（８））に対し、各周波数帯について設定された境界量としての音響インピーダンスZの値を適用し、各周波数帯での境界計算式を設定する。例えば、１２５Ｈｚ帯の音響インピーダンスZ_125が式中のZに代入された境界計算式（例えば上式（４）,（７）,（８））を、１２５Ｈｚ帯での境界計算式として用意する。

【0062】

また、周波数帯毎計算部１１１の離散化式設定部１１１ｃは、使用すべき離散化式に対し、必要なパラメータ値を適用（代入）し、離散化式を設定（用意）する。さらに境界・空間内計算部１１１ｄは、上記の周波数帯毎に、設定された当該周波数帯での境界計算式と、設定された離散化式とを用いて、室壁（境界）に隣接する室壁（境界）外の位置における当該周波数帯での音圧（計算対象量）と、これら室壁（境界）外の位置での音圧（計算対象量）を用いて計算される室壁（境界）での計算対象量を含む、対象室における当該周波数帯での音圧（計算対象量）と、を算出する。

【0063】

＜波動関連量算出手段＞
同じく図１の機能ブロック図において、波動関連量算出部１１２は、上記の周波数帯毎に算出された当該周波数帯での音圧（計算対象量）から、対象室におけるインパルス応答（波動関連量）を算出する。具体的に本実施形態では、波動関連量算出部１１２は、
（ａ）上記の周波数帯毎に、算出された当該周波数帯での音圧から、設定された音源位置と受音位置との間における当該周波数帯でのインパルス応答を算出し、
（ｂ）算出した各周波数帯でのインパルス応答を合算して、又は所定の周波数帯に係るインパルス応答に対しより大きい若しくは小さい重みを付与した上で合算して、対象室における（設定された音源位置と受音位置との間での）インパルス応答を算出する。

【0064】

ここで、上記（ａ）で用いられる音源位置及び受音位置は、本実施形態において音源・受音位置設定部１１２ａが、例えば入出力ＩＦ部１０４や外部ＰＣ２から（通信インタフェース部１０１や入出力制御部１２１を介し）取得し、例えば入力データ保存部１０２に保存したものとすることができる。また、上記（ａ）において波動関連量算出部１１２は、周波数帯毎計算部１１１において１つの周波数帯での音圧（分布）が算出されたタイミングで、この音圧計算結果を用いて、この１つの周波数帯でのインパルス応答を算出することも好ましい。勿論、全ての周波数帯での音圧計算結果を受け取って、全ての周波数帯でのインパルス応答を算出する設定にすることも可能である。

【0065】

さらに、上記（ｂ）において、合算される各周波数帯でのインパルス応答は、バンドパスフィルタ処理部１１２ｂによって、合算前にバンドパスフィルタ処理が施されたものとすることも好ましい。ここで、このバンドパスフィルタ処理は、この後示す具体例においても説明するが、該当する周波数帯に合わせて通過帯域が設定されたフィルタ処理とすることができる。これにより、ある周波数帯でのインパルス応答に含まれる、他の周波数帯に係る成分を低減又は除去することができ、最終的に算出されるインパルス応答の精度を向上させることが可能となる。

【0066】

さらに、上記（ｂ）において、算出した各周波数帯でのインパルス応答は、（例えば上述したバンドパスフィルタ処理が施された後、）所定の周波数帯に係るインパルス応答に対しより大きい若しくは小さい重みを付与した上で、合算されることも好ましい。また、この重み付け処理は、上述したバンドパスフィルタ処理における周波数帯毎のゲインの調整によって実現することも可能である。

【0067】

いずれにしてもこのような重み付け処理により、例えばＦＤＴＤ法で発生しがちな誤差を補正することも可能となる。実際、ＣＥ-ＦＤＴＤ法においては通常、サンプリング周波数あたりで音響が若干強めに推定されてしまう傾向がみられるところ、この重み付け処理を適切に実施することによって、このような誤差を低減又は解消することも可能となる。また、所定の周波数帯の音響を強める又は弱めることになるインパルス応答を、意図的に生成することもできる。これにより、例えばコンテンツ制作において、より自然な音響を演出すべく、所定の低周波数帯の音響を強めた音響信号を生成することに貢献することも可能となるのである。

【0068】

なお、この波動関連量算出部１１２でのインパルス応答算出結果や、上述した周波数帯毎計算部での音圧計算結果は、計算結果保存部１０３に保存され、この計算結果保存部１０３から適宜、読み出される設定とすることも好ましい。

【0069】

以下、以上に説明した周波数帯毎のインパルス応答の算出処理、及び目的である対象室におけるインパルス応答の算出処理の具体例１を説明する。なおこの具体例１は、各周波数帯でのインパルス応答を時間軸上で合算するものとなっている。

【0070】

（具体例１）
本具体例１では、境界条件に係る複数の周波数帯として、上述した１２５Ｈｚ帯、２５０Ｈｚ帯、・・・、１６０００Ｈｚ帯を採用し、また、サンプリング周波数fsを48kHzとする。サンプリング周波数は通常、使用する周波数の２倍（以上）に設定される。したがって本具体例１のサンプリング周波数fs（=48kHz）は、24kHzまでの周波数帯を採用可能な余裕のある値となっている。

【0071】

最初に、対象室・室壁の形状・サイズに係るデータや、その他ＣＥ-ＦＤＴＤ法を実施するのに必要となるパラメータ値、さらには対象室内における音源位置及び受音位置を設定する。ここで媒質は空気であり、また室壁（境界）は、岩綿吸音板で構成されているものとする。この岩綿吸音板の（垂直入射）吸音率αは、α_125＝0.60、α_250＝0.65、α_500＝0.66、α_1000＝0.70、α_2000＝0.79、α_4000＝0.95、・・・となっている。ここから、各周波数帯での境界の音響インピーダンス（Z_125, Z_250, Z_500, ・・・, Z_16000）を、上式（１６）を用いて決定するのである。なお一般に、境界の吸音率α、ひいては境界の音響インピーダンスZは、このように音波の周波数帯に大きく依存する量となっている。

【0072】

次いで、上記の周波数帯毎に、当該周波数帯での決定した音響インピーダンス（例えば当該周波数帯が１２５Ｈｚ帯ならばZ_125）を、境界計算式（例えば上式（４）,（７）,（８））に適用（代入）し、この当該周波数帯での境界計算式、及び離散化式（上式（２））を用いて、ＣＥ-ＦＤＴＤ法による音圧シミュレーションを行う。具体的には、音源位置にインパルス音圧を与えて、対象室内での音圧（分布）の時間変化を計算し、その際の各時点tでの受音位置での音圧値を、当該周波数帯での境界条件を適用したインパルス応答x(t)とする。これにより、各周波数帯でのインパルス応答（x_125(t), x_250(t), x_500(t), ・・・, x_16000(t)）が算出される。

【0073】

次いで、上記の周波数帯毎に、当該周波数帯に対応したバンドパスフィルタb(t)を設計し、当該周波数帯でのインパルス応答x(t)に対しこのバンドパスフィルタb(t)を畳み込み、フィルタ処理後のインパルス応答y(t)を算出する。本具体例１では、このバンドパスフィルタ（b_125(t), b_250(t), b_500(t), ・・・, b_16000(t)）として、上限遮断周波数が中心周波数の2^0.5（＝√2）倍であって下限遮断周波数が中心周波数の1／2^0.5倍であるオクターブバンドパスフィルタを採用する。この場合、例えば１２５Ｈｚ帯用のバンドパスフィルタb_125(t)は、通過帯域が90～177Hzのオクターブバンドパスフィルタとなる。なお、このようなバンドパスフィルタの設計にあたっては、後に行う合算処理において時間軸上でのずれを抑制するべく、全ての周波数帯におけるサンプル数を同一にすることが好ましい。

【0074】

いずれにしても、以上に説明したフィルタ処理において、例えば１２５Ｈｚ帯でのフィルタ処理後のインパルス応答y_125(t)は、＊を畳み込み積分演算子として次式
（１７） y_125(t)＝x_125(t)＊b_125(t)
によって算出される。これにより、目的である対象室におけるインパルス応答z(t)は、次式
（１８） z(t)＝y_125(t)＋y_500(t)＋・・・＋y_16000(t)
によって、すなわち各周波数帯でのインパルス応答y(t)を合算することによって算出されるのである。ちなみに上述したように、上式（１８）において、各周波数帯でのインパルス応答y(t)に所定の重み付けを行った上で合算することも可能である。

【0075】

同じく図１の機能ブロック図において、波動関連量算出部１１２は、以上説明した実施形態（及び具体例１）とは別の実施形態として、
（ａ）上記の周波数帯毎に、当該周波数帯での音圧（計算対象量）から当該周波数帯でのインパルス応答（波動関連量）を算出し、
（ｂ）算出した各周波数帯でのインパルス応答に対しフーリエ変換処理を施し、このフーリエ変換処理を施した各波数帯でのインパルス応答から、対象室におけるフーリエ変換されたインパルス応答を決定し、
（ｃ）決定したフーリエ変換されたインパルス応答に対し逆フーリエ変換処理を施して、対象室におけるインパルス応答を算出する
ことも好ましい。

【0076】

このように、フーリエ変換を用いて周波数軸上でインパルス応答を計算することによって、各周波数帯でのインパルス応答から適応したい周波数成分のみを取り出して、目的とする対象室におけるインパルス応答を精度良く算出することが可能となる。すなわち、このようなフーリエ変換を用いたインパルス応答算出処理によれば、実際には通過帯域の直前及び直後でフィルタ処理対象量が若干残留してしまうバンドパスフィルタを用いることなく、周波数成分の調整を行うことが可能となり、これにより、目的であるインパルス応答の算出精度を向上させることもできるのである。

【0077】

以下、このようなフーリエ変換を用いたインパルス応答算出処理の具体例２を説明する。なおこの具体例２は、各周波数帯でのインパルス応答を周波数軸上で合算するものとなっている。

【0078】

（具体例２）
本具体例２においては最初に、上記の具体例１と同様にして、各周波数帯でのインパルス応答（x_125(t), x_250(t), x_500(t), ・・・, x_16000(t)）を算出する。次いで、これらのインパルス応答の各々に対しフーリエ変換を施し、フーリエ変換後のインパルス応答（X_125(k), X_250(k), X_500(k), ・・・, x_16000(k)）を取得する。ここでkは、フーリエ変換のサイズとサンプリング周波数とで決まる周波数ポイントとしての周波数ビンである。

【0079】

次いで、これらのフーリエ変換後のインパルス応答から、対象室におけるフーリエ変換されたインパルス応答Z(k)を決定する。具体的には、次式
（１９） Z(k)＝X_125(k) （90＜k≦177）
＝X_250(k) （177＜k≦354）
・・・
＝X_16000(k) （11314＜k≦22627）
のように、インパルス応答Z(k)が決定される。最後に、このインパルス応答Z(k)に対し逆フーリエ変換を施すことによって、目的である対象室におけるインパルス応答z(t)を算出するのである。

【0080】

＜サンプリング周波数に関する他の実施形態＞
次に、ＦＤＴＤ法に係るサンプリング周波数に関する好適な他の実施形態の説明を行う。同じく図１の機能ブロック図において、周波数帯毎計算部１１１のサンプリング周波数設定部１１１ａは、本実施形態において、上記の周波数帯毎に、当該周波数帯に合わせて設定されたサンプリング周波数fsに対応した境界計算式及び離散化式を用いて、当該周波数帯での音圧（計算対象量）を算出する。

【0081】

例えば、各周波数帯において、サンプリング周波数fsを、次式
（２０） fs＝(当該周波数帯の上限遮断周波数)×2
(当該周波数帯の中心周波数)×2×2^0.5
を用いて設定してもよい。この場合、例えば１２５Ｈｚ帯でのサンプリング周波数fs_125は、354（≒125××2×2^0.5）Hzであって、サンプリングを行うにあたり十分に大きな値となっているのであり、例えば１６０００Ｈｚ帯に合わせた（上記の具体例１及び２で設定した）48kHz等にまで高くする必要はないのである。

【0082】

このように、各周波数帯でのサンプリング周波数（fs_125, fs_250, fs_500, ・・・, fs_16000）を、できるだけ小さく設定すれば、サンプリング周波数の逆数としての時間ステップΔtはより大きく設定される。ここでＦＤＴＤ法では通常、安定条件を表すＣＦＬ数χ（＝c×Δt／Δ，上式（３）参照）が所定値に固定されている。したがって、各周波数帯においてグリッド間隔Δも、時間ステップΔtに比例してより大きく設定され、これにより、グリッド数がより小さい値で済むことになる。その結果、波動計算の負担を減らし、一律のサンプリング周波数を用いる場合と比べて計算時間を短縮することも可能となるのである。

【0083】

ちなみに、以上に述べた周波数帯毎のサンプリング周波数の設定は、例えば具体例１で示した実施形態でも、具体例２で示した実施形態でも、さらにはＣＥ-ＦＤＴＤ法以外のＦＤＴＤ法を用いる場合においても実施することができる。また勿論、いずれの場合においても計算時間の短縮化が可能となるのである。

【0084】

同じく図１の機能ブロック図において、波動関連量算出部１１２は本実施形態において、対象室におけるインパルス応答の計算結果を、例えば計算結果保存部１０３に保存した上で、入出力制御部１２１を介して入出力ＩＦ部１０４へ出力してもよい。この場合、入出力ＩＦ部１０４は、ディスプレイを備えている場合とはなるが、このインパルス応答計算結果を表示することができる。また、このインパルス応答計算結果は、入出力制御部１２１及び通信インタフェース部１０１を介して外部ＰＣ２へ送信され、この外部ＰＣ２で利用されてもよい。

【0085】

以上詳細に説明したように、本発明によれば、ＦＤＴＤ法を用いた波動関連量算出処理において、周波数帯によって生じる、計算値の実測値からの乖離を抑制し、波動関連量算出処理の精度向上を図ることができる。例えば、所定の音場における波動関連量としてのインパルス応答を、より高い精度で算出することも可能となるのである。

【0086】

また本発明によれば、１つの応用例とはなるが、収音空間のインパルス応答をより高い精度で推定して、この収音空間での音響を、受聴者に対しより忠実に再現してみせることに貢献することもできる。

【0087】

さらに例えば、所定のコンサートホール、劇場や演説会場におけるより精度の高いインパルス応答を算出して、あたかもこれらのコンサートホール、劇場や演説会場において収音されたかのような豊かな音響特性を有する、演奏、演劇や演説に関する音楽・音声教材を作成することも可能となる。またこのように作成された音楽・音声教材を例えばオンライン授業の中で活用し、都市部だけでなく地方の子供達や受講者に対し、より質の高い教育を受ける機会を提供することも可能となる。すなわち本発明によれば、国連が主導する持続可能な開発目標（ＳＤＧｓ）の目標４「すべての人々に包摂的かつ公平で質の高い教育を提供し、生涯学習の機会を促進する」に貢献することも可能となるのである。

【0088】

以上に述べた本発明の種々の実施形態について、本発明の技術思想及び見地の範囲の種々の変更、修正及び省略は、当業者によれば容易に行うことができる。前述の説明はあくまで例であって、何ら制約しようとするものではない。本発明は、特許請求の範囲及びその均等物として限定するものにのみ制約される。

【符号の説明】

【0089】

１ 波動計算装置
１０１ 通信インタフェース部
１０２ 入力データ保存部
１０３ 計算結果保存部
１０４ 入出力インタフェース（ＩＦ）部
１１１ 周波数帯毎（計算対象量）計算部
１１１ａ サンプリング周波数設定部
１１１ｂ 境界計算式設定部
１１１ｃ 離散化式設定部
１１１ｄ 境界・空間内計算部
１１２ 波動関連量算出部
１１２ａ 音源・受音位置設定部
１１２ｂ バンドパスフィルタ処理部
１２１ 入出力制御部
２ 外部ＰＣ
３ 外部データベース（ＤＢ）

【図1】

【図2】

【図3】

【図4】