JP6288561B2 - Blind signal separation method and apparatus - Google Patents

Description

本発明は、未知の畳み込み混合系により混在した互いに統計的に独立な未知信号源信号を、観測信号のみから推定するブラインド信号分離方法に係わり、特に、最小２乗型同時対角化問題の解法を用いて高い精度で信号を分離することができるブラインド信号分離装置に関する。  The present invention relates to a blind signal separation method for estimating unknown signal source signals that are statistically independent from each other mixed by an unknown convolutional mixing system from observation signals only, and in particular, a method for solving a least-squares simultaneous diagonalization problem. The present invention relates to a blind signal separation device that can separate a signal with high accuracy by using the.

複数の未知信号源信号が未知の畳み込み混合系により混在されて観測されるとき、観測信号を分離して混在前の未知信号源信号を推定する処理をブラインド信号分離という。ブラインド信号分離では、未知信号源信号間の統計的独立性のみを条件として、観測信号から未知信号源信号を推定する方法であり、信号源の位置或いは観測信号の到来方向の推定を必ずしも必要としない方法である。  When a plurality of unknown signal source signals are mixed and observed by an unknown convolutional mixing system, the process of separating the observation signals and estimating the unknown signal source signals before mixing is called blind signal separation. Blind signal separation is a method of estimating an unknown signal source signal from an observation signal only on the condition of statistical independence between unknown signal source signals, and it is not always necessary to estimate the position of the signal source or the arrival direction of the observation signal. It is a way not to.

勾配法に基づく周波数領域の信号分離方法が非特許文献１で提案されている。勾配法に基づき周波数ビン毎に分離行列を更新することによって分離信号を推定する。  Non-Patent Document 1 proposes a frequency domain signal separation method based on the gradient method. Based on the gradient method, the separation signal is estimated by updating the separation matrix for each frequency bin.

最小２乗型同時対角化問題の解法を用いたブラインド信号分離方法が非特許文献２と非特許文献３で提案されている。これらの方法は、最小２乗型同時対角化問題の解法による混合行列の推定、最小２乗型一般化逆行列を用いた混合行列からの分離行列の推定、スケーリング問題の解法、パーミュテーション問題の解法の４つの手順から成る。  Non-Patent Document 2 and Non-Patent Document 3 propose a blind signal separation method using a method of solving the least square simultaneous diagonalization problem. These methods include estimation of the mixing matrix by solving the least-squares simultaneous diagonalization problem, estimation of the separation matrix from the mixing matrix using the least-squares generalized inverse matrix, solving the scaling problem, permutation It consists of four procedures for solving the problem.

周波数ビン間の電力比の相関に基づきパーミュテーション行列を推定する方法が非特許文献４で提案されている。  A method for estimating a permutation matrix based on the correlation of power ratios between frequency bins is proposed in Non-Patent Document 4.

Ｌ．Ｐａｒｒａ ａｎｄ Ｃ．Ｓｐｅｎｃｅ，”Ｃｏｎｖｏｌｕｔｉｖｅ ｂｌｉｎｄ ｓｅｐａｒａｔｉｏｎ ｏｆ ｎｏｎ−ｓｔａｔｉｏｎａｒｙ Ｓｏｕｒｃｅｓ”，ＩＥＥＥ Ｔｒａｎｓａｃｔｉｏｎｓ ｏｎ Ｓｐｅｅｃｈ ａｎｄ Ａｕｄｉｏ Ｐｒｏｃｅｓｓｉｎｇ，ｖｏｌ．８，ｎｏ．３，ｐｐ．３２０−３２７，Ｍａｙ，２０００  L. Parara and C.M. Spence, “Contributive blind separation of non-stationary Sources”, IEEE Transactions on Speech and Audio Processing, vol. 8, no. 3, pp. 320-327, May, 2000 Ａ．Ｙｅｒｅｄｏｒ，”Ｎｏｎ−ｏｒｔｈｏｇｏｎａｌ ｊｏｉｎｔ ｄｉａｇｏｎａｌｉｚａｔｉｏｎ ｉｎ ｔｈｅ ｌｅａｓｔ−ｓｑｕａｒｅｓ ｓｅｎｓｅ ｗｉｔｈ ａｐｐｌｉｃａｔｉｏｎ ｉｎ ｂｌｉｎｄ ｓｏｕｒｃｅ ｓｅｐａｒａｔｉｏｎ”，ＩＥＥＥ Ｔｒａｎｓａｃｔｉｏｎｓ ｏｎ Ｓｉｇｎａｌ Ｐｒｏｃｅｓｓｉｎｇ，ｖｏｌ．５０，ｎｏ．７，ｐｐ．１５４５−１５５３，Ｊｕｌｙ ２００２  A. Yedor, “Non-orthogonal joint diagonalization in the least-squares sense with application in blind source separation, IEEE Transactions on Signal. 50, no. 7, pp. 1545-1553, July 2002 Ｋ．Ｒａｈｂａｒ ａｎｄ Ｊ．Ｐ．Ｒｅｉｌｌｙ，”Ａ ｆｒｅｑｕｅｎｃｙ ｄｏｍａｉｎ ｍｅｔｈｏｄ ｆｏｒ ｂｌｉｎｄ ｓｏｕｒｃｅ ｓｅｐａｒａｔｉｏｎ ｏｆ ｃｏｎｖｏｌｕｔｉｖｅ ａｕｄｉｏ ｍｉｘｔｕｒｅｓ”，ＩＥＥＥ Ｔｒａｎｓａｃｔｉｏｎｓ ｏｎ Ｓｐｅｅｃｈ ａｎｄ Ａｕｄｉｏ Ｐｒｏｃｅｓｓｉｎｇ，ｖｏｌ．１３，ｎｏ．５，ｐｐ．８３２−８４４，Ｓｅｐｔ．２００５  K. Rahbar and J.M. P. Reilly, "A frequency domain method for blind source separation of convenient audio mixtures", IEEE Transactions on Audio and Proceeding and Audio. 13, no. 5, pp. 832-844, Sept. 2005 Ｈ．Ｓａｗａｄａ，Ｓ．Ａｒａｋｉ，ａｎｄ Ｓ．Ｍａｋｉｎｏ，”Ｍｅａｓｕｒｉｎｇ ｄｅｐｅｎｄｅｎｃｅ ｏｆ ｂｉｎ−ｗｉｓｅ ｓｅｐａｒａｔｅｄ ｓｉｇｎａｌｓ ｆｏｒ ｐｅｒｍｕｔａｔｉｏｎ ａｌｉｇｎｍｅｎｔ ｉｎ ｆｒｅｑｕｅｎｃｙ−ｄｏｍａｉｎ ＢＳＳ”，Ｐｒｏｃ．ＩＥＥＥ Ｉｎｔ．Ｓｙｍｐ．Ｃｉｒｃｕｉｔｓ Ｓｙｓｔ．，ｐｐ．３２４７−３２５０，Ｍａｙ ２００７  H. Sawada, S .; Araki, and S.R. Makino, “Measuring dependency of bin-wise separated signals for permutation alignment in frequency-domain BSS”, Proc. IEEE Int. Symp. Circuits Syst. , Pp. 3247-3250, May 2007

従来のブラインド信号分離方法で使用されている勾配法は演算量が少ないが、収束が遅く、実環境下で十分な信号分離精度を得ることができない。  Although the gradient method used in the conventional blind signal separation method has a small amount of calculation, convergence is slow and sufficient signal separation accuracy cannot be obtained in an actual environment.

従来のブラインド信号分離方法では、時間周波数領域で周波数ビンとエポック毎、即ち、ポイント毎に観測信号の空間相関行列をそのフロベニウスノルムで正規化しているため、未知信号源信号が音声の場合、音声の特徴を表す周波数特性が失われ、周波数ビンに対して振幅特性が一定になり、信号分離精度が劣化する。  In the conventional blind signal separation method, the spatial correlation matrix of the observed signal is normalized by the Frobenius norm for each frequency bin and epoch, that is, for each point in the time-frequency domain. The frequency characteristic representing the characteristic is lost, the amplitude characteristic becomes constant with respect to the frequency bin, and the signal separation accuracy deteriorates.

従来のブラインド信号分離方法では、時間周波数領域で観測信号の空間相関行列から最小２乗型同時対角化問題の解法を用いて混合行列を推定し、次いで最小２乗型一般化逆行列を用いて混合行列から分離行列を推定しているため、最小２乗型同時対角化問題に分離行列の推定が含まれず、仮に混合行列を精度良く推定できても分離行列が精度良く推定できるとは限らず、信号分離精度が劣化する場合がある。  In the conventional blind signal separation method, the mixing matrix is estimated from the spatial correlation matrix of the observed signal in the time-frequency domain using the method of solving the least-squares simultaneous diagonalization problem, and then the least-squares generalized inverse matrix is used. Because the separation matrix is estimated from the mixing matrix, the estimation of the separation matrix is not included in the least square simultaneous diagonalization problem, and even if the mixing matrix can be estimated with high accuracy, the separation matrix can be estimated with high accuracy. Not limited to this, signal separation accuracy may deteriorate.

