JP7792263B2

JP7792263B2 - Sound field reproduction device and program

Info

Publication number: JP7792263B2
Application number: JP2022024213A
Authority: JP
Inventors: 陽佐々木
Original assignee: Japan Broadcasting Corp
Current assignee: Japan Broadcasting Corp
Priority date: 2022-02-18
Filing date: 2022-02-18
Publication date: 2025-12-25
Anticipated expiration: 2042-02-18
Also published as: JP2023121021A

Description

本発明は、スピーカアレイを用いて音場を再現する音場再現装置及びプログラムに関する。 The present invention relates to a sound field reproduction device and program that reproduces a sound field using a speaker array.

従来、多数のスピーカを用いて、ある空間に任意の音場を形成する音場再現技術に関する研究が進められている。音場再現の主要な技術としては、波面合成法（ＷＦＳ：Wave Field Synthesis）、境界音場制御（ＢｏＳＣ：Boundary Surface Control）、高次アンビソニックス（ＨＯＡ：Higher Order Ambisonics）、スペクトル除算法（ＳＤＭ：Spectral Division Method）等が知られている。 Research has been ongoing into sound field reproduction technologies that use multiple speakers to create an arbitrary sound field in a given space. Key sound field reproduction technologies include Wave Field Synthesis (WFS), Boundary Surface Control (BoSC), Higher Order Ambisonics (HOA), and Spectral Division Method (SDM).

図７は、波面合成法（ＷＦＳ）におけるスピーカアレイの配置例を示す図である。ＷＦＳは、レイリー積分に基づいた手法であり、直線状または平面状に配置されたスピーカアレイ１００－１を用いて、所望音場の境界平面上の音圧または音圧勾配を再現することで、境界外部から到来する音源による波面を形成するものである。 Figure 7 shows an example of speaker array placement in wave field synthesis (WFS). WFS is a method based on the Rayleigh integral, and uses a speaker array 100-1 arranged linearly or planarly to reproduce the sound pressure or sound pressure gradient on the boundary plane of the desired sound field, thereby forming a wavefront due to a sound source arriving from outside the boundary.

図８は、境界音場制御（ＢｏＳＣ）におけるスピーカアレイの配置例を示す図である。ＢｏＳＣは、キルヒホッフ・ヘルムホルツ積分方程式に基づいた手法である。ＢｏＳＣは、音場を再現したい領域の外部に、領域内部に向けて聴取者を囲むようにスピーカアレイ１００－２を配置し、逆システムを用いてその領域の境界面上の音圧及び音圧勾配を制御することで、領域内部に所望音場を再現するものである。 Figure 8 shows an example of speaker array placement in boundary sound field control (BoSC). BoSC is a method based on the Kirchhoff-Helmholtz integral equation. BoSC places speaker array 100-2 outside the area where you want to reproduce the sound field, surrounding the listener toward the inside of the area, and uses an inverse system to control the sound pressure and sound pressure gradient on the boundary surface of that area, thereby reproducing the desired sound field inside the area.

図９は、高次アンビソニックス（ＨＯＡ）におけるスピーカアレイの配置例を示す図である。ＨＯＡは、球面調和関数を用いて所望音場を表現し、球面内部に向けて聴取者を囲むように配置した球面のスピーカアレイ１００－３を用いて、再現音場の球面調和係数を所望音場の球面調和係数とマッチングさせることで、音場を再現する手法である。 Figure 9 shows an example of a speaker array arrangement in Higher Order Ambisonics (HOA). HOA is a method of reproducing a sound field by expressing the desired sound field using spherical harmonic functions and using a spherical speaker array 100-3 arranged to surround the listener facing the inside of the sphere, matching the spherical harmonic coefficients of the reproduced sound field with the spherical harmonic coefficients of the desired sound field.

図１０は、スペクトル除算法（ＳＤＭ）におけるスピーカアレイの配置例を示す図である。ＳＤＭは、角度スペクトルを用いて所望音場を表現し、直線状または平面状に配置されたスピーカアレイ１００－４を用いて、再現音場の角度スペクトルを所望音場の角度スペクトルとマッチングさせることで、音場を再現する手法である。ＳＤＭの詳細については、例えば非特許文献１を参照されたい。 Figure 10 shows an example of speaker array placement in the spectral division method (SDM). SDM is a technique for expressing a desired sound field using an angular spectrum, and reproducing the sound field by matching the angular spectrum of the reproduced sound field with the angular spectrum of the desired sound field using a speaker array 100-4 arranged in a line or plane. For details on SDM, see, for example, Non-Patent Document 1.

〔ＳＤＭにおける駆動信号の算出方法〕
以下、一般的なＳＤＭにおける駆動信号の算出方法について説明する。図１１は、一般的なＳＤＭを説明する図である。 [Method of calculating drive signal in SDM]
A method for calculating a drive signal in a typical SDM will be described below with reference to FIG.

ｘｙｚ空間上のｙ＝ｙ₀，ｚ＝ｚ₀に配置された無限直線音源１０１上の座標ｒ0＝（ｘ₀, ｙ₀，ｚ₀）の点において、角周波数ωの駆動信号をＤ（ｒ0，ω）とする。この場合の座標ｒ＝（ｘ，ｙ_ref，ｚ₀）の点における音圧Ｐ（ｒ，ω）は、以下の式にて表される。
Let D(r0,ω) be the drive signal with angular frequency ω at the point of coordinate r0 = ( _x0 , _y0 , _z0 ) on the infinite linear sound source 101 placed at y = _y0 , z = _z0 in the xyz space. In this case, the sound pressure P(r,ω) at the point of coordinate r = (x, _yref , _z0 ) is expressed by the following equation.

ここで、Ｇ（ｒ－ｒ0，ω）は、ベクトル（ｒ－ｒ0）で表現される方向及び距離の伝達関数であり、再現音場として自由音場を想定した場合、３次元の自由音場グリーン関数で表現され、以下の式にて表される。
ｋは波数である。 Here, G(r-r0, ω) is a transfer function of direction and distance expressed by the vector (r-r0), and when a free sound field is assumed as the reproduced sound field, it is expressed by a three-dimensional free sound field Green's function and is expressed by the following equation.
k is the wave number.

音圧Ｐ（ｒ，ω）は、駆動信号Ｄ（ｒ0，ω）及び伝達関数Ｇ（ｒ－ｒ0，ω）の空間上の畳み込みによって得られるものと解釈でき、ｙ＝ｙ_refにおいてｘ軸の方向に沿ってフーリエ変換することにより、以下の角度スペクトル表現Ｐ^（ｋ_x，ｙ_ref，ｚ₀，ω）が得られる。
The sound pressure P(r,ω) can be interpreted as being obtained by spatially convolving the drive signal D(r0,ω) and the transfer function G(r-r0,ω), and by performing a Fourier transform along the x-axis at y = y _ref , the following angular spectrum representation P^(k _x , y _ref , z ₀ ,ω) is obtained.

Ｐ^（ｋ_x，ｙ_ref，ｚ₀，ω）は、音圧Ｐ（ｒ，ω）をｙ＝ｙ_refで空間フーリエ変換することにより得られる角度スペクトルである。また、Ｄ^（ｋ_x，ｙ₀，ｚ₀，ω）は、駆動信号Ｄ（ｒ0，ω）をｙ＝ｙ₀で空間フーリエ変換することにより得られる角度スペクトルである。さらに、Ｇ^（ｋ_x，ｙ_ref－ｙ₀，ｚ₀，ω）は、伝達関数Ｇ（ｒ－ｒ0，ω）をｙ＝ｙ_refで空間フーリエ変換することにより得られる角度スペクトルである。ｋ_xはｘ軸方向の波数である。 P^( _kx , _yref , _z0 , ω) is the angular spectrum obtained by spatial Fourier transforming the sound pressure P(r, ω) with y = _yref . Furthermore, D^( _kx , _y0 , _z0 , ω) is the angular spectrum obtained by spatial Fourier transforming the drive signal D(r0, ω) with y = _y0 . Furthermore, G^( _kx , _yref - _y0 , _z0 , ω) is the angular spectrum obtained by spatial Fourier transforming the transfer function G(r - r0, ω) with y = _yref . _kx is the wave number in the x-axis direction.

ここで、波面を再現したい所望音場のｙ＝ｙ_des（以下、所望境界という。）における音圧分布を空間フーリエ変換することにより得られる角度スペクトルをＰ^_d（ｋ_x，ｙ_des，ｚ₀，ω）とすると、所望音場の波面を再現するための駆動信号の角度スペクトルＤ^（ｋ_x，ｙ₀，ｚ₀，ω）は、以下の式にて表される。
再現境界ｙ_refは、スピーカアレイを設置した、実際に音場を再現する空間での音圧分布を再現する境界を示す。所望境界ｙ_desは、再現境界ｙ_refで再現したい実音源による音圧分布の境界を示す。 Here, if the angular spectrum obtained by spatial Fourier transforming the sound pressure distribution at y = y _des (hereinafter referred to as the desired boundary) of the desired sound field whose wavefront you want to reproduce is P^ _d ( _kx , _ydes , _z0 , ω), the angular spectrum D^( _kx , _y0 , _z0 , ω) of the drive signal for reproducing the wavefront of the desired sound field is expressed by the following equation.
The reproduction boundary _yref indicates the boundary for reproducing the sound pressure distribution in a space where a speaker array is installed and a sound field is actually reproduced. The desired boundary _ydes indicates the boundary of the sound pressure distribution due to the actual sound source that is desired to be reproduced at the reproduction boundary _yref .

座標ｒ_s＝（ｘ_s，ｙ_s，ｚ₀）に配置された点音源が角周波数ωにおいてＳ（ω）で駆動された場合、所望音場のｙ＝ｙ_desでの角度スペクトルＰ^_d（ｋ_x，ｙ_des，ｚ₀，ω）、すなわち所望音場角度スペクトルは、以下の式にて表される。
When a point sound source located at coordinates r _s = (x _s , y _s , z ₀ ) is driven by S(ω) at angular frequency ω, the angular spectrum P^ _d (k _x , y _des , z ₀ , ω) of the desired sound field at y = y _des , i.e., the desired sound field angular spectrum, is expressed by the following equation.

ここで、Ｈ₀ ⁽²⁾は０次第２種ハンケル関数であり、Ｋ₀は０次変形ベッセル関数である。また、０≦ｋ²－ｋ_x ²の場合の音波は、遠方まで伝搬する音波を示し、０＞ｋ²－ｋ_x ²の場合の音波はエバネッセント波と呼ばれ、振幅が急速に減衰する遠方まで伝搬しない音波を示す。 Here, H ₀ ⁽²⁾ is the zero-order Hankel function of the second kind, and K ₀ is the zeroth-order modified Bessel function. When 0≦k ² −k _x ² , the sound wave propagates over long distances, while when 0>k ² −k _x ² , the sound wave is called an evanescent wave, and the amplitude rapidly decays, so the sound wave does not propagate over long distances.

同様に、角度スペクトルＧ^（ｋ_x，ｙ_ref－ｙ₀，ｚ₀，ω）、すなわち再現音場角度スペクトルは、以下の式にて表される。
ただし、０＞ｋ²－ｋ_x ²のときはいずれの場合も、０次変形ベッセル関数Ｋ₀は、音源近傍のみに伝搬するエバネッセント波を表しており、エバネッセント波の再現は数値的に不安定になり易いため、一般には省かれることが多い。 Similarly, the angular spectrum G^(k _x , y _ref −y ₀ , z ₀ , ω), that is, the reproduced sound field angular spectrum, is expressed by the following equation.
However, in any case where 0>k ² −k _x ² , the zeroth-order modified Bessel function K ₀ represents an evanescent wave that propagates only near the sound source, and since the reproduction of evanescent waves tends to be numerically unstable, it is generally omitted.

したがって、駆動信号の角度スペクトルＤ^（ｋ_x，ｙ₀，ｚ₀，ω）は、以下の式にて表す。
以下、特に断りなく０＞ｋ²－ｋ_x ²の場合はＤ^（ｋ_x，ｙ₀，ｚ₀，ω）＝０とし、０≦ｋ²－ｋ_x ²の場合についてのみ説明するものとする。 Therefore, the angular spectrum D^(k _x , y ₀ , z ₀ , ω) of the drive signal is expressed by the following equation.
Unless otherwise specified, in the following, when 0>k ² -k _x ² , it is assumed that D^(k _x , y ₀ , z ₀ , ω)=0, and only the case when 0≦k ² -k _x ² will be explained.

図１２は、再現領域及び再現境界を説明する図である。直線状に配置されたスピーカアレイ１００が無限直線音源である場合、当該無限直線音源を、前記式（７）の駆動信号を用いて駆動することにより、ｙ≧ｙ_refの領域において音場を再現することができる。このため、ここではｙ＝ｙ_refを再現境界、ｙ≧ｙ_refの領域を再現領域と呼ぶこととする。 12 is a diagram illustrating the reproduction area and reproduction boundary. When the speaker array 100 arranged in a line is an infinite linear sound source, a sound field can be reproduced in the area y≧ _yref by driving the infinite linear sound source using the drive signal of the above-mentioned equation (7). For this reason, y= _yref will be called the reproduction boundary, and the area y≧ _yref will be called the reproduction area.

