ホーム > 特許ランキング > 日本電信電話株式会社
知財求人 - 知財ポータルサイト「IP Force」
▶ 日本電信電話株式会社の特許一覧 ▶ 国立大学法人電気通信大学の特許一覧
(19)【発行国】日本国特許庁(JP)
(12)【公報種別】公開特許公報(A)
(11)【公開番号】P2025037291
(43)【公開日】2025-03-18
(54)【発明の名称】音圧推定装置、音圧推定方法、プログラム
(51)【国際特許分類】
H04R 3/00 20060101AFI20250311BHJP
G10K 15/00 20060101ALI20250311BHJP
G01H 17/00 20060101ALI20250311BHJP
【ＦＩ】
H04R3/00 320
G10K15/00 L
G01H17/00 Z
【審査請求】未請求
【請求項の数】7
【出願形態】ＯＬ
(21)【出願番号】P 2023144105
(22)【出願日】2023-09-06
(71)【出願人】
【識別番号】000004226
【氏名又は名称】日本電信電話株式会社
(71)【出願人】
【識別番号】504133110
【氏名又は名称】国立大学法人電気通信大学
(74)【代理人】
【識別番号】100121706
【弁理士】
【氏名又は名称】中尾 直樹
(74)【代理人】
【識別番号】100128705
【弁理士】
【氏名又は名称】中村 幸雄
(74)【代理人】
【識別番号】100147773
【弁理士】
【氏名又は名称】義村 宗洋
(72)【発明者】
【氏名】伊藤 弘章
(72)【発明者】
【氏名】小塚 詩穂里
(72)【発明者】
【氏名】鎌土 記良
(72)【発明者】
【氏名】信夫 直樹
(72)【発明者】
【氏名】吉松 亨真
(72)【発明者】
【氏名】羽田 陽一
【テーマコード（参考）】
2G064
5D220
【Ｆターム（参考）】
2G064AB01
2G064AB02
2G064AB13
2G064BA02
2G064BD02
2G064CC02
2G064CC41
2G064DD02
5D220BA01
5D220BA30
(57)【要約】
【課題】周波数領域への変換を介することなく、頭部近傍に配置した複数のマイクを用いて観測された音圧から耳元の音圧を推定する技術を提供する。
【解決手段】L個のマイクが頭部近傍の球面上に配置される点における音圧p1, …, pLから次数0からNまでの音圧p1, …, pLの重み付け和q0, …, qNを計算する重み付け加算部と、次数0からNまでの重み付け和q0, …, qNから次数0からNまでのIIRフィルタη0, …, ηNを用いて耳元の音圧^pを計算する音圧計算部とを含む音圧推定装置であって、(re, Ω1), …, (re, ΩL)（ただし、reは頭部近傍の球面の半径）をL個のマイクが配置される点、(a, Ωa)を耳元の位置を表す点（ただし、aは耳元が位置する球面の半径）とし、次数nのIIRフィルタηnは、周波数領域における所定の関数Hn(kr)に基づいて設計された時間領域フィルタである。
【選択図】図７
【解決手段】L個のマイクが頭部近傍の球面上に配置される点における音圧p1, …, pLから次数0からNまでの音圧p1, …, pLの重み付け和q0, …, qNを計算する重み付け加算部と、次数0からNまでの重み付け和q0, …, qNから次数0からNまでのIIRフィルタη0, …, ηNを用いて耳元の音圧^pを計算する音圧計算部とを含む音圧推定装置であって、(re, Ω1), …, (re, ΩL)（ただし、reは頭部近傍の球面の半径）をL個のマイクが配置される点、(a, Ωa)を耳元の位置を表す点（ただし、aは耳元が位置する球面の半径）とし、次数nのIIRフィルタηnは、周波数領域における所定の関数Hn(kr)に基づいて設計された時間領域フィルタである。
【選択図】図７
【特許請求の範囲】
【請求項1】
Lを頭部近傍の球面上に配置されるマイクの数（ただし、Lは2以上の整数）、p1, …, pLをL個のマイクが配置される点における音圧とし、
音圧p1, …, pLから、次数0からNまでの音圧p1, …, pLの重み付け和q0, …, qNを計算する重み付け加算部と、
次数0からNまでの重み付け和q0, …, qNから、次数0からNまでのIIRフィルタη0, …, ηNを用いて耳元の音圧^pを計算する音圧計算部と
を含む音圧推定装置であって、
(re, Ω1), …, (re, ΩL)（ただし、reは頭部近傍の球面の半径）をL個のマイクが配置される点、(a, Ωa)を耳元の位置を表す点（ただし、aは耳元が位置する球面の半径）とし、
次数nのIIRフィルタηn(n=0, …, N)は、音圧の内挿関係を表す周波数領域における関数Hn(kr)(n=0, …, N, ただし、kは波数)に基づいて設計された時間領域フィルタである
音圧推定装置。
音圧p1, …, pLから、次数0からNまでの音圧p1, …, pLの重み付け和q0, …, qNを計算する重み付け加算部と、
次数0からNまでの重み付け和q0, …, qNから、次数0からNまでのIIRフィルタη0, …, ηNを用いて耳元の音圧^pを計算する音圧計算部と
を含む音圧推定装置であって、
(re, Ω1), …, (re, ΩL)（ただし、reは頭部近傍の球面の半径）をL個のマイクが配置される点、(a, Ωa)を耳元の位置を表す点（ただし、aは耳元が位置する球面の半径）とし、
次数nのIIRフィルタηn(n=0, …, N)は、音圧の内挿関係を表す周波数領域における関数Hn(kr)(n=0, …, N, ただし、kは波数)に基づいて設計された時間領域フィルタである
音圧推定装置。
【請求項2】
請求項１に記載の音圧推定装置であって、
前記重み付け加算部は、次式により次数nの重み付け和qn(n=0, …, N)を計算するものであり、
（ただし、Yn mは球面調和関数を表す。）
前記音圧計算部は、次式により音圧^pを計算するものである
ことを特徴とする音圧推定装置。
前記重み付け加算部は、次式により次数nの重み付け和qn(n=0, …, N)を計算するものであり、
【数40】
【数41】
（ただし、Yn mは球面調和関数を表す。）
前記音圧計算部は、次式により音圧^pを計算するものである
【数42】
ことを特徴とする音圧推定装置。
【請求項3】
請求項１に記載の音圧推定装置であって、
関数Hn(kr)(n=0, …, N)は、次式で表される
（ただし、jnはn次の球ベッセル関数、hn (2)はn次の第2種球ハンケル関数を表す。）
ことを特徴とする音圧推定装置。
関数Hn(kr)(n=0, …, N)は、次式で表される
【数43】
【数44】
（ただし、jnはn次の球ベッセル関数、hn (2)はn次の第2種球ハンケル関数を表す。）
ことを特徴とする音圧推定装置。
【請求項4】
請求項３に記載の音圧推定装置であって、
次数0のIIRフィルタη0は、関数H0(kr)のs領域上の表現である関数H0(s)をs-z変換することにより得られる次式の関数H0(z)の係数am(m=1, 2, 3, 2u, 2u+1, 2u+2, 2u+3), bm(m=0, 1, 2, 3)を用いて設計された時間領域フィルタである
ことを特徴とする音圧推定装置。
次数0のIIRフィルタη0は、関数H0(kr)のs領域上の表現である関数H0(s)をs-z変換することにより得られる次式の関数H0(z)の係数am(m=1, 2, 3, 2u, 2u+1, 2u+2, 2u+3), bm(m=0, 1, 2, 3)を用いて設計された時間領域フィルタである
【数45】
ことを特徴とする音圧推定装置。
【請求項5】
請求項３に記載の音圧推定装置であって、
次数nのIIRフィルタηn(n=0, …, N)は、関数Hn(kr)の極となる周波数fw(w=1, …, W、Wは関数Hn(kr)の極の数)を用いて次式により表される関数Hn(z)に基づいて得られるIIRフィルタ係数ηn(w) (w=1, …, 2W)を用いて設計された時間領域フィルタである
（ただし、vw=r×exp(jθw), θw=2πfw/fsであり、rは1に略等しく、fsはサンプリング周波数を表す。）
ことを特徴とする音圧推定装置。
次数nのIIRフィルタηn(n=0, …, N)は、関数Hn(kr)の極となる周波数fw(w=1, …, W、Wは関数Hn(kr)の極の数)を用いて次式により表される関数Hn(z)に基づいて得られるIIRフィルタ係数ηn(w) (w=1, …, 2W)を用いて設計された時間領域フィルタである
【数46】
（ただし、vw=r×exp(jθw), θw=2πfw/fsであり、rは1に略等しく、fsはサンプリング周波数を表す。）
ことを特徴とする音圧推定装置。
【請求項6】
Lを頭部近傍の球面上に配置されるマイクの数（ただし、Lは2以上の整数）、p1, …, pLをL個のマイクが配置される点における音圧とし、
音圧推定装置が、音圧p1, …, pLから、次数0からNまでの音圧p1, …, pLの重み付け和q0, …, qNを計算する重み付け加算ステップと、
前記音圧推定装置が、次数0からNまでの重み付け和q0, …, qNから、次数0からNまでのIIRフィルタη0, …, ηNを用いて耳元の音圧^pを計算する音圧計算ステップと
を含む音圧推定方法であって、
(re, Ω1), …, (re, ΩL)（ただし、reは頭部近傍の球面の半径）をL個のマイクが配置される点、(a, Ωa)を耳元の位置を表す点（ただし、aは耳元が位置する球面の半径）とし、
次数nのIIRフィルタηn(n=0, …, N)は、音圧の内挿関係を表す周波数領域における関数Hn(kr)(n=0, …, N, ただし、kは波数)に基づいて設計された時間領域フィルタである
音圧推定方法。
音圧推定装置が、音圧p1, …, pLから、次数0からNまでの音圧p1, …, pLの重み付け和q0, …, qNを計算する重み付け加算ステップと、
前記音圧推定装置が、次数0からNまでの重み付け和q0, …, qNから、次数0からNまでのIIRフィルタη0, …, ηNを用いて耳元の音圧^pを計算する音圧計算ステップと
を含む音圧推定方法であって、
(re, Ω1), …, (re, ΩL)（ただし、reは頭部近傍の球面の半径）をL個のマイクが配置される点、(a, Ωa)を耳元の位置を表す点（ただし、aは耳元が位置する球面の半径）とし、
次数nのIIRフィルタηn(n=0, …, N)は、音圧の内挿関係を表す周波数領域における関数Hn(kr)(n=0, …, N, ただし、kは波数)に基づいて設計された時間領域フィルタである
音圧推定方法。
【請求項7】
請求項１ないし５のいずれか１項に記載の音圧推定装置としてコンピュータを機能させるためのプログラム。
【発明の詳細な説明】
【技術分野】
【0001】
本発明は、頭部近傍の能動騒音制御（アクティブノイズコントロール）技術に関し、特に複数のマイクを用いて所望の制御点の音圧を推定する技術に関する。
【背景技術】
【0002】
参照マイク、エラーマイク、キャンセルスピーカを用いて、特定の位置での騒音を抑圧する能動騒音制御技術がある（非特許文献１参照）。
【0003】
座席上の能動騒音制御システムにおいて、耳元における騒音抑圧性能を担保するため、頭部近傍に配置したマイクを用いて耳元の音圧を推定する手法が検討されている。例えば、剛球であると仮定した頭部を含むような球面上に配置したマイクを用い、球面調和関数展開および音圧の内挿関係を利用することで耳元の音圧を推定することが考えられる。
【先行技術文献】
【非特許文献】
【0004】
【非特許文献1】梶川嘉延, “アクティブノイズコントロールの最近の話題と応用,” 情報処理学会研究報告, Vol.2015-MUS-107, No.3, pp.1-6, 2015.
【発明の概要】
【発明が解決しようとする課題】
【0005】
しかし、上述した方法では周波数領域における演算を用いるため、時間領域の信号を周波数領域の信号に変換する処理が必要となる。したがって、上述した方法を即時性の高い処理が必要となる能動騒音制御システムに適用しようとすると、要求される処理性能を満たすことができない場合がある。
【0006】
そこで本発明では、周波数領域への変換を介することなく、頭部近傍に配置した複数のマイクを用いて観測された音圧から耳元の音圧を推定する技術を提供することを目的とする。
【課題を解決するための手段】
【0007】
本発明の一態様は、Lを頭部近傍の球面上に配置されるマイクの数（ただし、Lは2以上の整数）、p1, …, pLをL個のマイクが配置される点における音圧とし、音圧p1, …, pLから、次数0からNまでの音圧p1, …, pLの重み付け和q0, …, qNを計算する重み付け加算部と、次数0からNまでの重み付け和q0, …, qNから、次数0からNまでのIIRフィルタη0, …, ηNを用いて耳元の音圧^pを計算する音圧計算部とを含む音圧推定装置であって、(re, Ω1), …, (re, ΩL)（ただし、reは頭部近傍の球面の半径）をL個のマイクが配置される点、(a, Ωa)を耳元の位置を表す点（ただし、aは耳元が位置する球面の半径）とし、次数nのIIRフィルタηn(n=0, …, N)は、音圧の内挿関係を表す周波数領域における関数Hn(kr)(n=0, …, N, ただし、kは波数)に基づいて設計された時間領域フィルタである。
【発明の効果】
【0008】
本発明によれば、周波数領域における演算を用いることなく、耳元の音圧を推定することが可能となる。
【図面の簡単な説明】
【0009】
【図1】複数のマイクが頭部近傍の球面上に配置される様子を示す図である。
【図2】関数H0(kr)の概形を示す図である。
【図3】関数Hn(kr)の極となる周波数fkの様子を示す図である。
【図4】帯域制限パルスに対する推定波形と観測波形を示す図（補正前）である。
【図5】帯域制限パルスに対する推定波形と観測波形を示す図（補正後）である。
【図6】白色雑音に対する推定波形と観測波形を示す図（補正後）である。
【図7】音圧推定装置１００の構成を示すブロック図である。
【図8】音圧推定装置１００の動作を示すフローチャートである。
【図9】本発明の実施形態における各装置を実現するコンピュータの機能構成の一例を示す図である。
【発明を実施するための形態】
【0010】
以下、本発明の実施形態について、詳細に説明する。なお、同じ機能を有する構成部には同じ番号を付し、重複説明を省略する。
【0011】
各実施形態の説明に先立って、この明細書における表記方法について説明する。
【0012】
^（キャレット）は上付き添字を表す。例えば、xy^zはyzがxに対する上付き添字であり、xy^zはyzがxに対する下付き添字であることを表す。また、_（アンダースコア）は下付き添字を表す。例えば、xy_zはyzがxに対する上付き添字であり、xy_zはyzがxに対する下付き添字であることを表す。
【0013】
ある文字xに対する^xや~xのような上付き添え字の”^”や”~”は、本来”x”の真上に記載されるべきであるが、明細書の記載表記の制約上、^xや~xと記載しているものである。
【0014】
＜技術的背景＞
本発明の実施形態では、頭部を剛球であると仮定し、頭部近傍の球面上に配置した複数のマイクを用いて観測される信号の音圧の重み付け和に対して、耳元の音圧を推定するために設計された時間領域フィルタであるIIRフィルタを畳み込むことにより、耳元の音圧を推定する。ここで、頭部近傍の球面とは、頭部を含むことができる球面のことである。図１は、複数のマイクが頭部近傍の球面上に配置される様子を示す図である。以下、本発明の実施形態における音圧推定処理と当該処理で用いるIIRフィルタの導出方法について説明する。
本発明の実施形態では、頭部を剛球であると仮定し、頭部近傍の球面上に配置した複数のマイクを用いて観測される信号の音圧の重み付け和に対して、耳元の音圧を推定するために設計された時間領域フィルタであるIIRフィルタを畳み込むことにより、耳元の音圧を推定する。ここで、頭部近傍の球面とは、頭部を含むことができる球面のことである。図１は、複数のマイクが頭部近傍の球面上に配置される様子を示す図である。以下、本発明の実施形態における音圧推定処理と当該処理で用いるIIRフィルタの導出方法について説明する。
【0015】
＜＜１：時間領域における音圧推定処理＞＞
まず、頭部近傍の球面上に配置した2以上のマイクを用いて観測される信号の音圧から、周波数領域において耳元の音圧を推定することを考える。以下、空間の点を極座標(r, Ω)（ただし、Ωは2つの偏角θ, φを表す）を用いて表すこととする。Lをマイクの数とし、(re, Ω1), …, (re, ΩL)をL個のマイクが配置される点（ただし、reはL個のマイクが配置される頭部近傍の球面の半径）、(a, Ωa)を耳元の位置を表す点（ただし、aは耳元が位置する球面の半径）とする。また、p1, …, pLをL個のマイクが配置される点における音圧とする。
まず、頭部近傍の球面上に配置した2以上のマイクを用いて観測される信号の音圧から、周波数領域において耳元の音圧を推定することを考える。以下、空間の点を極座標(r, Ω)（ただし、Ωは2つの偏角θ, φを表す）を用いて表すこととする。Lをマイクの数とし、(re, Ω1), …, (re, ΩL)をL個のマイクが配置される点（ただし、reはL個のマイクが配置される頭部近傍の球面の半径）、(a, Ωa)を耳元の位置を表す点（ただし、aは耳元が位置する球面の半径）とする。また、p1, …, pLをL個のマイクが配置される点における音圧とする。
【0016】
耳元の位置を表す点(a, Ωa)における音圧^p(a, Ωa)は、球面調和関数展開を用いると、次式により表される。
ここで、Yn mは球面調和関数を表す。また、Pnm(r)は球面調和関数Yn mに対する半径r上の展開係数を表す。
【数1】
ここで、Yn mは球面調和関数を表す。また、Pnm(r)は球面調和関数Yn mに対する半径r上の展開係数を表す。
【0017】
式(1)は、音圧の内挿関係を用いると、以下のように変形できる。
ここで、
（ただし、kは波数であり、jnはn次の球ベッセル関数、hn (2)はn次の第2種球ハンケル関数を表す）であり、右肩の*, ’はそれぞれ複素共役、微分を表す。式(3)及び式(4)が音圧の内挿関係を表す式である。また、fを周波数、cを音速とすると、波数kはk=2πf/cと表される。なお、式(2)において周波数に依存する項はHnのみであることに留意する。
【数2】
ここで、
【数3】
【数4】
（ただし、kは波数であり、jnはn次の球ベッセル関数、hn (2)はn次の第2種球ハンケル関数を表す）であり、右肩の*, ’はそれぞれ複素共役、微分を表す。式(3)及び式(4)が音圧の内挿関係を表す式である。また、fを周波数、cを音速とすると、波数kはk=2πf/cと表される。なお、式(2)において周波数に依存する項はHnのみであることに留意する。
【0018】
式(2)からわかるように、項Hnに対応するIIRフィルタηnを用いることができれば、時間領域において音圧を推定することが可能となる。このときの時間領域における音圧推定処理は、以下の手順となる。
【0019】
（１）L個のマイクが配置される点における音圧p1, …, pLを入力とし、次数nごとに音圧p1, …, pLの重み付け和qn(n=0, …, N)を計算する。
ここで、γλnが音圧pλに対する次数nの重みであり、重みγλnは予め計算しておくとよい。
【数5】
【数6】
ここで、γλnが音圧pλに対する次数nの重みであり、重みγλnは予め計算しておくとよい。
【0020】
（２）次数0からNまでの重み付け和q0, …, qNを入力とし、次数0からNまでのIIRフィルタη0, …, ηNを用いて次式により耳元の音圧^pを計算する。
つまり、耳元の音圧^pは、IIRフィルタηnを用いて重み付け和qnをフィルタリングして得られるηnqnの和として得られる。
【数7】
つまり、耳元の音圧^pは、IIRフィルタηnを用いて重み付け和qnをフィルタリングして得られるηnqnの和として得られる。
【0021】
＜＜２：IIRフィルタの導出方法＞＞
ここでは、項Hnに対応するIIRフィルタηnを導出する方法について説明する。具体的には、以下の２つの方法について説明する。
ここでは、項Hnに対応するIIRフィルタηnを導出する方法について説明する。具体的には、以下の２つの方法について説明する。
【0022】
（１）s-z変換を用いる方法
（２）関数Hn(kr)の極の数値計算を利用する方法
方法（１）、（２）はいずれもz平面を用いて導出する方法である。そこで、方法（１）、（２）の説明を始める前に、z平面上の伝達関数H(z)とIIRフィルタのフィルタ係数の関係について考察する。一般に、IIRフィルタの出力信号y(i)（ただし、iはサンプル時刻を表す）は、x(i)を入力信号として、次式により表現することができる。
ここで、am(m=1, …, M), bm(m=0, …, W)はフィルタ係数である。
（２）関数Hn(kr)の極の数値計算を利用する方法
方法（１）、（２）はいずれもz平面を用いて導出する方法である。そこで、方法（１）、（２）の説明を始める前に、z平面上の伝達関数H(z)とIIRフィルタのフィルタ係数の関係について考察する。一般に、IIRフィルタの出力信号y(i)（ただし、iはサンプル時刻を表す）は、x(i)を入力信号として、次式により表現することができる。
【数8】
ここで、am(m=1, …, M), bm(m=0, …, W)はフィルタ係数である。
【0023】
上式の両辺をz変換し、z平面上の伝達関数H(z)を考えると、次式が得られる。
これより、フィルタ係数am(m=1, …, M), bm(m=0, …, W)は関数H(z)を表す分数多項式の係数として現れることがわかる。したがって、音圧の内挿関係を表す式(3)に対応する周波数領域における関数Hn(kr)をz平面上の関数Hn(z)に変換することができれば、IIRフィルタηnを導出することができる。
【数9】
これより、フィルタ係数am(m=1, …, M), bm(m=0, …, W)は関数H(z)を表す分数多項式の係数として現れることがわかる。したがって、音圧の内挿関係を表す式(3)に対応する周波数領域における関数Hn(kr)をz平面上の関数Hn(z)に変換することができれば、IIRフィルタηnを導出することができる。
【0024】
＜＜＜２－１：s-z変換を用いる方法＞＞＞
周波数領域における関数Hn(kr)をz平面上の関数Hn(z)に変換する手順は、以下の通りである。
周波数領域における関数Hn(kr)をz平面上の関数Hn(z)に変換する手順は、以下の通りである。
【0025】
（１）関数Hn(kr)をs領域表現にすることで、s平面上の関数Hn(s)を得る。そのために、周波数領域がs平面の虚軸上の集合に対応することを利用する。
【0026】
（２）s平面上の関数Hn(s)をz平面上の表現とすることで、z平面上の関数Hn(z)を得る。そのために、s平面上の関数をz平面に対応付ける変換手法であるs-z変換を利用する。
【0027】
まず、手順（１）について説明する。関数Hn(kr)をs領域上の表現とするために、関数Hn(kr)をs=jω（ただし、jは虚数単位、ω(=2πf)は角周波数である。）の関数とみる。周波数領域における関数をs=jωの関数とみる一例として、0次の第2種球ハンケル関数h0
(2)(kr)から関数h0
(2)(s)を得る例を示す。まず、関数h0
(2)(kr)を以下のように変形する。
これを用いて、関数h0 (2)(s)を得ることができる。
関数Hn(s)を得るためには、式(4)からわかるように、球ベッセル関数jn(kr)、第2種球ハンケル関数hn (2)(kr)、球ベッセル関数の導関数jn’(kr)、第2種球ハンケル関数の導関数hn’(2)(kr)のs領域上の表現が求まれば十分である。ここで、これら４種類の関数のs領域上の表現を求めるために、いくつかの式に注目する。球ベッセル関数jn(kr)と2種類の球ハンケル関数hn (1), hn (2)の間には以下の式が成り立つ。
また、m=1, 2として、球ハンケル関数の導関数hn’(m)と球ハンケル関数hn (m)の間に以下の式が成り立つ。
これらの式から、関数jn(kr), hn (2)(kr), jn’(kr), hn’(2)(kr)のs領域上の表現である関数jn(s)、hn (2)(s)、jn’(s)、hn’(2)(s)が求まるためには、第1種球ハンケル関数hn (1)(s)、第2種球ハンケル関数hn (2)(s)が求まればよいことがわかる。
【数10】
これを用いて、関数h0 (2)(s)を得ることができる。
【数11】
関数Hn(s)を得るためには、式(4)からわかるように、球ベッセル関数jn(kr)、第2種球ハンケル関数hn (2)(kr)、球ベッセル関数の導関数jn’(kr)、第2種球ハンケル関数の導関数hn’(2)(kr)のs領域上の表現が求まれば十分である。ここで、これら４種類の関数のs領域上の表現を求めるために、いくつかの式に注目する。球ベッセル関数jn(kr)と2種類の球ハンケル関数hn (1), hn (2)の間には以下の式が成り立つ。
【数12】
また、m=1, 2として、球ハンケル関数の導関数hn’(m)と球ハンケル関数hn (m)の間に以下の式が成り立つ。
【数13】
【数14】
これらの式から、関数jn(kr), hn (2)(kr), jn’(kr), hn’(2)(kr)のs領域上の表現である関数jn(s)、hn (2)(s)、jn’(s)、hn’(2)(s)が求まるためには、第1種球ハンケル関数hn (1)(s)、第2種球ハンケル関数hn (2)(s)が求まればよいことがわかる。
【0028】
第1種球ハンケル関数hn
(1)(s)、第2種球ハンケル関数hn
(2)(s)をn=0, 1, 2に対して求めると、以下のようになる。
ここで、
である。
【数15】
【数16】
【数17】
【数18】
【数19】
【数20】
ここで、
【数21】
である。
【0029】
以上の結果を用いて、関数H0(s)を求めると次式が得られる。
次に、手順（２）について説明する。H0(s)からH0(z)を得るため、s-z変換を施す。その結果、以下の式が得られる。
ここで、Tはサンプル周期を表す。なお、rがreに近い値であるものとして、|α|<1, u=floor((re-a)/c×1/T)（ただし、uは2以上の整数であることが好ましい）とするとき、以下の式が成り立つことを利用している。
なお、am, b’mはいずれも実数となる。
【数22】
次に、手順（２）について説明する。H0(s)からH0(z)を得るため、s-z変換を施す。その結果、以下の式が得られる。
【数23】
ここで、Tはサンプル周期を表す。なお、rがreに近い値であるものとして、|α|<1, u=floor((re-a)/c×1/T)（ただし、uは2以上の整数であることが好ましい）とするとき、以下の式が成り立つことを利用している。
【数24】
なお、am, b’mはいずれも実数となる。
【0030】
したがって、r=reとすると、u=7の場合にはフィルタ係数am, bmとして以下の表に示す値が算出される。
ただし、
である。なお、αは(re-a)/cTの非整数部分が小さくなるほど1に近い値になるように調整する。例えば、re=0.15, a=0.10の場合、(re-a)/cTはおよそ7.059となるので、αを1とする。
【表1】
ただし、
【数25】
【数26】
である。なお、αは(re-a)/cTの非整数部分が小さくなるほど1に近い値になるように調整する。例えば、re=0.15, a=0.10の場合、(re-a)/cTはおよそ7.059となるので、αを1とする。
【0031】
特に、(re-a)/cTが整数となるように半径reの値を選択すると、式(6)とフィルタ係数am, bmは以下のように簡単になる。
以上より、式(6)とそのフィルタ係数、式(6)’とそのフィルタ係数を用いることで、IIRフィルタη0を求めることができる。
【数27】
【表2】
以上より、式(6)とそのフィルタ係数、式(6)’とそのフィルタ係数を用いることで、IIRフィルタη0を求めることができる。
【0032】
＜＜＜２－２：関数Hn(kr)の極の数値計算を利用する方法＞＞＞
ここでは、関数Hn(kr)が実数値関数になることに着目する。図２は、関数H0(kr)の概形を示す図である。図２は周波数を横軸として表したグラフになっている。図２からわかるように、関数H0(kr)の周波数特性は周期的に発散しており、極が支配的になっている。そこで、IIRフィルタを関数Hn(kr)が単位円上に極を持つフィルタとして導出する。そのために、関数Hn(kr)の極の位置を数値計算により求める。IIRフィルタを導出する手順は、以下の通りである。
ここでは、関数Hn(kr)が実数値関数になることに着目する。図２は、関数H0(kr)の概形を示す図である。図２は周波数を横軸として表したグラフになっている。図２からわかるように、関数H0(kr)の周波数特性は周期的に発散しており、極が支配的になっている。そこで、IIRフィルタを関数Hn(kr)が単位円上に極を持つフィルタとして導出する。そのために、関数Hn(kr)の極の位置を数値計算により求める。IIRフィルタを導出する手順は、以下の通りである。
【0033】
【0034】
（２）得られた極となる周波数fwのすべてを配置した伝達関数Hn(z)を求める。極となる周波数fwに対して、z平面上の偏角はθw=2πfw/fs, 座標はvw=r×exp(jθw)（ただし、rは1に近い値であり、例えば、0.999とするとよい。）となる。したがって、Wを得られた極となる周波数の数とすると、伝達関数Hn(z)は次式により与えられる。
ここで、式(7)はvwの複素共役であるvw *も用いていることに留意する。
【数29】
ここで、式(7)はvwの複素共役であるvw *も用いていることに留意する。
【0035】
（３）伝達関数Hn(z)を展開し、IIRフィルタ係数ηn(w)(w=1, …, 2W)を求める。具体的には、次式のようにIIRフィルタ係数ηn(w)(w=1, …, 2W)を用いて表した伝達関数Hn(z)と式(7)を用いてIIRフィルタ係数ηn(w)(w=1, …, 2W)を求めることができる。
上記手順で求めたIIRフィルタの性能を計算機シミュレーションにより確認する。具体的には、導出したIIRフィルタを用いて得られる耳元での波形（以下、推定波形という）と耳元で観測される波形（以下、観測波形という）とを比較する。耳元での音圧^p(i)はIIRフィルタ係数ηn(w)(w=1, …, 2W)を用いた次式により計算される。
ただし、iはサンプル時刻、nは球面調和関数の次数、xn(i)はサンプル時刻iにおけるL個のマイクで観測された信号の音圧の重み付け和、yn(i)はサンプル時刻iにおけるIIRフィルタηnの出力信号である。
【数30】
上記手順で求めたIIRフィルタの性能を計算機シミュレーションにより確認する。具体的には、導出したIIRフィルタを用いて得られる耳元での波形（以下、推定波形という）と耳元で観測される波形（以下、観測波形という）とを比較する。耳元での音圧^p(i)はIIRフィルタ係数ηn(w)(w=1, …, 2W)を用いた次式により計算される。
【数31】
ただし、iはサンプル時刻、nは球面調和関数の次数、xn(i)はサンプル時刻iにおけるL個のマイクで観測された信号の音圧の重み付け和、yn(i)はサンプル時刻iにおけるIIRフィルタηnの出力信号である。
【0036】
音圧推定対象とする騒音は、帯域制限パルス(100Hz～700Hz)と白色雑音(100Hz～700Hz)とする。また、実験条件は以下の通りとする。
【0037】
（１）騒音が左耳の左側から到来するものとする。したがって、剛球の左表面での観測波形を推定することとなる。
【0038】
（２）頭部に相当する剛球の半径を0.10m、12個のマイクを設置する球の半径を0.15mとする。
【0039】
（３）球面調和関数の次数nは2以下、IIRフィルタ長は14とする。なお、IIRフィルタ長は、極の個数が7であることに依存した結果である。
【0040】
上記実験条件のもと、推定波形と観測波形とを取得する。図４は帯域制限パルスに対する推定波形と観測波形を示す図である。図４から、式(8)を用いて計算すると、推定波形と観測波形との間に数サンプル分の時刻のずれが発生すること、推定波形と観測波形の振幅も一致しないことがわかる。このことから、反射騒音による先進の考慮やフィルタゲインの調整が必要なことがわかる。そこで、式(8)を以下のように補正する。
式(8)を補正した式(8)’では、新たにdとAの２つのパラメータを導入している。dは観測値の遅延に対応するパラメータ、Aは観測値の重み付けを調整するパラメータである。パラメータd, Aの決定方法として、例えば、計算機シミュレーションにより推定波形を取得し、ピークの時間差分やピーク値の振幅比を読み取って設定する方法や相互相関を用いて設計する方法がある。
【数32】
式(8)を補正した式(8)’では、新たにdとAの２つのパラメータを導入している。dは観測値の遅延に対応するパラメータ、Aは観測値の重み付けを調整するパラメータである。パラメータd, Aの決定方法として、例えば、計算機シミュレーションにより推定波形を取得し、ピークの時間差分やピーク値の振幅比を読み取って設定する方法や相互相関を用いて設計する方法がある。
【0041】
図５は補正した式(8)’を用いたときの帯域制限パルスに対する推定波形と観測波形を示す図である。図５から、ずれが解消し、推定が安定していることがわかる。また、図６は補正した式(8)’を用いたときの白色雑音に対する推定波形と観測波形を示す図である。図６から、白色雑音に対しても推定が安定していることがわかる。
【0042】
＜第１実施形態＞
音圧推定装置１００は、＜技術的背景＞で説明した通り、頭部近傍の球面上に配置した2以上のマイクを用いて観測される信号の音圧から、耳元の音圧を推定する。Lを頭部近傍の球面上に配置されるマイクの数（ただし、Lは2以上の整数）、p1, …, pLをL個のマイクが配置される点における音圧、(re, Ω1), …, (re, ΩL)（ただし、reは頭部近傍の球面の半径）をL個のマイクが配置される点、(a, Ωa)を耳元の位置を表す点（ただし、aは耳元が位置する球面の半径）とする。
音圧推定装置１００は、＜技術的背景＞で説明した通り、頭部近傍の球面上に配置した2以上のマイクを用いて観測される信号の音圧から、耳元の音圧を推定する。Lを頭部近傍の球面上に配置されるマイクの数（ただし、Lは2以上の整数）、p1, …, pLをL個のマイクが配置される点における音圧、(re, Ω1), …, (re, ΩL)（ただし、reは頭部近傍の球面の半径）をL個のマイクが配置される点、(a, Ωa)を耳元の位置を表す点（ただし、aは耳元が位置する球面の半径）とする。
【0043】
以下、図７～図８を参照して音圧推定装置１００を説明する。図７は、音圧推定装置１００の構成を示すブロック図である。図８は、音圧推定装置１００の動作を示すフローチャートである。図７に示すように音圧推定装置１００は、重み付け加算部１１０と、音圧計算部１２０と、記録部１９０を含む。記録部１９０は、音圧推定装置１００の処理に必要な情報を適宜記録する構成部である。記録部１９０は、例えば、音圧pλに対する次数nの重みγλn(λ=1, …, L, n=0, …, N)を記録しておく。
【0044】
図８に従い音圧推定装置１００の動作について説明する。
【0045】
Ｓ１１０において、重み付け加算部１１０は、音圧p1, …, pLを入力とし、音圧p1, …, pLから、次数0からNまでの音圧p1, …, pLの重み付け和q0, …, qNを計算し、出力する。具体的には、重み付け加算部１１０は、次式により次数nの重み付け和qn(n=0, …, N)を計算する。
（ただし、Yn mは球面調和関数を表す。）
Ｓ１２０において、音圧計算部１２０は、Ｓ１１０で計算した次数0からNまでの重み付け和q0, …, qNを入力とし、重み付け和q0, …, qNから、次数0からNまでのIIRフィルタη0, …, ηNを用いて耳元の音圧^pを計算し、出力する。具体的には、音圧計算部１２０は、次式により音圧^pを計算する。
ここで、次数nのIIRフィルタηn(n=0, …, N)は、次式で表される関数Hn(kr)(n=0, …, N)に基づいて設計された時間領域フィルタである。
（ただし、kは波数であり、jnはn次の球ベッセル関数、hn (2)はn次の第2種球ハンケル関数を表す。）
したがって、関数Hn(kr)(n=0, …, N)は、音圧の内挿関係を表す周波数領域における関数であるといえる。
【数33】
【数34】
（ただし、Yn mは球面調和関数を表す。）
Ｓ１２０において、音圧計算部１２０は、Ｓ１１０で計算した次数0からNまでの重み付け和q0, …, qNを入力とし、重み付け和q0, …, qNから、次数0からNまでのIIRフィルタη0, …, ηNを用いて耳元の音圧^pを計算し、出力する。具体的には、音圧計算部１２０は、次式により音圧^pを計算する。
【数35】
ここで、次数nのIIRフィルタηn(n=0, …, N)は、次式で表される関数Hn(kr)(n=0, …, N)に基づいて設計された時間領域フィルタである。
【数36】
【数37】
（ただし、kは波数であり、jnはn次の球ベッセル関数、hn (2)はn次の第2種球ハンケル関数を表す。）
したがって、関数Hn(kr)(n=0, …, N)は、音圧の内挿関係を表す周波数領域における関数であるといえる。
【0046】
具体的な設計例を以下に示す。
（例１）
次数0のIIRフィルタη0は、関数H0(kr)のs領域上の表現である関数H0(s)をs-z変換することにより得られる次式の関数H0(z)の係数am(m=1, 2, 3, 2u, 2u+1, 2u+2, 2u+3、ただし、u=floor((re-a)/c×1/T)), bm(m=0, 1, 2, 3)を用いて設計された時間領域フィルタである。なお、uは2以上の整数であることが好ましい。
（例２）
次数nのIIRフィルタηn(n=0, …, N)は、関数Hn(kr)の極となる周波数fw(w=1, …, W、Wは関数Hn(kr)の極の数)を用いて次式により表される関数Hn(z)に基づいて得られるIIRフィルタ係数ηn(w) (w=1, …, 2W)を用いて設計された時間領域フィルタである。
（ただし、vw=r×exp(jθw), θw=2πfw/fsであり、rは1に略等しく、fsはサンプリング周波数を表す。）
本発明の実施形態によれば、周波数領域における演算を用いることなく（換言すると、時間領域における演算のみにより）、耳元の音圧を推定することが可能となる。これにより、即時性の高い処理が必要となる能動騒音制御システムに適用することが可能となる。
（例１）
次数0のIIRフィルタη0は、関数H0(kr)のs領域上の表現である関数H0(s)をs-z変換することにより得られる次式の関数H0(z)の係数am(m=1, 2, 3, 2u, 2u+1, 2u+2, 2u+3、ただし、u=floor((re-a)/c×1/T)), bm(m=0, 1, 2, 3)を用いて設計された時間領域フィルタである。なお、uは2以上の整数であることが好ましい。
【数38】
（例２）
次数nのIIRフィルタηn(n=0, …, N)は、関数Hn(kr)の極となる周波数fw(w=1, …, W、Wは関数Hn(kr)の極の数)を用いて次式により表される関数Hn(z)に基づいて得られるIIRフィルタ係数ηn(w) (w=1, …, 2W)を用いて設計された時間領域フィルタである。
【数39】
（ただし、vw=r×exp(jθw), θw=2πfw/fsであり、rは1に略等しく、fsはサンプリング周波数を表す。）
本発明の実施形態によれば、周波数領域における演算を用いることなく（換言すると、時間領域における演算のみにより）、耳元の音圧を推定することが可能となる。これにより、即時性の高い処理が必要となる能動騒音制御システムに適用することが可能となる。
【0047】
＜補記＞
上述した各装置の各部の処理をコンピュータにより実現してもよく、この場合は各装置が有すべき機能の処理内容はプログラムによって記述される。そして、このプログラムを図９に示すコンピュータ２０００の記録部２０２０に読み込ませ、演算処理部２０１０、入力部２０３０、出力部２０４０、補助記録部２０２５などを動作させることにより、上記各装置における処理機能がコンピュータ上で実現される。
上述した各装置の各部の処理をコンピュータにより実現してもよく、この場合は各装置が有すべき機能の処理内容はプログラムによって記述される。そして、このプログラムを図９に示すコンピュータ２０００の記録部２０２０に読み込ませ、演算処理部２０１０、入力部２０３０、出力部２０４０、補助記録部２０２５などを動作させることにより、上記各装置における処理機能がコンピュータ上で実現される。
【0048】
本発明の装置は、例えば単一のハードウェアエンティティとして、ハードウェアエンティティの外部から信号を入力可能な入力部、ハードウェアエンティティの外部に信号を出力可能な出力部、ハードウェアエンティティの外部に通信可能な通信装置（例えば通信ケーブル）が接続可能な通信部、演算処理部であるCPU（Central Processing Unit、キャッシュメモリやレジスタなどを備えていてもよい）、メモリであるRAMやROM、ハードディスクである外部記憶装置並びにこれらの入力部、出力部、通信部、CPU、RAM、ROM、外部記憶装置の間のデータのやり取りが可能なように接続するバスを有している。また必要に応じて、ハードウェアエンティティに、CD-ROMなどの記録媒体を読み書きできる装置（ドライブ）などを設けることとしてもよい。このようなハードウェア資源を備えた物理的実体としては、汎用コンピュータなどがある。
【0049】
ハードウェアエンティティの外部記憶装置には、上述の機能を実現するために必要となるプログラムおよびこのプログラムの処理において必要となるデータなどが記憶されている（外部記憶装置に限らず、例えばプログラムを読み出し専用記憶装置であるROMに記憶させておくこととしてもよい）。また、これらのプログラムの処理によって得られるデータなどは、RAMや外部記憶装置などに適宜に記憶される。
【0050】
ハードウェアエンティティでは、外部記憶装置（あるいはROMなど）に記憶された各プログラムとこの各プログラムの処理に必要なデータが必要に応じてメモリに読み込まれて、適宜にCPUで解釈実行、処理される。その結果、CPUが所定の機能（上記、…部、…手段などと表した各構成部）を実現する。つまり、本発明の実施形態の各構成部は、処理回路(Processing Circuitry)により構成されてもよい。
【0051】
既述のように、上記実施形態において説明したハードウェアエンティティ（本発明の装置）における処理機能をコンピュータによって実現する場合、ハードウェアエンティティが有すべき機能の処理内容はプログラムによって記述される。そして、このプログラムをコンピュータで実行することにより、上記ハードウェアエンティティにおける処理機能がコンピュータ上で実現される。
【0052】
この処理内容を記述したプログラムは、コンピュータで読み取り可能な記録媒体に記録しておくことができる。コンピュータで読み取り可能な記録媒体は、例えば、非一時的な記録媒体であり、具体的には、磁気記録装置、光ディスク等である。
【0053】
また、このプログラムの流通は、例えば、そのプログラムを記録したDVD、CD-ROM等の可搬型記録媒体を販売、譲渡、貸与等することによって行う。さらに、このプログラムをサーバコンピュータの記憶装置に格納しておき、ネットワークを介して、サーバコンピュータから他のコンピュータにそのプログラムを転送することにより、このプログラムを流通させる構成としてもよい。
【0054】
このようなプログラムを実行するコンピュータは、例えば、まず、可搬型記録媒体に記録されたプログラムもしくはサーバコンピュータから転送されたプログラムを、一旦、自己の非一時的な記憶装置である補助記録部２０２５に格納する。そして、処理の実行時、このコンピュータは、自己の非一時的な記憶装置である補助記録部２０２５に格納されたプログラムを記録部２０２０に読み込み、読み込んだプログラムに従った処理を実行する。また、このプログラムの別の実行形態として、コンピュータが可搬型記録媒体から直接プログラムを記録部２０２０に読み込み、そのプログラムに従った処理を実行することとしてもよく、さらに、このコンピュータにサーバコンピュータからプログラムが転送されるたびに、逐次、受け取ったプログラムに従った処理を実行することとしてもよい。また、サーバコンピュータから、このコンピュータへのプログラムの転送は行わず、その実行指示と結果取得のみによって処理機能を実現する、いわゆるASP（Application Service Provider）型のサービスによって、上述の処理を実行する構成としてもよい。なお、本形態におけるプログラムには、電子計算機による処理の用に供する情報であってプログラムに準ずるもの（コンピュータに対する直接の指令ではないがコンピュータの処理を規定する性質を有するデータ等）を含むものとする。
【0055】
また、この形態では、コンピュータ上で所定のプログラムを実行させることにより、本装置を構成することとしたが、これらの処理内容の少なくとも一部をハードウェア的に実現することとしてもよい。
【0056】
本発明は上述の実施形態に限定されるものではなく、本発明の趣旨を逸脱しない範囲で適宜変更が可能である。
【図1】
【図2】
【図3】
【図4】
【図5】
【図6】
【図7】
【図8】
【図9】