従来のブラインド信号分離方法では、時間周波数領域で観測信号の空間相関行列から制約条件を課すことなく最小２乗型同時対角化問題の解法を用いて混合行列を推定し、最小２乗型一般化逆行列を用いて混合行列から分離行列を推定する前に、塁乗法により制約条件を課しているので、手順が収束に要する反復回数は大幅に増加する。  In the conventional blind signal separation method, the mixture matrix is estimated using the solution of the least-squares simultaneous diagonalization problem without imposing constraints from the spatial correlation matrix of the observed signal in the time-frequency domain. Since the constraint is imposed by the power method before estimating the separation matrix from the mixing matrix using the generalized inverse matrix, the number of iterations required for the procedure to converge is greatly increased.

更に、従来のブラインド信号分離方法において混合行列の最小２乗型一般化逆行列を単に分離行列として取り扱うと、加法性雑音が存在する場合、雑音の影響が大きくなって信号分離性能の劣化につながる。  Furthermore, when the least square generalized inverse matrix of the mixing matrix is simply treated as a separation matrix in the conventional blind signal separation method, if additive noise is present, the influence of noise increases and leads to degradation of signal separation performance. .

同一信号源から発生した信号の隣接周波数ビンに相関があることを利用した従来のパーミュテーション問題の解法では、１度誤りが生じると、これ以降の解法が誤り続ける確率が非常に高くなる。  In the conventional solution of the permutation problem using the fact that there is a correlation between adjacent frequency bins of signals generated from the same signal source, once an error occurs, the probability that the subsequent solution will continue to be very high becomes very high.

本発明はこのような事情を鑑みてなされたものであり、信号源とブラインド信号分離装置出力間を表した混合行列と分離行列の縦続接続モデルを最小２乗型同時対角化問題に取り込み、更にモデルに因果性を与えることによって信号分離性能の向上と反復手順の収束高速化の両立、更に、同一信号源から発生した信号の周波数ビン間の電力比相関性と多数決を組み合わせによってパーミュテーション問題を正確に解法することを目的とする。  The present invention has been made in view of such circumstances, and incorporates a mixed matrix and a cascaded model of separation matrices representing between the signal source and the output of the blind signal separation device into the least square simultaneous diagonalization problem, Furthermore, by providing causality to the model, it is possible to improve signal separation performance and speed up the convergence of the iterative procedure. In addition, permutation is achieved by combining the power ratio correlation between the frequency bins of signals generated from the same signal source and majority vote. The goal is to solve the problem accurately.

このような目的に応えるために本発明（請求項１記載の発明）に係るブラインド信号分離方法は、未知の畳み込み混合系により混在した互いに統計的に独立な未知信号源信号を、観測信号のみからブラインドで推定する方法であって、時間周波数領域においてエポック時刻毎に観測信号ベクトルから空間相関行列を求め、各エポック時刻において周波数ビンの中からフロベニウスノルムが最大となる空間相関行列を求め、そのノルムで全ての同エポック時刻の空間相関行列を正規化した後、正規化された空間相関行列から混合行列を推定するために、制約条件付きフォワードモデル型最小２乗型同時対角化問題をラグランジュの未定乗数法と反復法の組み合わせにより解法する手順と、手順の収束後、累乗法を１回用いて最小２乗型同時対角化問題の解となる混合行列を求め、混合行列のランクが落ちていた場合は、混合行列がフルランクになるように特異値分解を用いて基底を補い、次いで、ランク落ちした混合行列の特異値とその総和が変わらないように特異値を補い、最小２乗型一般化逆行列を用いて混合行列から分離行列を推定し、バックワード型最小２乗型同時対角化問題を最小２乗法で解いて対角行列を求め、同一信号源から発生した信号の周波数ビンに相関があることを利用してパーミュテーション問題を、多数決を利用して解法することによって信号分離精度の高い分離行列を求めることを特徴とする。  In order to meet such an object, the blind signal separation method according to the present invention (the invention according to claim 1) is configured to extract unknown signal source signals that are statistically independent from each other mixed only by an unknown convolutional mixture system from only observation signals. A blind estimation method that obtains a spatial correlation matrix from the observed signal vector at each epoch time in the time-frequency domain, and obtains a spatial correlation matrix that maximizes the Frobenius norm from the frequency bins at each epoch time. In order to estimate the mixture matrix from the normalized spatial correlation matrix after normalizing the spatial correlation matrix for all the same epochs in, Lagrange's forward model least squares simultaneous diagonalization problem is A procedure to solve by a combination of the undetermined multiplier method and the iterative method, and after convergence of the procedure, a least square simultaneous diagonal using the power method once Find the mixing matrix that solves the problem, and if the rank of the mixing matrix has dropped, use the singular value decomposition to make the mixing matrix full rank, and then use the singular value of the mixing matrix that has fallen in rank Singular values are supplemented so that the sum does not change, the separation matrix is estimated from the mixing matrix using the least square generalized inverse matrix, and the backward least square simultaneous diagonalization problem is solved by the least square method. Solving the diagonal matrix and solving the permutation problem using the fact that there is a correlation in the frequency bins of signals generated from the same signal source, and obtaining a separation matrix with high signal separation accuracy by using majority voting It is characterized by seeking.

本発明（請求項２記載の発明）に係るブラインド信号分離方法は、信号源とブラインド信号分離装置出力間の伝達関数を表した混合行列と分離行列の縦続接続モデルにラグランジュの未定乗数を導入した制約条件付きフォワードモデル型最小２乗型同時対角化問題とバックワード型最小２乗型同時対角化問題を導入し、更に分離行列に遅延を与えることによって信号源とブラインド信号分離装置出力間の因果的なモデルを推定することを特徴とする。  In the blind signal separation method according to the present invention (the invention described in claim 2), a Lagrange's undetermined multiplier is introduced into a cascade connection model of a mixing matrix and a separation matrix representing a transfer function between the signal source and the output of the blind signal separation device. Introducing a constrained forward model least-squares simultaneous diagonalization problem and backward-type least squares simultaneous diagonalization problem, and further delaying the separation matrix between the signal source and the blind signal separator output It is characterized by estimating a causal model.

本発明（請求項３記載の発明）に係るブラインド信号分離方法は、基準周波数ビンを複数選択し、基準周波数ビン間において最も電力比の相関が大きいパーミュテーション行列を推定し、複数の基準周波数ビンから基準周波数ビンを１つ選択し、全ての周波数ビン間で電力比の相関が最も大きいパーミュテーション行列を推定する手順を、全ての基準周波数ビンが１度選択されるまで繰り返した後、各周波数ビンに複数割り当てられたパーミュテーション行列から多数決によって周波数ビンに割り当てられるパーミュテーション行列を決定し、多数決によってパーミュテーション行列が決定できない場合には、最も相関値が大きいパーミュテーション行列を採用することを特徴とする。  A blind signal separation method according to the present invention (invention of claim 3) selects a plurality of reference frequency bins, estimates a permutation matrix having the largest correlation of power ratios among the reference frequency bins, and generates a plurality of reference frequencies. After selecting one reference frequency bin from the bin and estimating the permutation matrix having the largest correlation of the power ratio among all frequency bins, repeating until all the reference frequency bins are selected once, If the permutation matrix assigned to the frequency bin is determined by majority from the permutation matrix assigned to each frequency bin, and the permutation matrix cannot be determined by majority, the permutation matrix with the largest correlation value It is characterized by adopting.

すなわち、本発明によれば、観測信号の空間相関行列を制約条件付き最小２乗型同時対角化問題の対象に、ラグランジュの未定乗数法と反復法を組み合わせることによって最小２乗型同時対角化問題の近似解、即ち、混合行列を求めることができ、推定性能の向上と反復回数の低減を両立させる。  That is, according to the present invention, the least square simultaneous diagonalization is performed by combining the Lagrange's undetermined multiplier method and the iterative method with the spatial correlation matrix of the observed signal as a target of the constrained least square simultaneous diagonalization problem. An approximate solution of the optimization problem, that is, a mixing matrix can be obtained, and both improvement in estimation performance and reduction in the number of iterations are achieved.

また、本発明によれば、各周波数ビンに複数割り当てられたパーミュテーション行列から多数決によって当該周波数ビンに割り当てられるパーミュテーション行列を決定し、全体の信号分離性能を向上させる。  Further, according to the present invention, the permutation matrix assigned to the frequency bin is determined by majority from the permutation matrix assigned to each frequency bin, thereby improving the overall signal separation performance.

本発明によれば、時間周波数領域においてエポック時刻毎に観測信号ベクトルから空間相関行列を求め、各エポック時刻において周波数ビンの中からフロベニウスノルムが最大となる空間相関行列を求め、そのノルムで全ての同エポック時刻の空間相関行列を正規化した後、これに制約条件付き最小２乗型同時対角化問題を適用する。制約条件付き最小２乗型同時対角化問題では、ラグランジュの未定乗数法と反復法を組み合わせることによって近似解、即ち、混合行列を求める。本発明に係るブラインド信号分離方法及びブラインド信号分離装置では、混合行列のランクが落ちていた場合は、混合行列がフルランクになるように特異値分解を用いて基底と特異値を補い、最小２乗型一般化逆行列を用いて混合行列から分離行列を推定し、バックワード型最小２乗型同時対角化問題を最小２乗法で解いて対角行列を求めるため、少ない反復回数で近似解に収束するという効果がある。  According to the present invention, the spatial correlation matrix is obtained from the observation signal vector at each epoch time in the time-frequency domain, and the spatial correlation matrix that maximizes the Frobenius norm is obtained from the frequency bins at each epoch time. After normalizing the spatial correlation matrix of the same epoch time, the least square simultaneous diagonalization problem with constraints is applied thereto. In the constrained least-squares simultaneous diagonalization problem, an approximate solution, that is, a mixing matrix is obtained by combining Lagrange's undetermined multiplier method and iterative method. In the blind signal separation method and the blind signal separation apparatus according to the present invention, when the rank of the mixing matrix is lowered, the basis and the singular value are supplemented by using singular value decomposition so that the mixing matrix becomes the full rank, and the minimum 2 Estimating the separation matrix from the mixed matrix using the squared generalized inverse matrix and solving the backward least square simultaneous diagonalization problem by the least square method to obtain the diagonal matrix, the approximate solution with a small number of iterations Has the effect of converging.