前記式（７）の駆動信号の角度スペクトルＤ^（ｋ_x，ｙ₀，ｚ₀，ω）は、無限直線音源を用いて波面合成を行うことを想定しているが、実際には、直線状に配置した有限個のスピーカユニットからなるスピーカアレイ１００を用いて実現する。 The angular spectrum D^( _kx , _y0 , _z0 , ω) of the drive signal in equation (7) is assumed to be produced by wave field synthesis using an infinite linear sound source, but in reality it is realized using a speaker array 100 consisting of a finite number of speaker units arranged in a line.

このため、前記式（７）の駆動信号の角度スペクトルＤ^（ｋ_x，ｙ₀，ｚ₀，ω）において、周波数、角度スペクトルを離散化すると共に、有限長で打ち切り、さらに離散逆フーリエ変換（ＩＤＦＴ：Inverse Discrete Fourier Transform）を行うことで、スピーカユニット毎の時間周波数領域の駆動信号を算出する必要がある。 For this reason, in the angular spectrum D^( _kx , _y0 , _z0 , ω) of the drive signal in equation (7), it is necessary to discretize the frequency and angular spectrum, truncate it to a finite length, and then perform an inverse discrete Fourier transform (IDFT) to calculate the drive signal in the time-frequency domain for each speaker unit.

スピーカアレイ１００を構成するスピーカユニット数をＭ（自然数）、スピーカユニットの間隔をΔｘ、スピーカアレイ長をＬ（＝（Ｍ－１）Δｘ）とする。前記式（７）の駆動信号の角度スペクトルＤ^（ｋ_x，ｙ₀，ｚ₀，ω）に対し、離散逆フーリエ変換を行うことで、ｘに配置されたスピーカユニットの時間周波数領域の駆動信号Ｄ~（ｘ，ω）は、以下の式で表される。
ここで、ｍは、０～Ｍ－１（ｍは整数）であり、ｋ^m _x＝２πｍ／（ΔｘＭ）は、ｘ軸方向の波数である。尚、スピーカアレイ１００は、ｘｙｚ空間において、予め設定されたｙ＝ｙ₀，ｚ＝ｚ₀のｘ軸上に配置される。 Let M (a natural number) be the number of speaker units constituting the speaker array 100, Δx be the spacing between the speaker units, and L (=(M-1)Δx) be the length of the speaker array. By performing an inverse discrete Fourier transform on the angular spectrum D^(k _x , y ₀ , z ₀ , ω) of the drive signal in the time-frequency domain of the speaker unit located at x, the drive signal D~(x, ω) can be expressed by the following equation:
Here, m is an integer between 0 and M-1, and k ^m _x = 2πm/(ΔxM) is the wave number in the x-axis direction. The speaker array 100 is placed on the x-axis at a preset y=y ₀ and z=z ₀ in the xyz space.

そして、スピーカアレイ１００を構成するスピーカユニット毎に、時間領域の駆動信号が算出される。 Then, a time-domain drive signal is calculated for each speaker unit that makes up the speaker array 100.

このようなＳＤＭ等の音場再現技術を用いることにより、ユーザが指定する座標に配置される点音源によって形成される音場を再現することができる。 By using sound field reproduction technologies such as SDM, it is possible to reproduce a sound field formed by a point sound source placed at coordinates specified by the user.

従来のスピーカ駆動信号算出方法によれば、波面を再現したい所望音場において、音源が自由空間に設置してあることが想定されており、壁等の遮蔽物の存在は考慮されていない。 Conventional speaker drive signal calculation methods assume that the sound source is installed in free space in the desired sound field where the wavefront is to be reproduced, and do not take into account the presence of obstructions such as walls.

J.Ahrens, and S.Spors,“Sound field reproduction using planar and linear arrays of loudspeakers”, IEEE Trans. Audio, Speech, Lang. Process., vol.18, no.8, pp.2038-2050, Nov.2010.J. Ahrens, and S. Spors, “Sound field reproduction using planar and linear arrays of loudspeakers”, IEEE Trans. Audio, Speech, Lang. Process., vol.18, no.8, pp.2038-2050, Nov.2010.

しかしながら、ＳＤＭを用いた音響コンテンツの制作においては、壁等の遮蔽物によるオクルージョンまたは音色の変化を考慮したいという演出も考えられる。 However, when producing audio content using SDM, it may be necessary to take into account occlusions caused by walls or other obstructions, or changes in tone.

そこで、本発明は前記課題を解決するためになされたものであり、その目的は、ＳＤＭによりスピーカアレイの駆動信号を生成する際に、所望音場の音源と所望境界との間に開口部を有する壁を設置した場合において、当該壁の影響を考慮した波面を再現可能な音場再現装置及びプログラムを提供することにある。 The present invention was conceived to solve the above-mentioned problems, and its purpose is to provide a sound field reproduction device and program that can reproduce a wavefront that takes into account the influence of a wall with an opening placed between the sound source and the desired boundary of the desired sound field when generating a drive signal for a speaker array using SDM.

前記課題を解決するために、請求項１の音場再現装置は、複数のスピーカユニットからなるスピーカアレイを用いて、所望音源により形成される音場を再現するための駆動信号を生成する音場再現装置において、前記スピーカアレイが、ｘｙｚ空間上のｚ＝ｚ₀のｙ＝ｙ₀上に配置されているものとし、ｋを波数、ｋ_xをｘ軸方向の波数、Δｋ_xを前記波数ｋ_xの離散化間隔、ωを角周波数、ｃを音速、（ｘ_s，ｙ_s）を前記所望音源の座標である仮想音源座標、ｓ（ｎ）を時間サンプルｎ毎の前記所望音源の音源信号、ｙ_wを壁位置、Ｈ（ｋ_x）を前記壁位置ｙ_wのｙ軸方向粒子速度透過率分布角度スペクトル、ｙ_desを所望境界、ｙ_refを再現境界、Ｍをスピーカユニット数、Δｘをスピーカユニット間隔、ｆ_sをサンプリング周波数、ＦをＤＦＴ点数、Ｈ₀ ⁽²⁾を０次第２種ハンケル関数、ｉを虚数単位として、予め設定された前記スピーカユニット数Ｍ及び前記スピーカユニット間隔Δｘ、並びにインデックスｍ＝０～Ｍ－１（ｍは整数）を用いて、式：ｋ_x＝２πｍ／（ΔｘＭ）により、前記波数ｋ_x（＝ｋ¹ _x，ｋ² _x，・・・，ｋ^M _x）を算出し、前記スピーカアレイが等間隔の前記複数のスピーカユニットにより構成された場合の前記離散化間隔Δｋ_xを算出するｘ軸方向波数算出部と、予め設定された前記サンプリング周波数ｆ_s及び前記ＤＦＴ点数Ｆ、並びに周波数インデックスｌ＝０～Ｆ－１（ｌは整数）を用いて、式：ω＝２πｆ_sｌ／Ｆにより、前記角周波数ωを算出し、前記角周波数ω及び音速ｃを用いて、式：ｋ＝ω／ｃにより、前記波数ｋを算出する波数算出部と、前記ｘ軸方向波数算出部により算出された前記波数ｋ_x、前記波数算出部により算出された前記波数ｋ、及び予め設定された前記再現境界ｙ_ref及び前記値ｙ₀を用いて、式：
により、再現音場角度スペクトルを算出する再現音場角度スペクトル算出部と、所望音場角度スペクトルを算出する所望音場角度スペクトル算出部と、前記所望音場角度スペクトル算出部により算出された前記所望音場角度スペクトルを、前記再現音場角度スペクトル算出部により算出された前記再現音場角度スペクトルで除算することで、前記駆動信号の角度スペクトルを求める除算部と、を備え、前記所望音場角度スペクトル算出部が、所定バッファ長のバッファリング区間における前記仮想音源座標（ｘ_s，ｙ_s）及び前記音源信号ｓ（ｎ）を用いて、前記波数ｋ_x毎の前記壁位置ｙ_wの音圧角度スペクトルＰ^_d（ｋ_x，ｙ_w）を算出し、前記波数算出部により算出された前記波数ｋ、前記ｘ軸方向波数算出部により算出された前記波数ｋ_x、及び前記音圧角度スペクトルＰ^_d（ｋ_x，ｙ_w）を用いて、式：
により、ベクトルＰ^を算出し、前記波数ｋ、前記波数ｋ_x、予め設定された前記ｙ軸方向粒子速度透過率分布角度スペクトルＨ（ｋ_x）、前記所望境界ｙ_des及び前記壁位置ｙ_wを用いて、式：
（ｋⁿ _x及びｋ^m _xにおいてｎ，ｍ＝１～Ｍ（ｎ，ｍは整数））により、角度スペクトルＨ￣（ｋⁿ _x，ｋ^m _x）を算出し、式：
により、行列Ｈを求め、前記ｘ軸方向波数算出部により算出された前記離散化間隔Δｋ_x、前記ベクトルＰ^及び前記行列Ｈを乗算することで、前記所望音場角度スペクトルを求める、ことを特徴とする。 In order to solve the above problem, the sound field reproduction device of claim 1 is a sound field reproduction device that generates a drive signal for reproducing a sound field formed by a desired sound source using a speaker array consisting of a plurality of speaker units, wherein the speaker array is arranged on y= _y0 of z= _z0 in xyz space, k is the wave number, _kx is the wave number in the x-axis direction, _Δkx is the discretization interval of the wave number _kx , ω is the angular frequency, c is the speed of sound, ( _xs , _ys ) are virtual sound source coordinates that are the coordinates of the desired sound source, s(n) is the sound source signal of the desired sound source for each time sample n, _yw is the wall position, H( _kx ) is the particle velocity transmittance distribution angular spectrum in the y-axis direction at the wall position _yw , _ydes is the desired boundary, _yref is the reproduction boundary, M is the number of speaker units, Δx is the speaker unit interval, _fs is the sampling frequency, F is the number of DFT points, _H0 an x-axis direction wave number calculation unit that calculates the wave number k x (= k 1 x , k ² _x , _... , k M x ) by the formula k _x = 2πm/(ΔxM) using the number M of speaker units and the speaker unit interval Δx, and an index m = ⁰ to ^M _- ₁ (m is an integer), where ( ² ) is a zero-order second kind Hankel function, i is an imaginary unit, and calculates the discretization interval Δk _x when the speaker array is composed of the plurality of equally spaced speaker units; a wave number calculation unit that calculates the angular frequency ω by the formula ω = 2πf _s l/F using the preset sampling frequency f _s , the number of DFT points F, and a frequency index l = 0 to F-1 (l is an integer), and calculates the wave number k by the formula k = ω/c using the angular frequency ω and the sound speed _c ; Using the wave number k calculated by the wave number calculation unit, and the preset reproduction boundary y _ref and value y ₀ , the following equation is calculated:
_a desired sound field angular spectrum calculation unit _that calculates a desired _sound _field _angular _spectrum _according _to _the _formula _:
A vector P^ is calculated by the following equation using the wave number k, the wave number _kx , the preset y-axis direction particle velocity transmittance distribution angular spectrum H( _kx ), the desired boundary _ydes , and the wall position _yw :
(where n and m = 1 to M (n and m are integers) in k ⁿ _x and ^km _x ), the angular spectrum H (k ⁿ _x , ^km _x ) is calculated using the formula:
and multiplying the discretization interval Δk _x calculated by the x-axis direction wave number calculation unit, the vector P^, and the matrix H to calculate the desired sound field angular spectrum.

