JP2010164863A - Sound field measurement system - Google Patents
Sound field measurement system Download PDFInfo
- Publication number
- JP2010164863A JP2010164863A JP2009008340A JP2009008340A JP2010164863A JP 2010164863 A JP2010164863 A JP 2010164863A JP 2009008340 A JP2009008340 A JP 2009008340A JP 2009008340 A JP2009008340 A JP 2009008340A JP 2010164863 A JP2010164863 A JP 2010164863A
- Authority
- JP
- Japan
- Prior art keywords
- transfer function
- sound
- matrix
- calculation
- evaluation
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Pending
Links
Images
Landscapes
- Stereophonic System (AREA)
Abstract
Description
本発明は、音場測定システムに関し、より詳細には、測定音場において複数のスピーカ(音源)から一斉に評価音を出力すると共に、複数のマイクで複数のスピーカより出力された評価音を収録することにより、各スピーカから各マイクに対応するそれぞれの伝達関数を一度の測定で算出することが可能な音場測定システムに関する。 The present invention relates to a sound field measurement system, and more specifically, the evaluation sound is simultaneously output from a plurality of speakers (sound sources) in the measurement sound field, and the evaluation sound output from the plurality of speakers is recorded by a plurality of microphones. Thus, the present invention relates to a sound field measurement system capable of calculating each transfer function corresponding to each microphone from each speaker by a single measurement.
従来より、ある空間(原音場)の音場特性を収録(測定)した後、その空間とは異なる空間(再生音場)において、測定した原音場の音場環境を仮想的に再現する音場再現方法が提案されている(例えば、特許文献1参照)。 Conventionally, after recording (measuring) the sound field characteristics of a certain space (original sound field), a sound field that virtually reproduces the sound field environment of the measured original sound field in a different space (reproduced sound field). A reproduction method has been proposed (see, for example, Patent Document 1).
上述した音場再現方法では、まず、ある空間における各スピーカから各マイクへの全ての伝達関数を原音場の音響特性として求めると共に、異なる空間(再生音場)における各スピーカから各マイクへの全ての伝達関数を再生音場の音響特性として求め、この2つの音響特性に基づいて音場再現を行うための音場再現フィルタを算出する。次に、算出された音場再現フィルタを用いて、再生音場において出力される音響信号に対して音響処理(例えば、FIR(Finite Impulse Response)フィルタを用いた畳み込み演算処理)を行うことによって、再生音場へ出力される出力音に音響処理を施す。このようにして音響処理が施された音響信号を再生音場で出力することによって、聴取者はあたかも原音場で聴取しているかのような感覚を再生音場で得ることが可能となる。 In the sound field reproduction method described above, first, all transfer functions from each speaker to each microphone in a certain space are obtained as the acoustic characteristics of the original sound field, and all from each speaker to each microphone in a different space (reproduction sound field). Is obtained as the acoustic characteristics of the reproduced sound field, and a sound field reproduction filter for reproducing the sound field is calculated based on these two acoustic characteristics. Next, by performing acoustic processing (for example, convolution calculation processing using a FIR (Finite Impulse Response) filter) on the acoustic signal output in the reproduced sound field using the calculated sound field reproduction filter, Acoustic processing is performed on the output sound output to the reproduction sound field. By outputting the acoustic signal subjected to the acoustic processing in this way in the reproduction sound field, the listener can obtain a feeling as if listening in the original sound field in the reproduction sound field.
具体的に、既知の入力信号(入力波形)x(t)を入力(スピーカから出力)させ、この入力に対応する出力信号(出力波形)y(t)を測定する(マイクで収録する)ことにより原音場あるいは再生音場における伝達関数h(t)を求める場合について説明する。 Specifically, a known input signal (input waveform) x (t) is input (output from a speaker), and an output signal (output waveform) y (t) corresponding to this input is measured (recorded with a microphone). The case where the transfer function h (t) in the original sound field or the reproduced sound field is obtained will be described.
上述した入力信号x(t)、伝達関数h(t)、出力信号y(t)のそれぞれをZ変換(離散関数のフーリエ変換)することにより、X(z)、H(z)、Y(z)が求められるとすると、X(z)、H(z)、Y(z)のそれぞれは、
このとき、システムの伝達関数H(z)は既知の入力信号X(z)と出力信号Y(z)とを用いて
式4に示す式は、具体的に以下の2通りの計算方法を用いることにより、求めることができる。 The expression shown in Expression 4 can be obtained by using the following two calculation methods.
(1)時間領域で計算する方法
時間領域で計算する方法とは、入力信号X(z)の逆特性である1/X(z)を求め、この1/X(z)の特性を持つフィルタを出力信号Y(z)に畳み込むことにより(畳み込み演算を行うことによって)、伝達関数H(z)を計算する方法を意味する。この方法を用いて伝達関数h(t)を求めるためには、1/X(z)の特性を持つフィルタが既知であることが必要とされる。例えば、入力信号が時間伸張パルス(TSP:Time stretched pulse)の場合は、入力信号の逆特性を持つフィルタの時間応答が既知となるため、この方法を用いることができる。
(1) Method of calculating in the time domain The method of calculating in the time domain is to obtain 1 / X (z), which is the inverse characteristic of the input signal X (z), and to have a filter having the 1 / X (z) characteristic. Is transferred to the output signal Y (z) (by performing a convolution operation) to calculate the transfer function H (z). In order to obtain the transfer function h (t) using this method, it is necessary that the filter having the characteristic of 1 / X (z) is known. For example, when the input signal is a time stretched pulse (TSP), this method can be used because the time response of the filter having the inverse characteristics of the input signal is known.
(2)周波数領域で計算する方法
周波数領域で計算する方法とは、式4をフーリエ変換し、周波数領域上で周波数成分毎に、
上述した(1)の方法および(2)の方法のいずれの方法を用いて伝達関数h(t)を求める場合であっても、S個の音源(スピーカ)より出力される音響信号をX(z)(音源1から音源Sまでの各音源より出力される信号をX1(z)〜XS(z)とする)、M個のマイクにより録音される音響信号をY(z)(マイク1からマイクMまでの各マイクにおいて録音される信号をY1(z)〜YM(z)とする)とし、それぞれの音響信号X(z)および音響信号Y(z)に基づいて求められる伝達関数をH(z)(音源jより出力される音響信号に対するマイクiの伝達関数をHij(z)とする)とすると、音響信号の入出力関係は、次の式で表現することができる。
Even when the transfer function h (t) is obtained using any one of the above-described methods (1) and (2), the acoustic signals output from the S sound sources (speakers) are represented by X ( z) (signals output from each sound source from the
Y(z)=H(z)X(z) ・・・式6
入力データとして、ある音源jについてのみ入力信号を加えた状態を測定すると、音響信号の入出力関係は、次の式で表現することができる。
各音源(j=1からj=SまでのS個の音源)に対応する伝達関数Hij(z)を順番に測定することによって、音場全体の伝達関数h(t)を測定することが可能となる。 It is possible to measure the transfer function h (t) of the entire sound field by sequentially measuring the transfer function H ij (z) corresponding to each sound source (S sound sources from j = 1 to j = S). It becomes possible.
ところで、従来の伝達関数の測定方法では、1つの音源より出力された音響信号を、複数のマイクで同時に収録する作業に基づいて、スピーカ(音源)毎の伝達関数を求める。このように、スピーカ(音源)毎に順番に伝達関数を求めることによって、全ての音源に対応する伝達関数を求める方法を採用している。 By the way, in the conventional transfer function measurement method, the transfer function for each speaker (sound source) is obtained based on the operation of simultaneously recording the acoustic signals output from one sound source with a plurality of microphones. In this way, a method is used in which transfer functions corresponding to all sound sources are obtained by obtaining transfer functions in order for each speaker (sound source).
このため、例えば、測定用のマイクを人間の耳位置に取り付けて、複数のスピーカ(音源)とマイクとの間の伝達関数を測定する場合、聴取者(人間)が姿勢を完全に静止させようとしても、姿勢変化を完全に抑制することが現実には困難であった。従って、スピーカ毎に1つずつ音響信号の出力を行っている間に、マイクの位置が変化してしまう恐れが高く、複数の音源に対して完全に同じ測定条件(この場合はマイク位置が変動しないこと)で測定データを得ることが極めて困難であるという問題があった。 For this reason, for example, when a measurement microphone is attached to the position of a person's ear and a transfer function between a plurality of speakers (sound sources) and the microphone is measured, the listener (human) will make the posture completely stationary. Even so, it was actually difficult to completely suppress the posture change. Therefore, there is a high possibility that the position of the microphone will change while outputting one acoustic signal for each speaker, and the same measurement conditions (in this case, the microphone position fluctuates) for a plurality of sound sources. There is a problem that it is extremely difficult to obtain measurement data.
さらに、非常に多くのスピーカを使用してスピーカ毎の伝達関数を測定する場合には、全てのスピーカの伝達関数を測定するための時間が、スピーカの設置数に比例して極めて長くなってしまうという問題があった。例えば、1つのスピーカについての測定時間が数秒の時間を要するとしても、スピーカの設置数が60以上の場合には、総測定時間が数分単位で必要になる。 Furthermore, when measuring the transfer function for each speaker using a very large number of speakers, the time for measuring the transfer function of all the speakers becomes extremely long in proportion to the number of speakers installed. There was a problem. For example, even if the measurement time for one speaker requires several seconds, if the number of speakers installed is 60 or more, the total measurement time is required in units of several minutes.
また、多数のスピーカ(音源)より同時に多数の音響信号を出力させ、出力された音響信号を、複数のマイクで同時に録音することにより、全てのスピーカ(音源)およびマイクに対応する伝達関数を同時に求める方法も考えられるが、従来の伝達関数の測定方法では、複数のスピーカ(音源)で出力され、複数のマイクにおいて収録された音響信号に基づいて、スピーカ(音源)およびマイクに対応する伝達関数成分を、スピーカ(音源)毎に分解する方法が確立されていなかった。このため、多数のスピーカ(音源)より同時に多数の音響信号を出力させ、出力された音響信号を、複数のマイクで同時に録音することにより、全てのスピーカ(音源)およびマイクに対応する伝達関数を同時に高精度に測定する方法を用いることができないという問題があった。 In addition, by outputting a large number of acoustic signals simultaneously from a large number of speakers (sound sources) and recording the output acoustic signals simultaneously with a plurality of microphones, transfer functions corresponding to all the speakers (sound sources) and microphones can be simultaneously performed. The transfer function measurement method can be considered, but in the conventional transfer function measurement method, the transfer function corresponding to the speaker (sound source) and the microphone is output based on the acoustic signals output from the plurality of speakers (sound source) and recorded in the plurality of microphones. A method for decomposing components for each speaker (sound source) has not been established. For this reason, by outputting a large number of acoustic signals simultaneously from a large number of speakers (sound sources) and recording the output acoustic signals simultaneously with a plurality of microphones, transfer functions corresponding to all the speakers (sound sources) and microphones can be obtained. At the same time, there is a problem that a method for measuring with high accuracy cannot be used.
本発明は、上記問題に鑑みて成されたものであり、多数のスピーカ(音源)より同時に音響信号を出力させ、出力された音響信号を、複数のマイクで同時に収録することにより、全てのスピーカ(音源)およびマイクに対応する伝達関数を同時かつ高精度に測定することが可能な音場測定システムを提供することを課題とする。 The present invention has been made in view of the above-described problem, and by outputting sound signals simultaneously from a large number of speakers (sound sources) and recording the output sound signals simultaneously with a plurality of microphones, It is an object of the present invention to provide a sound field measurement system that can simultaneously and accurately measure a transfer function corresponding to a (sound source) and a microphone.
上記課題を解決するために、本発明に係る音場測定システムは、特定の音響空間に設定される複数のスピーカより評価音を同時に出力し、前記音響空間に設置される複数のマイクより前記評価音を同時に収録することによって、各スピーカから各マイクに対応する伝達関数を算出するための音場測定システムであって、前記スピーカより出力するための評価音の生成を行う評価音生成手段と、前記マイクを用いて収録された収録音と前記評価音生成手段により生成された評価音とに基づいて、各スピーカから各マイクに対応する伝達関数を算出する伝達関数算出手段とを備え、前記伝達関数算出手段は、前記収録音と前記評価音とに基づいて、時間領域における反復解法を適用することにより、前記伝達関数の算出処理を実行し、前記反復解法の適用において用いる相関演算処理および畳み込み演算処理のみ、周波数領域への変換処理によって演算処理を行った後に前記時間領域に演算結果を逆変換させる演算方法を用い、前記評価音生成手段は、前記評価音におけるスピーカ毎の畳み込み行列を行列要素とし、各スピーカに対応する前記行列要素を備えた第1行列と、該第1行列の転置行列を成す第2行列との積の逆行列が存在し得るように評価音を生成することを特徴とする。 In order to solve the above-described problem, the sound field measurement system according to the present invention simultaneously outputs evaluation sounds from a plurality of speakers set in a specific acoustic space, and the evaluations from a plurality of microphones installed in the acoustic space. A sound field measuring system for calculating a transfer function corresponding to each microphone from each speaker by simultaneously recording sound, and an evaluation sound generating means for generating an evaluation sound for output from the speaker; Transfer function calculating means for calculating a transfer function corresponding to each microphone from each speaker based on the recorded sound recorded using the microphone and the evaluation sound generated by the evaluation sound generating means, The function calculation means executes the transfer function calculation process by applying an iterative solution in the time domain based on the recorded sound and the evaluation sound, and the iterative solution Only the correlation calculation process and the convolution calculation process used in the application of the calculation method using the calculation method of inversely converting the calculation result to the time domain after performing the calculation process by the conversion process to the frequency domain, the evaluation sound generating means, There may be an inverse matrix of a product of a first matrix having the matrix element corresponding to each speaker and a second matrix constituting a transposed matrix of the first matrix, with the convolution matrix for each speaker in the sound as a matrix element. As described above, the evaluation sound is generated.
このように、本発明に係る音場測定システムでは、反復解法を用いて伝達関数を算出する場合において、伝達関数算出手段が、相関演算処理と畳み込み演算処理とを行う場合にのみ、周波数領域における演算結果を利用して時間領域における解を求めるので、演算処理の負担軽減と演算処理に利用されるメモリ容量の低減とを図ることが可能となる。 As described above, in the sound field measurement system according to the present invention, when the transfer function is calculated using the iterative solution method, the transfer function calculating unit performs the correlation calculation process and the convolution calculation process only in the frequency domain. Since the solution in the time domain is obtained using the calculation result, it is possible to reduce the burden of the calculation process and reduce the memory capacity used for the calculation process.
さらに、評価音生成手段が、評価音におけるスピーカ毎の畳み込み行列を行列要素とし、各スピーカに対応する行列要素を備えた第1行列と、第1行列の転置行列を成す第2行列との積の逆行列が存在し得るように評価音を生成するため、伝達関数算出手段が、複数のスピーカより出力された評価音からなる複数要素の入力行列と、複数のマイクより収録された収録音からなる複数要素の出力行列とに基づいて各伝達関数を求める場合において、相関演算処理と畳み込み演算処理とを行う際に対応する伝達関数の値を一意的に求めることが可能となる。 Further, the evaluation sound generating means uses a convolution matrix for each speaker in the evaluation sound as a matrix element, and a product of a first matrix having a matrix element corresponding to each speaker and a second matrix forming a transposed matrix of the first matrix. In order to generate an evaluation sound so that there can be an inverse matrix of the transfer function, the transfer function calculating means uses a multi-element input matrix composed of evaluation sounds output from a plurality of speakers and a recording sound recorded from a plurality of microphones. When obtaining each transfer function based on the output matrix of a plurality of elements, it is possible to uniquely obtain the value of the corresponding transfer function when performing the correlation calculation process and the convolution calculation process.
このため、複数のスピーカより評価音を同時に出力し、前記音響空間に設置される複数のマイクで評価音を同時に収録する場合であっても、伝達関数算出手段において、各スピーカから各マイクに対応する伝達関数をそれぞれ算出することが可能となる。 For this reason, even when the evaluation sound is simultaneously output from a plurality of speakers and the evaluation sound is simultaneously recorded by a plurality of microphones installed in the acoustic space, the transfer function calculation means supports each microphone from each speaker. Each transfer function to be calculated can be calculated.
また、上述した音場測定システムにおいて、前記伝達関数算出手段は、前記伝達関数の算出処理において、前記伝達関数を求めるための計算式に対して時間窓を構成する対角行列を掛け合わせ、該対角行列の要素を変更することにより、前記複数のスピーカから連続して評価音を出力すると共に前記複数のマイクから連続して収録音の収録を行った場合における任意の時間範囲の伝達関数を算出するものであってもよい。 Further, in the sound field measurement system described above, the transfer function calculation means multiplies a calculation formula for obtaining the transfer function by a diagonal matrix constituting a time window in the transfer function calculation process, By changing the elements of the diagonal matrix, the evaluation function is output continuously from the plurality of speakers, and the transfer function in an arbitrary time range when recording the recorded sound is continuously recorded from the plurality of microphones. It may be calculated.
このように、本発明に係る音場測定システムでは、長時間連続して同時に複数のスピーカから評価音を出力すると共に、スピーカから出力された評価音を連続してマイクで同時に測定した場合であっても、対角行列要素による時間重みを適用して任意の時間範囲における伝達関数を算出することができる。このため、連続的に収録音の収録位置が音響空間において変更される場合であっても、時間変化に応じた適切な伝達関数を求めることが可能となる。 As described above, in the sound field measurement system according to the present invention, the evaluation sound is output from a plurality of speakers simultaneously for a long time, and the evaluation sound output from the speakers is continuously measured simultaneously by a microphone. However, a transfer function in an arbitrary time range can be calculated by applying time weights based on diagonal matrix elements. For this reason, even when the recording position of the recorded sound is continuously changed in the acoustic space, it is possible to obtain an appropriate transfer function corresponding to the time change.
さらに、上述した音場測定システムにおいて、前記評価音生成手段は、前記評価音における低域周波数帯域の出力値が高域周波数帯域の出力値よりも高くなるように、評価音の周波数特性を変更することにより、前記マイクにより収録される収録音の低域周波数帯域におけるS/Nの改善を図るものであってもよい。 Furthermore, in the sound field measurement system described above, the evaluation sound generation means changes the frequency characteristic of the evaluation sound so that the output value of the low frequency band in the evaluation sound is higher than the output value of the high frequency band. Thus, the S / N in the low frequency band of the recorded sound recorded by the microphone may be improved.
一般的に、特定の音響空間においてスピーカより出力される音をマイクを用いて収録する場合には、マイクで収録された収録音に対してノイズ成分が含まれる場合がある。特に、マイクを使用して音響特性を測定する場合には、低域周波数であるほど、レベルの高い環境ノイズが含まれる傾向がある。 Generally, when a sound output from a speaker in a specific acoustic space is recorded using a microphone, a noise component may be included in the recorded sound recorded by the microphone. In particular, when measuring acoustic characteristics using a microphone, the lower the frequency, the higher the level of environmental noise tends to be included.
このため、本発明に係る音場測定システムでは、評価音における低域周波数帯域の出力値が高域周波数帯域の出力値よりも高くなるように、評価音の周波数特性を変更することにより、低域周波数帯域に含まれるノイズの出力値を評価音の出力値に対して相対的に低くすることができ、収録音の低域周波数帯域におけるS/Nの改善を図ることが可能となる。 For this reason, in the sound field measurement system according to the present invention, by changing the frequency characteristic of the evaluation sound so that the output value of the low frequency band in the evaluation sound is higher than the output value of the high frequency band, The output value of noise included in the frequency band can be made relatively lower than the output value of the evaluation sound, and the S / N in the low frequency band of the recorded sound can be improved.
また、上述した音場測定システムにおいて、前記伝達関数算出手段は、周波数領域への変換処理によって演算処理を行った後に前記時間領域に演算結果を逆変換させて求められた前記相関演算処理または前記畳み込み演算処理の演算結果のうち、前記時間領域のみで前記相関演算処理または前記畳み込み演算処理を行った場合に得られる演算結果との誤差が生じ得ない範囲のデータのみを用いて、前記時間領域における前記伝達関数の算出処理を行うものであってもよい。 Further, in the sound field measurement system described above, the transfer function calculating means performs the calculation process by the conversion process to the frequency domain and then reversely converts the calculation result to the time domain, or the correlation calculation process or the Of the calculation results of the convolution calculation process, only the data in a range in which no error can occur with the calculation result obtained when the correlation calculation process or the convolution calculation process is performed only in the time domain, the time domain is used. The transfer function may be calculated in the above.
このように、上述した音場測定システムにおいて、伝達関数算出手段が、周波数領域から時間領域に逆変換させた相関演算処理または前記畳み込み演算処理の演算結果のうち、時間領域のみで演算処理を行った場合に得られる演算結果との誤差が生じ得ない範囲のデータのみを用いて、伝達関数の算出処理を実行することによって、最終的に算出される伝達関数の値に誤差等が含まれてしまうことを効果的に防止することが可能となる。 Thus, in the sound field measurement system described above, the transfer function calculation means performs the calculation process only in the time domain among the calculation results of the correlation calculation process or the convolution calculation process in which the frequency domain is inversely transformed into the time domain. When the transfer function calculation process is executed using only data in a range that does not cause an error from the calculation result obtained in the case of the It can be effectively prevented.
また、上述した音場測定システムにおいて、前記誤差が生じない範囲は、前記評価音の残響時間に応じて予め設定される有限の応答時間長さに基づいて決定されるものであってもよい。 In the sound field measurement system described above, the range in which the error does not occur may be determined based on a finite response time length set in advance according to the reverberation time of the evaluation sound.
このように、誤差が生じない範囲が、評価音の残響時間に応じて予め設定される有限の応答時間長さに基づいて決定される場合には、算出された伝達関数の性能(品質)を十分に確保しつつ、最終的に算出される伝達関数の値に誤差等が含まれてしまうことを抑制することが可能となる。 As described above, when the range in which no error occurs is determined based on the finite response time length set in advance according to the reverberation time of the evaluation sound, the performance (quality) of the calculated transfer function is It is possible to suppress an error or the like from being included in the finally calculated transfer function value while ensuring sufficiently.
本発明に係る音場測定システムによれば、評価音生成手段が、評価音におけるスピーカ毎の畳み込み行列を行列要素とし、各スピーカに対応する前記行列要素を備えた第1行列と、該第1行列の転置行列を成す第2行列との積の逆行列が存在し得るように評価音を生成するため、伝達関数算出手段が、複数のスピーカより出力された評価音からなる複数要素の入力行列と、複数のマイクより収録された収録音からなる複数要素の出力行列とに基づいて各伝達関数を求める場合において、相関演算処理と畳み込み演算処理とを行う際に対応する伝達関数の値を一意的に求めることが可能となる。 According to the sound field measurement system of the present invention, the evaluation sound generation means uses the convolution matrix for each speaker in the evaluation sound as a matrix element, the first matrix including the matrix element corresponding to each speaker, and the first matrix In order to generate an evaluation sound so that an inverse matrix of a product with the second matrix constituting the transpose matrix of the matrix may exist, the transfer function calculation means includes an input matrix of a plurality of elements composed of the evaluation sounds output from the plurality of speakers. And the transfer function values based on the multi-element output matrix consisting of recorded sounds from multiple microphones, the values of the corresponding transfer functions must be unique when performing correlation calculation processing and convolution calculation processing. Can be obtained.
このため、複数のスピーカより評価音を同時に出力し、前記音響空間に設置される複数のマイクより前記評価音を同時に収録する場合であっても、伝達関数算出手段において、各スピーカから各マイクに対応する伝達関数をそれぞれ算出することが可能となる。 For this reason, even when the evaluation sound is simultaneously output from a plurality of speakers and the evaluation sound is simultaneously recorded from a plurality of microphones installed in the acoustic space, the transfer function calculation means causes each speaker to each microphone. Corresponding transfer functions can be calculated respectively.
以下、本発明に係る音場測定システムについて、図面を用いて詳細に説明を行う。 Hereinafter, a sound field measurement system according to the present invention will be described in detail with reference to the drawings.
[実施の形態1]
図1は、本実施の形態1に係る音場測定システムの概略構成を示した図である。音場測定システム1は、再現を望む音響空間(原音場)あるいは、原音場の音響環境の再現を行う音響空間(再生音場)において設置されるシステムである。
[Embodiment 1]
FIG. 1 is a diagram showing a schematic configuration of a sound field measurement system according to the first embodiment. The sound
音場測定システム1は、図1に示すように、音響空間2に設置されるスピーカ部3と、スピーカ部3に対して音声出力を行う音声出力装置4と、音声出力装置4に対して評価音を出力するコンピュータ5と、音響空間2に設置されるマイクアレイ7と、マイクアレイ7により収録された収録音の入力を行う音声入力装置8とによって概略構成されている。
As shown in FIG. 1, the sound
スピーカ部3は、図1に示すように、マイクアレイ7を中心として等距離を保つように円形に配置された複数(例えばS個)のスピーカ9で構成されている。このスピーカ9は、音声出力装置4より出力された評価音(音響信号)を出力することが可能となっている。
As shown in FIG. 1, the
なお、図1に示すスピーカ部3の配置は、再生音場における一般的なスピーカの配置状態として実験室を用いた場合を一例として示している。例えば、本実施の形態1に係る音場測定システム1を原音場に設置する場合には、その原音場の状況に応じてスピーカの設置位置およびスピーカの設置個数が決定される。例えば、原音場として車両の室内空間を想定する場合には、車室の四隅に対して4個のスピーカが均等に設置される場合などが考えられる。
In addition, arrangement | positioning of the
音声出力装置4は、コンピュータ5より出力された評価音を、スピーカ部3のS個のスピーカ9に対して同時に出力する役割を有している。ここで評価音とは、マイクアレイ7および音声入力装置8を介して収録された収録音に基づいて音響空間2の伝達関数(音響特性)を算出するために用いられる出力音であり、一般的には、収録された収録音をインパルス応答として測定することに適した音(例えば、M系列信号や時間伸張パルス(TSP:Time Stretched Pulse))が該当する。なお、本実施の形態において用いられる評価音については、後述する。
The audio output device 4 has a role of simultaneously outputting the evaluation sound output from the
マイクアレイ7は、円形に配設された複数(例えばM個)のマイク10の一群により構成されている。マイクアレイ7は、各スピーカ9から等距離となる音響空間2の中心位置に設置されている。音声入力装置8は、マイクアレイ7によって測定されたスピーカ部3の出力音をコンピュータ5に出力する役割を有している。
The
コンピュータ5は、上述した評価音を生成すると共に、マイクアレイ7により収録された収録音に基づいて音響空間2における伝達関数(伝達関数マトリックス、音響特性)を算出する機能を有している。
The
コンピュータ5は、評価音生成部(評価音生成手段)12と、収録音記録部13と、伝達関数算出部(伝達関数算出手段)14と、伝達関数記録部15と、制御部16とによって概略構成されている。評価音生成部12は、上述したインパルス応答を収録するための評価音を生成する機能を有している。評価音生成部12は、制御部16の制御命令に応じて評価音を生成し、制御部16を介して音声出力装置4に評価音を出力する。
The
収録音記録部13は、マイクアレイ7および音声入力装置8を介して収録されたマイク毎の収録音を記録する機能を有している。伝達関数算出部14は、音声出力装置4により出力された評価音と、収録音記録部13に記録された収録音とに基づいて、スピーカ部3のS個のスピーカ9とマイクアレイ7のM個のマイク10との間におけるインパルス応答をそれぞれ求めることにより、各スピーカと各マイクとに対応した伝達関数を算出する。なお、各スピーカと各マイクとに対応した伝達関数の要素の集合により、音響空間2における全体の音響特性が求められる。
The recorded
伝達関数算出部14は、伝達関数の算出処理を行うCPU(Central Processing Unit)と、伝達関数を算出するためのプログラム等が記録されるROM(Read Only Memory)、ワークエリア等に使用されるRAM(Random Access Memory)等の一般的な演算処理器によって構成されている。
The transfer
伝達関数記録部15は、伝達関数算出部14により求められた伝達関数を記録する機能を有している。伝達関数記録部15に記録される伝達関数は、各スピーカと各マイクとに対応付けてそれぞれ記録される。
The transfer
制御部16は、音響空間2における伝達関数(音響特性)を求めるために必要な処理を行う役割を有している。制御部16は、上述したように評価音生成部12に対して評価音を生成させると共に、生成された評価音を音声出力装置4に出力する役割を有している。また、制御部16は、音声入力装置8より入力される収録音を収録音記録部13に記録させると共に、伝達関数算出部14に伝達関数を算出させて、伝達関数記録部15に記録させる役割を有している。
The
なお、本実施の形態1に係る音場測定システム1では、評価音生成部12において生成された評価音が、スピーカ部3の全て(S個)のスピーカ9より出力される。マイクアレイ7では、全て(S個)のスピーカより出力された出力音を全て(M個)のマイク10で同時に収録し、収録された収録音は、収録音記録部13においてマイク毎に記録される。
In the sound
次に、音声出力装置4より出力される評価音と、収録音記録部13に記録される収録音とに基づいて、伝達関数算出部14で音響空間2の伝達関数を算出する方法について説明する。
Next, a method for calculating the transfer function of the
まず、スピーカ部3より出力される評価音を測定時の既知の入力信号波形(評価音信号、入力信号)をxj(t)(jはスピーカを示し、j=1,2,・・・,S)とし、マイクアレイ7において収録される収録音を既知の出力信号波形(収録音信号、出力信号)をyi(t)(iはマイクを示し、i=1,2,・・・,M)とし、それぞれをz変換で表現したものを、X(z)、Y(z)とする。また、入力信号波形xj(t)および出力信号波形yi(t)により求められる伝達関数をhij(t)とし、z変換で表現されたX(z)およびY(z)に基づいて求められる未知の伝達関数マトリックスをH(z)とする。
First, a known input signal waveform (evaluation sound signal, input signal) at the time of measuring the evaluation sound output from the
音場測定システム1における測定時間を有限とし、入力信号として有限の時間長さNxの信号を用いると、入力信号波形xj(t)の値が0でない時間範囲は、0≦t≦Nx−1となる。
When the measurement time in the sound
また、伝達関数hij(t)の時間応答を有限の応答時間長さLで近似できると仮定し、hij(t)の値が0でない時間範囲を、0≦t≦L−1であるとする。このとき、出力信号波形yi(t)の値が入力信号波形xj(t)に依存する時間応答も有限となり、その時間は、0≦t≦Nx+L−2(入力信号波形xj(t)の0でない時間範囲(Nx−1)とhij(t)の値が0でない時間範囲(L−1)とを加えた時間範囲)の範囲となる。
Further, assuming that the time response of the transfer function h ij (t) can be approximated by a finite response time length L, the time range where the value of h ij (t) is not 0 is 0 ≦ t
従って、音場測定システム1において測定する出力信号波形yi(t)の測定サンプル数をNyとすると、測定サンプル数Nyは、出力信号波形yi(t)の応答時間に基づいてNx+L−1となるように選べばよくなる。測定サンプル数NyとしてNy=Nx+L−1を選ぶことにより、出力信号波形yi(t)の応答時間に対応する値が0以上となる測定サンプルを求めることが可能となる。
Therefore, when the number of measurement samples of the output signal waveform y i (t) measured in the sound
複数のスピーカ9と複数のマイク10とを備えた音場測定システム1における測定系の入出力関係は、以下の関係式で表すことができる。
Y(z)=H(z)X(z) ・・・式12
The input / output relationship of the measurement system in the sound
Y (z) = H (z) X (z)
但し、X(z)は、
また、Y(z)は、
また、H(z)は、
式12の両辺の転置をとると、
Yt(z)=Xt(z)Ht(z) ・・・式19
となり、Xt(z)は、
Xt(z)=(X1(z) X2(z) ・・・ XS(z))
・・・ 式20
Yt(z)は、
Yt(z)=(Y1(z) Y2(z) ・・・ YM(z))
・・・ 式21
Ht(z)は、
Y t (z) = X t (z) H t (z) Expression 19
And X t (z) is
X t (z) = (X 1 (z) X 2 (z)... X S (z))
... Equation 20
Y t (z) is
Y t (z) = (Y 1 (z) Y 2 (z)... Y M (z))
... Formula 21
H t (z) is
ここで、z変換同士の積は、畳み込み演算マトリックスとベクトルの積として表すことができ、式19は、等価な連立一次方程式である
D=AQ ・・・式23
で表すことができる。但し、
従って、伝達関数H(z)を求めることは、上述した連立一次方程式D=AQを解く問題に帰着することができる。
Here, the product of z-transforms can be expressed as a product of a convolution matrix and a vector, and Equation 19 is an equivalent simultaneous linear equation D = AQ Equation 23
Can be expressed as However,
Therefore, obtaining the transfer function H (z) can be reduced to the problem of solving the simultaneous linear equation D = AQ described above.
音場測定システム1において、図2に示すように、測定対象である伝達関数H(z)が線形の伝達関数モデルである限り、D=AQを満たす解は、必ず存在する(但し、等号が成立するのは、測定データ中にノイズがない場合となる)。このため、D=AQの解を求めることによって伝達関数H(z)を算出できるようにするためには、D=AQを満たす解が1つだけ存在するように方程式を構築する必要がある。
In the sound
具体的に、D=AQを満たす解が1つだけ存在するように方程式を構築するためには、AtAの逆行列が存在するように入力信号波形x(t)の組を与えることが必要とされる。つまり、入力信号波形(評価音信号、入力信号)xj(t)の畳み込み行列(式26)を行列要素とし、各スピーカに対応する行列要素(式26)を備えた行列A(式25:第1行列)と、この行列A(第1行列)の転置行列を成す行列At(第2行列)との積であるAtAの逆行列が存在し得るように入力信号波形(評価信号、評価音)xj(t)が設定されればよい。 Specifically, in order to build the equations as solutions satisfying D = AQ there is only one, is to provide a set of input signal waveform x (t) so that there is an inverse matrix of A t A Needed. In other words, the matrix A (Equation 25: Eq. 25) having a matrix element (Equation 26) corresponding to each speaker is a convolution matrix (Equation 26) of the input signal waveform (evaluation sound signal, input signal) x j (t). Input signal waveform (evaluation signal) so that an inverse matrix of A t A which is the product of the first matrix) and a matrix A t (second matrix) that forms a transposed matrix of the matrix A (first matrix) can exist. , Evaluation sound) x j (t) may be set.
AtAの逆行列が存在する場合、Qは、D=AQの最小二乗解として求められ、次の連立一次方程式
(AtD)=(AtA)Q ・・・式28
の解として、次式で求めることができる。
Q=(AtA)−1(AtD) ・・・式29
When an inverse matrix of A t A exists, Q is obtained as a least squares solution of D = AQ, and the following simultaneous linear equation (A t D) = (A t A) Q Equation 28
Can be obtained by the following equation.
Q = (A t A) −1 (A t D) Equation 29
上述した計算原理を用いることにより、複数の音源を同時に再生した場合の測定データを用いて伝達関数マトリックスのすべての成分を正確に計算することが可能となる。このため、1回の測定ですべてのスピーカ・マイク間の伝達関数を求めることが可能になる。 By using the calculation principle described above, it is possible to accurately calculate all the components of the transfer function matrix using the measurement data when a plurality of sound sources are reproduced simultaneously. For this reason, it becomes possible to obtain | require the transfer function between all the speakers and microphones by one measurement.
具体的な数値解を得るためには、式23および式28を具体的に解く数値計算アルゴリズムおよび入力信号の与え方が重要になる。一般に式23および式28で示される方程式は、計算規模(演算量)が極めて大きく、汎用の行列解法を用いる場合には、極めて処理能力の高い計算環境を用意する必要があるため、方程式の演算を実行することが困難である。 In order to obtain a specific numerical solution, it is important to provide a numerical calculation algorithm that specifically solves Equations 23 and 28 and to provide an input signal. In general, the equations shown in Expression 23 and Expression 28 have a very large calculation scale (computation amount), and when a general-purpose matrix solution is used, it is necessary to prepare a calculation environment with extremely high processing capacity. Is difficult to perform.
具体的な数値を挙げて説明すると、式28の係数行列AtAの計算規模は、音源数(スピーカ9数)をS、伝達関数の時間応答長さをLとした場合、S×L行、S×L列の行列を計算するための計算規模が必要となる。伝達関数を求める場合、時間応答長さLは、通常、数千〜数万のオーダーで設定されるため、仮に時間応答長さLを215=32768に設定し、データ精度を単精度実数(4byte)とすると、係数行列の演算処理に用いられる記憶容量だけでも、S×S×4Gbyte程度のメモリ容量が必要となる。従って、16音源を一度に測定したい場合には、1000Gbyte程度のメモリ容量が係数記憶のためだけに必要になる。さらに、このような大規模の方程式を解くための演算コストや計算時間も膨大なものとなってしまう。 Explaining with specific numerical values, the calculation scale of the coefficient matrix A t A in Equation 28 is S × L rows when the number of sound sources (the number of speakers 9) is S and the time response length of the transfer function is L. , A calculation scale for calculating a matrix of S × L columns is required. When obtaining the transfer function, the time response length L is usually set in the order of several thousand to several tens of thousands. Therefore, the time response length L is set to 2 15 = 32768, and the data accuracy is a single precision real number ( 4 bytes), a memory capacity of about S × S × 4 Gbytes is required only for the storage capacity used for the coefficient matrix calculation process. Therefore, if it is desired to measure 16 sound sources at once, a memory capacity of about 1000 Gbytes is required only for coefficient storage. Furthermore, the calculation cost and calculation time for solving such a large-scale equation become enormous.
そのため、本実施の形態1に係る伝達関数算出部14では、畳み込み問題の性質に特化して効率的な計算を行う計算アルゴリズムを用いている。図3〜図6は、本実施の形態1における計算アルゴリズムを示したフローチャートである。
For this reason, the transfer
図3は、D=AQを解くための数学計算方式の概要を示した図である。図3に示すように、まず、伝達関数算出部14は、方程式D=AQを最小二乗法で解くために、P=RQに基づいて
P=AtA
R=AtD
とおいて、PおよびRの計算を行い(ステップS.1)、この計算式に基づいて連立一次方程式P=RQを解くことにより、時間領域における伝達関数マトリックス(行列Q)(Q=R−1P)を求める(ステップS.2)
FIG. 3 is a diagram showing an outline of a mathematical calculation method for solving D = AQ. As shown in FIG. 3, first, the transfer
R = A t D
Then, P and R are calculated (step S.1), and a simultaneous linear equation P = RQ is solved based on this calculation formula, whereby a transfer function matrix (matrix Q) in the time domain (Q = R −1 P) is obtained (step S.2).
以下、伝達関数算出部14が、時間領域において伝達関数マトリックス(行列Q)を求める処理を説明する。伝達関数算出部14は、時間領域において
Q=(AtA)−1(AtD) ・・・式29
の最適解を求めるために、反復解法を用いるものとし、特に本実施の形態1では反復解法の解法例として共役勾配法を用いて演算を行うものとする。
Hereinafter, a process in which the transfer
It is assumed that an iterative solution method is used in order to obtain the optimal solution, and in particular, in
伝達関数マトリックス(行列Q)を算出するためには、上述したように、
Q=(AtA)−1(AtD) ・・・式29
を解く必要がある。ここで、反復解法を用いて伝達関数マトリックス(行列Q)を算出する場合、反復計算の処理において必要とされる行列形式は、規則的な畳み込み行列になるという特徴がある。このため、本実施の形態1に係る伝達関数算出部14では、この畳み込み行列を利用した特徴的な演算処理を行うことによって、計算過程に必要なメモリ量をフィルタ段数Lに比例するオーダーまで削減することが可能となる。このように、連立一次方程式を解く大きな枠組みとして、反復法(共役勾配法)を用いることによって、係数行列の規則性を保ち、メモリの使用量を抑えたままでの計算を可能にすることができる。
In order to calculate the transfer function matrix (matrix Q), as described above,
Q = (A t A) −1 (A t D) Equation 29
Need to be solved. Here, when the transfer function matrix (matrix Q) is calculated using the iterative solution method, the matrix format required in the iterative calculation process is a regular convolution matrix. For this reason, the transfer
また、本実施の形態1に係る計算アルゴリズムでは、係数行列AtAが相関行列となるため、方程式の組み立て処理の段階で相関行列の規則性を利用することにより、計算に必要な記憶容量をL×Lに比例するオーダーから、Lに比例するオーダーへと大幅に削減することができる。 Further, the calculation algorithm according to the first embodiment, since the coefficient matrix A t A is the correlation matrix, by utilizing the regularity of the correlation matrix at the stage of the assembly process of the equation, the storage capacity required for the calculation It is possible to greatly reduce the order proportional to L × L to the order proportional to L.
従って、方程式の組み立て処理と反復計算の処理とに用いられる主要な行列演算を、高速フーリエ変換(FFT)を利用した相関演算処理と畳み込み演算処理とで演算することによって、計算量を大幅に削減することができ、音場測定システムの規模が大きくなっても、一般的なコンピュータを用いて演算処理を行うことが可能となる。 Therefore, the amount of calculation is greatly reduced by calculating the main matrix operations used in the equation assembly process and the iterative calculation process by the correlation calculation process using fast Fourier transform (FFT) and the convolution calculation process. Therefore, even if the scale of the sound field measurement system increases, it becomes possible to perform arithmetic processing using a general computer.
なお、上述する計算アルゴリズムにおいて入力信号波形として入力される評価音(入力信号)には、連立一次方程式
(AtD)=(AtA)Q ・・・式28
の係数行列AtAの逆行列が存在するような信号系列を選ぶことが必要である。
Note that the evaluation sound (input signal) input as the input signal waveform in the calculation algorithm described above includes simultaneous linear equations (A t D) = (A t A) Q (Equation 28)
It is necessary to select a signal sequence in which an inverse matrix of the coefficient matrix A t A exists.
このように、係数行列AtAの逆行列が存在するような信号系列を選ぶことによって、原理上、計算アルゴリズムにおいて解の精度が劣化することを防止できる。なお、有限桁数の数値計算で有限の計算時間で解を求める関係上、係数行列が単位行列に近くなるように入力信号を選ぶことが望ましい。係数行列が単位行列に近いほど、連立一次方程式を反復解法で解く場合の収束回数が少なくなるので、有限の反復回数で計算処理を打ち切った場合に得られる解の精度を向上させることができる。 Thus, by selecting a signal sequence in which an inverse matrix of the coefficient matrix A t A exists, in principle, it is possible to prevent the accuracy of the solution from being deteriorated in the calculation algorithm. Note that it is desirable to select an input signal so that the coefficient matrix is close to the unit matrix because a solution is obtained in a finite calculation time by numerical calculation with a finite number of digits. The closer the coefficient matrix is to the unit matrix, the smaller the number of times of convergence when solving simultaneous linear equations by the iterative method, so that the accuracy of the solution obtained when the calculation process is terminated with a finite number of iterations can be improved.
このような入力信号として、正規分布や一様分布などの自己相関特性が鋭くて(自己相関波形が時刻0でのみ大きな値を持ち、その他の時刻では0に近い値を持つもの)、相互相関の低い信号を用いることが好ましい。従って、評価音生成部12では、このような条件を満たし得る入力信号(評価音)の生成を行う。
As such an input signal, autocorrelation characteristics such as normal distribution and uniform distribution are sharp (the autocorrelation waveform has a large value only at
図4は、図3に示したステップS.1の処理内容である方程式を組み立てるための計算処理を示したフローチャートである。また、図5および図6は、図3に示したステップS.2の処理内容である方程式を解く計算処理を示したフローチャートであり、共役勾配法を使用して伝達関数マトリックス(行列Q)を算出する処理を示している。 FIG. 4 shows the steps S.E. It is the flowchart which showed the calculation process for assembling the equation which is the processing content of 1. FIG. 5 and 6 are the same as those shown in FIG. 2 is a flowchart showing calculation processing for solving an equation which is the processing content of 2, and shows processing for calculating a transfer function matrix (matrix Q) using a conjugate gradient method.
フローチャートに従って時間領域だけで演算処理を行う場合には、次述する(1)〜(4)の演算処理において処理負担が増大するという問題があった。
(1)行列R:R=AtAの演算処理
(2)行列P:P=AtDの演算処理
(3)初期値計算におけるRとベクトルとの積Rq0の演算処理
(4)反復計算におけるRとベクトルとの積Rσμの演算処理
また、行列Rの演算に要する記憶容量も増大するという問題があった。
When the arithmetic processing is performed only in the time domain according to the flowchart, there is a problem that the processing load increases in the arithmetic processing of (1) to (4) described below.
(1) Matrix R: R = A t A arithmetic processing (2) Matrix P: P = A t D arithmetic processing (3) Initial value calculation R and vector product Rq 0 arithmetic processing (4) iteration There is a problem that the calculation processing of the product Rσ μ of R and the vector in the calculation and the storage capacity required for the calculation of the matrix R increase.
本実施の形態1に係る計算アルゴリズムでは、上述したように、(1)〜(4)の演算を行う場合に、高速フーリエ変換(FFT)を用いて周波数領域における演算処理を行い、演算処理された演算結果を逆高速フーリエ変換(IFFT)によって時間領域の演算結果に逆変換する(時間領域の演算結果に戻す)ことによって、演算処理負担の軽減を図っている。 In the calculation algorithm according to the first embodiment, as described above, when performing the calculations (1) to (4), the calculation process is performed in the frequency domain using the fast Fourier transform (FFT). The calculation processing load is reduced by inversely transforming the calculation result into the time domain calculation result by inverse fast Fourier transform (IFFT) (returning to the time domain calculation result).
まず、図4に示すフローチャートに基づいて、初期演算として行列Rと行列Pとを計算する場合について説明する。 First, a case where the matrix R and the matrix P are calculated as an initial calculation will be described based on the flowchart shown in FIG.
行列Rは、上述したようにR=AtAで示される。
ここで行列Rの部分行列Ri,jは、相関行列であって規則正しいパターンを持ち、斜め方向の要素の値が同一の値になるという特徴を有している。この部分行列Ri,jでは、2L−1個の異なる値を備えている。このため、係数行列Rの演算を行う場合には、重複する値をメモリ等の記憶領域に記憶せず、異なる値のみをベクトルとして記憶することによって、伝達関数算出部14の演算処理で使用されるメモリ使用量をLのオーダーで削減することが可能となる。
The matrix R is represented by R = A t A as described above.
Here, the partial matrix R i, j of the matrix R is a correlation matrix, has a regular pattern, and has the characteristic that the values of the elements in the oblique direction are the same value. This submatrix R i, j has 2L−1 different values. For this reason, when calculating the coefficient matrix R, the overlapping values are not stored in a storage area such as a memory, but only different values are stored as vectors. The amount of memory used can be reduced on the order of L.
また、行列Ri,jの各要素ri,jは、伝達関数同士の相関演算の形を含んでいる。このため、相関演算処理を行う場合、高速フーリエ変換(FFT)を用いてri,jの値を周波数領域の値に変換して解を求めた後に、求められた解を逆フーリエ変換(IFFT)により逆変換して、時間領域におけるri,jの相関演算処理の解を求める。このように周波数領域において演算処理を行うことによって、処理負担を軽減させることができる。 Each element r i, j of the matrix R i, j includes a form of correlation calculation between transfer functions. For this reason, when performing the correlation calculation process , the value of r i, j is converted into a frequency domain value using fast Fourier transform (FFT), and then a solution is obtained. Then, the obtained solution is subjected to inverse Fourier transform (IFFT). ) To obtain a solution of r i, j correlation calculation processing in the time domain. By performing the arithmetic processing in the frequency domain in this way, the processing load can be reduced.
また、行列Pは、上述したようにP=AtDで示され、行列Pi,jの各要素pi,jも、行列Rと同様に、伝達関数同士の相関演算の形を含んでいるため、pi,jに対して高速フーリエ変換を適用し、pi,jの値を周波数領域の値に変換して解を求め、求められた解を逆フーリエ変換して時間領域の解に逆変換することにより値を求めることができる。このように周波数領域の解を時間領域の解に逆変換する処理を用いることにより、演算処理の負担を軽減させることができる。 Further, the matrix P is represented by P = A t D as described above, and each element p i, j of the matrix P i, j includes the form of the correlation operation between the transfer functions similarly to the matrix R. Therefore , fast Fourier transform is applied to p i, j , the value of p i, j is converted to a frequency domain value to obtain a solution, and the obtained solution is inverse Fourier transformed to obtain a time domain solution. The value can be obtained by inversely transforming to. Thus, by using the process of inversely transforming the frequency domain solution into the time domain solution, it is possible to reduce the burden of the arithmetic processing.
ここで、伝達関数算出部14は、RおよびPの高速フーリエ変換を行う場合における段数の数Ncorrとして、Ncorr≧Ny+L−1を満たす値、例えば、Ncorr=Ny+Lを設定する(ステップS.21)。
Here, the transfer
そして、伝達関数算出部14は、入力信号波形であるxj(t)(但しj=1,2,・・・,S)および出力信号波形であるyi(t)(但しi=1,2,・・・,M)に対してNcorr点の高速フーリエ変換を適用する(ステップS.22)。
Then, the transfer
なお、xj(t)およびyi(t)における時系列のデータ数がNcorrに満たない場合、伝達関数算出部14は、各データの後ろに0を追加することにより高速フーリエ変換を実行する。このように0を追加してデータ長を増やしてから高速フーリエ変換を行うことにより、後で逆フーリエ変換(IFFT)を適用することにより求められた時間領域における演算結果が、求めたい範囲に対応する値となるように調整することができる。このように調整を行うことにより、共役勾配法を用いて求められる伝達関数に誤差が生じてしまうことを防止することができる。
When the number of time series data in x j (t) and y i (t) is less than Ncorr, the transfer
高速フーリエ変換を適用することにより、
A(ω)≡Xt(ω)=(X1(ω) X2(ω) …XS(ω))
・・・式30
D(ω)≡Yt(ω)=(Y1(ω) Y2(ω) …YM(ω))
・・・式31
となる。
By applying the fast Fourier transform,
A (ω) ≡X t (ω) = (X 1 (ω) X 2 (ω)... X S (ω))
... Formula 30
D (ω) ≡Y t (ω) = (Y 1 (ω) Y 2 (ω)... Y M (ω))
... Formula 31
It becomes.
次に、伝達関数算出部14は、求められたA(ω)およびD(ω)を用いて、R(ω)およびP(ω)の周波数領域における相関演算を実施する(ステップS.23)。
式30および式31より、R(ω)およびP(ω)は、
R(ω)=A*(ω)A(ω) ・・・式32
P(ω)=A*(ω)D(ω) ・・・式33
で示すことができる。ここで、A*(ω)は、A(ω)の共役転置行列を示しており、
A*(ω)は、
From Equation 30 and Equation 31, R (ω) and P (ω) are
R (ω) = A * (ω) A (ω) Equation 32
P (ω) = A * (ω) D (ω) Equation 33
Can be shown. Where A * (ω) represents the conjugate transpose matrix of A (ω),
A * (ω) is
また、R(ω)の要素ri,j(ω)は、
ri,j(ω)=conj(Xi(ω))Xj(ω)
で表すことができ、P(ω)の要素pi,j(ω)は、
pi,j(ω)=conj(Xi(ω))Yj(ω)
で表すことができる。
なお、conj(X)はXの複素共役を示し、conj(Y)はYの複素共役を示している。
The element r i, j (ω) of R (ω) is
r i, j (ω) = conj (X i (ω)) X j (ω)
The element p i, j (ω) of P (ω) is
p i, j (ω) = conj (X i (ω)) Y j (ω)
Can be expressed as
Note that conj (X) represents the complex conjugate of X, and conj (Y) represents the complex conjugate of Y.
このように、フーリエ変換されたri,j(ω)およびpi,j(ω)の値を周波数成分の掛け算を用いて求めた後に、求めた解を逆フーリエ変換して時間領域の時系列データに変換することによって、時間領域において演算を行う場合に比べて迅速かつ容易に解を求めることができる。 As described above, after the values of r i, j (ω) and p i, j (ω) subjected to Fourier transformation are obtained by multiplication of frequency components, the obtained solution is subjected to inverse Fourier transformation to obtain the time domain time. By converting to series data, a solution can be obtained quickly and easily as compared with the case of performing computation in the time domain.
伝達関数算出部14は、求められたri,j(ω)の値を逆フーリエ変換(IFFT)することにより、時系列データであるri,j(t)に変換する(ステップS.24)。この場合、伝達関数算出部14は、逆フーリエ変換によって求められた解が正確な値となるように(時間領域のみの演算により求められる解とに誤差が生じないように)、−L+1≦t≦L−1の範囲のデータのみを用いて解を求める。つまり、誤差の生じない範囲は、出力音の残響時間等に応じて予め設定される有限の応答時間長さLに応じて決定される。
The transfer
伝達関数算出部14は、−L+1≦t≦L−1の範囲のデータのうち、t<0のときに、ri,j(t+Ncorr)を負の時刻のri,j(t)として所定の記憶手段に記憶し、t≧0のときに、ri,j(t)をそのまま記憶する。
Predetermined transfer
一方で、伝達関数算出部14は、求められたpi,j(ω)の値を逆フーリエ変換(IFFT)することにより、時系列データであるpi,j(t)に変換する(ステップS.25)。この場合、伝達関数算出部14は、逆フーリエ変換によって求められた解が正確な値となるように、0≦t≦L−1の範囲のデータのみを用いて解を求め、記憶する。
On the other hand, the transfer
次に、図5および図6に示すフローチャートに基づいて、伝達関数算出部14の具体的な伝達関数の演算処理について説明する。
Next, a specific transfer function calculation process of the transfer
伝達関数算出部14は、伝達関数の具体的な演算処理においても、行列Ri,jの各要素ri,jに含まれる伝達関数同士の相関演算の形を利用する。伝達関数算出部14は、相関演算処理を行う場合において、高速フーリエ変換(FFT)を用いてri,jの値を周波数領域の値に変換して解を求めた後に、求められたri,jを周波数領域において演算処理し、演算処理された演算結果を逆フーリエ変換(IFFT)を用いて逆変換して、時間領域におけるri,jの相関演算処理の解を求める。このように周波数領域において演算処理を行うことによって、処理負担を軽減させることができる。
The transfer
まず、伝達関数算出部14は、Rに対して高速フーリエ変換を行う場合の段数をNCGとして、NCG≧2L−1を満たす値、例えば、NCG=2Lを設定する(ステップS.31)。
First, the transfer
そして、伝達関数算出部14は、図4のステップS.24において求められたri,j(t)を並べ直したベクトルr’i,j(t) 0≦t≦NCG−1を定義する(ステップS.32)。具体的にr’i,j(t)は、tの値に応じて、下記のように示される。
伝達関数算出部14は、r’i,j(t)のNCG点における高速フーリエ変換(FFT)を行うことにより、周波数領域におけるRCG(ω)のij成分を計算する(ステップS.33)。
The transfer
次に、伝達関数算出部14は、初期値としてjに1を代入し、μに0を代入する(ステップS.34)、ここでjは、音響空間2におけるスピーカ9の数を示している。またμは、行列Qのj列目の値を算出するために実行された演算回数(反復回数)を示すパラメータであり、具体的には、ステップS.36〜ステップS.38の処理を繰り返し実行した回数を示している。
Next, the transfer
次に伝達関数算出部14は、初期設定を行う(ステップS.35)。初期設定では、P,Qのj列目の値をp,qとおき、適当な初期値q0を用いて、
ε0=p−Rq0 ・・・式37
を求める。なお、εμは、反復回数μでの暫定的な解をqμとしたときの、誤差ベクトルを示したものである。
Next, the transfer
Ask for. Note that ε μ represents an error vector when a provisional solution at the number of iterations μ is q μ .
ここで。ε0の算出に用いられる行列Rと任意のベクトルVとの演算処理(例えば、Rq0の演算)について、図7に示すフローチャートに沿って説明する。 here. A calculation process (for example, calculation of Rq 0 ) between the matrix R used for calculating ε 0 and an arbitrary vector V will be described with reference to the flowchart shown in FIG.
係数行列Rに対する積算対象となる任意のベクトルVを
このようにして示されるベクトルVのそれぞれを、ステップS.31において設定したNCG点で高速フーリエ変換(FFT)すると、
このようにして求められたV(ω)とステップS.33において求められたRCG(ω)との積を、周波数領域において算出すると、
例えば、Rq0の演算を行う場合には、積算対象となるベクトルq0をベクトルVと置き換えて、図7に示すフローチャートの計算処理を行うことにより、周波数領域における演算処理を用いてRq0の演算を行うことができる。 For example, when performing the calculation of Rq 0 is a vector q 0 as the integration target is replaced with the vector V, by performing the calculation processing of the flowchart shown in FIG. 7, the Rq 0 by using an arithmetic processing in the frequency domain Arithmetic can be performed.
具体的に伝達関数算出部14は、q0をNCG点で高速フーリエ変換(FFT)してq0(ω)を計算する。伝達関数算出部14では、ステップS.33において計算されたRCG(ω)と、q0(ω)とを周波数領域において積算することにより、RCG(ω)q0(ω)を算出し、算出されたRCG(ω)q0(ω)の算出結果を時間領域に変換(逆フーリエ変換)することにより(Rq0)’である
(Rq0)’=IFFT(RCG(ω)q0(ω)) ・・・式42
を求める。そして、伝達関数算出部14は、求められた(Rq0)’の内のS個の逆フーリエ変換(IFFT)結果に対してそれぞれ0≦t≦L−1の時間窓をかけることによりRq0を求めることができる。このように、伝達関数算出部14において、周波数成分における演算処理を用いて計算を行うことにより、伝達関数算出部14における処理負担を軽減させることができる。
Specifically, the transfer
Ask for. Then, the transfer
また、伝達関数算出部14は、
さらに、伝達関数算出部14は、共役勾配法の収束判定に用いられるパラメータをηとして設定する(ステップS.35)。ηは収束を判断するための閾値パラメータとしての定数であって、このパラメータに基づいて誤差が収束判定条件を満たす値となった場合(次式45を満たす場合)、伝達関数算出部14は、解が収束したと判断して反復計算を終了する。
Furthermore, the transfer
そして、伝達関数算出部14は、収束条件の判断処理(ステップS.36)を行う。収束条件の判断式は、
収束条件を満たさなかった場合(ステップS.36においてNoの場合)、伝達関数算出部14は、Rσμの計算を行う(ステップS.37)
When the convergence condition is not satisfied (No in step S.36), the transfer
Rσμの計算を行う場合にも、図7に示すフローチャートを用いて説明したように、行列Rとベクトルとの積を求める計算に該当する畳み込み演算処理を用いて計算を行うことが可能である。 Even when Rσ μ is calculated, as described with reference to the flowchart shown in FIG. 7, the calculation can be performed using the convolution operation processing corresponding to the calculation for obtaining the product of the matrix R and the vector. .
具体的に伝達関数算出部14は、σμをNCG点で高速フーリエ変換(FFT)してσμ(ω)を計算する。伝達関数算出部14では、ステップS.33において計算されたRCG(ω)と、σμ(ω)とを周波数領域において積算することにより、RCG(ω)σμ(ω)を算出し、算出されたRCG(ω)σμ(ω)の算出結果を時間領域に変換(逆フーリエ変換)することにより、
(Rσμ)’=IFFT(RCG(ω)σμ(ω)) ・・・式46
を求める。
Specifically, the transfer
(Rσ μ ) ′ = IFFT (R CG (ω) σ μ (ω)) Equation 46
Ask for.
この場合においても、伝達関数算出部14は、時間領域において解に誤差が生じないように、0≦t≦L−1の範囲のデータのみを用いて解を求める。従って、伝達関数算出部14は、求められた(Rσμ)’の内のS個の逆フーリエ変換(IFFT)結果に対してそれぞれ0≦t≦L−1の時間窓をかけることによりRq0を求める。このように、周波数成分における演算処理を用いて計算を行うことにより、伝達関数算出部14における処理負担を軽減させることができる。
Also in this case, the transfer
そして、伝達関数算出部14は、求められたRσμの値と、εμ、σμとを用いて、
そして、伝達関数算出部14は、ステップS.38において算出されたεμとPとを用いて、収束条件の判断処理(ステップS.36)を繰り返し実行する。
Then, the transfer
収束条件を満たす場合(ステップS.36においてYesの場合)、伝達関数算出部14は、ステップS.38において算出されたqμの値を、行列Qにおけるj列目の解(伝達関数)とする(ステップS.39)。
When the convergence condition is satisfied (Yes in step S.36), the transfer
その後、伝達関数算出部14は、全てのjについて、つまり全てのスピーカ9に対応する各マイク10の伝達関数マトリックス(行列Q)における全ての伝達関数を求めたか否かを判断し(ステップS.40)、全ての列において伝達関数を求めていない場合(ステップS.40においてNoの場合)には、jに1を加え、μに0を代入した後に(ステップS.41)、処理をステップS.35に移行して上述した処理を繰り返し実行する。
Thereafter, the transfer
全ての列において伝達関数を求めた場合(ステップS.40においてYesの場合)、伝達関数算出部14は、全ての列の伝達関数の値を求めたものと判断して、方程式を解く計算処理(伝達関数マトリックス(行列Q)の算出処理)を終了する。
When transfer functions are obtained for all columns (Yes in step S.40), the transfer
このようにして伝達関数算出部14が、複数のスピーカ9より出力される入力信号波形に関するデータxj(t)と、複数のマイクにより収録される出力信号波形に関するデータyi(t)とに基づいて伝達関数の算出を行うことによって、同時に複数のスピーカ9から評価音を出力すると共に、スピーカ9から出力された評価音をマイク10で同時に測定することにより、各スピーカ9から各マイク10に対応する伝達関数を同時に算出することが可能となる。このため、従来のように、評価音をそれぞれのスピーカ9より順番に出力することにより、出力先となるスピーカ9と収録先となるマイク10とを対応させて伝達関数を算出する必要がなくなり、伝達関数の測定時間の短縮化を図ることが可能となる。
In this way, the transfer
また、一度に全てのスピーカから評価音を出力させて、同時に全てのマイクで収録を行うことにより伝達関数の算出を行うため、従来のように集音時におけるマイク位置の変動などによる影響を受けにくくなり、測定精度の向上を図ることが可能となる。 In addition, since the evaluation function outputs the evaluation sound from all the speakers at the same time and records with all the microphones at the same time, the transfer function is calculated, so that it is affected by fluctuations in the microphone position during sound collection as before. It becomes difficult to improve measurement accuracy.
[実施の形態2]
次に、実施の形態2に係る音場測定システムについて説明を行う。
[Embodiment 2]
Next, the sound field measurement system according to
実施の形態1に係る音場測定システム1では、複数のスピーカ9から同時に評価音を出力すると共に、複数のマイク10から同時に評価音の収録を行うことにより、一度の測定でそれぞれのマイク10とスピーカ9とに対応させた伝達関数を求める方法について説明を行った。実施の形態1に示した方法を用いることにより、一度の測定で、各スピーカ9から各マイク10に対応する伝達関数を求めることができるので、マイク位置の変動(聴取者の耳位置のぐらつき等)があっても高い検出精度で伝達関数の測定を行うことが可能となる。
In the sound
一方で、スピーカ9からマイク10までの伝達関数は、そのスピーカ9に対するマイク10の設置位置に応じて変化する。従って、特定の空間における伝達関数を算出する場合、通常は、聴取者の聴取位置を予め決定し、決定された聴取位置で最適な音響効果を得ることができるように、決定位置にマイク10を設置して収録作業を行うことが一般的である。
On the other hand, the transfer function from the speaker 9 to the
一方で、再現音場において聴取者の聴取位置(つまりマイク10の設置位置)が固定されておらず、移動(変化)する場合などにおいては、その変動位置に応じて、個別に伝達関数の測定を行う必要があった。特に、従来の方法では、測定位置において各スピーカ9より順番に評価音を出力させることによって、各スピーカ9から各マイク10に対応する伝達関数を算出する必要が有るため、連続的に測定位置を変更させて伝達関数の測定を行うことが困難であった。
On the other hand, when the listener's listening position (that is, the installation position of the microphone 10) is not fixed in the reproduced sound field and moves (changes), the transfer function is measured individually according to the fluctuation position. Had to do. In particular, in the conventional method, it is necessary to calculate the transfer function corresponding to each
しかしながら、実施の形態1に示した方法を用いることにより、一度の測定で、各スピーカから各マイクに対応する全ての伝達関数を求めることができるため、測定位置を連続的に移動させることにより、より具体的には、長時間の連続した測定データの中から、一部のデータを用いて伝達関数を求めることにより、伝達関数の時間的変化を求めることが可能となる。実施の形態2では、このような時間的な伝達関数の変化を求めるために、時間重みを適用した場合における伝達関数の算出方法について説明する。 However, by using the method shown in the first embodiment, it is possible to obtain all transfer functions corresponding to each microphone from each speaker in one measurement, so by continuously moving the measurement position, More specifically, it is possible to obtain a temporal change of the transfer function by obtaining a transfer function using a part of data from long-time continuous measurement data. In the second embodiment, a method for calculating a transfer function in the case where a time weight is applied in order to obtain such a change in temporal transfer function will be described.
なお、実施の形態2に係る音場測定システムは、実施の形態1において図1を用いて説明した音場測定システム1と同一の構成であり、同一の機能を有するシステムを用いることができるため、実施の形態2において同一符号を用いると共に、その詳細な説明は省略する。また、時間重みを適用した具体的な伝達関数の算出は、実施の形態1に示した伝達関数算出部14で行われるものとし、その詳細な構成についての説明は省略する。以下、伝達関数算出部14における処理内容について説明を行う。
The sound field measurement system according to the second embodiment has the same configuration as the sound
測定データの全てを用いた関係式は、式19に示すように
Yt(z)=Xt(z)Ht(z)
で表すことができ、式23に示すように、この式と等価な連立一次方程式を示す
D=AQ
を解くことにより、伝達関数Ht(z)を求めることができる。
The relational expression using all of the measurement data is as follows: Y t (z) = X t (z) H t (z)
As shown in Equation 23, a simultaneous linear equation equivalent to this equation is shown. D = AQ
, The transfer function H t (z) can be obtained.
このようにして求められる伝達関数が、長時間の連続した測定で得られた測定データの一部を用いることにより求められる場合には、長時間の測定により得られた測定データのうち、任意の時間範囲を解析する処理が必要とされる。このため、本実施の形態2に係る伝達関数算出部14では、時間窓を表す対角行列Wを用いることにより、式23に示す方程式を求めて、任意の時間範囲の伝達関数を求める。
When the transfer function obtained in this way is obtained by using a part of the measurement data obtained by continuous measurement for a long time, any measurement data obtained by the long-time measurement is arbitrarily selected. Processing to analyze the time range is required. For this reason, in the transfer
具体的には、式23に示す連立一次方程式の両辺に対して対角行列W
WD=(WA)Q ・・・式51
を解くことにより、所定の時間範囲に対応する伝達関数を求めることができる。
なお、この対角行列Wは、対角成分が0である行は、伝達関数を求めるための計算において無視され、その行の測定データ(その時刻の測定データ)が使用されないことを意味する。
Specifically, the diagonal matrix W is applied to both sides of the simultaneous linear equations shown in Equation 23.
WD = (WA) Q Formula 51
By solving, a transfer function corresponding to a predetermined time range can be obtained.
The diagonal matrix W means that a row having a diagonal component of 0 is ignored in the calculation for obtaining the transfer function, and the measurement data of that row (measurement data at that time) is not used.
WD=(WA)Q ・・・式51
の最小二乗解は、
AtW2D=(AtW2A)Q ・・・式52
の解として、
Q=(AtW2A)−1(AtW2D) ・・・式53
と表すことができる。
この式53に示す連立方程式の数値解を得るための計算手順を図8〜図12に示す。
WD = (WA) Q Formula 51
The least squares solution of
A t W 2 D = (A t W 2 A) Q ···
As a solution of
Q = (A t W 2 A) −1 (A t W 2 D) Expression 53
It can be expressed as.
A calculation procedure for obtaining a numerical solution of the simultaneous equations shown in Equation 53 is shown in FIGS.
図8〜図12の計算手順は、連立一次方程式を解く大きな枠組みとして共役勾配法を採用し、個々の計算過程で必要な演算を畳み込み問題に特化して効率化する方法を用いたものとなっている。共役勾配法の実施で問題となる係数行列とベクトルの演算については、図8に示す畳み込み行列Aと任意のベクトルVとの積の計算と、図9に示す行列Aの転置行列Atと任意のベクトルVとの積の計算を用いることにより、高速フーリエ変換(FFT)を利用した畳み込み演算または相関演算で効率的に計算する方法を採用している。
まず、図8に示す、行列AとベクトルVとの積の計算について説明する。
The calculation procedures in FIGS. 8 to 12 employ a conjugate gradient method as a large framework for solving simultaneous linear equations, and use a method for improving efficiency by specializing a convolution problem in each calculation process. ing. Regarding the calculation of the coefficient matrix and the vector, which is a problem in the implementation of the conjugate gradient method, the calculation of the product of the convolution matrix A and the arbitrary vector V shown in FIG. 8, the transposed matrix At of the matrix A shown in FIG. By using the calculation of the product with the vector V, a method of efficiently calculating by a convolution operation or a correlation operation using fast Fourier transform (FFT) is adopted.
First, calculation of the product of the matrix A and the vector V shown in FIG. 8 will be described.
行列Aに対する積算対象となる任意のベクトルVを
このようにして示されるVのそれぞれに対して、Nconv≧Nx+L−1の条件を満たすNconvを決定する(ステップS.62)。本実施の形態2では、例えば、Nconv=Nx+Lを用いる。 For each V indicated in this way, Nconv satisfying the condition of Nconv ≧ N x + L−1 is determined (step S.62). In the second embodiment, for example, Nconv = N x + L is used.
次に、伝達関数算出部14は、入力信号波形であるxj(t)(但しj=1,2,・・・,S)に対してNconv点の高速フーリエ変換を適用する(ステップS.63)。
Next, the transfer
なお、xj(t)における時系列のデータ数がNconvに満たない場合、伝達関数算出部14は、データの後ろに0を追加することにより高速フーリエ変換を実行する。このように0を追加してデータ長を増やしてから高速フーリエ変換を行うことにより、後で逆フーリエ変換(IFFT)を適用することにより求められる時間領域の演算結果が、求めたい範囲に対応する値となるように調整することができる。
When the number of time-series data in x j (t) is less than Nconv, the transfer
高速フーリエ変換を適用することにより、
A(ω)≡Xt(ω)=(X1(ω) X2(ω) …XS(ω))
・・・式55
となる。
By applying the fast Fourier transform,
A (ω) ≡X t (ω) = (X 1 (ω) X 2 (ω)... X S (ω))
... Formula 55
It becomes.
一方で、伝達関数算出部14は、ステップ61において定義されたベクトルVに対して、Nconv点の高速フーリエ変換を適用する(ステップS.64)。
On the other hand, the transfer
なお、ベクトルVの構成要素であるvj(t)(j=1,2,・・・,S)における時系列のデータ数がNconvに満たない場合、伝達関数算出部14は、データの後ろに0を追加することにより高速フーリエ変換を実行する。このように0を追加することによってデータ長を増やしてから高速フーリエ変換を行うことにより、後で逆フーリエ変換(IFFT)を適用して求められる時間領域の演算結果が、求めたい範囲に対応する値となるように調整することができる。
When the number of time-series data in v j (t) (j = 1, 2,..., S) that is a component of the vector V is less than Nconv, the transfer
ステップS.64においてNconv点で高速フーリエ変換(FFT)すると、
このようにして求められたV(ω)とステップS.63において求められたA(ω)との積を、周波数領域において算出する(周波数領域において畳み込み算出を実施する)と、
そして、伝達関数算出部14は、このようにしても求められた(AV)(ω)の演算結果に対して逆フーリエ変換を適用することにより、時間領域におけるデータ(時系列データ)AVに変換を行う(ステップS.66)。
Then, the transfer
このように、行列AおよびベクトルVをそれぞれ高速フーリエ変換することにより周波数領域におけるデータの変換を行い、変換された周波数領域において畳み込み演算を行った後に、演算結果を逆フーリエ変換により時間領域に戻して演算を行うことにより、伝達関数算出部14における演算処理の負担を軽減させることができる。
In this way, the data in the frequency domain is converted by performing fast Fourier transform on each of the matrix A and the vector V. After performing the convolution operation in the converted frequency domain, the calculation result is returned to the time domain by inverse Fourier transform. By performing the calculation in this way, it is possible to reduce the burden of calculation processing in the transfer
次に、図9に示す、転置行列AtとベクトルVとの積の計算について説明する。 Next, FIG. 9, illustrating the calculation of a product of a transposed matrix A t and the vector V.
転置行列Atに対する積算対象となる任意のベクトルVを
このようにして示されるVのそれぞれに対して、Ncorr≧Nx+L−1の条件を満たすNcorrを決定する(ステップS.72)。本実施の形態2では、例えば、Ncorr=Nx+Lを用いる。 Ncorr satisfying the condition of Ncorr ≧ N x + L−1 is determined for each of V shown in this way (step S.72). In the second embodiment, for example, Ncorr = N x + L is used.
次に、伝達関数算出部14は、入力信号波形であるxj(t)(但しj=1,2,・・・,S)に対してNcorr点の高速フーリエ変換を適用する(ステップS.73)。
Next, the transfer
なお、xj(t)における時系列のデータ数がNcorrに満たない場合、伝達関数算出部14は、データの後ろに0を追加することにより高速フーリエ変換を実行する。このように0を追加することによってデータ長を増やしてから高速フーリエ変換を行うことにより、後で逆フーリエ変換(IFFT)を適用することにより求められた時間領域における演算結果が、求めたい範囲に対応する値となるように調整することができる。
When the number of time-series data in x j (t) is less than Ncorr, the transfer
高速フーリエ変換を適用することにより、
A(ω)≡Xt(ω)=(X1(ω) X2(ω) …XS(ω))
・・・式59
となる。
By applying the fast Fourier transform,
A (ω) ≡X t (ω) = (X 1 (ω) X 2 (ω)... X S (ω))
... Formula 59
It becomes.
一方で、伝達関数算出部14は、ステップ71において定義されたベクトルVに対して、Ncorr点の高速フーリエ変換を適用する(ステップS.74)。
なお、ベクトルVの構成要素であるvj(t)(j=1,2,・・・,S)における時系列のデータ数がNcorrに満たない場合、伝達関数算出部14は、データの後ろに0を追加することにより高速フーリエ変換を実行する。このように0を追加することによってデータ長を増やしてから高速フーリエ変換を行うことにより、後で逆フーリエ変換(IFFT)を適用することにより求められた時間領域における演算結果が、求めたい範囲に対応する値となるように調整することができる。
On the other hand, the transfer
When the number of time-series data in v j (t) (j = 1, 2,..., S) that is a component of the vector V is less than Ncorr, the transfer
ステップS.74においてNcorr点で高速フーリエ変換(FFT)すると、
このようにして求められたV(ω)と転置行列Atとの積を、周波数領域において算出する(周波数領域において畳み込み算出を実施する)と、
ここで、A*(ω)は、A(ω)の共役転置行列を意味しており、
そして、伝達関数算出部14は、上述した式により算出された(AtV)(ω)の演算結果に対して逆フーリエ変換を適用することにより、時間領域におけるデータ(時系列データ)(AtV)’に変換を行う(ステップS.76)。
Then, the transfer
伝達関数算出部14は、この(AtV)’の内のS個の逆フーリエ変換(IFFT)結果に対してそれぞれ0≦t≦L−1の時間窓をかけることによりAtVを求める(ステップS.76)。このように、0≦t≦L−1の範囲の時間窓を利用して、対応するデータのみを用いて計算を行うことによって、逆フーリエ変換によって求められた解が正確な値となるように(時間領域のみの演算により求められる解とに誤差が生じないように)することが可能となる。
The transfer
上述した図8および図9に示す処理手順に基づいて演算処理を行うことによって、伝達関数算出部14における演算処理負担の低減を図ることが可能となる。
By performing the arithmetic processing based on the processing procedure shown in FIG. 8 and FIG. 9 described above, it is possible to reduce the arithmetic processing burden in the transfer
次に、図10〜図12を用いて、長時間の連続した測定で得られた測定データの一部を用いることにより伝達関数を求める処理内容について説明する。 Next, processing contents for obtaining a transfer function by using a part of measurement data obtained by continuous measurement for a long time will be described with reference to FIGS.
上述したように、長時間の連続した測定で得られた測定データの一部を用いることにより伝達関数を求めるためには、
AtW2D=(AtW2A)Q ・・・式52
の解である、
Q=(AtW2A)−1(AtW2D) ・・・式53
を求めることが必要となる。
As described above, in order to obtain a transfer function by using a part of measurement data obtained by continuous measurement for a long time,
A t W 2 D = (A t W 2 A) Q ···
Is the solution of
Q = (A t W 2 A) −1 (A t W 2 D) Expression 53
Is required.
伝達関数算出部14は、式52に示す方程式を最小二乗法で解くために、式52の左辺をP、右辺の係数行列をRとし、P=RQに基づいて
P=At(W2D)
R=At(W2A)
とおいて、PおよびRの計算を行い、この計算式に基づいて連立一次方程式P=RQを解くことにより、時間領域における行列Q:Q=R−1Pを求める。
In order to solve the equation shown in
R = A t (W 2 A)
Then, P and R are calculated, and the simultaneous linear equation P = RQ is solved based on this calculation formula to obtain a matrix Q: Q = R −1 P in the time domain.
なお、
Q=(AtW2A)−1(AtW2D) ・・・式53
の連立一次方程式を求める場合、Wが単位行列でない場合には、式53の係数行列であるAtW2Aの規則性が崩れてしまう。このため、伝達関数算出部14では、式53の演算を行う際に、AtW2Aを相関関数のベクトルとして記憶させて演算を行うことができない。
In addition,
Q = (A t W 2 A) −1 (A t W 2 D) Expression 53
When W is not a unit matrix, the regularity of A t W 2 A, which is the coefficient matrix of Expression 53, is lost. For this reason, the transfer
簡単な数値例を挙げて説明すると、
x1(t)=1,2,3,4,0・・・
x2(t)=11,12,13,14,0・・・
w(t)=0,0,0,1,1,1,0・・・
とし、
x 1 (t) = 1, 2, 3, 4, 0.
x 2 (t) = 11, 12, 13, 14, 0...
w (t) = 0,0,0,1,1,1,0 ...
age,
しかしながら、時間重みのある場合には、
このため、本実施の形態2に係る伝達関数算出部14では、伝達関数の計算処理において、係数行列AtW2Aそのものを求めることではなく、(AtW2A)vの形式の行列のベクトル積の演算結果を求める方法を用いている。(AtW2A)vをAt(W2(Av))とする計算は、後述する図12の計算手順と、図8および図9に示した一連の計算手順を用いることにより、畳み込み演算、時間重み、相関演算を順番に実施することで効率的に行うことが可能になる。従って、実施の形態2に示す数値計算方法を用いることによって、係数行列の記憶に関する問題を回避することができる。以下、上述した計算方法を用いて演算を行う手順について説明する。
For this reason, the transfer
まず、伝達関数算出部14は、図10に示すように、P=At(W2D)の計算を行う(ステップS.81)。
First, as shown in FIG. 10, the transfer
具体的に、伝達関数算出部14は、図11に示すように、初期値としてiに1を代入する(ステップS.91)、ここでiは、音響空間2におけるマイク10の数を示しており、1に設定されたiがマイク10の設置数Mに該当するまでの間、伝達関数算出部14はステップS.92〜ステップS.95の処理を繰り返し実行する。
Specifically, as shown in FIG. 11, the transfer
次に伝達関数算出部14は、行列Dのi列目の要素diをとして抽出する(ステップS.92)。そして、伝達関数算出部14は、対角行列Wとdiとを用いて、W2diの計算を行う(ステップS.93)。具体的に、伝達関数算出部14は、w(t)yi(t)の計算を行う。
Next, the transfer
その後、伝達関数算出部14は、行列Dのi列目の要素diを用いて求められたW2diの計算結果をベクトルVとして、既に図9を用いて説明したように、転置行列AtとベクトルVとの演算方法を用いることにより、At(W2di)の計算を行う(ステップS.94)。
Thereafter, the transfer
そして、伝達関数算出部14は、全てのiに対してAt(W2di)の計算を行ったか否か(つまり、i=Mであるか否か)の判断を行う(ステップS.95)。全てのiに対して計算を行っていないと判断した場合(ステップS.95においてNoの場合)、伝達関数算出部14は、iに1を加えた後に、上述したステップS.92〜ステップS.95の処理を繰り返し実行する(ステップS.96)。
Then, the transfer
一方で、全てのiに対して計算を行ったと判断した場合(ステップS.95においてYesの場合)、伝達関数算出部14は、Dの全ての列に対応するiのAt(W2di)を算出したものと判断し、つまり、At(W2D)(=P)を算出したものと判断して、処理を図10のステップS.81に戻す。
On the other hand, if it is determined that were calculated for all i (in the case of Yes in step S.95), the transfer
次に、伝達関数算出部14は、初期値としてjに1を代入し、μに0を代入する(ステップS.82)、ここでjは、音響空間2におけるスピーカ9の数を示している。またμは、行列Qのj列目の値を算出するために実行された演算回数(反復回数)を示すパラメータであり、具体的には、ステップS.83〜ステップS.88の処理を繰り返し実行した回数を示している。
Next, the transfer
次に伝達関数算出部14は、初期設定を行う(ステップS.83)。初期設定では、P,Qのj列目の値をp,qとおき、適当な初期値q0を用いて、
ε0=p−At(W2(Aq0)) ・・・式67
を求める。なお、εμは、反復回数μでの暫定的な解をqμとしたときの、誤差ベクトルを示したものである。
Next, the transfer
ε 0 = p−A t (W 2 (Aq 0 )) Equation 67
Ask for. Note that ε μ represents an error vector when a provisional solution at the number of iterations μ is q μ .
ε0=p−At(W2(Aq0)) ・・・式67
におけるAt(W2(Aq0))の計算は、後述するステップS.85に示す
Rσμ=At(W2(Aσμ))
のσμに対して初期値としてq0を設定することによって求めることができる。
Rσμ=At(W2(Aσμ))
の計算方法については、後述する。
なお、p、q、σμは、長さLのベクトルをS個集めた長さSLのベクトルとして定義される。
ε 0 = p−A t (W 2 (Aq 0 )) Equation 67
Calculation of A t (W 2 (Aq 0 )) in the later-described step S. Rσ μ = A t shown in 85 (W 2 (Aσ μ ))
Can be obtained by setting q 0 as an initial value for σ μ .
Rσ μ = A t (W 2 (Aσ μ ))
The calculation method will be described later.
Note that p, q, and σ μ are defined as vectors of length SL obtained by collecting S vectors of length L.
また、伝達関数算出部14は、求められたε0を用いて
さらに、伝達関数算出部14は、共役勾配法の収束判定に用いられるパラメータをηとして設定する。ηは収束を判断するための閾値パラメータとしての定数であって、このパラメータに基づいて誤差が収束判定条件を満たす値となった場合(次述するステップS.84においてYesの場合)、伝達関数算出部14は、解が収束したと判断して反復計算を終了する。
Furthermore, the transfer
そして、伝達関数算出部14は、収束条件の判断処理(ステップS.84)を行う。収束条件の判断式は、
収束条件を満たさなかった場合(ステップS.84においてNoの場合)、伝達関数算出部14は、Rσμの計算を行う(ステップS.85)。
When the convergence condition is not satisfied (No in step S.84), the transfer
図12は、Rσμの計算方法、つまり、At(W2(Aσμ))を計算する方法を示したフローチャートである。 Figure 12 is a method of calculating Arushiguma mu, that is, a flowchart illustrating a method of calculating the A t (W 2 (Aσ μ )).
伝達関数算出部14は、まず、図8に示した行列AとベクトルVとの積の計算方法において、ベクトルVにベクトルσμを用いることにより、Aσμを計算する(ステップS.101)。次に、伝達関数算出部14は、対角行列Wの積であるW2にAσμを積算したW2(Aσμ)の計算を行う(ステップS.102)。
First, in the method for calculating the product of the matrix A and the vector V shown in FIG. 8, the transfer
そして、伝達関数算出部14は、図9に示した転置行列AtとベクトルVとの積の計算方法において、ベクトルVにベクトルW2(Aσμ)を用いることにより、At(W2(Aσμ))を計算する(ステップS.103)。
Then, the transfer
このように、At(W2(Aσμ))の計算を行う場合にも、伝達関数算出部14において、周波数成分における演算処理を用いて計算を行うことにより、伝達関数算出部14における処理負担を軽減させることができる。
Thus, even when the calculation of A t (W 2 (Aσ μ )), the transfer
そして、伝達関数算出部14は、求められたRσμの値と、εμ、σμとを用いて、
そして、伝達関数算出部14は、ステップS.86において算出されたεμとPとを用いて、収束条件の判断処理(ステップS.84)を繰り返し実行する。
Then, the transfer
収束条件を満たす場合(ステップS.84においてYesの場合)、伝達関数算出部14は、ステップS.86において算出されたqμの値を、行列Qにおけるj列目の解(伝達関数)とする(ステップS.87)。
When the convergence condition is satisfied (Yes in step S.84), the transfer
その後、伝達関数算出部14は、全ての列jについて、つまり全てのスピーカに対応する各マイクの伝達関数(伝達関数マトリックス(行列Q)における全ての伝達関数)を求めたか否かを判断し(ステップS.88)、全ての列において伝達関数を求めていない場合(ステップS.88においてNoの場合)には、jに1を加え、μに0を代入した後に(ステップS.89)、処理をステップS.83に移行して上述した処理を繰り返し実行する。
Thereafter, the transfer
全ての列において伝達関数を求めた場合(ステップS.88においてYesの場合)、伝達関数算出部14は、全ての列の伝達関数の値を求めたものと判断して、伝達関数の計算処理(伝達関数マトリックス(行列Q)の算出処理)を終了する。
When transfer functions are obtained for all columns (Yes in step S.88), the transfer
このように、実施の形態2に係る伝達関数算出部14では、長時間連続して同時に複数のスピーカ9から評価音を出力すると共に、スピーカ9から出力された評価音を連続してマイク10で同時に測定した場合であっても、時間重みを適用して伝達関数を算出することができる。このため、連続的に測定位置が変更される場合であっても時間変化に応じた適切な伝達関数を求めることが可能となる。
As described above, in the transfer
なお、実施の形態2では、上述したように、長時間連続して同時に複数のスピーカ9から評価音を同時に出力すると共に、スピーカ9から出力された評価音を連続してマイク10で同時に測定した場合について説明を行ったが、実施の形態2に係る発明は、複数のスピーカ9から同時に評価音を出力する場合に限定されるものではない。例えば、従来の方法を用いて伝達関数を求める場合のように、各スピーカから順番に評価音を出力することにより該当するスピーカとマイクとに対応する伝達関数を算出する場合であっても、各評価音の出力を長時間連続して行うことにより、該当するスピーカと各マイクとに対応する伝達関数の算出を、時間重みを適用して求めることが可能となる。
In the second embodiment, as described above, the evaluation sounds are simultaneously output from the plurality of speakers 9 simultaneously for a long time, and the evaluation sounds output from the speakers 9 are continuously measured by the
[実施の形態3]
実施の形態1に示す方法を用いることにより、複数のスピーカから同時に評価音を出力すると共に、複数のマイクから同時に評価音の収録を行うことによって、一度の測定でそれぞれのマイクとスピーカとに対応させた伝達関数を求めることが可能となる。また、実施の形態2に示す方法を用いることにより、時間重みを適用して伝達関数を算出することができるので、連続的に測定位置が変更される場合であっても時間変化に応じた適切な伝達関数を求めることが可能となる。
[Embodiment 3]
By using the method shown in the first embodiment, the evaluation sound is simultaneously output from a plurality of speakers, and the evaluation sound is simultaneously recorded from the plurality of microphones, thereby corresponding to each microphone and the speaker in one measurement. It is possible to obtain the transferred transfer function. In addition, since the transfer function can be calculated by applying the time weight by using the method shown in the second embodiment, even if the measurement position is continuously changed, it is appropriate to respond to the time change. It is possible to obtain a simple transfer function.
しかしながら、音響空間2にノイズが存在する場合、マイク10で測定された収録音(出力信号、測定データ)Yには、評価音(入力信号)Xの成分(寄与成分)に加えて、ノイズ成分が含まれることになる。
However, when noise is present in the
例えば、図13に示すように、測定データY(z)を、入力信号X(z)と、伝達関数H(z)と、ノイズ成分N(z)とを用いて表すと、
このように、測定データY(z)に、入力信号X(z)の寄与成分と、ノイズ成分N(z)とが含まれている場合において、ノイズの影響を低減させるためには、ノイズ信号のパワーに対して入力信号のパワーを相対的に大きくすることにより、測定S/N(Signal to Noise ratio)を大きくすることが必要となる。 As described above, when the measurement data Y (z) includes the contribution component of the input signal X (z) and the noise component N (z), in order to reduce the influence of noise, the noise signal It is necessary to increase the measurement S / N (Signal to Noise ratio) by increasing the power of the input signal relative to the power of the signal.
ところで、マイク10を使用して、音響空間2における音響特性を測定する場合には、低域周波数であるほど、レベルの高い環境ノイズが含まれる傾向があり、高い周波数に比例して低域周波数での測定S/Nが劣化する傾向がある。このため、スピーカ9より出力される評価音(入力信号)の低域周波数成分のレベルを大きくすることにより、収録音(測定データ)に含まれる低い周波数成分のS/Nを改善することが考えられる。
By the way, when measuring the acoustic characteristics in the
従って、本実施の形態3に係る伝達関数算出部14では、入力信号として低域周波数成分のレベルの高い信号を用いることによって、測定データに含まれる低い周波数成分のS/Nを改善する方法を採用する。
Therefore, the transfer
図14(a)(b)は、入力信号の周波数特性をそれぞれ変えた場合における入力信号の信号パワー値とマイクで集音されるノイズの信号パワーとの対応を示した図である。(a)は、入力信号の周波数特性が平坦な場合(つまり、低域周波数成分の信号パワーの補強を行わない場合)を示しており、(b)は、入力信号の周波数特性を変えることにより、低域周波数成分の信号パワーが大きくなるようにした場合を示している。また、(a)(b)には、周波数帯域に応じて信号パワーが変化するノイズの状態を対比して示している。 14A and 14B are diagrams showing the correspondence between the signal power value of the input signal and the signal power of noise collected by the microphone when the frequency characteristics of the input signal are changed. (A) shows the case where the frequency characteristic of the input signal is flat (that is, the case where the signal power of the low frequency component is not reinforced), and (b) shows that the frequency characteristic of the input signal is changed. This shows a case where the signal power of the low frequency component is increased. Also, (a) and (b) show a comparison of noise states in which the signal power changes according to the frequency band.
図14(a)に示す入力信号とノイズとの信号のパワー比をとると、低域周波数帯域において、ノイズのパワー値と入力信号のパワー値とが近い値となるため、S/Nが小さな値になってしまう。一方で、図14(b)に示す入力信号とノイズとの信号のパワー比をとると、入力信号の信号パワーが周波数の減少にともない増加する状態を示しているため、低域周波数帯域であっても、ノイズのパワー値と入力信号のパワー値との間隔が一定に保たれることになり、S/Nの減少を抑制することが可能となる。 When the power ratio between the input signal and the noise shown in FIG. 14A is taken, the noise power value and the input signal power value are close to each other in the low frequency band, so the S / N is small. It becomes a value. On the other hand, when the power ratio between the input signal and noise shown in FIG. 14B is taken, the signal power of the input signal increases as the frequency decreases. However, the interval between the power value of the noise and the power value of the input signal is kept constant, and it is possible to suppress a decrease in S / N.
しかしながら、低域周波数成分の信号パワーが大きな入力信号を用いると、自己相関特性が劣化して、連立一次方程式の係数行列(AtA)の性質が悪くなるため、連立方程式を反復解法で解くときの収束性が劣化してしまう。このように、解の収束性が悪くなると、反復解法に必要な収束回数が多く必要になり、有限の反復回数で得られる解の精度が劣化してしまう。そこで、実施の形態3に係る伝達関数算出部14では、収束性を改善するために、入力信号(評価音)に対して前処理を行うことにより、収束性の劣化を回避して解の精度を保つ方法が用いられている。
However, if an input signal having a large signal power of low frequency components is used, the autocorrelation characteristic is deteriorated and the property of the coefficient matrix (A t A) of the simultaneous linear equations is deteriorated. Therefore, the simultaneous equations are solved by an iterative method. Convergence at the time will deteriorate. Thus, when the convergence of the solution is deteriorated, the number of convergence required for the iterative solution is increased, and the accuracy of the solution obtained with a finite number of iterations is deteriorated. Therefore, the transfer
下記の式75は、低域成分が高いパワーとなる試験信号を用いる場合において、収束性の改善を実現することが可能な数式を示している。実施の形態3の入力信号Xには、予めマイク10に加えられるノイズの周波数特性を考慮して周波数特性が平坦ではない信号を使用しているものとする。このため、
Yt(z)=Xt(z)Ht(z)
を解く代わりに、この式の両辺にg(z)を畳み込む前処理を実施して、
(g(z)Yt(z))=(g(z)Xt(z))Ht(z) ・・・式75
を解く。
The following Expression 75 shows an expression that can realize improvement in convergence when using a test signal whose low frequency component has high power. It is assumed that a signal having a non-flat frequency characteristic is used as the input signal X of the third embodiment in consideration of the frequency characteristic of noise applied to the
Y t (z) = X t (z) H t (z)
Instead of solving, perform a preprocessing that convolves g (z) on both sides of this equation,
(G (z) Y t (z)) = (g (z) X t (z)) H t (z) Expression 75
Solve.
ここで、
g(z)Xt(z)
=(g(z)X1(z) g(z)X2(z) … g(z)XS(z))
・・・ 式76
g(z)Yt(z)
=(g(z)Y1(z) g(z)Y2(z) … g(z)YM(z))
・・・ 式77
である。
here,
g (z) X t (z)
= (G (z) X 1 (z) g (z) X 2 (z)... G (z) X S (z))
... 76
g (z) Y t (z)
= (G (z) Y 1 (z) g (z) Y 2 (z)... G (z) Y M (z))
... Formula 77
It is.
式75を解くことは、これを代数方程式で表した
GD=(GA)Q ・・・ 式78
を解いて、
Q=((GA)tGA)−1((GA)tD)・・・ 式79
を得ることと同じである。
Solving Equation 75 is expressed as an algebraic equation: GD = (GA) Q Equation 78
Solve
Q = ((GA) t GA) −1 ((GA) t D) Equation 79
Is the same as getting
ここで、Gはg(z)の畳み込みマトリックスを示しており、
式79に示す式の係数行列(GA)t(GA)は、
(GA)t(GA)=At(GtG)A
が単位行列に近づくようなGを与えることで、解の収束性を改善することができる。例えば、前処理に用いるフィルタとして、入力信号Xの周波数特性とほぼ逆の周波数特性を持つ最小位相フィルタを利用することができる。
The coefficient matrix (GA) t (GA) of the equation shown in Equation 79 is
(GA) t (GA) = A t (G t G) A
By giving G such that approaches the unit matrix, the convergence of the solution can be improved. For example, a minimum phase filter having a frequency characteristic almost opposite to the frequency characteristic of the input signal X can be used as a filter used for preprocessing.
一般的に、周波数特性がフラットな信号の相関特性は、周波数特性がフラットでない信号より鋭い相関ピークを持ち、相関特性がよい。また最小位相フィルタを用いることにより前処理フィルタの段数を必要最小限にし、前処理に必要な計算コストの低減を図ることが可能となる。 In general, a correlation characteristic of a signal with a flat frequency characteristic has a sharper correlation peak than a signal with a non-flat frequency characteristic, and the correlation characteristic is good. Further, by using the minimum phase filter, it is possible to minimize the number of stages of the preprocessing filter and to reduce the calculation cost required for the preprocessing.
以上、本発明に係る音場測定システムについて、図面を用いて詳細に説明を行ったが、本発明に係る音場測定システムは、上述したものに限定されるわけではない。当業者であれば、特許請求の範囲に記載された範疇内において、各種の変更例または修正例に想到しうることは明らかであり、それらについても当然に本発明の技術的範囲に属するものと了解される。 The sound field measurement system according to the present invention has been described in detail with reference to the drawings. However, the sound field measurement system according to the present invention is not limited to the above-described one. It will be apparent to those skilled in the art that various changes and modifications can be made within the scope of the claims, and these are naturally within the technical scope of the present invention. Understood.
なお、上述した音場測定システムを用いて、多数のスピーカから一度に評価音を出力し、多数のマイクにおいて同時に評価音の収録を行うことによって、各スピーカから各マイクに対応する伝達関数を一度に算出する方法が、従来のように各スピーカより順番に評価音を出力し、多数のマイクで評価音の収録を行うことにより伝達関数を算出する場合に比べて、測定時間の短縮化を図ることが可能であることを、具体的に説明する。 In addition, by using the above-described sound field measurement system, evaluation sounds are output from a large number of speakers at a time, and evaluation sounds are simultaneously recorded by a large number of microphones, whereby a transfer function corresponding to each microphone is once transmitted from each speaker. Compared to the case where the transfer function is calculated by outputting the evaluation sound in order from each speaker and recording the evaluation sound with a large number of microphones as in the conventional method, the calculation time is shortened. The fact that it is possible will be specifically described.
スピーカ(音源)の数をSとし、伝達関数を測定するための応答時間がLで打ち切られると仮定する。従来のように、各スピーカより順番に評価音を出力し、多数のマイクにおいて同時に評価音の収録を行うことにより伝達関数を算出する場合において、1つのスピーカ(音源)より出力される評価音をマイクで収録して測定を行う時間は、概ね
α×L
程度と見積もることができる。
Assume that the number of speakers (sound sources) is S, and the response time for measuring the transfer function is cut off at L. When the transfer function is calculated by outputting the evaluation sound sequentially from each speaker and recording the evaluation sound simultaneously in a number of microphones as in the past, the evaluation sound output from one speaker (sound source) The time to record and measure with a microphone is approximately α × L
The degree can be estimated.
ここでαは測定条件(積分処理に要する時間や応答時間等)を考慮した定数であるが、測定データの積分時間や平均回数を大きくしてノイズ環境下での測定S/Nを確保するために、αは1より大きな値とする。全てのスピーカから出力される評価音を測定するためには、
α×L×S
程度の時間が必要とされる。
Here, α is a constant that takes into account measurement conditions (time required for integration processing, response time, etc.), but in order to ensure measurement S / N in a noise environment by increasing the integration time and average number of measurement data. In addition, α is a value larger than 1. To measure the evaluation sound output from all speakers,
α × L × S
About time is required.
一方で、本実施の形態1〜実施の形態3に示した音場測定システムを用いる場合における測定時間を求めるためには、
1)ノイズ環境下において測定S/Nを確保するために必要な測定時間
2)方程式の係数行列が逆行列を持つために原理上必要となる入力信号の長さ
を考慮する必要がある。
On the other hand, in order to obtain the measurement time in the case of using the sound field measurement system shown in the first to third embodiments,
1) Measurement time required to ensure measurement S / N in a noisy environment 2) Since the coefficient matrix of the equation has an inverse matrix, it is necessary to consider the length of the input signal that is necessary in principle.
まず1つ目の「ノイズ環境下での測定S/Nを確保するために必要な測定時間」について考える。本実施の形態に示した測定方式は、他のスピーカ(音源)から再生される音が測定ノイズとして加算されることがない方式である。このため、測定S/Nを確保するために必要な測定時間は、本実施の形態に示したように複数のスピーカから同時に評価音を出力する方式であっても、従来の方式により1つのスピーカで出力される音を測定する方式であっても、一度の測定で必要とされる測定時間は同等である。従って、測定時間は音源の数に依存せず、
α×L
と見積もることができる。
First, consider the first “measurement time required to ensure measurement S / N in a noisy environment”. The measurement method shown in this embodiment is a method in which sounds reproduced from other speakers (sound sources) are not added as measurement noise. For this reason, the measurement time required to secure the measurement S / N is one speaker according to the conventional method even if the evaluation sound is output simultaneously from a plurality of speakers as shown in the present embodiment. The measurement time required for a single measurement is the same even if the sound output by the method is measured. Therefore, the measurement time does not depend on the number of sound sources,
α × L
Can be estimated.
次に2つ目の「方程式の係数行列が逆行列を持つために原理上必要となる入力信号の長さ」は
β×L×S
と見積もることができる。
なお、βの値は、入力信号(評価音)の長さとして演算処理上必要な値を示しており、応答時間Lに比例した値となる。このβの値は、測定環境のノイズの大きさには無関係であり、β=2程度で十分であると判断できる。
Next, the second “length of the input signal required in principle because the coefficient matrix of the equation has an inverse matrix” is β × L × S
Can be estimated.
Note that the value of β indicates a value necessary for calculation processing as the length of the input signal (evaluation sound), and is a value proportional to the response time L. The value of β is irrelevant to the magnitude of noise in the measurement environment, and it can be determined that approximately β = 2 is sufficient.
以上2つの条件を考慮すると、本実施の形態に示した測定方法を用いる際に必要とされる測定時間は、
MAX(α×L,β×L×S)
(:(α×L)と(β×L×S)とを比較して大きな値を有効とする)
と見積もることができる。
Considering the above two conditions, the measurement time required when using the measurement method shown in this embodiment is:
MAX (α × L, β × L × S)
(: (Α × L) and (β × L × S) are compared and a larger value is valid)
Can be estimated.
従って、従来方式を用いて測定する時間と本提案方式を用いて測定する時間との比は以下のようになる。
1 …音場測定システム
2 …音響空間
3 …スピーカ部
4 …音声出力装置
5 …コンピュータ
7 …マイクアレイ
8 …音声入力装置
9 …スピーカ
10 …マイク
12 …評価音生成部(評価音生成手段)
13 …収録音記録部
14 …伝達関数算出部(伝達関数算出手段)
15 …伝達関数記録部
16 …制御部
DESCRIPTION OF
13 ... Recorded
15 ... transfer
Claims (5)
前記スピーカより出力するための評価音の生成を行う評価音生成手段と、
前記マイクを用いて収録された収録音と前記評価音生成手段により生成された評価音とに基づいて、各スピーカから各マイクに対応する伝達関数を算出する伝達関数算出手段と
を備え、
前記伝達関数算出手段は、前記収録音と前記評価音とに基づいて、時間領域における反復解法を適用することにより、前記伝達関数の算出処理を実行し、前記反復解法の適用において用いる相関演算処理および畳み込み演算処理のみ、周波数領域への変換処理によって演算処理を行った後に前記時間領域に演算結果を逆変換させる演算方法を用い、
前記評価音生成手段は、前記評価音におけるスピーカ毎の畳み込み行列を行列要素とし、各スピーカに対応する前記行列要素を備えた第1行列と、該第1行列の転置行列を成す第2行列との積の逆行列が存在し得るように評価音を生成する
ことを特徴とする音場測定システム。 A transfer function corresponding to each microphone from each speaker by simultaneously outputting evaluation sounds from a plurality of speakers set in a specific acoustic space and simultaneously recording the evaluation sounds from a plurality of microphones installed in the acoustic space. A sound field measuring system for calculating
Evaluation sound generating means for generating evaluation sound for output from the speaker;
Transfer function calculating means for calculating a transfer function corresponding to each microphone from each speaker based on the recorded sound recorded using the microphone and the evaluation sound generated by the evaluation sound generating means, and
The transfer function calculation means executes the transfer function calculation process by applying an iterative solution in the time domain based on the recorded sound and the evaluation sound, and a correlation calculation process used in the application of the iterative solution And only the convolution calculation process, after performing the calculation process by the conversion process to the frequency domain, using the calculation method to inversely convert the calculation result in the time domain,
The evaluation sound generation means uses a convolution matrix for each speaker in the evaluation sound as a matrix element, a first matrix having the matrix element corresponding to each speaker, and a second matrix that forms a transposed matrix of the first matrix; A sound field measurement system, characterized in that an evaluation sound is generated so that an inverse matrix of the product of can exist.
を特徴とする請求項1に記載の音場測定システム。 In the transfer function calculation process, the transfer function calculation means multiplies a calculation formula for obtaining the transfer function by a diagonal matrix that constitutes a time window, and changes elements of the diagonal matrix. An evaluation sound is continuously output from the plurality of speakers, and a transfer function in an arbitrary time range is calculated when recording sound is continuously recorded from the plurality of microphones. The sound field measurement system according to 1.
を特徴とする請求項1又は請求項2に記載の音場測定システム。 The evaluation sound generating means is recorded by the microphone by changing the frequency characteristic of the evaluation sound so that the output value of the low frequency band in the evaluation sound is higher than the output value of the high frequency band. The sound field measurement system according to claim 1 or 2, wherein the S / N is improved in the low frequency band of the recorded sound.
ことを特徴とする請求項1乃至請求項3のいずれか1項に記載の音場測定システム。 The transfer function calculating means, among the calculation results of the correlation calculation process or the convolution calculation process obtained by performing inverse conversion of the calculation result in the time domain after performing the calculation process by the conversion process to the frequency domain, The transfer function calculation process in the time domain is performed using only data in a range in which an error from the calculation result obtained when the correlation calculation process or the convolution calculation process is performed only in the time domain. The sound field measurement system according to claim 1, wherein the sound field measurement system is a sound field measurement system.
を特徴とする請求項4に記載の音場測定システム。 The sound field measurement system according to claim 4, wherein the range in which the error does not occur is determined based on a finite response time length set in advance according to the reverberation time of the evaluation sound.
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
JP2009008340A JP2010164863A (en) | 2009-01-19 | 2009-01-19 | Sound field measurement system |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
JP2009008340A JP2010164863A (en) | 2009-01-19 | 2009-01-19 | Sound field measurement system |
Publications (1)
Publication Number | Publication Date |
---|---|
JP2010164863A true JP2010164863A (en) | 2010-07-29 |
Family
ID=42581056
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
JP2009008340A Pending JP2010164863A (en) | 2009-01-19 | 2009-01-19 | Sound field measurement system |
Country Status (1)
Country | Link |
---|---|
JP (1) | JP2010164863A (en) |
Cited By (1)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
WO2020022347A1 (en) * | 2018-07-24 | 2020-01-30 | ソニー株式会社 | Measurement device and measurement system |
-
2009
- 2009-01-19 JP JP2009008340A patent/JP2010164863A/en active Pending
Cited By (4)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
WO2020022347A1 (en) * | 2018-07-24 | 2020-01-30 | ソニー株式会社 | Measurement device and measurement system |
JPWO2020022347A1 (en) * | 2018-07-24 | 2021-08-02 | ソニーグループ株式会社 | Measuring device and measuring system |
JP7347424B2 (en) | 2018-07-24 | 2023-09-20 | ソニーグループ株式会社 | Measuring device and measuring system |
US11805377B2 (en) | 2018-07-24 | 2023-10-31 | Sony Corporation | Measurement device and measurement system |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
EP3232686B1 (en) | Neural network-based loudspeaker modeling with a deconvolution filter | |
CN104041074B (en) | Method and apparatus for processing signals of a spherical microphone array on a rigid sphere used for generating an ambisonics representation of the sound field | |
JP4499920B2 (en) | Generation of nonlinear model and generation of driving signal for simulation test using it | |
EP1772713A1 (en) | Impulse-responsive measurement method and device | |
CN103634726A (en) | Automatic loudspeaker equalization method | |
CN111031463B (en) | Microphone array performance evaluation method, device, equipment and medium | |
CN106255027B (en) | A kind of the sound quality Small Enclosure appraisal procedure and system of non-linear audio system | |
CN106168942B (en) | A kind of fluctuation types dynamic data reconstructing method based on singular boundary method | |
Takeuchi et al. | Source directivity approximation for finite-difference time-domain simulation by estimating initial value | |
CN103916810B (en) | A kind of time domain acoustic energy compared with control method and system | |
CN1659926B (en) | Method and system of representing a sound field | |
JPH10320008A (en) | Efficient composition for composite system to be driven | |
CN115438589A (en) | Fishing rod adjustability prediction model based on BP neural network and optimization method thereof | |
Antonello et al. | Source localization and signal reconstruction in a reverberant field using the FDTD method | |
JP5791081B2 (en) | Sound source separation localization apparatus, method, and program | |
US9668075B2 (en) | Estimating parameter values for a lumped parameter model of a loudspeaker | |
JP5986966B2 (en) | Sound field recording / reproducing apparatus, method, and program | |
CN110489800B (en) | Structural dynamic load sparse identification method based on matrix regularization | |
JP2010164863A (en) | Sound field measurement system | |
CN112016235A (en) | Impact load identification method and system of flexible antenna structure | |
JP5106936B2 (en) | Sound field reproduction filter calculation device and sound field reproduction system | |
JP2019054344A (en) | Filter coefficient calculation device, sound pickup device, method thereof, and program | |
JP4653674B2 (en) | Signal separation device, signal separation method, program thereof, and recording medium | |
JP2022500710A (en) | Combined sound source localization and separation method for acoustic sources | |
JP2018077139A (en) | Sound field estimation device, sound field estimation method and program |