また、本発明に係るブラインド信号分離方法及びブラインド信号分離装置では、信号源とブラインド信号分離装置出力間の伝達関数を表した混合行列と分離行列の縦続接続モデルにラグランジュの未定乗数を導入した制約条件付きフォワードモデル型最小２乗型同時対角化問題とバックワード型最小２乗型同時対角化問題を導入し、更に分離行列に遅延を与えることによって信号源とブラインド信号分離装置出力間の因果的なモデルを推定することによって雑音の影響を最小にする分離行列を推定できるという効果がある。  Further, in the blind signal separation method and the blind signal separation device according to the present invention, a constraint in which a Lagrange's undetermined multiplier is introduced to a cascade connection model of a mixing matrix and a separation matrix representing a transfer function between a signal source and a blind signal separation device output. Introducing conditional forward model least-squares simultaneous diagonalization problem and backward least-squares simultaneous diagonalization problem, and further delaying the separation matrix between signal source and blind signal separator output By estimating the causal model, the separation matrix that minimizes the influence of noise can be estimated.

更に、パーミュテーション問題の解法では、基準周波数ビンを複数選択し、基準周波数ビン間において最も電力比の相関が大きいパーミュテーション行列を推定し、複数の基準周波数ビンから基準周波数ビンを１つ選択し、全ての周波数ビン間で電力比の相関が最も大きいパーミュテーション行列を推定する手順を全ての基準周波数ビンが１度選択されるまで繰り返した後、各周波数ビンに複数割り当てられたパーミュテーション行列から多数決によって周波数ビンに割り当てられるパーミュテーション行列を決定し、多数決によってパーミュテーション行列が決定できない場合には、最も相関値が大きいパーミュテーション行列を採用することを特徴とすることによって本発明に係るブラインド信号分離方法及びブラインド信号分離装置には、信号分離性能を高めることができるという効果がある。  Further, in solving the permutation problem, a plurality of reference frequency bins are selected, a permutation matrix having the largest correlation of the power ratio among the reference frequency bins is estimated, and one reference frequency bin is selected from the plurality of reference frequency bins. After selecting and repeating the procedure for estimating the permutation matrix having the highest power ratio correlation among all frequency bins until all the reference frequency bins are selected once, a plurality of permutations assigned to each frequency bin are assigned. The permutation matrix assigned to the frequency bin is determined from the muting matrix by majority voting, and if the permutation matrix cannot be determined by voting, the permutation matrix having the largest correlation value is adopted. By the blind signal separation method and blind signal separation device according to the present invention , There is an effect that it is possible to enhance the signal separation performance.

本発明に係るブラインド信号分離方法の実施の形態を示す図である。  It is a figure which shows embodiment of the blind signal separation method which concerns on this invention. 本発明に係るブラインド信号分離方法において最小２乗型同時対角化問題の対角化行列のラグランジェ未定乗数法の適用について説明するためのフローチャートである。  6 is a flowchart for explaining application of a Lagrangian undetermined multiplier method of a diagonalization matrix of a least square simultaneous diagonalization problem in the blind signal separation method according to the present invention. 本発明に係るブラインド信号分離方法において最小２乗型同時対角化問題の対角行列の解法について説明するためのフローチャートである。  6 is a flowchart for explaining a method of solving a diagonal matrix of the least square simultaneous diagonalization problem in the blind signal separation method according to the present invention. 本発明に係るブラインド信号分離方法のパーミュテーション問題の解法について説明するためのフローチャートである。  5 is a flowchart for explaining a solution of a permutation problem of the blind signal separation method according to the present invention. 本発明に係る信号分離方法の実施例１における信号源（スピーカ）とマイクロホンの位置関係を表す平面図である。  It is a top view showing the positional relationship of the signal source (speaker) and microphone in Example 1 of the signal separation method which concerns on this invention. 本発明に係る信号分離方法の実施例１における部屋の残響時間と信号分離性能の関係を示す図である。  It is a figure which shows the relationship between the reverberation time of a room and signal separation performance in Example 1 of the signal separation method which concerns on this invention. 本発明に係る信号分離方法の実施例１におけるＳＮＲと信号分離性能の関係を示す図である。  It is a figure which shows the relationship between SNR and signal separation performance in Example 1 of the signal separation method based on this invention. 実施例２における従来のパーミュテーション問題の解法によるパーミュテーション行列の割り当て結果を示す図である。  It is a figure which shows the allocation result of the permutation matrix by the solution of the conventional permutation problem in Example 2. FIG. 実施例２における本発明に係るブラインド信号分離方法によるパーミュテーション行列の割り当て結果を示す図である。  It is a figure which shows the allocation result of the permutation matrix by the blind signal separation method which concerns on this invention in Example 2. FIG. 分離行列の推定に本発明に係るブラインド信号分離方法を使用して、パーミュテーション問題の解法には非特許文献４を使用した場合と、本発明に係るブラインド信号分離方法の実施例２におけるパーミュテーション行列の割り当て正答率と信号分離性能を比較する図である。  When the blind signal separation method according to the present invention is used for estimation of the separation matrix and Non-Patent Document 4 is used for the solution of the permutation problem, and in the second embodiment of the blind signal separation method according to the present invention, It is a figure which compares the allocation correct answer rate of a mutation matrix, and signal separation performance.

本発明に係るブラインド信号分離方法の実施の形態について図面を参照して説明する。  An embodiment of a blind signal separation method according to the present invention will be described with reference to the drawings.

１．畳み込み混合モデル
図１に示すように、時刻ｔにおいてＮ個の信号源１１、１２、…、１Ｎから発せられた信号源信号ｓ_ｊ（ｔ）が畳み込み混合されてｘ_ｉ（ｔ）として観測される。ｓ_ｊ（ｔ）は平均０で互いに統計的独立な非定常信号である。ｈ_ｉｊ（ｔ）は信号源１ｊからマイクロホン２ｉまでの経路の時不変なインパルス応答で、因巣的で非最小位相系である。また、ｎ_ｉ（ｔ）はマイクロホン２ｉに加わる平均０、分散σ^２のガウス性白色雑音で、ｓ_ｊ（ｔ）と統計的独立である。時刻ｔにおいてＪ個のマイクロホン２１、２２、…、２Ｊで観測される観測信号ｘ_ｉ（ｔ）は式（１）で表される。ここで、Ｊ≧Ｎ≧２とする。1. As shown in FIG. 1, the signal source signals s _j (t) emitted from N signal sources 11, 12,..., 1N at time t are convolved and observed as x _i (t). The s _j (t) are non-stationary signals that are 0 on average and statistically independent of each other. h _ij (t) is a time-invariant impulse response of a path from the signal source 1j to the microphone 2i, and is a focal non-minimum phase system. Further, n _i (t) is a Gaussian white noise having an average of 0 and a variance σ ² applied to the microphone 2i, and is statistically independent of s _j (t). An observation signal x _i (t) observed by J microphones 21, 22,..., 2J at time t is expressed by Expression (1). Here, J ≧ N ≧ 2.

ここで、＊は畳み込み演算を表す。

Here, * represents a convolution operation.