また、請求項２のプログラムは、複数のスピーカユニットからなるスピーカアレイを用いて、所望音源により形成される音場を再現するための駆動信号を生成する音場再現装置を構成するコンピュータを、前記スピーカアレイが、ｘｙｚ空間上のｚ＝ｚ₀のｙ＝ｙ₀上に配置されているものとし、ｋを波数、ｋ_xをｘ軸方向の波数、Δｋ_xを前記波数ｋ_xの離散化間隔、ωを角周波数、ｃを音速、（ｘ_s，ｙ_s）を前記所望音源の座標である仮想音源座標、ｓ（ｎ）を時間サンプルｎ毎の前記所望音源の音源信号、ｙ_wを壁位置、Ｈ（ｋ_x）を前記壁位置ｙ_wのｙ軸方向粒子速度透過率分布角度スペクトル、ｙ_desを所望境界、ｙ_refを再現境界、Ｍをスピーカユニット数、Δｘをスピーカユニット間隔、ｆ_sをサンプリング周波数、ＦをＤＦＴ点数、Ｈ₀ ⁽²⁾を０次第２種ハンケル関数、ｉを虚数単位として、予め設定された前記スピーカユニット数Ｍ及び前記スピーカユニット間隔Δｘ、並びにインデックスｍ＝０～Ｍ－１（ｍは整数）を用いて、式：ｋ_x＝２πｍ／（ΔｘＭ）により、前記波数ｋ_x（＝ｋ¹ _x，ｋ² _x，・・・，ｋ^M _x）を算出し、前記スピーカアレイが等間隔の前記複数のスピーカユニットにより構成された場合の前記離散化間隔Δｋ_xを算出するｘ軸方向波数算出部、予め設定された前記サンプリング周波数ｆ_s及び前記ＤＦＴ点数Ｆ、並びに周波数インデックスｌ＝０～Ｆ－１（ｌは整数）を用いて、式：ω＝２πｆ_sｌ／Ｆにより、前記角周波数ωを算出し、前記角周波数ω及び音速ｃを用いて、式：ｋ＝ω／ｃにより、前記波数ｋを算出する波数算出部、前記ｘ軸方向波数算出部により算出された前記波数ｋ_x、前記波数算出部により算出された前記波数ｋ、及び予め設定された前記再現境界ｙ_ref及び前記値ｙ₀を用いて、式：
により、再現音場角度スペクトルを算出する再現音場角度スペクトル算出部、所望音場角度スペクトルを算出する所望音場角度スペクトル算出部、及び、前記所望音場角度スペクトル算出部により算出された前記所望音場角度スペクトルを、前記再現音場角度スペクトル算出部により算出された前記再現音場角度スペクトルで除算することで、前記駆動信号の角度スペクトルを求める除算部として機能させるためのプログラムであって、前記所望音場角度スペクトル算出部が、所定バッファ長のバッファリング区間における前記仮想音源座標（ｘ_s，ｙ_s）及び前記音源信号ｓ（ｎ）を用いて、前記波数ｋ_x毎の前記壁位置ｙ_wの音圧角度スペクトルＰ^_d（ｋ_x，ｙ_w）を算出し、前記波数算出部により算出された前記波数ｋ、前記ｘ軸方向波数算出部により算出された前記波数ｋ_x、及び前記音圧角度スペクトルＰ^_d（ｋ_x，ｙ_w）を用いて、式：
により、ベクトルＰ^を算出し、前記波数ｋ、前記波数ｋ_x、予め設定された前記ｙ軸方向粒子速度透過率分布角度スペクトルＨ（ｋ_x）、前記所望境界ｙ_des及び前記壁位置ｙ_wを用いて、式：
（ｋⁿ _x及びｋ^m _xにおいてｎ，ｍ＝１～Ｍ（ｎ，ｍは整数））により、角度スペクトルＨ￣（ｋⁿ _x，ｋ^m _x）を算出し、式：
により、行列Ｈを求め、前記ｘ軸方向波数算出部により算出された前記離散化間隔Δｋ_x、前記ベクトルＰ^及び前記行列Ｈを乗算することで、前記所望音場角度スペクトルを求める、ことを特徴とする。 The program of claim 2 is a computer program that configures a sound field reproduction device that generates a drive signal for reproducing a sound field formed by a desired sound source using a speaker array consisting of a plurality of speaker units, the speaker array being arranged on y= _y0 of z= _z0 in an xyz space, and that configures the computer program as follows: k is the wave number, _kx is the wave number in the x-axis direction, Δkx is the discretization interval of the _wave number _kx , ω is the angular frequency, c is the speed of sound, ( _xs , _ys ) are virtual sound source coordinates that are the coordinates of the desired sound source, s(n) is the sound source signal of the desired sound source for each time sample n, _yw is the wall position, H( _kx ) is the particle velocity transmittance distribution angular spectrum in the y-axis direction at the wall position _yw , _ydes is the desired boundary, _yref is the reproduction boundary, M is the number of speaker units, Δx is the speaker unit interval, _fs is the sampling frequency, F is the number of DFT points, _H0 an x-axis direction wave number calculation unit that calculates the wave number k x (= k 1 x , k ² _x , ..., _k M x ) by the formula k _x = 2πm/(ΔxM) using the preset number of speaker units M and the speaker unit interval Δx, and an index m = ⁰ to M- ¹ ( _m is an _integer ), where ( ² ) is a zero-order second kind Hankel function, i is an imaginary unit, and calculates the discretization interval Δk _x when the speaker array is composed of the plurality of equally spaced speaker units; a wave number calculation unit that calculates the angular frequency ω by the formula ω = 2πf _s l/F using the preset sampling frequency f _s and the number of DFT points F, and a frequency index l = 0 to F-1 (l is an integer), and calculates the wave number k by the formula k = ω/c using the angular frequency ω and the sound speed _c; Using the wave number k calculated by the wave number calculation unit, and the preset reproduction boundary y _ref and value y ₀ , the following equation is calculated:
a reproduction sound field angular spectrum calculation unit that calculates a reproduction sound field angular spectrum, a desired sound field angular spectrum calculation unit that calculates a desired sound field angular spectrum, and a division unit that determines the angular spectrum of the drive signal by dividing the desired sound field angular spectrum calculated by the desired sound field angular spectrum calculation unit by the reproduction sound field angular spectrum calculated by the reproduction sound field angular spectrum calculation unit, wherein the desired sound field angular spectrum calculation unit calculates a sound pressure angular spectrum P^d(kx, yw) of the wall position yw for each wave number _kx using the virtual sound source coordinates ( _xs , _ys ) and the sound source signal s(n) in a buffering section of a predetermined buffer length, and calculates a sound pressure angular spectrum P^ _d ( _kx , _yw ) at the wall position _yw for each wave number kx using the wave number k calculated by the wave number calculation unit, the wave number _kx calculated by the x-axis direction wave number calculation unit, and the sound pressure angular spectrum P^ _d ( _kx , _yw ) using the equation:
A vector P^ is calculated by the following equation using the wave number k, the wave number _kx , the preset y-axis direction particle velocity transmittance distribution angular spectrum H( _kx ), the desired boundary _ydes , and the wall position _yw :
(where n and m = 1 to M (n and m are integers) in k ⁿ _x and ^km _x ), the angular spectrum H (k ⁿ _x , ^km _x ) is calculated using the formula:
and multiplying the discretization interval Δk _x calculated by the x-axis direction wave number calculation unit, the vector P^, and the matrix H to calculate the desired sound field angular spectrum.

以上のように、本発明によれば、ＳＤＭによるスピーカアレイの駆動信号を生成する際に、所望音場の音源と所望境界との間に開口部を有する壁を設置した場合において、当該壁の影響を考慮した波面を再現することができる。その結果、壁の存在によるオクルージョンまたは音色の変化も再現することができる。 As described above, according to the present invention, when generating a drive signal for a speaker array using SDM, if a wall with an opening is placed between the sound source and the desired boundary of the desired sound field, it is possible to reproduce a wavefront that takes into account the influence of the wall. As a result, it is also possible to reproduce occlusion or changes in tone due to the presence of the wall.

本発明の実施形態において想定する音場の構成例を説明する図である。1A and 1B are diagrams illustrating an example of the configuration of a sound field assumed in an embodiment of the present invention. 本発明の実施形態による音場再現装置の構成例を示すブロック図である。1 is a block diagram illustrating an example of the configuration of a sound field reproduction device according to an embodiment of the present invention. 駆動信号算出部の構成例を示すブロック図である。FIG. 3 is a block diagram showing an example of the configuration of a drive signal calculation unit. 駆動信号算出部の処理例を示すフローチャートである。10 is a flowchart illustrating an example of processing by a drive signal calculation unit. 再現音場角度スペクトル算出部の構成例を示すブロック図である。10 is a block diagram showing an example of the configuration of a reproduced sound field angular spectrum calculation unit. FIG. 所望音場角度スペクトル算出部の構成例を示すブロック図である。10 is a block diagram showing an example of the configuration of a desired sound field angular spectrum calculation unit. FIG. 波面合成法（ＷＦＳ）におけるスピーカアレイの配置例を示す図である。FIG. 1 is a diagram showing an example of speaker array arrangement in wave field synthesis (WFS). 境界音場制御（ＢｏＳＣ）におけるスピーカアレイの配置例を示す図である。FIG. 1 is a diagram showing an example of speaker array arrangement in boundary sound field control (BoSC). 高次アンビソニックス（ＨＯＡ）におけるスピーカアレイの配置例を示す図である。FIG. 1 is a diagram illustrating an example of a speaker array arrangement in Higher Order Ambisonics (HOA). スペクトル除算法（ＳＤＭ）におけるスピーカアレイの配置例を示す図である。FIG. 1 is a diagram showing an example of speaker array arrangement in the spectral division method (SDM). 一般的なＳＤＭを説明する図である。FIG. 1 is a diagram illustrating a typical SDM. 再現領域及び再現境界を説明する図である。FIG. 2 is a diagram illustrating a reproduction area and a reproduction boundary.

以下、本発明を実施するための形態について図面を用いて詳細に説明する。本発明は、ＳＤＭによるスピーカアレイの駆動信号を生成する際に、壁の影響を考慮した音場を定式化することを特徴とする。 The following describes in detail the embodiments of the present invention, using the accompanying drawings. The present invention is characterized by formulating a sound field that takes into account the influence of walls when generating drive signals for a speaker array using SDM.

〔スピーカアレイの駆動信号〕
まず、本発明の実施形態による音場再現装置が生成する、スピーカアレイの駆動信号について説明する。 [Speaker array drive signal]
First, a speaker array driving signal generated by a sound field reproduction device according to an embodiment of the present invention will be described.

図１は、本発明の実施形態において想定する音場の構成例を説明する図である。ｘｙｚ座標空間上の点音源の座標を（ｘ_s，ｙ_s，ｚ₀）とし、以下に示す式において、ｚ₀及びωは省略する。 1 is a diagram illustrating an example of the configuration of a sound field assumed in an embodiment of the present invention. The coordinates of a point sound source in the xyz coordinate space are ( _xs , _ys , _z0 ), and _z0 and ω are omitted in the following equations.

ｙ＝ｙ_w＞ｙ_sに厚みのない遮蔽版（開口部を有する壁５０）が設置され、例えば点音源によるｙ＜ｙ_w方向の到来音波の音圧角度スペクトルＰ^_d（ｋ_x，ｙ_w）は既知であるとする。 It is assumed that a thin shielding plate (wall 50 with an opening) is placed at y = y _w > y _s , and the sound pressure angular spectrum P^ _d (k _x , y _w ) of the sound wave arriving in the direction y < y _w from a point sound source is known.

このとき、この直線状の壁５０におけるｙ＜ｙ_w方向の（壁５０が存在しないときの）ｙ軸方向粒子速度角度スペクトルＷ（ｋ_x，ｙ_w）は、以下で与えられる。
ここで、ρは空気の密度、ｃは音速である。 In this case, the particle velocity angular spectrum W(k _x , y _w ) in the y<y _w direction (when the wall 50 does not exist) in the straight wall 50 is given by:
where ρ is the density of air and c is the speed of sound.

また、壁５０の形状及び特性等により定められるｙ＝ｙ_wにおけるｙ軸方向粒子速度の等価率分布をｈ（ｘ）とし、その角度スペクトル表現をｙ軸方向粒子速度透過率分布角度スペクトルＨ（ｋ_x）とする。そうすると、直線状の壁５０におけるｙ＞ｙ_w方向の（壁５０が存在するときの）ｙ軸方向粒子速度角度スペクトルＷ^（ｋ_x，ｙ_w）は、以下のように波数領域（空間周波数領域）での畳み込みで表現することができる。
Furthermore, the equivalent rate distribution of the y-axis direction particle velocity at y = _yw , which is determined by the shape and characteristics of the wall 50, is defined as h(x), and its angular spectrum expression is defined as the y-axis direction particle velocity transmittance distribution angular spectrum H( _kx ). Then, the y-axis direction particle velocity angular spectrum W^( _kx , _yw ) in the y > _yw direction (when the wall 50 exists) for the linear wall 50 can be expressed by convolution in the wave number domain (spatial frequency domain) as follows:

さらに、波数領域における波動場の外挿を行い、ｙ軸方向粒子速度角度スペクトルＷ^（ｋ_x，ｙ_w）から音圧角度スペクトルＰ~（ｋ_x，ｙ_des）への変換を行うことで、ｙ＝ｙ_des＞ｙ_wにおける所望音場の音圧角度スペクトルＰ~（ｋ_x，ｙ_des）は、以下の式となる。
Furthermore, by extrapolating the wave field in the wave number domain and converting the y-axis particle velocity angular spectrum W^( _kx , _yw ) to the sound pressure angular spectrum P~( _kx , _ydes ), the sound pressure angular spectrum P~( _kx , _ydes ) of the desired sound field when y = _ydes > _yw is given by the following formula.

実際に計算機で上記の積分を算出するにあたっては、ｋ_xについて、離散化及び有限長での打ち切りを行う必要がある。ここで、Ｐ^_d（ｋ_x，ｙ_w）及びＰ~（ｋ_x，ｙ_des）について、|ｋ_x|≦ｋにおいては遠方まで伝搬する平面波成分を表しており、|ｋ_x|＞ｋにおいては急激に減衰するエバネッセント波成分を表している。 When actually calculating the above integral on a computer, _kx must be discretized and truncated to a finite length. Here, P^ _d ( _kx , _yw ) and P~( _kx , _ydes ) represent plane wave components that propagate far when | _kx |≦k, and represent evanescent wave components that rapidly attenuate when | _kx |>k.

したがって、Ｐ^_d（ｋ_x，ｙ_w）＝０，Ｐ~（ｋ_x，ｙ_des）＝０（|ｋ_x|＞ｋ）として計算しても、実用上差し支えがない。つまり、前記式（１１）の積分範囲を有限で打ち切ることができ、さらに離散近似することにより、前記式（１１）を、行列Ｈ及びベクトルＰ^を用いた以下の式で表現することができる。
Therefore, there is no practical problem in performing the calculations assuming P^ _d ( _kx , _yw ) = 0, P~( _kx , _ydes ) = 0 (| _kx | > k). In other words, the integral range of the above equation (11) can be truncated to a finite value, and further, by performing discrete approximation, the above equation (11) can be expressed by the following equation using the matrix H and the vector P^.