観測信号ｘ_ｉ（ｔ）を３１で短時間フーリエ変換すると、フレーム時刻ｍにおける観測信号は式（２）により表される。式（２）において、ｗｉｎ（ｔ）は窓関数、Ｋは短時間フーリエ変換の点数、Ｔ_ｂは２つの重複窓間のシフトサイズ、ω_ｋ＝２πｋ／Ｋ、ｋ＝０，１，…，Ｋ−１、をそれぞれ表す。離散フーリエ変換の点数Ｋがｈ_ｉｊ（ｔ）のインパルス応答長より十分に大きいとき観測信号は式（３）により近似される。ここで、ｈ_ｉｊ（ｔ）のＫ点離散フーリエ変換をｈ_ｉｊ（ω_ｋ）、ｓ_ｊ（ｔ）に窓関数を乗算した後、Ｋ点短時間フーリエ変換で時間周波数領域に変換したフレーム時刻ｍの信号源信号をｓ_ｊ（ω_ｋ，ｍ）、同様に、ｎ_ｉ（ｔ）に窓関数を乗算した後、Ｋ点短時間フーリエ変換で時間周波数領域に変換したフレーム時刻ｍの雑音をｎ_ｉ（ω_ｋ，ｍ）とそれぞれ表記している。また、式（３）において、ｘ（ω_ｋ，ｍ）はフレーム時刻ｍに各マイクロホンでの観測信号ベクトル、ｓ（ω_ｋ，ｍ）はフレーム時刻ｍに各信号源信号ベクトル、混合行列Ｈ（ω_ｋ）はＮ個の信号源からＪ個のマイクロホンまでの混合行列、ｎ（ω_ｋ，ｍ）は雑音ベクトルでそれぞれ式（４）、（７）、（５）、（８）により定義される。入手可能なエポック時刻の総数をＭとすると、１≦ｍ≦Ｍとなる。信号源信号の共分散行列はＰ_ｓ（ω_ｋ，ｍ）＝Ｅ［ｓ（ω_ｋ，ｍ）ｓ（ω_ｋ，ｍ）^Ｈ］∈Ｒ^Ｎ×Ｎで、対角行列となる。Ｅ［・］と上付き添字^Ｈは期待値と複素共役転置をそれぞれ表す。また、上付き添字^ＴとＲ^Ｎ×Ｎは転置とＮ×Ｎの実数空間を表す。When the observation signal x _i (t) is Fourier-transformed for a short time by 31, the observation signal at the frame time m is expressed by Expression (2). In equation (2), win (t) is a window function, K is the number of short-time Fourier transform points, T _b is a shift size between two overlapping windows, ω _k = 2πk / K, k = 0, 1,. K-1, respectively. When the discrete Fourier transform score K is sufficiently larger than the impulse response length of h _ij (t), the observed signal is approximated by equation (3). Here, a frame time obtained by multiplying _Kij discrete Fourier transform of _hij (t) by _hij ([omega] _k ), _sj (t) by a window function and then transforming it into the time-frequency domain by K-point short-time Fourier transform. The noise of the frame time m obtained by multiplying the signal source signal of m by s _j (ω _k , m) and similarly by multiplying n _i (t) by the window function and then converting it to the time-frequency domain by K-point short-time Fourier transform. It is written as n _i (ω _k , m). In Equation (3), x (ω _k , m) is an observation signal vector at each microphone at frame time m, and s (ω _k , m) is each signal source signal vector and mixing matrix H (at frame time m. ω _k ) is a mixing matrix from N signal sources to J microphones, and n (ω _k , m) is a noise vector defined by equations (4), (7), (5), and (8), respectively. The If the total number of available epoch times is M, 1 ≦ m ≦ M. The covariance matrix of the signal source signal is P _s (ω _k , m) = E [s (ω _k , m) s (ω _k , m) ^H ] ∈R ^{N × N,} which is a diagonal matrix. E [•] and superscript ^H represent the expected value and the complex conjugate transpose, respectively. Superscript ^T and ^{RN × N} represent transpose and N × N real space.

信号を分離するには５１、５２、…、５Ｋで周波数ビン毎に式（９）を満足する分離行列Ｗ（ω_ｋ）を推定し、６０で信号源の割り当てを定めるパーミュテーション行列Π（ω_ｋ）∈Ｒ^Ｎ×Ｎを決定する。周波数ビン毎に独立にΠ（ω_ｋ）を決定しても信号が完全に分離する保証はなく、同一信号源から発生した信号の隣接または近接周波数ビンに相関があることを利用してパーミュテーション行列Π（ω_ｋ）を決定する。In order to separate the signals, a separation matrix W (ω _k ) satisfying Equation (9) is estimated for each frequency bin at 51, 52,..., 5K, and a permutation matrix 定 め る (60) that determines the signal source assignment at 60 Determine ω _k ) εR ^{N × N.} Even if Π (ω _k ) is determined independently for each frequency bin, there is no guarantee that the signals will be completely separated, and the permutation will be made by utilizing the correlation between adjacent or adjacent frequency bins of signals generated from the same signal source. Determine the rotation matrix Π (ω _k ).

ここで、Ｄ（ω_ｋ）∈Ｃ^Ｎ×Ｎは周波数ビン毎に異なる任意の対角行列である。

Here, D (ω _k ) εC ^{N × N} is an arbitrary diagonal matrix that differs for each frequency bin.

スケーリング問題とパーミュテーション問題を順に解法した後、７１、７２、…、７Ｋでｘ（ω_ｋ，ｍ）に左から分離行列Ｗ（ω_ｋ）を乗算すると、周波数ビンω_ｋにおける分離信号ｙ（ω_ｋ，ｍ）は式（１０）で表される。尚、スケーリング問題の解法については後述する。式（１０）を８０で短時間逆フーリエ変換と重複加算によって時間領域に変換すると分離信号ｙ_ｉ（ｔ）が求められる。雑音の分散σ^２が十分に小さいとき、ｙ_ｉ（ｔ）≒ｓ_ｉ（ｔ）になる。尚、分離信号ベクトルｙ（ω_ｋ，ｍ）は式（１１）により表される。After solving the scaling problem and the permutation problem in order, multiplying x (ω _k , m) by 71, 72,..., 7K from the left with the separation matrix W (ω _k ) _yields the separated signal y in the frequency bin ω _k . (Ω _k , m) is expressed by Equation (10). A method for solving the scaling problem will be described later. When Expression (10) is converted into the time domain by short-time inverse Fourier transform and overlap addition at 80, a separated signal y _i (t) is obtained. When the noise variance σ ² is sufficiently small, y _i (t) ≈s _i (t). Note that the separated signal vector y (ω _k , m) is expressed by Expression (11).

本発明に関するブラインド信号分離方法について、図１乃至図４を参照して詳細に説明する。図２乃至図４は、図１の３１における短時間フーリエ変換後、ブラインド信号分離システム４０において本発明により周波数ビン毎に推定される分離行列の算出手順を示したものである。図４は、分離行列の算出後、ブラインド信号分離システム４０において本発明によりパーミュテーション行列の算出手順を示したものである。  The blind signal separation method according to the present invention will be described in detail with reference to FIGS. 2 to 4 show the calculation procedure of the separation matrix estimated for each frequency bin in the blind signal separation system 40 after the short-time Fourier transform in 31 of FIG. 1 according to the present invention. FIG. 4 shows a procedure for calculating a permutation matrix according to the present invention in the blind signal separation system 40 after the separation matrix is calculated.

観測信号ｘ（ω_ｋ，ｍ）の共分散行列Ｐ_ｘ（ω_ｋ，ｍ）∈Ｃ^Ｊ×Ｊは式（１２）で与えられる。式（１３）の制約条件を課して式（１４）を満足する対角化行列Ｂ（ω_ｋ）と対角行列Λ（ω_ｋ，ｍ）を求めると、式（９）よりＢ（ω_ｋ）とＷ（ω_ｋ）の関係は式（１５）で与えられる。ただし、Ｃ^Ｊ×ＪはＪ×Ｊの複素空間を表す。The covariance matrix P _x (ω _k , m) εC ^{J × J} of the observation signal x (ω _k , m) is given by equation (12). When the diagonalization matrix B (ω _k ) and the diagonal matrix Λ (ω _k , m) satisfying the expression (14) with the constraint condition of the expression (13) are _obtained , B (ω _k ) and W (ω _k ) are given by equation (15). However, C ^{J × J} represents a J × J complex space.

ただし、Ｉは単位行列である。

Here, I is a unit matrix.

２．最小２乗型同時対角化問題とその解法
観測信号ｘ（ω_ｋ，ｍ）の共分散行列Ｐ_ｘ（ω_ｋ，ｍ）の推定値Ｐ_ｘ（ω_ｋ，ｍ）を正規化して、式（１７）を最小にする対角化行列Ｂ（ω_ｋ）と対角行列Λ（ω_ｋ，ｍ）を求める。式（１７）は最小２乗型同時対角化問題の解法として知られている。

2. Least-squares simultaneous diagonalization problem and its solution Normalize the estimated value P _x (ω _k , m) of the covariance matrix P _x (ω _k , m) of the observed signal x (ω _k , m) A diagonalization matrix B (ω _k ) and a diagonal matrix Λ (ω _k , m) that minimize (17) are _obtained . Expression (17) is known as a method for solving the least square simultaneous diagonalization problem.

本発明では、制約条件付き最小２乗型同時対角化問題を解法することによって混合行列Ｂ（ω_ｋ）を推定した後、分離行列Ｗ（ω_ｋ）を求める。次いで、対角行列Λ（ω_ｋ，ｍ）を推定するために、最小２乗法を用いて分離行列Ｗ（ω_ｋ）を用いた評価量を最小化する。本発明では、混合行列と対角行列の推定を交互に繰り返す。音声は低域周波数帯にフォルマントと呼ばれる振幅スペクトルのピークを有している。この音声波形の特徴を失うことなく、式（１９）によって観測信号ｘ（ω_ｋ，ｍ）の共分散行列Ｐ_ｘ（ω_ｋ，ｍ）を正規化することが、本発明の特徴の一つである。In the present invention, the separation matrix W (ω _k ) is _obtained after estimating the mixing matrix B (ω _k ) by solving the constrained least square simultaneous diagonalization problem. Next, in order to estimate the diagonal matrix Λ (ω _k , m), the evaluation quantity using the separation matrix W (ω _k ) is minimized using the least square method. In the present invention, estimation of the mixing matrix and the diagonal matrix is repeated alternately. Voice has a peak of an amplitude spectrum called formant in a low frequency band. One of the features of the present invention is to normalize the covariance matrix P _x (ω _k , m) of the observation signal x (ω _k , m) by the equation (19) without losing the characteristics of the speech waveform. It is.

２．１ 時間周波数領域における観測信号の共分散行列の正規化
時間領域観測信号は式（１８）の短時間フーリエ変換によって時間周波数領域に変換される。2.1 Normalization of observation signal covariance matrix in the time-frequency domain The time-domain observation signal is transformed into the time-frequency domain by the short-time Fourier transform of equation (18).