ここで、Δｋ_xは、ｘ軸方向の波数ｋ_xの離散化間隔である。離散化したｋ_xは、|ｋⁿ _x|≦ｋ（ｎは正の整数であり、離散化された波数のインデックスを示す。後述するｍも同様である。ｍ，ｎ＝１～Ｍ）とし、以下の式を定義する。Ｍはスピーカユニット数である。
Here, Δk _x is the discretization interval of the wave number k _x in the x-axis direction. The discretized k _x satisfies |k ⁿ _x |≦k (n is a positive integer and indicates the index of the discretized wave number. The same applies to m, which will be described later; m, n = 1 to M), and the following equation is defined, where M is the number of speaker units.

ただし、Ｈ￣（ｋⁿ _x，ｋ^m _x）は以下とする。
However, H(k ⁿ _x , ^km _x ) is as follows.

したがって、前記式（４）をｋ_xについて離散化し、前記式（５）のＰ^_d（ｋ_x，ｙ_des，ｚ₀，ω）の代わりに、前記式（１２）で得られたＰ~すなわちＰ~（ｋ_x，ｙ_des）（ｚ₀，ωを省略する。）を代入することにより、駆動信号を得ることができる。 Therefore, the drive signal can be obtained by discretizing the equation (4) with respect to _kx and substituting P~ obtained in the equation (12), i.e., P _~ ( _kx , _ydes ) ( _z0 , ω are omitted), in place of P^d( _kx , _ydes , _z0 , ω) in the equation (5).

すなわち、前記式（４）をｋ_xについて離散化したときの駆動信号の角度スペクトルＤ^（ｋ_x，y₀，ｚ₀，ω）は、所望音場の音圧角度スペクトルである所望音場角度スペクトルＰ~（ｋ_x，ｙ_des）を、再現音場角度スペクトルＧ^（ｋ_x，ｙ_ref－ｙ₀，ｚ₀，ω）で除算することで得られる。 That is, when the above equation (4) is discretized with respect to _kx , the angular spectrum D^( _kx , _y0 , _z0 , ω) of the drive signal is obtained by dividing the desired sound field angular spectrum P~( _kx , _ydes ), which is the sound pressure angular spectrum of the desired sound field, by the reproduced sound field angular spectrum G^( _kx , _yref _-y0 , _z0 , ω).

このように、本発明の実施形態による音場再現装置は、前記式（１２）を用いて、所望音場の音圧角度スペクトルＰ~（ｋ_x，ｙ_des）を算出する。そして、音場再現装置は、前記式（４）をｋ_xについて離散化し、前記式（５）のＰ^_d（ｋ_x，ｙ_des，ｚ₀，ω）の代わりに、前記式（１２）のＰ~（ｋ_x，ｙ_des）を代入した式を用いて、所望音場の音圧角度スペクトルＰ~（ｋ_x，ｙ_des）及び再現音場角度スペクトルＧ^（ｋ_x，ｙ_ref－ｙ₀，ｚ₀，ω）から駆動信号の角度スペクトルＤ^（ｋ_x，ｙ₀，ｚ₀，ω）を算出する。 In this way, the sound field reproduction device according to the embodiment of the present invention calculates the sound pressure angular spectrum P~( _kx , _ydes ) of the desired sound field using the above formula (12).The sound field reproduction device then discretizes the above formula (4) with respect to _kx , and calculates the angular spectrum _D ^( _kx , _y0 , _z0 , ω) of the drive signal from the sound pressure angular spectrum P~( _kx , _ydes ) of the desired sound field and the reproduced sound field angular spectrum G^( _kx , yref - _y0 , _z0 , ω) using an equation in which P~( _kx , _ydes ) of the above _formula (12) is substituted for P^d( _kx , _ydes , _z0 , ω) of the above formula (5).

〔音場再現装置〕
次に、本発明の実施形態による音場再現装置について説明する。図２は、本発明の実施形態による音場再現装置の構成例を示すブロック図である。この音場再現装置１は、駆動信号算出部１０、空間周波数領域離散逆フーリエ変換部１１及び時間周波数領域離散逆フーリエ変換部１２を備えている。 [Sound field reproduction device]
Next, a sound field reproduction device according to an embodiment of the present invention will be described. Fig. 2 is a block diagram showing an example of the configuration of a sound field reproduction device according to an embodiment of the present invention. This sound field reproduction device 1 includes a drive signal calculation unit 10, a spatial frequency domain discrete inverse Fourier transform unit 11, and a time frequency domain discrete inverse Fourier transform unit 12.

駆動信号算出部１０は、所定バッファ長のバッファリング区間についての点音源の座標である仮想音源座標（ｘ_s，ｙ_s）及び音源信号ｓ（ｎ）を入力すると共に、予め設定された壁５０の位置での（壁位置ｙ_wにおける）ｙ軸方向粒子速度透過率分布角度スペクトルＨ（ｋ_x）、壁位置ｙ_w、所望境界ｙ_des、再現境界ｙ_ref、スピーカユニット数Ｍ、スピーカユニット間隔Δｘ、サンプリング周波数ｆ_s、及びＤＦＴ（またはＦＦＴ）タップ数であるＤＦＴ点数Ｆを入力する。ｎは、時間サンプルである。 The drive signal calculation unit 10 inputs virtual sound source coordinates ( _xs , _ys ), which are the coordinates of a point sound source for a buffering section of a predetermined buffer length, and sound source signal s(n), as well as the y-axis particle velocity transmittance distribution angular spectrum H( _kx ) at a preset position of the wall 50 (at wall position _yw ), the wall position _yw , the desired boundary _ydes , the reproduction boundary _yref , the number of speaker units M, the speaker unit spacing Δx, the sampling frequency _fs , and the number of DFT points F, which is the number of DFT (or FFT) taps, where n is a time sample.

尚、仮想音源（点音源）は、固定の仮想音源座標（ｘ_s，ｙ_s）に位置する全指向性の静止した音源であってもよいし、指向性を有し、変化する仮想音源座標（ｘ_s，ｙ_s）にて任意の軌跡を描いて移動及び回転する音源であってもよい。 The virtual sound source (point sound source) may be an omnidirectional stationary sound source located at fixed virtual sound source coordinates ( _xs , _ys ), or may be a directional sound source that moves and rotates along an arbitrary trajectory at changing virtual sound source coordinates ( _xs , _ys ).