式（１７）において、ｗ（ｔ）は窓関数、Ｋは短時間フーリエ変換の点数、Ｔ_ｓは２つの重複窓間のシフトサイズ、Ｔ_ｓはエポックサイズ、Ｎ_ｓ＋１は各エポックにおける総重複フレーム数でＫ＋Ｎ_ｓＴ_ｓ≦Ｔ_ｂ、ｌ＝０，１，…，Ｎ_ｓをそれぞれ表す。

In equation (17), w (t) is a window function, K is the number of short-time Fourier transform points, T _s is the shift size between two overlapping windows, T _s is the epoch size, and N _s +1 is the total overlap in each epoch The number of frames represents K + N _s T _s ≦ T _b , l = 0, 1,..., N _s .

ステップＳ１０１においてエポック時刻ｍにおける時間周波数領域観測信号の共分散行列を式（１９）によって推定される。  In step S101, the covariance matrix of the time-frequency domain observation signal at the epoch time m is estimated by equation (19).

混合行列Ｂ（ω_ｋ）と分離行列Ｗ（ω_ｋ）を縦続に接続したとき、そのインパルス応答は式（１５）を最小にすることによって求められる。式（２１）の評価量ｅ（ω_ｋ）をＷ（ω_ｋ）によって微分すると、分離行列Ｗ（ω_ｋ）は式（２２）によって求められる。Ｂ（ω_ｋ）^ＨＢ（ω_ｋ）のランクがＮのときのみ、式（２１）の評価量ｅ（ω_ｋ）は零になる。一方、Ｂ（ω_ｋ）^ＨＢ（ω_ｋ）のランクがＮ未満のとき、ｅ（ω_ｋ）は零より大きくなる。そこで、制約条件‖ｂ_ｊ（ω_ｋ）‖_２＝１に制約条件ｒａｎｋ（Ｂ（ω_ｋ）^ＨＢ（ω_ｋ））＝Ｎを付け加え、混合行列Ｂ（ω_ｋ）の推定のための最小２乗型同時対角化問題を解法する。ここで、ｂ_ｊ（ω_ｋ）はＢ（ω_ｋ）のｊ番目の列ベクトル、‖・‖_２はユークリッドノルム、ｒａｎｋ（Ａ）は行列Ａのランクをそれぞれ表す。When the mixing matrix B (ω _k ) and the separation matrix W (ω _k ) are connected in cascade, the impulse response is obtained by minimizing the equation (15). When the evaluation quantity e (ω _k ) of Expression (21) is differentiated by W (ω _k ), the separation matrix W (ω _k ) is obtained by Expression (22). Only when the rank of B (ω _k ) ^H B (ω _k ) is N, the evaluation quantity e (ω _k ) of the equation (21) becomes zero. On the other hand, when the rank of B (ω _k ) ^H B (ω _k ) is less than N, e (ω _k ) is greater than zero. Therefore, the constraint condition rank (B (ω _k ) ^H B (ω _k )) = N is added to the constraint condition ‖ b _j (ω _k ) ‖ ₂ = 1, and the minimum for estimating the mixing matrix B (ω _k ) Solves the squared diagonalization problem. Here, b _j (ω _k ) represents the j-th column vector of B (ω _k ), ‖ · ‖ ₂ represents the Euclidean norm, and rank (A) represents the rank of the matrix A.

ここで、‖・‖_Ｆはフロベニウスノルムを表す。

Here, ‖ / ‖ _F represents the Frobenius norm.

ステップＳ１０２においてＰ_ｘ（ω_ｋ，ｍ）を式（２３）によって正規化する。In step S102, P _x (ω _k , m) is normalized by equation (23).

２．２ 対角化行列の解法
制約条件‖ｂ_ｊ（ω_ｋ）‖_２＝１を課したフォワードモデル型最小２乗型同時対角化問題を周波数ビンω_ｋ毎に解くことによって、対角化行列Ｂ（ω_ｋ）、即ち、混合行列を求める。
評価量を式（２４）に示す。式（２４）はフォワードモデル型最小２乗型同時対角化問題として知られている。2.2 Solving the diagonalization matrix By solving the forward model least-squares simultaneous diagonalization problem with constraints ‖b _j (ω _k ) ‖ ₂ = 1 for each frequency bin ω _k , the diagonal A quantization matrix B (ω _k ), that is, a mixing matrix is obtained.
The evaluation amount is shown in Formula (24). Equation (24) is known as a forward model type least squares type simultaneous diagonalization problem.

ここで、γ_ｉはラグランジェの未定乗数を表す。

Here, γ _i represents Lagrange's undetermined multiplier.

ベクトル表現を用いると、式（２４）の評価関数は式（２６）のように表現することができる。ここで、ｒ_ｘ（ω_ｋ，ｍ）、Ｇ（ω_ｋ）、ｄ（ω_ｋ，ｍ）、Ｇ（ω_ｋ）ｄ（ω_ｋ，ｍ）はそれぞれ式（２８）〜（３１）により表される。ただし、ｖｅｃ｛Ａ｝は行列Ａの列を積み重ね

れぞれ表す。λ_ｉは対角行列Ａのｉ番目の要素を表す。When the vector expression is used, the evaluation function of Expression (24) can be expressed as Expression (26). Here, r _x (ω _k , m), G (ω _k ), d (ω _k , m), and G (ω _k ) d (ω _k , m) are expressed by equations (28) to (31), respectively. Is done. Where vec {A} stacks the columns of matrix A

Represent each. λ _i represents the i-th element of the diagonal matrix A.

反復法によって混合行列Ｂ（ω_ｋ）を求めるためにｚ_ｉ（ω_ｋ）とＴ（ω_ｋ）を式（３２）、（３３）によって定義され、ステップＳ１０３において作成される。In order to obtain the mixing matrix B (ω _k ) by the iterative method, z _i (ω _k ) and T (ω _k ) are defined by the equations (32) and (33), and are created in step S103.

ステップＳ１０４においてｇ_ｊ（ω_ｋ）を求める際、ｇ_ｊ（ω_ｋ）（ｊ≠ｉ）を式（３４）のように、ｚ_ｊ（ω_ｋ）を式（３４）のようにそれぞれ定数に設定して、Ｆ_ｊ（ω_ｋ）を計算する。ここで、ｇ_ｊ（ω_ｋ）はＧ（ω_ｋ）のｊ番目の列ベクトルである。When obtaining g _j (ω _k ) in step S104, g _j (ω _k ) (j ≠ i) is set to a constant as shown in equation (34), and z _j (ω _k ) is set as a constant as shown in equation (34). Set and calculate F _j (ω _k ). Here, g _j (ω _k ) is the j-th column vector of G (ω _k ).

式（２６）の制約条件付きフォワードモデル型最小２乗型同時対角化問題を式（３５）のように書き直すことができ、その近似解ｇ_ｊ（ω_ｋ）はステップＳ１０５において式（３６）のラグランジェの未定乗数法によって求められる。The forward model type least squares simultaneous diagonalization problem with constraints in Equation (26) can be rewritten as in Equation (35), and its approximate solution g _j (ω _k ) is expressed in Equation (36) in step S105. The Lagrange's undetermined multiplier method.

ここで、ｕｎｖｅｃ｛Ａ｝は、Ｊ^２×１の列ベクトルＡをＪ×Ｊの行列に変換することを表す。

Here, unvec {A} represents converting the J ² × 1 column vector A into a J × J matrix.

ステップＳ１０７においてｊ＝１，２，…，Ｎについて誤差の限界がε_Ｇの近似解ｇ_ｊ（ω_ｋ）を反復法によって推定した後、ステップＳ１０８において累乗法を１回用いて式（３７）を最小にするｂ_ｊ（ω_ｋ）を算出する。次いで、ステップＳ１０９において式（３８）のようにＢ（ω_ｋ）を特異値分解する。ここで、ｔｒ［Ａ］は行列Ａのトレースを表す。In step S107, an approximate solution g _j (ω _k ) with an error limit of ε _G for j = 1, 2,..., N is estimated by an iterative method, and then in step S108, the power method is used once to obtain equation (37). B _j (ω _k ) that minimizes is calculated. Next, in step S109, singular value decomposition is performed on B (ω _k ) as shown in equation (38). Here, tr [A] represents a trace of the matrix A.

ステップＳ１１０においてＢ（ω_ｋ）^ＨＢ（ω_ｋ）のランクがＮ未満のとき、式（２１）の評価量ｅ（ω_ｋ）を零にするために、Ｂ（ω_ｋ）を式（４４）の行列によって置き換える。ここで、ステップＳ１１１において正規直交基底ｖ_１（ω_ｋ），ｖ_２（ω_ｋ），…，ｖ_ｒ（ω_ｋ）によって張られる空間に直交する空間の正規直交基底ｖ_ｒ＋１（ω_ｋ），ｖ_ｒ＋２（ω_ｋ），…，ｖ_Ｎ（ω_ｋ）、同様に、ステップＳ１１２において正規直交基底ｕ_１（ω_ｋ），ｕ_２（ω_ｋ），…，ｕ_ｒ（ω_ｋ）によって張られる空間に直交する空間の正規直交基底ｕ_ｒ＋１（ω_ｋ），ｕ_ｒ＋２（ω_ｋ），…，ｕ_Ｎ（ω_ｋ）はそれぞれ求められる。When the rank of B (ω _k ) ^H B (ω _k ) is less than N in step S110, B (ω _k ) is changed to formula (44) in order to make the evaluation quantity e (ω _k ) of formula (21) zero. ) Matrix. Here, orthogonal bases _{v 1} in step _{S111 (ω k), v 2} (ω k), ..., v r orthonormal basis of the space orthogonal to the space spanned by _{(ω k) v r + 1} (ω k), v _{r + 2} (ω _k ),..., v _N (ω _k ), and similarly, stretched by orthonormal basis u ₁ (ω _k ), u ₂ (ω _k ),..., u _r (ω _k ) in step S112. The orthonormal bases u _{r + 1} (ω _k ), u _{r + 2} (ω _k ),..., U _N (ω _k ) of the space orthogonal to the space are obtained.