駆動信号算出部１０は、仮想音源座標（ｘ_s，ｙ_s）、音源信号ｓ（ｎ）、壁位置ｙ_w等を用いて、壁位置ｙ_wの音圧角度スペクトルＰ^_d（ｋ_x，ｙ_w）を算出する。そして、駆動信号算出部１０は、壁位置ｙ_wの音圧角度スペクトルＰ^_d（ｋ_x，ｙ_w）、ｙ軸方向粒子速度透過率分布角度スペクトルＨ（ｋ_x）等を用いて、前記式（１２）により、壁５０の影響を考慮した所望音場角度スペクトルＰ~（ｋ_x，ｙ_des）を算出する。駆動信号算出部１０は、ＳＤＭの除算式（前記式（４）をｋ_xについて離散化し、前記式（５）の代わりに前記式（１２）を代入した前記式（４）に対応する式）により、ｘ軸方向の波数ｋ_x毎の駆動信号の角度スペクトルＤ^（ｋ_x，y₀，ｚ₀，ω）を算出し、これを空間周波数領域離散逆フーリエ変換部１１に出力する。駆動信号算出部１０の詳細については後述する。 The drive signal calculation unit 10 calculates the sound pressure angular spectrum P^ _d ( _kx , yw) at the wall position _yw using the virtual sound source coordinates ( _xs , _ys ), the sound source signal s( _n ), the wall position _yw , etc. Then, the drive signal calculation unit 10 calculates the desired sound field angular spectrum P~(kx, _ydes ) taking into account the influence of the wall ₅₀ according to the above formula ( ₁₂ ) using the sound pressure angular spectrum P^ _d ( _kx , _yw ) at the wall position yw, the y-axis particle velocity transmittance distribution angular spectrum H( _kx ), etc. The drive signal calculation unit 10 calculates the angular spectrum D^(kx, y0, z0, ω) of the drive signal for each wave _number _kx in the x-axis direction using an SDM division formula (a formula corresponding to formula (4) obtained by discretizing formula (4) with respect to _kx and substituting formula (12 ₎ instead of formula ( ₅ )), and outputs this to the spatial frequency domain discrete inverse Fourier transform unit 11. Details of the drive signal calculation unit 10 will be described later.

空間周波数領域離散逆フーリエ変換部１１は、駆動信号算出部１０から、ｘ軸方向の波数ｋ_x毎の駆動信号の角度スペクトルＤ^（ｋ_x，ｙ₀，ｚ₀，ω）を入力すると共に、予め設定されたスピーカユニット数Ｍ及びスピーカユニット間隔Δｘを入力する。 The spatial frequency domain discrete inverse Fourier transform unit 11 inputs the angular spectrum D^(k _x , y ₀ , z ₀ , ω) of the drive signal for each wave number k _x in the x-axis direction from the drive signal calculation unit 10, as well as the preset number M of speaker units and speaker unit spacing Δx.

空間周波数領域離散逆フーリエ変換部１１は、駆動信号の角度スペクトルＤ^（ｋ_x，ｙ₀，ｚ₀，ω）に対し離散逆フーリエ変換を行うことで、スピーカアレイ１００を構成するスピーカユニット毎の時間周波数領域の駆動信号Ｄ~（ｘ，ω）を求める。具体的には、空間周波数領域離散逆フーリエ変換部１１は、スピーカユニット数Ｍ及びスピーカユニット間隔Δｘからスピーカアレイ長Ｌ＝（Ｍ－１）Δｘを算出し、前記式（８）の演算を行う。 The spatial frequency domain discrete inverse Fourier transform unit 11 performs a discrete inverse Fourier transform on the angular spectrum D^(k _x , y ₀ , z ₀ , ω) of the drive signal to determine the drive signal D~(x, ω) in the time frequency domain for each speaker unit that makes up the speaker array 100. Specifically, the spatial frequency domain discrete inverse Fourier transform unit 11 calculates the speaker array length L = (M-1)Δx from the number of speaker units M and the speaker unit spacing Δx, and performs the calculation of equation (8).

空間周波数領域離散逆フーリエ変換部１１は、スピーカユニット毎の時間周波数領域の駆動信号Ｄ~（ｘ，ω）を時間周波数領域離散逆フーリエ変換部１２に出力する。 The spatial frequency domain discrete inverse Fourier transform unit 11 outputs the time frequency domain drive signal D~(x, ω) for each speaker unit to the time frequency domain discrete inverse Fourier transform unit 12.

時間周波数領域離散逆フーリエ変換部１２は、空間周波数領域離散逆フーリエ変換部１１から、スピーカユニット毎の時間周波数領域の駆動信号Ｄ~（ｘ，ω）を入力する。そして、時間周波数領域離散逆フーリエ変換部１２は、時間周波数領域の駆動信号Ｄ~（ｘ，ω）に対し離散逆フーリエ変換を行うことで、スピーカアレイ１００を構成するスピーカユニット毎の時間領域のスピーカ駆動信号を求める。これにより、壁５０の影響を考慮した音場を再現するためのスピーカ駆動信号が生成される。 The time-frequency domain discrete inverse Fourier transform unit 12 receives the time-frequency domain drive signal D~(x,ω) for each speaker unit from the spatial frequency domain discrete inverse Fourier transform unit 11. The time-frequency domain discrete inverse Fourier transform unit 12 then performs a discrete inverse Fourier transform on the time-frequency domain drive signal D~(x,ω) to obtain a time-domain speaker drive signal for each speaker unit that makes up the speaker array 100. This generates a speaker drive signal for reproducing a sound field that takes into account the influence of the wall 50.

時間周波数領域離散逆フーリエ変換部１２は、スピーカユニット毎の時間領域のスピーカ駆動信号を、対応するスピーカアレイ１００を構成するスピーカユニットへ出力する。 The time-frequency domain discrete inverse Fourier transform unit 12 outputs the time domain speaker drive signal for each speaker unit to the corresponding speaker unit that makes up the speaker array 100.

（駆動信号算出部１０）
次に、図２に示した駆動信号算出部１０について詳細に説明する。図３は、駆動信号算出部１０の構成例を示すブロック図であり、図４は、駆動信号算出部１０の処理例を示すフローチャートである。 (Drive signal calculation unit 10)
Next, a detailed description will be given of the drive signal calculation unit 10 shown in Fig. 2. Fig. 3 is a block diagram showing an example of the configuration of the drive signal calculation unit 10, and Fig. 4 is a flowchart showing an example of processing by the drive signal calculation unit 10.

この駆動信号算出部１０は、ｘ軸方向波数算出部２０、波数算出部２１、再現音場角度スペクトル算出部２２、所望音場角度スペクトル算出部２３及び除算部２４を備えている。これらの構成部について、図４のフローチャートを用いて説明する。 This drive signal calculation unit 10 includes an x-axis wave number calculation unit 20, a wave number calculation unit 21, a reproduced sound field angular spectrum calculation unit 22, a desired sound field angular spectrum calculation unit 23, and a division unit 24. These components will be explained using the flowchart in Figure 4.

駆動信号算出部１０は、予め設定されたスピーカユニット数Ｍ、スピーカユニット間隔Δｘ、サンプリング周波数ｆ_s、ＤＦＴ点数Ｆ、再現境界ｙ_ref、壁５０の位置でのｙ軸方向粒子速度透過率分布角度スペクトルＨ（ｋ_x）、壁位置ｙ_w及び所望境界ｙ_desを入力する（ステップＳ４０１）。 The drive signal calculation unit 10 inputs the preset number of speaker units M, the speaker unit spacing Δx, the sampling frequency f _s , the number of DFT points F, the reproduction boundary y _ref , the y-axis direction particle velocity transmittance distribution angular spectrum H(k _x ) at the position of the wall 50, the wall position y _w and the desired boundary y _des (step S401).

（ｘ軸方向波数算出部２０）
ｘ軸方向波数算出部２０は、予め設定されたスピーカユニット数Ｍ及びスピーカユニット間隔Δｘを入力する。そして、ｘ軸方向波数算出部２０は、スピーカユニット数Ｍ及びスピーカユニット間隔Δｘ、並びに所定のインデックスｍを用いて、ｘ軸方向の波数ｋ_x（＝２πｍ／（ΔｘＭ））を算出すると共に、ｘ軸方向の波数ｋ_xの離散化間隔Δｋ_xを算出する（ステップＳ４０２）。 (x-axis direction wave number calculation unit 20)
The x-axis direction wave number calculation unit 20 receives a preset number of speaker units M and a speaker unit interval Δx. Then, using the number of speaker units M, the speaker unit interval Δx, and a predetermined index m, the x-axis direction wave number calculation unit 20 calculates the wave number _kx (=2πm/(ΔxM)) in the x-axis direction and calculates the discretization interval _Δkx of the wave number _kx in the x-axis direction (step S402).

ここで、ｘ軸方向波数算出部２０により、ｍ＝０～Ｍ－１に対応するＭ個の波数ｋ_x（＝ｋ¹ _x，ｋ² _x，・・・，ｋ^M _x）が算出される。また、スピーカアレイ１００を構成する隣り合うスピーカユニットは等間隔であるとし、任意の隣り合うスピーカユニットに対応する２つのｘ軸方向の波数ｋ_ｘを用いて、ｘ軸方向の波数ｋ_ｘの離散化間隔Δｋ_xが算出される。 Here, M wave numbers k _x (= k ¹ _x , k ² _x , ..., k ^M _x ) corresponding to m = 0 to M-1 are calculated by the x-axis direction wave number calculation unit 20. Also, assuming that adjacent speaker units constituting the speaker array 100 are equally spaced, the discretization interval Δk _{x of the wave numbers k x} in the x-axis direction is calculated using two x-axis direction _{wave numbers k x} _{corresponding} to any adjacent speaker units.

ｘ軸方向波数算出部２０は、ｘ軸方向の波数ｋ_x（＝ｋ¹ _x，ｋ² _x，・・・，ｋ^M _x）を再現音場角度スペクトル算出部２２及び所望音場角度スペクトル算出部２３に出力し、ｘ軸方向の波数ｋ_ｘの離散化間隔Δｋ_xを所望音場角度スペクトル算出部２３に出力する。 The x-axis direction wave number calculation unit 20 ^outputs the wave numbers _kx (= ^k1x , _k2x , _... , ^kMx ) in the _x -axis direction to the reproduced sound field angular spectrum calculation unit 22 and the desired sound field angular spectrum calculation unit 23, and outputs the discretization interval _Δkx of the wave numbers _kx in the x-axis direction to the desired sound field angular spectrum calculation unit 23.

（波数算出部２１）
波数算出部２１は、予め設定されたサンプリング周波数ｆ_s及びＤＦＴ点数Ｆを入力する。そして、波数算出部２１は、サンプリング周波数ｆ_s及びＤＦＴ点数Ｆ、並びに所定の周波数インデックスｌを用いて、波数ｋ（＝ω／ｃ）及び角周波数ω（＝２πｆ_sｌ／Ｆ）を算出する（ステップＳ４０３）。ｃは音速であり、周波数インデックスｌは、ｌ＝０～Ｆ－１（ｌは整数）である。 (Wavenumber calculation unit 21)
The wave number calculation unit 21 receives a preset sampling frequency f _s and a number of DFT points F. Then, using the sampling frequency f _s , the number of DFT points F, and a predetermined frequency index l, the wave number calculation unit 21 calculates a wave number k (=ω/c) and an angular frequency ω (=2πf _s l/F) (step S403), where c is the speed of sound and the frequency index l is l=0 to F−1 (l is an integer).

波数算出部２１は、波数ｋ及び角周波数ωを再現音場角度スペクトル算出部２２及び所望音場角度スペクトル算出部２３に出力する。 The wave number calculation unit 21 outputs the wave number k and the angular frequency ω to the reproduced sound field angular spectrum calculation unit 22 and the desired sound field angular spectrum calculation unit 23.

（再現音場角度スペクトル算出部２２）
再現音場角度スペクトル算出部２２は、予め設定された再現境界ｙ_refを入力すると共に、ｘ軸方向波数算出部２０からｘ軸方向の波数ｋ_ｘを入力し、波数算出部２１から波数ｋ及び角周波数ωを入力する。 (Reproduced sound field angular spectrum calculation unit 22)
The reproduced sound field angular spectrum calculation unit 22 receives a preset reproduction boundary _yref , receives the wave number _kx in the x-axis direction from the x-axis direction wave number calculation unit 20, and receives the wave number k and angular frequency ω from the wave number calculation unit 21.

再現音場角度スペクトル算出部２２は、再現境界ｙ_ref、ｘ軸方向の波数ｋ_ｘ、波数ｋ、角周波数ω及び予め設定された値ｙ₀を用いて、前記式（４）の分母に相当する再現音場角度スペクトルＧ^（ｋ_x，ｙ_ref－ｙ₀，ｚ₀，ω）を算出する（ステップＳ４０４）。 The reproduced sound field angular spectrum calculation unit 22 calculates the reproduced sound field angular spectrum G^(k _x , y ref - y ₀ , z ₀ , ω) corresponding to the denominator of equation (4) using the reproduction _boundary y _ref , the wave number k _x in the x _- axis direction, the wave number k, the angular frequency ω, and a preset value y 0 (step S404).

前述のとおり、スピーカアレイ１００は、ｘｙｚ空間において予め設定されたｚ＝ｚ₀のｙ＝ｙ₀上に配置されており、値ｙ₀は、スピーカアレイ１００が配置されているｘｙｚ空間のｙ値を示す。そして、再現音場角度スペクトル算出部２２は、再現音場角度スペクトルＧ^（ｋ_x，ｙ_ref－ｙ₀，ｚ₀，ω）を除算部２４に出力する。 As described above, speaker array 100 is arranged on y= _y0 of preset z= _z0 in xyz space, and the value _y0 indicates the y value in the xyz space where speaker array 100 is arranged. Then, reproduced sound field angular spectrum calculation unit 22 outputs reproduced sound field angular spectrum G^( _kx , _yref - _y0 , _z0 , ω) to division unit 24.

図５は、再現音場角度スペクトル算出部２２の構成例を示すブロック図である。この再現音場角度スペクトル算出部２２は、算出部３０及びメモリ３１を備えている。 Figure 5 is a block diagram showing an example configuration of the reproduced sound field angular spectrum calculation unit 22. This reproduced sound field angular spectrum calculation unit 22 includes a calculation unit 30 and a memory 31.

算出部３０は、予め設定された再現境界ｙ_refを入力すると共に、ｘ軸方向波数算出部２０から（スピーカユニット数Ｍ分の）ｘ軸方向の波数ｋ_ｘを、波数算出部２１から波数ｋ及び角周波数ωを入力する。 The calculation unit 30 inputs a preset reproduction boundary _yref , and also inputs the wave number _kx in the x-axis direction (for the number M of speaker units) from the x-axis direction wave number calculation unit 20, and the wave number k and angular frequency ω from the wave number calculation unit 21.

算出部３０は、以下の式にて、再現音場角度スペクトルを算出し、再現音場角度スペクトルをメモリ３１に格納する。
この式（１７）は、再現音場角度スペクトルを算出する前記式（６）における上式に相当する。 The calculation unit 30 calculates the reproduced sound field angular spectrum using the following formula, and stores the reproduced sound field angular spectrum in the memory 31.
This equation (17) corresponds to the upper equation in the above equation (6) for calculating the reproduced sound field angular spectrum.

具体的には、算出部３０は、波数ｋの２乗値からｘ軸方向の波数ｋ_ｘの２乗値を減算し、減算結果の平方根を求め、再現境界ｙ_refから値ｙ₀を減算し、減算結果の絶対値を求め、前記平方根及び前記絶対値を乗算する。そして、算出部３０は、乗算結果の０次第２種ハンケル関数Ｈ₀ ⁽²⁾を演算し、演算結果に（－ｉ／４）を乗算することで、再現音場角度スペクトルを求める。ｉは虚数単位である。この再現音場角度スペクトルは、再現音場の再現境界ｙ_refにおける角度スペクトルＧ^（ｋ_x，ｙ_ref－ｙ₀，ｚ₀，ω）に相当するデータである。 Specifically, the calculation unit 30 subtracts the square of the wave number _kx in the x-axis direction from the square of the wave number k, calculates the square root of the subtraction result, subtracts the value _y0 from the reproduction boundary _yref , calculates the absolute value of the subtraction result, and multiplies the square root and the absolute value.The calculation unit 30 then calculates the zero-order second-order Hankel function _H0 ⁽²⁾ of the multiplication result and multiplies the calculation result by (-i/4) to obtain the reproduced sound field angular spectrum, where i is the imaginary unit.This reproduced sound field angular spectrum is data equivalent to the angular spectrum G^( _kx , _yref _-y0 , _z0 , ω) at the reproduction boundary _yref of the reproduced sound field.

このように、再現音場角度スペクトルは、再現境界ｙ_ref等の予め設定されたパラメータのみを用いて算出される。このため、再現音場角度スペクトルは、当該駆動信号算出部１０による駆動信号の角度スペクトルＤ^（ｋ_x，ｙ₀，ｚ₀，ω）の算出処理に先立って予め算出され、メモリ３１に格納される。 In this way, the reproduced sound field angular spectrum is calculated using only preset parameters such as the reproduction boundary _yref , etc. Therefore, the reproduced sound field angular spectrum is calculated in advance prior to the calculation process of the angular spectrum D^( _kx , _y0 , _z0 , ω) of the drive signal by the drive signal calculation unit 10, and stored in the memory 31.

（所望音場角度スペクトル算出部２３）
図３及び図４に戻って、所望音場角度スペクトル算出部２３は、所定バッファ長のバッファリング区間についての仮想音源座標（ｘ_s，ｙ_s）及び音源信号ｓ（ｎ）を入力する（ステップＳ４０５）。また、所望音場角度スペクトル算出部２３は、予め設定された再現境界ｙ_ref、壁５０の位置でのｙ軸方向粒子速度透過率分布角度スペクトルＨ（ｋ_x）、壁位置ｙ_w及び所望境界ｙ_desを入力する。さらに、所望音場角度スペクトル算出部２３は、ｘ軸方向波数算出部２０からｘ軸方向の波数ｋ_ｘ及びｘ軸方向の波数ｋ_ｘの離散化間隔Δｋ_xを入力し、波数算出部２１から波数ｋ及び角周波数ωを入力する。 (Desired sound field angle spectrum calculation unit 23)
3 and 4 , the desired sound field angular spectrum calculation unit 23 receives input of the virtual sound source coordinates ( _xs , _ys ) and sound source signal s(n) for the buffering section of a predetermined buffer length (step S405). The desired sound field angular spectrum calculation unit 23 also receives input of the preset reproduction boundary _yref , the y-axis particle velocity transmittance distribution angular spectrum H( _kx ) at the position of the wall 50, the wall position _yw , and the desired boundary _ydes . The desired sound field angular spectrum calculation unit 23 also receives input of the wave number _kx in the x-axis direction and the discretization interval _Δkx of the wave number _kx in the x-axis direction from the x-axis direction wave number calculation unit 20, and inputs the wave number k and angular frequency ω from the wave number calculation unit 21.

所望音場角度スペクトル算出部２３は、所定バッファ長の仮想音源座標（ｘ_s，ｙ_s）及び音源信号ｓ（ｎ）、並びに壁位置ｙ_w等を用いて、例えば点音源が全指向性の静止した音源である場合には前記式（５）（後述する式（１８））により、壁位置ｙ_wの音圧角度スペクトルＰ^_d（ｋ_x，ｙ_w，ｚ₀，ω）を算出する（ステップＳ４０６）。 The desired sound field angular spectrum calculation unit 23 uses the virtual sound source coordinates ( _xs , _ys ) and sound source signal s(n) of a predetermined buffer length, as well as the wall position _yw , to calculate the sound pressure angular spectrum P^ _d ( _kx , _yw , _z0 , ω) at the wall position _yw , for example, using equation (5) (equation (18) described later) when the point sound source is an omnidirectional stationary sound source (step S406).

所望音場角度スペクトル算出部２３は、壁位置ｙ_wの音圧角度スペクトルＰ^_d（ｋ_x，ｙ_w，ｚ₀，ω）、壁５０の位置でのｙ軸方向粒子速度透過率分布角度スペクトルＨ（ｋ_x）等を用いて、壁５０の影響を考慮した前記式（１２）により、音圧角度スペクトルである所望音場角度スペクトルＰ~（ｋ_x，ｙ_des，ｚ₀，ω）を算出する（ステップＳ４０７）。そして、所望音場角度スペクトル算出部２３は、所望音場角度スペクトルＰ~（ｋ_x，ｙ_des，ｚ₀，ω）を除算部２４に出力する。 The desired sound field angular spectrum calculation unit 23 calculates the desired sound _field angular spectrum P~( _kx , _ydes , _z0 , ω), which is the sound pressure angular spectrum, by using the sound pressure angular spectrum P^ _d ( _kx , yw, _z0 , ω) at the wall position yw, the y-axis particle velocity transmittance distribution angular spectrum H( _kx ) at the position of the _wall 50, etc., according to the above formula (12) taking into account the influence of the wall 50 (step S407). Then, the desired sound field angular spectrum calculation unit 23 outputs the desired sound field angular spectrum P~( _kx , _ydes , _z0 , ω) to the division unit 24.

図６は、所望音場角度スペクトル算出部２３の構成例を示すブロック図である。この所望音場角度スペクトル算出部２３は、音圧角度スペクトル算出部４０、行列Ｈ算出部４１、ベクトルＰ^算出部４２及び乗算部４３を備えている。 Figure 6 is a block diagram showing an example configuration of the desired sound field angular spectrum calculation unit 23. This desired sound field angular spectrum calculation unit 23 includes a sound pressure angular spectrum calculation unit 40, a matrix H calculation unit 41, a vector P^ calculation unit 42, and a multiplication unit 43.

音圧角度スペクトル算出部４０は、所定バッファ長のバッファリング区間に対応する時間サンプル毎の仮想音源座標（ｘ_s，ｙ_s）及び音源信号ｓ（ｎ）を入力すると共に、予め設定された壁位置ｙ_wを入力する。また、音圧角度スペクトル算出部４０は、ｘ軸方向波数算出部２０からｘ軸方向の波数ｋ_ｘを入力し、波数算出部２１から波数ｋ及び角周波数ωを入力する。 The sound pressure angular spectrum calculation unit 40 receives the virtual sound source coordinates ( _xs , _ys ) and sound source signal s(n) for each time sample corresponding to a buffering section of a predetermined buffer length, as well as a preset wall position _yw . The sound pressure angular spectrum calculation unit 40 also receives the wave number _kx in the x-axis direction from the x-axis direction wave number calculation unit 20, and the wave number k and angular frequency ω from the wave number calculation unit 21.

音圧角度スペクトル算出部４０は、例えば点音源が全指向性の静止した音源である場合、前記式（５）において所望境界ｙ_desの代わりに壁位置ｙ_wを用いた以下の式により、所定バッファ長の仮想音源座標（ｘ_s，ｙ_s）及び音源信号ｓ（ｎ）、並びに波数ｋ、ｘ軸方向の波数ｋ_ｘ及び壁位置ｙ_wを用いて、壁位置ｙ_wの音圧角度スペクトルＰ^_d（ｋ_x，ｙ_w，ｚ₀，ω）を算出する。
Ｓ（ω）は、音源信号ｓ（ｎ）をフーリエ変換して得られる信号である。 For example, when the point sound source is an omnidirectional stationary sound source, the sound pressure angular spectrum calculation unit 40 calculates the sound pressure angular spectrum P^ _d (kx, yw, z0, ω) at the wall position _yw using the virtual sound source coordinates ( _xs , _ys ) of a predetermined buffer length, the sound source signal s(n), the wave number k, the wave number _kx in the x-axis direction, and the wall position _yw , using the following equation, which uses the wall position _yw instead _of the desired boundary _ydes in equation ₍ ₅ ):
S(ω) is a signal obtained by Fourier transforming the sound source signal s(n).