ここで、δ（ω_ｋ）＞０とする。ステップＳ１１３において追加される特異値δ（ω_ｋ）は、式（４４）の右辺の√Ｎ／（√Ｎ＋δ（ω_ｋ）Ｎ）によって条件ｔｒ［Σ（ω_ｋ）］＝√Ｎを満足するように設定される。

Here, δ (ω _k )> 0. The singular value δ (ω _k ) added in step S113 satisfies the condition tr [Σ (ω _k )] = √N according to √N / (√N + δ (ω _k ) N) on the right side of the equation (44). Is set as follows.

２．３ 対角行列の解法
ステップＳ１１４においてＢ（ω_ｋ）から分離行列Ｗ（ω_ｋ）を式（２２）の最小２乗型一般化逆行列によって求める。式（１４）、（１５）よりΛ（ω_ｋ，ｍ）の左からＷ（ω_ｋ）Ｂ（ω_ｋ）、右からＢ（ω_ｋ）^ＨＷ（ω_ｋ）^Ｈをそれぞれ乗算すると、式（４７）を得る。

2.3 Diagonal Matrix Solution In step S114, a separation matrix W (ω _k ) is _obtained from B (ω _k ) using the least square generalized inverse matrix of Equation (22). From equations (14) and (15), multiplying Λ (ω _k , m) from the left by W (ω _k ) B (ω _k ) and from the right by B (ω _k ) ^H W (ω _k ) ^H respectively (47) is obtained.

Ｐ_ｘ（ω_ｋ，ｍ）の推定値を使用して誤差Ψ（ω_ｋ，ｍ）を式（４８）で定義すると、式（４９）のバックワード型最小２乗型同時対角化問題に最小２乗法を適用すると、対角行列Λ（ω_ｋ，ｍ）はステップＳ１１５において式（５０）で推定される。式（５０）においてｄｉａｇ［Ａ］は行列Ａの対角行列を表す。Using the estimated value of P _x (ω _k , m) to define the error Ψ (ω _k , m) in equation (48), the backward least squares simultaneous diagonalization problem in equation (49) When the least square method is applied, the diagonal matrix Λ (ω _k , m) is estimated by Equation (50) in step S115. In the equation (50), diag [A] represents a diagonal matrix of the matrix A.

誤差の限界がε_Ｃの近似解ｇ_ｊ（ω_ｋ）と近似解Λ（ω_ｋ，ｍ）を推定するまで、上記のアルゴリズムはステップＳ１１６において繰り返される。The above algorithm is repeated in step S116 until an approximate solution g _j (ω _k ) and an approximate solution Λ (ω _k , m) with an error limit of ε _C are estimated.

３．パーミュテーション問題の解法
ステップＳ１１７において基準周波数ビンを複数選択し、基準周波数ビン間において電力比の相関に基づきパーミュテーション行列を推定するためにΞ（ω_ｋ）を式（５３）によりステップＳ１１８で算出する。3. Solution of Permutation Problem In step S117, a plurality of reference frequency bins are selected, and in order to estimate a permutation matrix based on the correlation of the power ratio between the reference frequency bins, Ξ (ω _k ) is expressed by equation (53) in step S118. Calculate with

ここで、Ｔｒ（・）は行列のトレースを表す。また、Ｑは行列の各行に１となる要素が１箇所、その他の要素は０で、１となる要素の位置が他の行と重複しない行列の集合である。

Here, Tr (•) represents a matrix trace. Q is a set of matrices in which one element is 1 in each row of the matrix, the other elements are 0, and the position of the element that is 1 does not overlap with other rows.

ステップＳ１１９において基準周波数ビン間で電力比の相関が最も大きいパーミュテーション行列を式（５４）によって推定する。  In step S119, the permutation matrix having the largest correlation of the power ratio between the reference frequency bins is estimated by Expression (54).

ステップＳ１２０において複数の基準周波数ビンから１つの基準周波数ビンを任意に選択し、ステップＳ１２１において全ての周波数ビン間で電力比の相関が最も大きいパーミュテーション行列を式（５４）によって推定する。  In step S120, one reference frequency bin is arbitrarily selected from the plurality of reference frequency bins, and in step S121, the permutation matrix having the largest correlation of the power ratio among all the frequency bins is estimated by equation (54).

ステップＳ１２２において選択された全ての基準周波数ビンが１度選択されるまで、上記のパーミュテーション行列の推定手順を繰り返す。この結果、基準周波数ビンを除き、各周波数ビンに複数のパーミュテーション行列が割り当てられることになる。  The permutation matrix estimation procedure is repeated until all the reference frequency bins selected in step S122 are selected once. As a result, a plurality of permutation matrices are assigned to each frequency bin except for the reference frequency bin.

ステップＳ１２３において多数決によって周波数ビンに割り当てられるパーミュテーション行列を決定する。ただし、多数決によってパーミュテーション行列が決定できない場合には、最も相関値が大きいパーミュテーション行列を採用する。式（５５）のように観測信号ｘ（ω_ｋ，ｍ）にΠ（ω_ｋ）Ｗ（ω_ｋ）を左から乗算して分離信号ｙ（ω_ｋ，ｍ）を得る。In step S123, the permutation matrix assigned to the frequency bin is determined by majority vote. However, when the permutation matrix cannot be determined by majority vote, the permutation matrix having the largest correlation value is adopted. As shown in the equation (55), the observation signal x (ω _k , m) is multiplied by Π (ω _k ) W (ω _k ) from the left to obtain the separated signal y (ω _k , m).

４．１ 評価データ
図５のように４．４５×３．５５×２．５メートルの部屋に３個の信号源（スピーカ）１１、１２、１３を半径１．２メートルの円の円周上に、円の中心に位置する一辺が２０センチメートルの正三方形の頂点に３個のマイクロホン２１、２２、２３をそれぞれ配置した。尚、図５は信号源（スピーカ）とマイクロホンの位置関係を示す平面図である。部屋の残響時間は１００ミリ秒から９００ミリ秒に設定し、標本化周波数８ｋＨｚ、量子化ビット数１６ビットで信号源とマイクロホンの間のインパルス応答は人工的に発生させた。実験条件は、１０００秒の音声データ、Ｋ＝８１９２点の短時間フーリエ変換、エポック当たり重複率９９％の２３個のフレームの使用、窓関数にはハニング窓を用いた。ＳＮＲは５ｄＢ間隔で０〜３０ｄＢの範囲で変化させた。マイクロホン２１、２２、２３のＳＮＲの設定方法については４．２で説明する。本発明に係るブラインド信号分離方法では、ε_Ｇ＝ε_Ｃ＝１０^−６、δ（ω_ｋ）＝σ_ｒ（ω_ｋ）を用いている。スケーリング問題は周波数ビン毎に分離行列の行ベクトルを正規化することによって解法した。Ｃ言語で作成したプログラムをインテル製コアｉ７−２６００ ３．４ＧＨｚプロセッサを用いて実行した。信号源信号からマイクロホンまでの経路は時不変のインパルス応答で、因果的で非最小位相系であるので、因果的な分離行列を実現するために、Π（ω_ｋ）^−１Ｄ（ω_ｋ）^−１Ｗ（ω_ｋ）にｅ^−ｊπｋを乗算した後、逆離散フーリエ変換をして分離フィルタのインパルス応答を得た。4.1 Evaluation Data As shown in Fig. 5, three signal sources (speakers) 11, 12, and 13 are placed on a circle with a radius of 1.2 meters in a 4.45 x 3.55 x 2.5 meter room. In addition, three microphones 21, 22, and 23 are respectively arranged at the apexes of a regular triangle having a side of 20 centimeters located at the center of the circle. FIG. 5 is a plan view showing the positional relationship between the signal source (speaker) and the microphone. The reverberation time of the room was set from 100 milliseconds to 900 milliseconds, the sampling frequency was 8 kHz, the number of quantization bits was 16 bits, and the impulse response between the signal source and the microphone was artificially generated. The experimental conditions were 1000 seconds of speech data, a short-time Fourier transform of K = 8192 points, the use of 23 frames with 99% overlap per epoch, and a Hanning window as the window function. The SNR was varied in the range of 0 to 30 dB at 5 dB intervals. A method for setting the SNR of the microphones 21, 22, and 23 will be described in 4.2. In the blind signal separation method according to the present invention, ε _G = ε _C = 10 ⁻⁶ and δ (ω _k ) = σ _r (ω _k ) are used. The scaling problem was solved by normalizing the row vector of the separation matrix for each frequency bin. A program created in C language was executed using an Intel Core i7-2600 3.4 GHz processor. Since the path from the source signal to the microphone is a time-invariant impulse response and is a causal and non-minimum phase system, in order to realize a causal separation matrix, Π (ω _k ) ⁻¹ D (ω _k ) After multiplying ⁻¹ W (ω _k ) by e ^−jπk , inverse discrete Fourier transform was performed to obtain an impulse response of the separation filter.