この壁位置ｙ_wの音圧角度スペクトルＰ^_d（ｋ_x，ｙ_w，ｚ₀，ω）は、ｘ軸方向の波数ｋ_x（＝ｋ¹ _x，ｋ² _x，・・・，ｋ^M _x）毎に算出され、前記式（１５）におけるＰ^_d（ｋ¹ _x，ｙ_w），Ｐ^_d（ｋ² _x，ｙ_w），・・・，Ｐ^_d（ｋ^M _x，ｙ_w）に相当する。 The sound pressure angular spectrum P^ _d ( _kx , _yw , _z0 , ω) at this wall position _yw is calculated for each wave number _kx (= ^k1x , ^k2x , _... , ^kMx ) in the x _{-axis direction, and corresponds to P^d(k1x, yw), P^d(k2x} _, _yw ⁾ _, _... _, ^P _^ _d ₍ ^kMx _, _yw ) in equation (15).

音圧角度スペクトル算出部４０は、ｘ軸方向の波数ｋ_x毎の壁位置ｙ_wの音圧角度スペクトルＰ^_d（ｋ_x，ｙ_w，ｚ₀，ω）をベクトルＰ^算出部４２に出力する。 The sound pressure angular spectrum calculation unit 40 outputs the sound pressure angular spectrum P^ _d ( _kx , _yw , _z0 , ω) at the wall position _yw for _each wave number kx in the x-axis direction to the vector P^ calculation unit 42.

尚、前述の例では、音圧角度スペクトル算出部４０は、点音源が全指向性の静止した音源である場合を示す前記式（１８）により、音圧角度スペクトルＰ^_d（ｋ_x，ｙ_w，ｚ₀，ω）を算出するようにした。これに対し、音圧角度スペクトル算出部４０は、点音源が指向性を有し、任意の軌跡を描いて移動及び回転する音源である場合において、音圧角度スペクトルＰ^_d（ｋ_x，ｙ_w，ｚ₀，ω）を算出するようにしてもよい。 In the above example, the sound pressure angular spectrum calculation unit 40 calculates the sound pressure angular spectrum P^ _d ( _kx , _yw , _z0 , ω) using the above formula (18), which indicates the case where the point sound source is an omnidirectional, stationary sound source. However, the sound pressure angular spectrum calculation unit 40 may also calculate the sound pressure angular spectrum P^ _d ( _kx , _yw , _z0 , ω) when the point sound source has directionality and moves and rotates along an arbitrary trajectory.

この場合の音圧角度スペクトルＰ^_d（ｋ_x，ｙ_w，ｚ₀，ω）の算出処理については、本件特許出願の同一の出願人及び発明者によりなされた、本件特許出願時に未公開の特願２０２１－１３４９４２号公報を参照されたい。音圧角度スペクトルＰ^_d（ｋ_x，ｙ_w，ｚ₀，ω）は、点音源の静止または動作状態、指向性等の様々な条件により異なるものとなるが、本発明は、音圧角度スペクトルＰ^_d（ｋ_x，ｙ_w，ｚ₀，ω）の算出処理を限定するものではない。 For the calculation process of the sound pressure angular spectrum P^ _d ( _kx , _yw , _z0 , ω) in this case, please refer to Japanese Patent Application No. 2021-134942, which was made by the same applicant and inventor of the present patent application and was unpublished at the time of filing the present patent application. The sound pressure angular spectrum P^ _d ( _kx , _yw , _z0 , ω) will vary depending on various conditions such as the stationary or moving state and directivity of the point sound source, but the present invention does not limit the calculation process of the sound pressure angular spectrum P^ _d ( _kx , _yw , _z0 , ω).

行列Ｈ算出部４１は、予め設定された壁５０の位置でのｙ軸方向粒子速度透過率分布角度スペクトルＨ（ｋ_x）、壁位置ｙ_w及び所望境界ｙ_desを入力すると共に、ｘ軸方向波数算出部２０からｘ軸方向の波数ｋ_ｘを入力し、波数算出部２１から波数ｋ及び角周波数ωを入力する。 The matrix H calculation unit 41 inputs the y-axis particle velocity transmittance distribution angular spectrum H( _kx ) at a preset position of the wall 50, the wall position _yw , and the desired boundary _ydes , as well as the wave number _kx in the x-axis direction from the x-axis wave number calculation unit 20 and the wave number k and angular frequency ω from the wave number calculation unit 21.

行列Ｈ算出部４１は、ｘ軸方向の波数ｋ_ｘにおける２つの組み合わせ（ｋⁿ _x，ｋ^m _x）毎に、ｘ軸方向の波数ｋⁿ _x，ｋ^m _x、波数ｋ、壁位置ｙ_w及び所望境界ｙ_desを用いて、前記式（１６）により、角度スペクトルＨ￣（ｋⁿ _x，ｋ^m _x）を算出し、これらを要素として構成される前記式（１４）に示した行列Ｈを求める。ｎ，ｍは、１からＭまでの整数である。 The matrix H calculation unit 41 calculates the _angular spectrum H( _knx , _kmx ) for each combination ( ^knx , ^kmx ) of the wave number _kx in the x-axis direction using the wave _numbers ^knx , ^kmx in the _x - ^axis direction, the wave number ^k , the wall position _yw , and the desired boundary _ydes according to equation ( ₁₆ ), and determines the matrix H shown in equation (14) that is configured using these as elements, where n and m are integers from 1 to M.

具体的には、行列Ｈ算出部４１は、前記式（１６）のとおり、ｘ軸方向の波数ｋ_ｘにおける２つの組み合わせ（ｋⁿ _x，ｋ^m _x）について、波数ｋの２乗値からｘ軸方向の波数ｋⁿ _ｘの２乗値を減算し、減算結果の平方根を求め、所望境界ｙ_desから壁位置ｙ_wを減算して減算結果を求め、（－ｉ）に前記平方根及び前記減算結果を乗算して乗算結果を求め、自然数ｅを底とし、前記乗算結果を指数とする指数関数を求める。 Specifically, as in equation (16), for two combinations ( ^knx , ^kmx ) of wave numbers _kx in the _x -axis direction, the matrix _H calculation unit 41 subtracts the squared value of the wave number knx in the x-axis direction from the squared value of the wave _number ^k , calculates the square root of the subtraction result, subtracts the wall position _yw from the desired _boundary ydes to calculate the subtraction result, multiplies the square root and the subtraction result by (-i) to calculate the multiplication result, and calculates an exponential function with the natural number e as the base and the multiplication result as the exponent.