４．２ 評価指標
ブラインド信号分離方法の信号分離性能を次の方法で評価した．式（５６）によって観測信号における所望信号源信号と干渉信号の電力の比、式（５７）によって出力信号における所望信号源信号と干渉信号の電力の比をそれぞれ計算し、ブラインド信号分離装置の各出力の信号分離性能を求める。各出力の平均を信号分離性能とした。γ_ｉｊ（ｔ）は式（５８）のΓ（ω_ｋ）のｉ行ｊ列の要素を、ｗ_ｉｊ（ｔ）はＷ（ω_ｋ）の要素をそれぞれ離散逆フーリエ変換したものである。また、分離行列の推定アルゴリズムにおいて収束に要した反復回数と計算時間も評価指標とする。ＳＮＲは、最適な分離行列ｅ^−ｊπｋＤ（ω_ｋ）^−１（Ｈ（ω_ｋ）^ＨＨ（ω_ｋ））^−１Ｈ（ω_ｋ）^Ｈとパーミュテーション行列Π_ｏｐｔ（ω_ｋ）を使用して観測信号から信号源信号を分離した後、分離信号ｙ_ｉ（ｔ）に含まれる雑音と干渉信号の電力と所望信号源信号の電力の比によって計算した。最適なパーミュテーション行列Π_ｏｐｔ（ω_ｋ）は式（５９）によって求めた。また、非ブラインド法は、受信信号を使用して分離行列を計算した後、混合行列が入手可能であるとして、式（６０）によってパーミュテーション行列を求めた。即ち、推定した分離行列に最適なパーミュテーション行列を求めることになり、ブラインド信号分離装置の性能の上限を与えることになる。4.2 Evaluation index The signal separation performance of the blind signal separation method was evaluated by the following method. The ratio of the power of the desired signal source signal and the interference signal in the observation signal is calculated by Expression (56), and the ratio of the power of the desired signal source signal and the interference signal in the output signal is calculated by Expression (57). Obtain the output signal separation performance. The average of each output was defined as the signal separation performance. γ _ij (t) is a discrete inverse Fourier transform of an element of i rows and j columns of Γ (ω _k ) of equation (58), and w _ij (t) is an element of W (ω _k ). In addition, the number of iterations required for convergence and the calculation time in the separation matrix estimation algorithm are also used as evaluation indexes. The SNR is obtained by ^calculating an optimal separation matrix e ^−jπk D (ω _k ) ⁻¹ (H (ω _k ) ^H H (ω _k )) ⁻¹ H (ω _k ) ^H and permutation matrix Π _opt (ω _k ). After the signal source signal was separated from the observation signal by use, calculation was performed by the ratio of the noise contained in the separated signal y _i (t), the power of the interference signal, and the power of the desired signal source signal. The optimal permutation matrix Π _opt (ω _k ) was obtained by Equation (59). In the non-blind method, a permutation matrix is obtained by using Equation (60), assuming that a mixing matrix is available after calculating a separation matrix using a received signal. That is, an optimum permutation matrix is obtained for the estimated separation matrix, and an upper limit of the performance of the blind signal separation device is given.

ここで、Ｃ_ｏｐｔ（ω_ｋ）＝ｅ^−ｊπｋＤ（ω_ｋ）^−１（Ｈ（ω_ｋ）^ＨＨ（ω_ｋ））^−１Ｈ（ω_ｋ）^ＨＨ（ω_ｋ）、Ｃ（ω_ｋ）＝Ｗ（ω_ｋ）Ｈ（ω_ｋ）である。

Here, C _opt (ω _k ) = e ^−jπk D (ω _k ) ⁻¹ (H (ω _k ) ^H H (ω _k )) ⁻¹ H (ω _k ) ^H H (ω _k ), C (ω _k ) = W (ω _k ) H (ω _k ).

４．３ 評価対象
勾配法を用いたバックワードモデル型ブラインド信号分離方法（非特許文献１）、最小２乗型同時対角化問題の解法を用いた２種類のフォワードモデル型ブラインド信号分離方法（非特許文献２、非特許文献３）を比較対象とする。従来のブラインド信号分離方法（非特許文献１、非特許文献２、非特許文献３）と本発明に係るブラインド信号分離方法における分離行列の推定精度を比較するため、パーミュテーション行列の推定法は共通の手法（非特許文献４）を使用した。尚、基準周波数ビンの番号には６１４を用いた。4.3 Evaluation target Backward model type blind signal separation method using gradient method (Non-Patent Document 1), two types of forward model type blind signal separation method using least square type simultaneous diagonalization problem solution ( Non-patent literature 2 and non-patent literature 3) are to be compared. In order to compare the estimation accuracy of the separation matrix in the conventional blind signal separation method (Non-Patent Document 1, Non-Patent Document 2, Non-Patent Document 3) and the blind signal separation method according to the present invention, the permutation matrix estimation method is: A common method (Non-Patent Document 4) was used. Note that 614 is used as the reference frequency bin number.

４．４ 評価結果
部屋の残響時間と信号分離性能の関係を図６に、ＳＮＲと信号分離性能の関係を図７にそれぞれ示す。尚、図６においてＳＮＲは２０ｄＢに設定している。両図において太字の数字が最も優れた性能を表している。両図の信号分離性能から明らかなように、計算時間では非特許文献１より劣るものの、本発明に係るブラインド信号分離方法が従来のブラインド信号分離方法よりも最も高い信号分離性能（高い出力ＳＩＲ）を最も少ない反復回数で得ることができた。この要因はラグランジェの未定乗数法を最小２乗型同時対角化問題に導入したこと、推定した混合行列がランク落ちしていた場合、補空間を補い分離行列を推定したことが高い信号分離性能の実現に貢献したと考えられる。また、非ブラインド法の出力ＳＩＲ、即ち、ブラインド信号分離装置の上限に近い値を、本発明に係るブラインド信号分離装置が実現できることが分かる。4.4 Evaluation Results FIG. 6 shows the relationship between the reverberation time of the room and the signal separation performance, and FIG. 7 shows the relationship between the SNR and the signal separation performance. In FIG. 6, the SNR is set to 20 dB. In both figures, bold numbers represent the best performance. As is apparent from the signal separation performance in both figures, the blind signal separation method according to the present invention has the highest signal separation performance (higher output SIR) than the conventional blind signal separation method, although the calculation time is inferior to that of Non-Patent Document 1. Was obtained with the least number of iterations. This is due to the fact that Lagrange's undetermined multiplier method has been introduced to the least-squares simultaneous diagonalization problem, and that if the estimated mixing matrix has dropped in rank, the supplementary space has been compensated and the separation matrix has been estimated. It is thought that it contributed to realization of performance. It can also be seen that the blind signal separation device according to the present invention can realize the output SIR of the non-blind method, that is, a value close to the upper limit of the blind signal separation device.

５．１ 評価データ
図５のように４．４５×３．５５×２．５メートルの部屋に３個の信号源（スピーカ）１１、１２、１３を半径１．２メートルの円の円周上に、円の中心に位置する一辺が２０センチメートルの正三方形の頂点に３個のマイクロホン２１、２２、２３をそれぞれ配置した。部屋の残響時間は７００ミリ秒に設定し、標本化周波数８ｋＨｚ、量子化ビット数１６ビットで信号源とマイクロホンの間のインパルス応答は人工的に発生させた。実験条件は、１０００秒の音声データ、Ｋ＝８１９２点の短時間フーリエ変換、エポック当たり重複率８０％の２個のフレームの使用、窓関数にはハニング窓を用いた。ＳＮＲは２０ｄＢに設定した。本発明に係るブラインド信号分離方法では、ε_Ｇ＝ε_Ｃ＝１０^−６、δ（ω_ｋ）＝σ_ｒ（ω_ｋ）、基準周波数ビンの番号は６１６、６１７、６１８を用いている。スケーリング問題は周波数ビン毎に分離行列の行ベクトルを正規化することによって解法した。信号源信号からマイクロホンまでの経路は時不変のインパルス応答で、因果的で非最小位相系であるので、因果的な分離行列を実現するために、Π（ω_ｋ）^−１Ｄ（ω_ｋ）^−１Ｗ（ω_ｋ）にｅ^−ｊπｋ乗算した後、逆離散フーリエ変換をして分離フィルタのインパルス応答を得た。5.1 Evaluation Data As shown in Fig. 5, three signal sources (speakers) 11, 12, and 13 are placed on a circle with a radius of 1.2 meters in a 4.45 x 3.55 x 2.5 meter room. In addition, three microphones 21, 22, and 23 are respectively arranged at the apexes of a regular triangle having a side of 20 centimeters located at the center of the circle. The room reverberation time was set to 700 milliseconds, the sampling frequency was 8 kHz, the number of quantization bits was 16 bits, and the impulse response between the signal source and the microphone was artificially generated. The experimental conditions were 1000 seconds of speech data, a short-time Fourier transform of K = 8192 points, the use of two frames with 80% overlap per epoch, and a Hanning window as the window function. The SNR was set to 20 dB. In the blind signal separation method according to the present invention, ε _G = ε _C = 10 ⁻⁶ , δ (ω _k ) = σ _r (ω _k ), and reference frequency bin numbers 616, 617, and 618 are used. The scaling problem was solved by normalizing the row vector of the separation matrix for each frequency bin. Since the path from the source signal to the microphone is a time-invariant impulse response and is a causal and non-minimum phase system, in order to realize a causal separation matrix, Π (ω _k ) ⁻¹ D (ω _k ) After multiplying ⁻¹ W (ω _k ) by e ^−jπk , inverse discrete Fourier transform was performed to obtain an impulse response of the separation filter.