行列Ｈ算出部４１は、前記指数関数を前記平方根で除算して除算結果を求め、ｘ軸方向の波数ｋⁿ _ｘからｘ軸方向の波数ｋ^m _ｘを減算して得られた減算結果をｋ_x”とした場合のｙ軸方向粒子速度透過率分布角度スペクトルＨ（ｋ_x”）を求める。そして、行列Ｈ算出部４１は、前記除算結果に前記ｙ軸方向粒子速度透過率分布角度スペクトルＨ（ｋ_x”）を乗算することで、角度スペクトルＨ￣（ｋⁿ _x，ｋ^m _x）を求める。 The matrix H calculation unit 41 divides the exponential function by the square root to obtain the division result, and calculates the y-axis direction particle velocity transmittance distribution angular spectrum H( ^kx " ₎ when the subtraction result obtained by subtracting the x-axis direction wavenumber _kmx from the x-axis direction wavenumber ^knx is set to _kx ". Then, the matrix H calculation unit 41 multiplies the y-axis direction particle velocity transmittance distribution angular spectrum H( _kx ") by the y-axis _direction _particle velocity transmittance distribution angular spectrum H( _kx ") to obtain the angular spectrum H( ^knx , ^kmx ).

このように、行列Ｈ算出部４１は、ｎ＝１，・・・，Ｍ及びｍ＝１，・・・，Ｍの組み合わせ（ｋ¹ _x，ｋ¹ _x），・・・，（ｋ¹ _x，ｋ^M _x），（ｋ² _x，ｋ¹ _x），・・・，（ｋ² _x，ｋ^M _x），・・・，（ｋ^M _x，ｋ¹ _x），・・・，（ｋ^M _x，ｋ^M _x）のそれぞれについて前述の処理を行うことで、前記式（１４）に示した行列Ｈにおける各要素である角度スペクトルＨ￣（ｋ¹ _x，ｋ¹ _x），・・・，Ｈ￣（ｋ¹ _x，ｋ^M _x），Ｈ￣（ｋ² _x，ｋ¹ _x），・・・，Ｈ￣（ｋ² _x，ｋ^M _x），・・・，Ｈ￣（ｋ^M _x，ｋ¹ _x），・・・，Ｈ￣（ｋ^M _x，ｋ^M _x）を求める。これにより、前記式（１４）に示した行列Ｈが算出される。 In this way, the matrix H calculation unit 41 performs the above-mentioned processing for each of the combinations ( ^k1x , ^k1x ), ..., ( _k1x , _kMx ), ( _k2x , _k1x ), _... , ( ^k2x , ^kMx ), ..., ( ^kMx , ^k1x ₎ , ..., ( ^kMx , ^kMx ), where ⁿ ₌ ¹ , _... , M _and _m ₌ 1, ..., _M , to obtain the angular spectra H( ^k1x , ^k1x ), _... , H ₍ ^k1x , ^kMx ), H( ^k2x , _k1x ), ..., H( _k2x , ^kMx ), ^... , _H ( ^kMx , ^k1x ), _... , H( ^kMx , ^kMx ), which are the ^elements of ^the _matrix H _shown _in _equation ₍ ¹⁴ ₎ . As a result, the matrix H shown in the above equation (14) is calculated.

行列Ｈ算出部４１は、前記式（１４）に示した行列Ｈを乗算部４３に出力する。 The matrix H calculation unit 41 outputs the matrix H shown in equation (14) to the multiplication unit 43.

ベクトルＰ^算出部４２は、音圧角度スペクトル算出部４０からｘ軸方向の波数ｋ_x毎の音圧角度スペクトルＰ^_d（ｋ_x，ｙ_w，ｚ₀，ω）を入力すると共に、予め設定された壁位置ｙ_wを入力する。また、ベクトルＰ^算出部４２は、ｘ軸方向波数算出部２０からｘ軸方向の波数ｋ_ｘを入力し、波数算出部２１から波数ｋ及び角周波数ωを入力する。 The vector P^ calculation unit 42 receives the sound pressure angular spectrum P^ _d ( _kx , _yw , _z0 , ω) for _each wave number kx in the x-axis direction from the sound pressure angular spectrum calculation unit 40, as well as the preset wall position _yw . The vector P^ calculation unit 42 also receives the wave number _kx in the x-axis direction from the x-axis direction wave number calculation unit 20, and the wave number k and angular frequency ω from the wave number calculation unit 21.

ベクトルＰ^算出部４２は、ｘ軸方向の波数ｋ_x毎に、波数ｋの２乗値から当該波数ｋ_xの２乗値を減算し、減算結果の平方根を求める。これにより、前記式（１５）に示したベクトルＰ^において、各要素に含まれる√（ｋ²－ｋ¹ _ｘ ²），√（ｋ²－ｋ² _ｘ ²），・・・，√（ｋ²－ｋ^M _ｘ ²）が算出される。 The vector P^ calculation unit 42 subtracts the square of the wave number _kx from the square of the wave number k for each wave number _kx in the x-axis direction, and calculates the square root of the subtraction result. As a result, √(k ² - ^{k 1} _x ² ), √(k 2 - k ² _x ² ), ..., √(k ² ^- ^{k M} _x ² ) contained in each element of the vector P^ shown in equation (15) are calculated.

ベクトルＰ^算出部４２は、ｘ軸方向の波数ｋ_x毎に、前記平方根に、音圧角度スペクトルＰ^_d（ｋ_x，ｙ_w，ｚ₀，ω）であるＰ^_d（ｋ¹ _ｘ，ｙ_w），Ｐ^_d（ｋ² _ｘ，ｙ_w），・・・，Ｐ^_d（ｋ^M _ｘ，ｙ_w）を乗算することで、前記式（１５）に示したベクトルＰ^における各要素を求める。これにより、前記式（１５）に示したベクトルＰ^が算出される。 The vector P^ calculation unit 42 _multiplies the square root by P^ _d (k1x,yw), P^d( _k2x , _yw ), ..., P^ _d ( _kMx , _yw ), which are the sound pressure angular spectrum P^ _d ( ^kx , _yw , _z0 , _ω ), for ^each wave number _kx in _{the x} ^- _axis direction, to find each element of the vector P^ shown in equation (15). In this way, the vector P^ shown in equation (15) is calculated.

ベクトルＰ^算出部４２は、前記式（１５）に示したベクトルＰ^を乗算部４３に出力する。 The vector P^ calculation unit 42 outputs the vector P^ shown in equation (15) to the multiplication unit 43.

乗算部４３は、行列Ｈ算出部４１から行列Ｈを入力すると共に、ベクトルＰ^算出部４２からベクトルＰ^を入力し、さらに、ｘ軸方向波数算出部２０からｘ軸方向の波数ｋ_ｘの離散化間隔Δｋ_xを入力する。 The multiplication unit 43 inputs the matrix H from the matrix H calculation unit 41, the vector P^ from the vector P^ calculation unit 42, and further inputs the discretization interval Δk _x of the wave number k _x in the x-axis direction from the x-axis direction wave number calculation unit 20.

乗算部４３は、ｘ軸方向の波数ｋ_ｘの離散化間隔Δｋ_x、行列Ｈ及びベクトルＰ^を乗算することで、所望音場角度スペクトルＰ~（ｋ_x，ｙ_des，ｚ₀，ω）を求める。これにより、前記式（１２）に示した音圧角度スペクトルである所望音場角度スペクトルＰ~（ｋ_x，ｙ_des，ｚ₀，ω）が算出される。 The multiplication unit 43 obtains the desired sound field angular spectrum P~( _kx , _ydes , _z0 , ω) by multiplying the discretization interval _Δkx of the wave number _kx in the x-axis direction, the matrix H, and the vector P^. This calculates the desired sound field angular spectrum P~( _kx , _ydes , _z0 , ω), which is the sound pressure angular spectrum shown in equation (12).

乗算部４３は、所望音場角度スペクトルＰ~（ｋ_x，ｙ_des，ｚ₀，ω）を除算部２４に出力する。 The multiplication unit 43 outputs the desired sound field angular spectrum P (k _x , y _des , z ₀ , ω) to the division unit 24 .

（除算部２４）
図３及び図４に戻って、除算部２４は、再現音場角度スペクトル算出部２２から再現音場角度スペクトルＧ^（ｋ_x，ｙ_ref－ｙ₀，ｚ₀，ω）を入力すると共に、所望音場角度スペクトル算出部２３から所定バッファ長毎の所望音場角度スペクトルＰ~（ｋ_x，ｙ_des，ｚ₀，ω）を入力する。 (Division section 24)
Returning to Figures 3 and 4, the division unit 24 inputs the reproduced sound field angular spectrum G^(k _x , y _ref - y ₀ , z ₀ , ω) from the reproduced sound field angular spectrum calculation unit 22, and also inputs the desired sound field angular spectrum P~(k _x , y _des , z ₀ , ω) for each specified buffer length from the desired sound field angular spectrum calculation unit 23.

除算部２４は、壁５０の影響を考慮したＳＤＭの除算式にて、所定バッファ長毎の所望音場角度スペクトルＰ~（ｋ_x，ｙ_des，ｚ₀，ω）を再現音場角度スペクトルＧ^（ｋ_x，ｙ_ref－ｙ₀，ｚ₀，ω）で除算することで、駆動信号の角度スペクトルＤ^（ｋ_x，ｙ₀，ｚ₀，ω）を求める。そして、除算部２４は、駆動信号の角度スペクトルＤ^（ｋ_x，ｙ₀，ｚ₀，ω）を空間周波数領域離散逆フーリエ変換部１１に出力する（ステップＳ４０８）。 The divider 24 obtains the angular spectrum D^(kx, y0, _z0 , ω) of the drive signal by dividing the desired sound field angular spectrum P~( _kx , _ydes , z0, ω) for each predetermined buffer length by the reproduced sound field angular spectrum G^( _kx , _yref _-y0 , _z0 , ω) using the SDM division formula that takes into account _the influence of the wall _50. The divider 24 then outputs the angular spectrum D^( _kx , _y0 , _z0 , ω) of the drive signal to _the spatial frequency domain discrete inverse Fourier transform unit 11 (step S408).

駆動信号算出部１０は、所定の終了条件を満たさない限り（ステップＳ４０９：Ｎ）、ステップＳ４０５へ移行する。そして、駆動信号算出部１０は、次の所定バッファ長のバッファリング区間について駆動信号の角度スペクトルＤ^（ｋ_x，ｙ₀，ｚ₀，ω）を算出するためのステップＳ４０５～Ｓ４０８の処理を繰り返す。一方、駆動信号算出部１０は、所定の終了条件を満たす場合（ステップＳ４０９：Ｙ）、処理を終了する。 Unless a predetermined termination condition is met (step S409: N), the drive signal calculation unit 10 proceeds to step S405. The drive signal calculation unit 10 then repeats the processes of steps S405 to S408 to calculate the angular spectrum D^(k _x , y ₀ , z ₀ , ω) of the drive signal for the next buffering section of the predetermined buffer length. On the other hand, if the predetermined termination condition is met (step S409: Y), the drive signal calculation unit 10 terminates the process.

以上のように、本発明の実施形態の音場再現装置１によれば、駆動信号算出部１０の再現音場角度スペクトル算出部２２は、前記式（１７）にて、ＳＤＭの除算式の分母に相当する再現音場角度スペクトルＧ^（ｋ_x，ｙ_ref－ｙ₀，ｚ₀，ω）を算出する。 As described above, according to the sound field reproduction device 1 of the embodiment of the present invention, the reproduced sound field angular spectrum calculation unit 22 of the drive signal calculation unit 10 calculates the reproduced sound field angular spectrum G^(k _x , y _ref - y ₀ , z ₀ , ω) corresponding to the denominator of the SDM division formula in equation (17).

駆動信号算出部１０の所望音場角度スペクトル算出部２３は、所定バッファ長の仮想音源座標（ｘ_s，ｙ_s）及び音源信号ｓ（ｎ）、並びに波数ｋ、ｘ軸方向の波数ｋ_ｘ及び壁位置ｙ_wを用いて、前記式（１８）により、壁位置ｙ_wの音圧角度スペクトルＰ^_d（ｋ_x，ｙ_w，ｚ₀，ω）を算出する。 The desired sound field angular spectrum calculation unit 23 of the drive signal calculation unit 10 uses the virtual sound source coordinates ( _xs , _ys ) and sound source signal s(n) of a predetermined buffer length, as well as the wave number k, the wave number kx in the _x -axis direction, and the wall position _yw to calculate the sound pressure angular spectrum P^ _d ( _kx , _yw , _z0 , ω) at the wall position _yw using equation (18).

駆動信号算出部１０の所望音場角度スペクトル算出部２３は、ｘ軸方向の波数ｋ_ｘにおける２つの組み合わせ（ｋⁿ _x，ｋ^m _x）毎に、ｘ軸方向の波数ｋⁿ _x，ｋ^m _x、波数ｋ、壁位置ｙ_w及び所望境界ｙ_desを用いて、前記式（１６）により、角度スペクトルＨ￣（ｋⁿ _x，ｋ^m _x）を算出し、前記式（１４）に示した行列Ｈを求める。 The desired sound field angular spectrum calculation unit 23 of the drive signal calculation unit 10 calculates the angular spectrum _H ( ^knx , ^kmx ) for each two combinations ( _knx , _kmx ) of wave number _kx in the x-axis direction using the wave numbers ^knx , _kmx in ^the ^x _- ^axis direction, wave number k, wall position _yw , and desired boundary _ydes according to equation ( ₁₆ ), and finds the matrix H shown in equation (14).

駆動信号算出部１０の所望音場角度スペクトル算出部２３は、波数ｋ、波数ｋ_x及び音圧角度スペクトルＰ^_d（ｋ_x，ｙ_w，ｚ₀，ω）であるＰ^_d（ｋ¹ _ｘ，ｙ_w），Ｐ^_d（ｋ² _ｘ，ｙ_w），・・・，Ｐ^_d（ｋ^M _ｘ，ｙ_w）を用いて、前記式（１５）により、ベクトルＰ^を算出する。 The desired sound field angular spectrum calculation unit 23 of the drive signal calculation unit 10 calculates the vector P _^ using the wave number k, wave number _kx and sound pressure angular spectrum P^ _d ( _kx , _yw _, _z0 , _ω ), i.e., P^ _d ( ^k1x _, _yw ), P^ _d ( ^k2x , yw), ... ^, P^d( _kMx , _yw ) according to the above equation (15).

駆動信号算出部１０の所望音場角度スペクトル算出部２３は、ｘ軸方向の波数ｋ_ｘの離散化間隔Δｋ_x、行列Ｈ及びベクトルＰ^を用いて、前記式（１２）により、所望音場角度スペクトルＰ~（ｋ_x，ｙ_des，ｚ₀，ω）を算出する。 The desired sound field angular spectrum calculation unit 23 of the drive signal calculation unit 10 calculates the desired sound field angular spectrum P~( _kx , _ydes , _z0 , ω) using the discretization interval _Δkx of the wave number _kx in the x-axis direction, the matrix H, and the vector P^ according to the above equation (12).

駆動信号算出部１０の除算部２４は、所望音場角度スペクトルＰ~（ｋ_x，ｙ_des，ｚ₀，ω）及び再現音場角度スペクトルＧ^（ｋ_x，ｙ_ref－ｙ₀，ｚ₀，ω）を用いて、ＳＤＭの除算式により、ｘ軸方向の波数ｋ_x毎の駆動信号の角度スペクトルＤ^（ｋ_x，ｙ₀，ｚ₀，ω）を算出する。 The division unit 24 of the drive signal calculation unit 10 uses the desired sound field angular spectrum P~(k _x , y _des , z ₀ , ω) and the reproduced sound field angular spectrum G^(k _x , y _ref - y ₀ , z ₀ , ω) to calculate the angular spectrum D^(k _x , y ₀ , z ₀ , ω) of the drive signal for each wave number k _x in the x-axis direction according to the SDM division formula.

空間周波数領域離散逆フーリエ変換部１１は、ｘ軸方向の波数ｋ_x毎の駆動信号の角度スペクトルＤ^（ｋ_x，ｙ₀，ｚ₀，ω）に対し空間周波数領域離散逆フーリエ変換を行うことで、スピーカユニット毎の時間周波数領域の駆動信号Ｄ~（ｘ，ω）を求める。 The spatial frequency domain discrete inverse Fourier transform unit 11 performs a spatial frequency domain discrete inverse Fourier transform on the angular spectrum D^(k _x , y ₀ , z ₀ , ω) of the drive signal for each wave number k _x in the x-axis direction, thereby obtaining the drive signal D~(x, ω) in the time frequency domain for each speaker unit.

時間周波数領域離散逆フーリエ変換部１２は、スピーカユニット毎の時間周波数領域の駆動信号Ｄ~（ｘ，ω）に対し時間周波数領域離散逆フーリエ変換を行うことで、スピーカユニット毎の時間領域のスピーカ駆動信号を求める。そして、時間周波数領域離散逆フーリエ変換部１２は、スピーカ駆動信号を、スピーカアレイ１００を構成するスピーカユニットへ出力する。 The time-frequency domain discrete inverse Fourier transform unit 12 performs a time-frequency domain discrete inverse Fourier transform on the time-frequency domain drive signal D~(x,ω) for each speaker unit to obtain a time-domain speaker drive signal for each speaker unit. The time-frequency domain discrete inverse Fourier transform unit 12 then outputs the speaker drive signals to the speaker units that make up the speaker array 100.

このように、本発明の実施形態では、ＳＤＭによるスピーカアレイ１００の駆動信号を生成する際に、所望音場の音源と所望境界との間に開口部を有する壁を設置した場合において、当該壁の影響を考慮した波面を再現することができる。その結果、壁の存在によるオクルージョンまたは音色の変化も再現することができる。 In this way, in an embodiment of the present invention, when generating a drive signal for the speaker array 100 using SDM, if a wall with an opening is placed between the sound source and the desired boundary of the desired sound field, it is possible to reproduce a wavefront that takes into account the influence of the wall. As a result, it is also possible to reproduce occlusions or changes in tone due to the presence of the wall.

以上、実施形態を挙げて本発明を説明したが、本発明は前記実施形態に限定されるものではなく、その技術思想を逸脱しない範囲で種々変形可能である。 The present invention has been described above using embodiments, but the present invention is not limited to these embodiments and can be modified in various ways without departing from the technical concept thereof.

尚、本発明の実施形態による音場再現装置１のハードウェア構成としては、通常のコンピュータを使用することができる。音場再現装置１は、ＣＰＵ、ＲＡＭ等の揮発性の記憶媒体、ＲＯＭ等の不揮発性の記憶媒体、及びインターフェース等を備えたコンピュータによって構成される。 The hardware configuration of the sound field reproduction device 1 according to an embodiment of the present invention can be a conventional computer. The sound field reproduction device 1 is composed of a computer equipped with a CPU, a volatile storage medium such as RAM, a non-volatile storage medium such as ROM, an interface, etc.

音場再現装置１に備えた駆動信号算出部１０、空間周波数領域離散逆フーリエ変換部１１及び時間周波数領域離散逆フーリエ変換部１２の各機能は、これらの機能を記述したプログラムをＣＰＵに実行させることによりそれぞれ実現される。 The functions of the drive signal calculation unit 10, spatial frequency domain discrete inverse Fourier transform unit 11, and time frequency domain discrete inverse Fourier transform unit 12 provided in the sound field reproduction device 1 are each realized by causing a CPU to execute a program that describes these functions.

これらのプログラムは、前記記憶媒体に格納されており、ＣＰＵに読み出されて実行される。また、これらのプログラムは、磁気ディスク（フロッピー（登録商標）ディスク、ハードディスク等）、光ディスク（ＣＤ－ＲＯＭ、ＤＶＤ等）、半導体メモリ等の記憶媒体に格納して頒布することもでき、ネットワークを介して送受信することもできる。 These programs are stored on the storage medium and are read and executed by the CPU. These programs can also be stored and distributed on storage media such as magnetic disks (floppy disks, hard disks, etc.), optical disks (CD-ROMs, DVDs, etc.), and semiconductor memory, and can also be sent and received over a network.

１音場再現装置
１０駆動信号算出部
１１空間周波数領域離散逆フーリエ変換部
１２時間周波数領域離散逆フーリエ変換部
２０ｘ軸方向波数算出部
２１波数算出部
２２再現音場角度スペクトル算出部
２３所望音場角度スペクトル算出部
２４除算部
３０算出部
３１メモリ
４０音圧角度スペクトル算出部
４１行列Ｈ算出部
４２ベクトルＰ^算出部
４３乗算部
５０壁
１００，１００－１，１００－２，１００－３，１００－４スピーカアレイ
１０１無限直線音源
（ｘ_s，ｙ_s）仮想音源座標
ｓ（ｎ）音源信号
ｙ_w 壁位置
ｙ_des 所望境界
ｙ_ref 再現境界
Ｈ（ｋ_x）ｙ軸方向粒子速度透過率分布角度スペクトル
Ｍスピーカユニット数
Δｘスピーカユニット間隔
ｆ_s サンプリング周波数
ＦＤＦＴ点数
ｋ波数
ｋ_x ｘ軸方向の波数
Δｋ_x ｋ_xの離散化間隔
ω 角周波数
Ｄ^（ｋ_x，ｙ₀，ｚ₀，ω）駆動信号の角度スペクトル
Ｐ~（ｋ_x，ｙ_des），Ｐ~（ｋ_x，ｙ_des，ｚ₀，ω）所望音場角度スペクトル
Ｇ^（ｋ_x，ｙ_ref－ｙ₀，ｚ₀，ω）再現音場角度スペクトル 1 Sound field reproduction device 10 Drive signal calculation unit 11 Spatial frequency domain discrete inverse Fourier transform unit 12 Time frequency domain discrete inverse Fourier transform unit 20 X-axis direction wave number calculation unit 21 Wave number calculation unit 22 Reproduced sound field angular spectrum calculation unit 23 Desired sound field angular spectrum calculation unit 24 Divider unit 30 Calculation unit 31 Memory 40 Sound pressure angular spectrum calculation unit 41 Matrix H calculation unit 42 Vector P^ calculation unit 43 Multiplication unit 50 Wall 100, 100-1, 100-2, 100-3, 100-4 Speaker array 101 Infinite linear sound source (x _s , y _s ) Virtual sound source coordinates s(n) Sound source signal y _w Wall position y _des Desired boundary y _ref Reproduced boundary H(k _x ) Y-axis direction particle velocity transmittance distribution angular spectrum M Number of speaker units Δx Speaker unit spacing f _s Sampling frequency F Number of DFT points k Wave number k _x Wave number Δk _x in the x-axis direction Discretization interval ω of k _x Angular frequency D^(k _x , y ₀ , z ₀ , ω) Angular spectrum of drive signal P~(k _x , y _des ), P~(k _x , y _des , z ₀ , ω) Desired sound field angular spectrum G^(k _x , y _ref - y ₀ , z ₀ , ω) Reproduced sound field angular spectrum

Claims

A sound field reproduction device that generates a drive signal for reproducing a sound field formed by a desired sound source using a speaker array consisting of a plurality of speaker units,
The speaker array is assumed to be arranged on y= _y0 at z= _z0 in the xyz space, k is the wave number, _kx is the wave number in the x-axis direction, _Δkx is the discretization interval of the wave number _kx , ω is the angular frequency, c is the speed of sound, ( _xs , _ys ) are virtual sound source coordinates which are the coordinates of the desired sound source, s(n) is the sound source signal of the desired sound source for each time sample n, _yw is the wall position, H( _kx ) is the particle velocity transmittance distribution angular spectrum in the y-axis direction at the wall position _yw , _ydes is the desired boundary, _yref is the reproduction boundary, M is the number of speaker units, Δx is the speaker unit interval, _fs is the sampling frequency, F is the number of DFT points, _H0 ⁽²⁾ is the zero-order second kind Hankel function, and i is the imaginary unit.
an x-axis wave number calculation unit that calculates the wave numbers k x (= k 1 x , k 2 x , ..., k M _x ) by the formula k _x = 2πm/(ΔxM) using a preset number _M of speaker units, a speaker unit interval _Δx , and an index _m = ⁰ to ^M - ¹ (m is an integer), and calculates the discretization interval Δk _x when the speaker array is composed of the plurality of speaker units that are equally spaced;
a wave number calculation unit that calculates the angular frequency ω by the formula: ω = 2πf _s l/F using the preset sampling frequency f _s , the number of DFT points F, and a frequency index l = 0 to F-1 (l is an integer), and calculates the wave number k by the formula: k = ω/c using the angular frequency ω and the sound speed c;
Using the wave number k _x calculated by the x-axis direction wave number calculation unit, the wave number k calculated by the wave number calculation unit, and the preset reproduction boundary y _ref and value y ₀ , the following equation is calculated:
a reproduced sound field angular spectrum calculation unit that calculates a reproduced sound field angular spectrum by
a desired sound field angular spectrum calculation unit that calculates a desired sound field angular spectrum;
a division unit that obtains the angular spectrum of the drive signal by dividing the desired sound field angular spectrum calculated by the desired sound field angular spectrum calculation unit by the reproduced sound field angular spectrum calculated by the reproduced sound field angular spectrum calculation unit,
The desired sound field angular spectrum calculation unit
calculating a sound pressure angular spectrum P^ _d (kx, yw) at the wall position _yw for each wave number _kx using the virtual sound source coordinates ( _xs , _ys ) and the sound source signal s(n) in a buffering section having _a _{predetermined} buffer length;
Using the wave number k calculated by the wave number calculation unit, the wave number k _x calculated by the x-axis direction wave number calculation unit, and the sound pressure angular spectrum P̂ _d (k _x , y _w ), the following equation is calculated:
The vector P^ is calculated by
Using the wave number k, the wave number _kx , the preset y-axis direction particle velocity transmittance distribution angular spectrum H( _kx ), the desired boundary _ydes , and the wall position _yw , the following equation is calculated:
(In k ⁿ _x and k ^m _x , n and m = 1 to M (n and m are integers))
The angular spectrum H(k ⁿ _x , ^km _x ) is calculated by
formula:
The matrix H is calculated by
a sound field reproducing device, characterized in that the desired sound field angular spectrum is obtained by multiplying the discretization interval Δk _x calculated by the x-axis direction wave number calculation unit, the vector P^, and the matrix H together.

A computer constituting a sound field reproduction device that generates a drive signal for reproducing a sound field formed by a desired sound source using a speaker array consisting of a plurality of speaker units,
The speaker array is assumed to be arranged on y= _y0 at z= _z0 in the xyz space, k is the wave number, _kx is the wave number in the x-axis direction, _Δkx is the discretization interval of the wave number _kx , ω is the angular frequency, c is the speed of sound, ( _xs , _ys ) are virtual sound source coordinates which are the coordinates of the desired sound source, s(n) is the sound source signal of the desired sound source for each time sample n, _yw is the wall position, H( _kx ) is the particle velocity transmittance distribution angular spectrum in the y-axis direction at the wall position _yw , _ydes is the desired boundary, _yref is the reproduction boundary, M is the number of speaker units, Δx is the speaker unit interval, _fs is the sampling frequency, F is the number of DFT points, _H0 ⁽²⁾ is the zero-order second kind Hankel function, and i is the imaginary unit.
an x-axis wave number calculation unit that calculates the wave numbers k x (= k 1 x , k 2 x , ..., k M _x ) by the formula k _x = 2πm/(ΔxM) using the preset number M of speaker units, the speaker unit interval _Δx , and an _index _m = ⁰ to ^M - ¹ (m is an integer), and calculates the discretization interval Δk _x when the speaker array is composed of the plurality of speaker units that are equally spaced;
a wave number calculation unit that calculates the angular frequency ω by the formula: ω = 2πf _s l/F using the preset sampling frequency f _s , the number of DFT points F, and a frequency index l = 0 to F-1 (l is an integer), and calculates the wave number k by the formula: k = ω/c using the angular frequency ω and the sound speed c;
Using the wave number k _x calculated by the x-axis direction wave number calculation unit, the wave number k calculated by the wave number calculation unit, and the preset reproduction boundary y _ref and value y ₀ , the following equation is calculated:
a reproduced sound field angular spectrum calculation unit that calculates a reproduced sound field angular spectrum by
a desired sound field angular spectrum calculation unit that calculates a desired sound field angular spectrum; and
a program for causing a division unit to function as a division unit that obtains an angular spectrum of the drive signal by dividing the desired sound field angular spectrum calculated by the desired sound field angular spectrum calculation unit by the reproduced sound field angular spectrum calculated by the reproduced sound field angular spectrum calculation unit,
The desired sound field angular spectrum calculation unit
calculating a sound pressure angular spectrum P^ _d (kx, yw) at the wall position _yw for each wave number _kx using the virtual sound source coordinates ( _xs , _ys ) and the sound source signal s(n) in a buffering section having _a _{predetermined} buffer length;
Using the wave number k calculated by the wave number calculation unit, the wave number k _x calculated by the x-axis direction wave number calculation unit, and the sound pressure angular spectrum P̂ _d (k _x , y _w ), the following equation is calculated:
The vector P^ is calculated by
Using the wave number k, the wave number _kx , the preset y-axis direction particle velocity transmittance distribution angular spectrum H( _kx ), the desired boundary _ydes , and the wall position _yw , the following equation is calculated:
(In k ⁿ _x and k ^m _x , n and m = 1 to M (n and m are integers))
The angular spectrum H(k ⁿ _x , ^km _x ) is calculated using the formula:
The matrix H is calculated by
a program for calculating the desired sound field angular spectrum by multiplying the discretization interval Δk _x calculated by the x-axis direction wave number calculation unit, the vector P^, and the matrix H together;