５．２ 評価指標
信号源とマイクロホンの個数が共に３である場合、式（６１）に示す６種類のパーミュテーション行列の何れか１つが各周波数ビンに割り当てられる。割り当てられたパーミュテーション行列が、任意のパーミュテーション行列に一致する割合と信号分離性能を計算する。5.2 Evaluation Index When both the number of signal sources and microphones are 3, any one of six types of permutation matrices shown in Expression (61) is assigned to each frequency bin. The ratio of the assigned permutation matrix to an arbitrary permutation matrix and the signal separation performance are calculated.

５．３ 評価対象
同一信号源から発生した信号の周波数ビン間の電力比に相関があることを利用したパーミュテーション問題の解法（非特許文献４）を比較対象とする。基準周波数ビンの番号には６１４を用いた。従来のブラインド信号分離方法（非特許文献１、非特許文献２、非特許文献３）と本発明に係るブラインド信号分離方法におけるパーミュテーション行列の推定精度を比較するため、分離行列の推定法は共通の手法（本発明に係る分離行列推定法）を使用した。5.3 Evaluation Target The solution to the permutation problem (non-patent document 4) using the fact that there is a correlation in the power ratio between frequency bins of signals generated from the same signal source is the comparison target. 614 was used as the reference frequency bin number. In order to compare the estimation accuracy of the permutation matrix in the conventional blind signal separation method (Non-Patent Document 1, Non-Patent Document 2, Non-Patent Document 3) and the blind signal separation method according to the present invention, the separation matrix estimation method is A common method (separation matrix estimation method according to the present invention) was used.

５．４ 評価結果
図８と図９に各周波数ビンに割り当てられたパーミュテーション行列の番号を×印で示す。パーミュテーション行列Π_ｉとパーミュテーション行列の番号ｉの関係を式（６１）に示している。図８と図９では、各周波数ビンでパーミュテーション行列の番号３に割り当てられると未知信号源への割り当てが揃うことになる。したがって、番号３を除く他の番号への割り当ては間違いになる。低周波数帯域（０〜２ｋＨｚ）と全周波数帯域におけるパーミュテーション行列の割り当て結果を図１０にまとめる。本発明に係るパーミュテーション行列の推定法が非特許文献４の方法に比べ正答率が向上していることが分かる。また、信号分離性能においても、本発明に係るパーミュテーション行列の推定法が高い出力ＳＩＲを達成することができた。5.4 Evaluation Results FIGS. 8 and 9 show the permutation matrix numbers assigned to the frequency bins with x marks. The relationship between the permutation matrix _{ｉ i} and the permutation matrix number i is shown in Equation (61). In FIG. 8 and FIG. 9, when the frequency bin is assigned to the permutation matrix number 3, the assignment to the unknown signal source is completed. Therefore, the assignment to other numbers except the number 3 is wrong. FIG. 10 summarizes the assignment results of the permutation matrix in the low frequency band (0 to 2 kHz) and the entire frequency band. It can be seen that the permutation matrix estimation method according to the present invention improves the correct answer rate compared to the method of Non-Patent Document 4. Also in terms of signal separation performance, the permutation matrix estimation method according to the present invention was able to achieve a high output SIR.

１１〜１Ｎ…信号源、２１〜２Ｊ…マイクロホン、３１…短時間フーリエ変換、４０…ブラインド信号分離システム、５１、５２、…、５Ｋ…最小２乗型同時対角化問題の解法、６０…パーミュテーション問題の解法、７１、７２、…、７Ｋ…畳み込み演算、８０…離散逆フーリエ変換と重複加算11 to 1N ... signal source, 21 to 2J ... microphone, 31 ... short-time Fourier transform, 40 ... blind signal separation system, 51, 52, ..., 5K ... solution of least squares simultaneous diagonalization problem, 60 ... par Solution of mutation problem, 71, 72,..., 7K ... convolution, 80 ... discrete inverse Fourier transform and overlap addition

Claims (4)

未知の畳み込み混合系により混在した互いに統計的に独立な未知信号源信号を、観測信号のみからブラインドで推定する方法であって、時間周波数領域においてエポック時刻毎に観測信号ベクトルから空間相関行列を求め、各エポック時刻において周波数ビンの中からフロベニウスノルムが最大となる空間相関行列を求め、そのノルムで全ての同エポック時刻の空間相関行列を正規化した後、正規化された空間相関行列から混合行列を推定するために、制約条件付きフォワードモデル型最小２乗型同時対角化問題をラグランジュの未定乗数法と反復法の組み合わせにより解法する手順と、手順の収束後、累乗法を１回用いて最小２乗型同時対角化問題の解となる混合行列を求め、混合行列のランクが落ちていた場合は、混合行列がフルランクになるように特異値分解を用いて基底を補い、次いで、ランク落ちした混合行列の特異値とその総和が変わらないように特異値を補い、最小２乗型一般化逆行列を用いて混合行列から分離行列を推定し、バックワード型最小２乗型同時対角化問題を最小２乗法で解いて対角行列を求め、同一信号源から発生した信号の周波数ビンの電力比に相関があることを利用してパーミュテーション問題を、多数決を利用して解法することによって信号分離精度の高い分離行列を求めることを特徴とするブラインド信号分離方法。  A method of blindly estimating unknown signal source signals that are statistically independent from each other mixed by an unknown convolutional mixing system from the observed signal alone, and obtains the spatial correlation matrix from the observed signal vector at each epoch time in the time-frequency domain. Find the spatial correlation matrix that maximizes the Frobenius norm from the frequency bins at each epoch time, normalize the spatial correlation matrix of all the same epochs with that norm, and then mix the mixed matrix from the normalized spatial correlation matrix To solve the constrained forward model least-squares simultaneous diagonalization problem using a combination of Lagrange's undetermined multiplier method and iterative method, and after convergence of the procedure, use the power method once Find the mixing matrix that solves the least-squares type diagonalization problem, and if the mixing matrix rank has dropped, the mixing matrix is full rank. The basis is supplemented using singular value decomposition, and then the singular value of the mixed matrix whose rank is reduced and the sum of the singular values are not changed. Then, the least square type generalized inverse matrix is used to calculate the singular value from the mixing matrix. Estimate the separation matrix, solve the backward least-squares simultaneous diagonalization problem using the least-squares method, find the diagonal matrix, and confirm that the power ratio of the frequency bins of signals generated from the same signal source is correlated. A blind signal separation method characterized by obtaining a separation matrix having high signal separation accuracy by solving a permutation problem using majority voting. 請求項１記載の最小２乗型同時対角化問題の解法を用いたブラインド信号分離方法において、信号源とブラインド信号分離装置出力間の伝達関数を表した混合行列と分離行列の縦続接続モデルにラグランジュの未定乗数を導入した制約条件付きフォワードモデル型最小２乗型同時対角化問題とバックワード型最小２乗型同時対角化問題を導入し、更に分離行列に遅延を与えることによって信号源とブラインド信号分離装置出力間の因果的なモデルを推定することを特徴とするブラインド信号分離方法。  A blind signal separation method using the method of solving the least-squares simultaneous diagonalization problem according to claim 1, wherein a mixed matrix and a separation matrix cascade model representing a transfer function between a signal source and a blind signal separation device output are used. Introducing Lagrange's undetermined multiplier with constraint forward model least squares simultaneous diagonalization problem and backward least squares simultaneous diagonalization problem, and further delaying the separation matrix And a blind signal separation device output, a causal model is estimated. 請求項１乃至請求項２のいずれか１項に記載の最小２乗型同時対角化問題の解法を用いたブラインド信号分離方法において、基準周波数ビンを複数選択し、基準周波数ビン間において最も電力比の相関が大きいパーミュテーション行列を推定し、複数の基準周波数ビンから基準周波数ビンを１つ選択し、全ての周波数ビン間で電力比の相関が最も大きいパーミュテーション行列を推定する手順を、全ての基準周波数ビンが１度選択されるまで繰り返した後、各周波数ビンに複数割り当てられたパーミュテーション行列から多数決によって周波数ビンに割り当てられるパーミュテーション行列を決定し、多数決によってパーミュテーション行列が決定できない場合には、最も相関値が大きいパーミュテーション行列を採用することを特徴とするブラインド信号分離方法。  The blind signal separation method using the method of solving the least square simultaneous diagonalization problem according to any one of claims 1 to 2, wherein a plurality of reference frequency bins are selected, and the most power among the reference frequency bins is selected. Estimating a permutation matrix having a large ratio correlation, selecting one reference frequency bin from a plurality of reference frequency bins, and estimating a permutation matrix having the largest power ratio correlation among all frequency bins. After all reference frequency bins are selected once, the permutation matrix assigned to the frequency bin is determined by majority from the permutation matrix assigned to each frequency bin, and permutation is performed by majority If the matrix cannot be determined, use the permutation matrix with the largest correlation value. Rind signal separation method. 請求項１乃至請求項３のいずれか１項に記載のブラインド信号分離方法を用いて信号源分離を行うように構成されていることを特徴とするブラインド信号分離方法を用いたブラインド信号分離装置。  A blind signal separation device using a blind signal separation method, wherein the blind signal separation method is configured to perform signal source separation using the blind signal separation method according to any one of claims 1 to 3.

Images