JP2022120750A - Method of separating signal waveform per incoming direction - Google Patents

Abstract

To enable a signal waveform to be separated per incoming direction.SOLUTION: When source signal space advance information is known, using an objective function having a constraint condition including information that a mixed matrix estimated from the source signal space advance information is an inverse matrix of a separation matrix, the separation matrix is updated per matrix block, and then a source signal is estimated from the obtained separation matrix and an observation signal. When the source signal space advance information is unknown, the source signal is estimated using advance information concerning statistical dispersion of the observation signal and the source signal, and then, the degree of attribution and a model parameter are alternately updated using a signal base vector as a feature quantity and assuming a von Mises-Fisher distribution as the mixed model, thereby solving a conversion indefiniteness problem of the estimated source signal.

Description

The present invention relates to a method for separating signal waveforms by direction of arrival, and in particular, a method for calculating separated signals based on source signal direction-of-arrival prior information (original signal space prior information) and a clustering method using von Mises-Fischer distribution as a mixture model. Regarding.

Problems to be Solved by the Invention

In a general road environment, the engine noise and tire rolling noise of vehicles, and in road bridges, the joint noise generated from the bridge structure due to the impact of a vehicle entering the bridge, and the low frequency vibration caused by the bending vibration of the girder caused by the running of the vehicle. Wave sound is a problem (Non-Patent Document 1). In particular, it has been reported that low-frequency sound generated from bridge structures poses a health hazard (Non-Patent Document 2), and it is important to identify the source and take countermeasures. Also, in residential areas, strange noises such as falling objects and operating noises of fans, which are difficult to identify, have become a problem. In order to investigate the source of noise and strange sounds, it is necessary to identify the sound in question from a large amount of observation data, which takes a lot of labor. In addition to this, in order to help predict construction noise, etc., construction machine sounds are generated from speakers, and multiple microphone arrays are used to specify the direction of arrival of the sounds, thereby performing noise prediction for each machine. There are some cases, but they are limited to direction identification only. Artificial intelligence is used to reduce labor when investigating the source of noise by calculating the frequency characteristics of sound for each direction, and when identifying problematic sounds from among multiple environmental sounds. In this case, the feature quantity of the signal waveform is required. At this time, if waveforms in the time domain or frequency domain for each direction of arrival of the source signal can be obtained, it will be possible to extract the feature amount of the target sound with higher accuracy, leading to highly accurate judgment. Conceivable.

In the field of anomaly detection, a method has been proposed to determine anomalies by extracting power spectra and waveform shape features from observed sound time waveforms and inputting them to artificial intelligence (Non-Patent Document 3). ). This abnormal sound detection method is used for various purposes, including machine failure prediction, crime detection by detecting gunshots, screams, and explosion sounds, and accident detection by detecting car crashes and slip sounds on roads. (Non-Patent Documents 4, 5, 6). If it is possible to separate waveforms by direction in advance, it is possible to have artificial intelligence learn the target waveform feature quantity and its direction frequency with higher accuracy, and to judge abnormalities based on the learning data. is considered to be

Subspace methods such as the music method (Non-Patent Document 7) have often been used in conventional noise and abnormal noise source investigations. For example, the music method is a method of obtaining the MUSIC spectrum for each direction based on the steering vector information determined from the correlation matrix of signals obtained at a plurality of observation points and the sensor arrangement. However, since the MUSIC spectrum expresses the ratio of relative sound levels in each direction, it is possible to evaluate the intensity ratio of incoming waves in each direction. A time waveform cannot be obtained. Using the music method to obtain time-series data of separated waveforms from observed signals in which multiple signal waveforms are mixed, it is necessary to estimate the mixing matrix from information on the direction of arrival of the source signal and use its inverse matrix as the separation matrix. be. However, in this method, the mixing matrix is ill-conditioned when the source signals exist close to each other, so the separation signal calculation accuracy is poor. Even when obtaining separated waveforms of wideband signals, the mixing matrix becomes ill-conditioned on the low frequency side. Therefore, in the method of calculating the inverse matrix of the mixing matrix estimated by the music method, separated signals can only be obtained in a situation where all the source signals are sufficiently separated and the frequency to be analyzed is limited to a narrow band. .

In the separation of music signals and speech signals, a technique for extracting a target sound from a plurality of mixed sounds and a technique for separating mixed sound signals by factors have been developed. Typical methods include independent component analysis and independent low-rank matrix analysis. These source signal separation methods are methods of separating signals by assuming statistical independence of the source signals. Even when the separation matrix is calculated by using independent component analysis or independent low-rank matrix analysis in the frequency domain, the estimation accuracy of the separated signals on the low frequency side of the wideband signal deteriorates, but the degree of deterioration depends on the mixing matrix. less than when using the inverse of Therefore, it can be applied to the case where the source signals exist close to each other or to calculate the separated signal from the broadband signal.

In the independent component analysis, separated signals are obtained by assuming the probability distribution shape of the source signal. However, since the probability distribution cannot be strictly defined, the signal separation results always contain a certain level of error. Further, when the number of source signals is larger than the number of observation points, there is a property that signals cannot be separated correctly, and when a plurality of noises other than the target sound are mixed, there is a problem that the signal separation accuracy deteriorates.

Independent component analysis and independent low-rank matrix analysis have a partition problem that the order in which the estimated values of the source signal appear cannot be determined. There is a problem that If the signal separation result contains an error, the directional feature obtained from the separation result is also affected by the error, making it difficult to solve the partition problem using the directional feature.

For the direction of arrival of the source signal, the tensor component of the mixing matrix, which is the inverse matrix of the separation matrix, is used to calculate the frequency transfer function from the source signal position to the sensor position, and the inter-sensor phase difference of the source signal is estimated from the frequency transfer function. Calculated by Therefore, if the calculation result of the separation matrix contains an error, the accuracy of the direction-of-arrival estimation result of the source signal deteriorates due to its influence.

In this way, the calculation error of the separation matrix is not limited to the error of the separated signals for each frequency bin, but also leads to the error of the cluster analysis of the signal basis vectors and the calculation error of the direction of arrival of the source signal.

When using independent component analysis, if the probability distribution shape can be determined strictly, all the above problems can be solved. The result of signal separation contains a certain level of error.

As methods for solving partitions, there have been methods that use the directivity pattern of the separation matrix as a clue (Non-Patent Documents 8 and 9) and a method that uses the continuity between frequencies of the inverse matrix of the separation matrix (Non-Patent Document 10). is proposed. However, these methods have different signal basis vectors for different frequencies, making it difficult to separate broadband signals. To solve this problem, a method has been proposed to construct consistent source signal labels at all frequencies by approximating the inverse matrix of the separation matrix with the frequency transfer function between the source signal and the observation point (non-patent Reference 11). In the above method, when a calculation error occurs in the separation matrix (separation signal) due to factors such as noise in the observation data, the signal basis vectors for each cluster are also affected by the calculation error. There is a problem that the optimum solution cannot be reached. In addition, it is necessary to determine the initial value of the average direction for each cluster of signal basis vectors at the time of cluster calculation. be. This initial value dependence appears conspicuously when a calculation error occurs in the separation matrix (separated signal).

Until now, independent component analysis has been utilized for signal separation, but a method for estimating the direction of arrival of the source signal with high accuracy after signal separation has not been established. A possible cause of this is that, in the signal arrival direction estimation method based on the separation matrix, the accuracy of the source signal arrival direction estimation result is degraded due to the influence of the error included in the separation matrix calculation result. Therefore, until now, the subspace method such as the music method has often been used for the direction of arrival of the source signal, and the independent component analysis has been used for the signal separation.

When the subspace method such as the music method and the camera image are used together, the direction of arrival of the source signal can be estimated with high accuracy. When the spatial information of the source signal is obtained in advance by such a method, the frequency transfer function between the sensor and the source signal can be obtained by actively using this spatial prior information. A mixing matrix is obtained using this frequency transfer function, and a separating matrix can be calculated from its inverse matrix. Theoretically, if this separation matrix is used and spatial prior information is used for signal separation, signals can be separated with high accuracy. However, when the analysis target frequencies of the source signals cover a wide band, or when the source signal positions are close to each other, the mixing matrix becomes a bad condition, and the separation matrix calculation accuracy deteriorates. At present, no method has been established to solve this problem and to perform signal separation with high accuracy using spatial prior information of signal waveforms.

Shinya Kimura, Kazuyuki Ono, Shun Kouji, Tetsusuke Kin, Mitsuro Kawatani: Study on Low-Frequency Sound Propagation Accompanying Traffic Vibration on Road Bridges, Structural Engineering Journal, Vol. 64A, pp307-314, 2018. Shunichi Nakamura, Yasuo Tokita, Atsushi Oda: Fundamental Research on Sensation and Evaluation of Low-Frequency Sound, 1979 Ministry of Education Scientific Research Fund “Environmental Science” Special Research, 1979. A. Ito, A. Aiba, M.; Ito and S.I. Makino, "Detection of abnormal sound using multi-stage GMM for surveillance microphone," Proc. Int. Conf. Information Assurance and Security, Vol. 1, pp. 733-736 (2009). C. N. Doukas and I. Maglogiannis, "Emergency fall incidents detection in assisted living environments utilizing motion, sound, and visual perceptual components," IEEE Trans. Inf. Technol. Biomed. , 15, 277-289 (2011). S. Lecomte, R.; Lengell'e, C.; Richard, F.; Capman andB. Ravera, "Abnormal events detection using unsupervised one-class SVM: Application to audio surveillance and evaluation," Proc. 8th IEEE Int. Conf. Advanced Video and Signal Based Surveillance, pp. 124-129 (2011). E. Marchi, F.; Vesperini, F.; Eyben, S.; Squartini and B. Schuller, "A novel approach for automatic acoustic novelty detection using a denoising autoencoder with bidirectional LSTM neural networks," Proc. ICASSP 2015, pp. 1996-2000 (2015). P. Stoica and A. Nehorai, "MUSIC, maximum likelihood, and Cramer-Rao bound," IEEE Trans. Acoust. Speech Signal Process. vol. 37, no. 5, pp. 720-741, May 1989. H. Saruwatari, S.; Kurita, K.; Takeda, F.; Itakura, T.; Nishikawa, KiyohiroShikano: Blind Source Separation Combining Independent Component Analysis and Beamforming, EURASIP Journal on Applied Signal Processing, 1135-1146, Nov. Jan. 2003. M. Z. Ikram and D. R. Morgan: Permutation inconsistency in blind speech separation: investigation and solutions, IEEE Trans. Speech and Audio Processing, Vol. 13, no. 1, pp. 1-13, 2005. F. Asano, S.; Ikeda, M.; Ogawa, H.; Asoh, and N.L. Kitawaki: Combined approach of array processing and independent component analysis for blind separation of acoustic signals, IEEE Trans. Speech and Audio Processing, Vol. 11, no. 3, pp. 204-215, May. 2003. H. Sawada, S.; Araki, R.; Mukai, S.; Makino: Blind extraction of a dominant source signal from mixtures of many sources, "International Conference on (Acoust Speech Signal Process), Vol. 3, 61-64) (2000). Hiroshi Sawada, Ryo Mukai, Akiko Araki, Shoji Makino: Frequency Domain Blind Source Separation for Multiple Sound Sources, The Japanese Society for Artificial Intelligence, JSAI Technical Report, SIG-CHallege-0522-3 (10/14), 2005. H. Sawada, H.; Kameoka, S.; Araki, andN. Ueda, "Multichannel extensions of non-negative matrix factorization with complex-valued data," IEEE Trans. ASLP, vol. 21, No. 5, pp. 971-982, 2013. MNMF1 Daichi Kitamura, Nobutaka Ono, Hiroshi Sawada, Hirokazu Kameoka, Hiroshi Saruwatari, "Blind source separation based on independent low-rank matrix analysis," IEICE Technical Report, EA2017-56, vol. 117, no. 255, pp. 73-80, Oct. 2017. N. Q. K. Duong, E. Vincent, and R. Gribonval, "Under determined reverberant audio source separation using a full rank spatial covariance model," IEEE Trans. ASLP, vol. 18, no. 7, pp. 1830-1840, 2010. MNMF2 C. F.e.votte, N.E. Bertin, andJ. -L. Durrieu, “Nonnegative matrix factorization with the Itakura-Saito divergence. With application to music analysis,” Neural Computation, vol. 21, no. 3, pp. 793-830, 2009. D. Kitamura, N.; Ono, H. Sawada, H.; Kameoka, and H. Saruwatari, “Determined blind source separation unifying independent vector analysis and nonnegative matrix factorization,” IEEE/ACM Trans. ASLP, vol. 24, no. 9, pp. 1626-1641, 2016. N. Ono, "Stable and fast update rules for independent vector analysis based on auxiliary function technique," Proc. WASPAA, pp. 189-192, 2011. Daichi Kitamura, Hayato Kadono, Norigen Takamune, Shinnosuke Takamichi, Hiroshi Saruwatari, Nobutaka Ono, "Experimental evaluation of multi-channel source separation based on independent deep learning matrix analysis," IEICE Technical Report, EA2017-104, vol. 117, no. 515, pp. 13-20, Mar. 2018. Hironori Mukumoto, Kazunori Hayashi, Megumi Kaneko, “Arrival angle estimation method by compression sensing using Khatri-Rao product extended array,” IEICE Technical Report, RCC2014-11, MICT2014-11, pp. 51-56, May 2014. J. Chen and X. Huo, "Theoretical results on sparse representations of multiple measurement vectors," IEEE Trans. Signal Process. , vol. 54, no. 12, pp. 4634-4643, 2006. W. -K. Ma, T. -H. Hsieh, and C.I. -Y. Chi, “DOA estimation of quasi-stationary signals with less sensors than sources and unknown spatial noise covariance: a KhatriRao substrate approach,” IEEE Trans. Signal Process. , vol. 58, no. 4, pp. 2168-2218, 2010. J. M. Kim, O. K. Lee, andJ. C. Ye, "Compressive MUSIC: revisiting the link between compressive sensing and array signal processing," IEEE Trans. Inf. Theory, vol. 58, no. 1, pp. 278-301, Jan. 2012.

SUMMARY OF THE INVENTION The present invention has been made in view of these points, and its object is to provide a method for separating signal waveforms by direction of arrival for calculating the directions of arrival of separated signals and source signals.

The purpose is to separate the signal waveforms according to the direction of arrival in an environment where multiple source signals coexist, and to estimate the direction of arrival of the source signals.

Means to solve problems

In order to solve the above-described problems and achieve highly accurate signal separation and source signal arrival direction estimation, the present invention uses source signal space prior information (original signal arrival direction prior information) to create a separation matrix with high accuracy. A separated signal calculation section that can be calculated, a clustering section that is not easily affected by separation matrix calculation errors, and a source signal arrival direction calculation section that estimates the direction of arrival of the source signal from the separation matrix calculation result, and corrects separated signals with low calculation accuracy. It is an object of the present invention to provide a calculation method characterized by having a separation signal corrector.

The direction-specific signal waveform separation method of the present invention includes a separation signal calculation unit capable of improving the separation matrix calculation accuracy by using the direction of arrival as prior information in frequency bins where the direction of arrival of the source signal is known, Using the von Mises-Fisher distribution, with the goal of maximizing its probability density function, we evaluated the accuracy of the clustering unit that clusters the signal basis vectors and the calculation accuracy of the separation matrix (separation signal) for each frequency bin and source signal. , the value obtained by multiplying the importance index of each frequency bin such as the power spectrum of the observed signal for each frequency bin by the calculation accuracy evaluation value is used as a weight value, and the weight of the direction of arrival of the source signal for each frequency bin is obtained for all frequency bins. It is characterized by having a source signal direction-of-arrival direction calculation unit capable of obtaining a highly accurate direction estimation result by taking an average.

In the signal separation calculation of this patent, the calculation performed by the separation signal calculation unit is selectively changed within the system depending on whether the spatial prior information of the source signal is obtained or not, and the spatial prior information is not obtained, the calculations of the clustering section, source signal direction-of-arrival calculating section, and separated signal correcting section are performed.

FIG. 2 is a schematic diagram of an apparatus for performing the operation of the directional signal waveform separation method; 4 is a flow chart showing the flow of arithmetic processing of the direction-specific signal waveform separation method. 4 is a flowchart of arithmetic processing for roughly calculating a source signal waveform using source signal arrival direction prior information; FIG. 10 is a flowchart of arithmetic processing for calculating the source signal waveform in detail by updating the source signal arrival direction a priori information, the separation matrix, and the estimated value of the source signal statistical variance; FIG. 4 is a flow chart of arithmetic processing performed by an initial clustering unit in base vector clustering operations performed by the clustering unit. FIG. 10 is a flowchart of arithmetic processing performed by a cluster determining unit in base vector clustering operations performed by the clustering unit; FIG. 4 is a flowchart of arithmetic processing performed by a source signal arrival direction calculating unit;

An embodiment of the direction-specific signal separation calculation according to the present invention will be described below with reference to the drawings. FIG. 1 is a schematic diagram of an apparatus for performing the direction-of-arrival signal waveform separation method. As shown in FIG. 1, in this direction-of-arrival signal waveform separation apparatus for executing the direction-specific signal separation method of the present embodiment, a signal is acquired using a microphone or the like, and signal data is transmitted to an observed signal separation arithmetic unit. . The observed signal separation calculation unit performs calculations for estimating the source signal waveforms for each direction of arrival and their directions of arrival from the observed signals. After obtaining the direction of arrival of the source signal and the corresponding source signal waveform by the calculations performed by the observed signal separation calculation device composed of the direction of arrival calculator and the separated signal corrector, the calculation results are transferred to a general-purpose computer or a dedicated device. Display using software.

FIG. 2 is a calculation flow executed by the observed signal separation calculation device in the direction-specific signal separation calculation in this embodiment. As shown in FIG. 2, the observation signal separation calculation in this embodiment is divided into cases when the direction of arrival of the source signal is known in advance and when it is not known. When the direction of arrival of the source signal is unknown, the separation signal calculator separates the observed signal into source signal waveforms for each direction of arrival. The signal corrector corrects the estimation error to estimate the source signal waveform, and at the same time, the source signal arrival direction calculator estimates their directions of arrival. When the spatial information of the source signals is obtained in advance, the separated signal calculating section calculates the source signals for each direction of arrival with higher accuracy using the prior information of the direction of arrival of the source signals.

[Signal Waveform Acquisition Unit]
The signal waveform acquisition unit is composed of two or more microphones or the like. It receives signal waveforms and sends signal data to an observed signal separation arithmetic unit.

[Waveform Separation Calculation in Observation Signal Separation Calculator]
An overview of the waveform separation calculation performed by the observation signal separation calculation device is given below. Details will be shown later.

[Separation signal calculator]
The waveform data acquired by the receiving sensor is converted into the frequency domain, and then the separation matrix is obtained by the separation signal calculator. The separation matrix is generated using independent component analysis, independent low-rank matrix analysis, or the like. In these methods, the KL information amount of the product of the joint distribution and the marginal distribution of the output of the separation operation is set as the objective function, and the separation matrix is updated and calculated so as to minimize the objective function.

In some cases, the direction of arrival of the source signal is known in advance, such as when the position of the source signal can be visually confirmed, or when the direction of arrival of the source signal is obtained by applying the music method, etc. to a specific frequency bin in advance. be. Thus, when the source signal arrival direction is known in advance, the source signal waveform can be estimated with high accuracy by using the prior information. FIG. 3 shows the calculation flow of the separation matrix and separated signals when the spatial a priori information of the source signal is obtained. First, using the information obtained from the frequency transfer function between the source signal and the observed sensor, the separation matrix and the source signal waveform are roughly calculated. FIG. 4 summarizes an example of a method for roughly estimating the separation matrix and the source signal waveform, which corresponds to the matter described in claim 2 of this patent. Using the roughly calculated separation matrix and source signal waveform, the separation matrix and source signal waveform are estimated in detail. A detailed estimation method for the separation matrix and the like corresponds to the matters described in claims 3 and 4. First, a detailed estimation method for the separation matrix and the like will be described below.

This paper summarizes the method of estimating the source signal in detail when the direction of arrival of the source signal is known in advance. To the objective function for estimating the source signal waveform from the mixed signal, a value obtained by adding prior information of the source signal arrival direction as a constraint condition is newly set as the objective function, and this objective function is minimized by: By updating the separation matrix, highly accurate signal separation can be achieved. The objective function for estimating the source signal waveform from the mixed signal includes, for example, the KL information amount of the product of the joint distribution and marginal distribution of the output of the separation operation. Claim 3 summarizes a method of setting a constraint condition in which the direction of arrival of the source signal is used as a priori information. Equation (9) can be used to calculate the frequency transfer function from each sensor to the source signal when the direction of arrival of the source signal and the location of the sensor placement are known. A matrix in which the frequency transfer functions are arranged such that the sensor number is a row vector and the source signal number is a column vector is called a mixing matrix, and is defined by Equation (10). A feature of the constraint condition according to claim 3 is that the information that the mixing matrix is the inverse matrix of the separation matrix is set as the constraint condition. When the mixing matrix is the inverse of the separation matrix, the product of the separation matrix and the mixing matrix is the identity matrix. This means that the inner product of the row vector of the separating matrix and the column vector of the mixing matrix having the same source signal number is 1, and the inner product of the row vector of the separating matrix and the column vector of the mixing matrix having different source signal numbers is 0. means to become Therefore, the constraint condition based on the direction of arrival of the source signal shown in Equation (41) imposes a constraint on the value of the inner product. Claim 3 describes the constraint on the inner product of the row vector of the separation matrix and the column vector of the mixing matrix. In order to keep the grating lobes out of the visible region, the sensor is placed at less than half the minimum wavelength. When this condition is met, the mixing matrix tends to be ill-conditioned at lower frequencies. Become. Therefore, the strength of the constraint based on the information that the inner product of the row vector of the separating matrix and the column vector of the mixing matrix with different source signal numbers is zero is adjusted. That is, the coefficients are adjusted for each frequency bin such that the more ill-conditioned the mixing matrix is, the weaker the strength of the constraint condition based on the direction of arrival of the source signal is. This adjustment of the coefficient emphasizes the constraint condition based on the source signal arrival direction that the product of the separation matrix and the mixing matrix becomes a unit matrix on the high frequency side where the reliability of the mixing matrix is high. Since the reliability of this mixing matrix varies depending on the sensor arrangement and the analysis symmetric frequency band of the signal, it cannot be generally concluded that the reliability drops at low frequencies. Therefore, it is desirable to evaluate the reliability of the mixing matrix and change the coefficient setting method of Equation (42) depending on the situation.

In order to use the constraint of claim 3, it is necessary to roughly estimate the magnitude of the variance of the source signal. The method of estimating the variance is summarized in Equations (11) to (36). To estimate the variance for each frequency bin and time frame, the Wiener mask is used to estimate the statistical variance of the source signal after roughly obtaining the separated signal waveforms. This method of roughly calculating the separated signal waveform corresponds to claim 2, and the specific calculation procedure is as follows. When the direction of arrival of the source signal is known, a mixing matrix in which the frequency transfer functions are arranged in a column vector can be obtained from equation (10). Under the condition that the mixing matrix and the observed signal are known, the logarithmic likelihood of the posterior probability of separated signals is represented by Equation (28). The separated signal waveform that maximizes this posterior probability is calculated by multiplying the observed signal by the inverse matrix of the mixing matrix plus the regularization term, as shown in Equation (29). In claim 2, the dispersion of separated signals is defined for each frequency _bin f and the update formula is obtained. It is also possible to change the notation to σ _n (τ), which is a constant between σ n (τ), and such a case is also covered by the scope of the present patent. The separation signal calculation method described in claim 2 can also calculate the separation matrix by itself and output the calculation result in the frequency domain or the time domain, and such cases are also covered by the scope of the present patent.

When the separation matrix is calculated using the direction of arrival of the source signal as a priori information, how to obtain the initial value is important. When the direction of arrival of the source signal is known, the frequency transfer function can be calculated using equation (9) from the positional relationship between the sensor and the source signal. In Equation (9), the phase delay at the reference sensor is set to 0 and the reference sensor is used as the viewpoint, but the viewpoint may be set at any position. By arranging this frequency transfer function so that the sensor number is a row vector and the source signal number is a column vector, the mixing matrix shown in Equation (10) is obtained. By using this mixing matrix as the initial value of the separation matrix and performing a solution search, the optimal solution can be reached efficiently. In addition, when the mixing matrix is ill-conditioned, such as when the source signals are close to each other, adding a regularization term to the objective function and then calculating the Moore-Penrose pseudo-inverse matrix of the mixing matrix stabilizes the A solution for the separated signals is obtained. This method is the content of claim 2, and can be realized by using equation (32). The separation matrix is repeatedly updated in the gradient descent direction of the objective function obtained by partially differentiating the objective function with the separation matrix to search for the optimum solution. In this way, the inverse of the mixing matrix or the Moore-Penrose pseudo-inverse of the mixing matrix obtained by adding the regularization term to the objective function is used as the initial value of the separation matrix, and the separation matrix is repeatedly updated by the gradient method. It is described in

When the separation matrix is calculated using the direction of arrival of the source signal as a priori information, how to obtain the initial value is important. When the direction of arrival of the source signal is known, the frequency transfer function can be calculated using equation (9) from the positional relationship between the sensor and the source signal. In Equation (9), the phase delay at the reference sensor is set to 0 and the reference sensor is used as the viewpoint, but the viewpoint may be set at any position. By arranging this frequency transfer function so that the sensor number is a row vector and the source signal number is a column vector, the mixing matrix shown in Equation (10) is obtained. By using this mixing matrix as the initial value of the separation matrix and performing a solution search, the optimal solution can be reached efficiently. In addition, when the mixing matrix is ill-conditioned, such as when the source signals are close to each other, a stable separation signal solution can be obtained by adding a regularization term to the square of the mixing matrix and then calculating the inverse matrix. is obtained. This method is the content of claim 2, and can be realized by using equation (30). A search for a solution is performed by repeatedly updating the separation matrix in the descending direction of the gradient obtained by partially differentiating the objective function with the separation matrix. In this way, the inverse matrix of the mixing matrix or the inverse matrix of the matrix obtained by adding the regularization term to the square of the mixing matrix is used as the initial value of the separation matrix, and the separation matrix is repeatedly updated by the gradient method. is doing.

In this way, the solution space of the objective function to which the constraint is added tends to be complicated, and the update of the separation matrix tends to be unstable. For this reason, we propose a method of stably updating the separation matrix when the solution space becomes complicated. First, the separation matrix is divided into row vectors or arbitrary ranges, and each range is defined as a matrix block. In normal independent component analysis, the update amount in the gradient descent direction of the objective function of the separation matrix is determined in advance by setting a constant parameter called learning rate. On the other hand, in the proposed method, the learning rate is updated for each matrix block in order to adjust the update amount of the separation matrix in the gradient descent direction. Therefore, using equations (47-1) and (47-4), for each matrix block, the cumulative value of the sum of squares of the gradients of the objective function up to the current time step and the update amount of the separation matrix are A cumulative value of the sum of squares is calculated respectively. Calculate the learning rate for each matrix block using the cumulative value of the sum of squares of the gradient of the objective function at the current time step and the cumulative value of the sum of squares of the update amount of the separation matrix at the previous time step. do. Formula (47-3 ) computes the update amount of the separation matrix at the current time step. In this way, the moving average is taken for the square of the gradient of the objective function of Equations (47-1) and (47-4) and the square of the separation matrix update amount, and based on these values, the gradient The optimum separation matrix is obtained by repeating a series of calculation flows for adjusting the update amount in the descending direction.

[Base vector calculator]
A separation matrix is used to separate the waveform data. Also, a basis vector is calculated from a separation matrix or a feature amount such as direction in a signal.

[Clustering part]
The basis vectors are clustered in the same direction or for each signal feature amount, and the signal number and the cluster number are associated with each other. Time-series data of the separated signals are obtained by collecting the separated signals for each time-frequency point for each cluster.

First, the concepts, terms, and variables used in cluster calculation are defined. Assume a mixed model in which the basis vector distribution is represented by a linear sum of C independent probability density functions. In this embodiment, the von Mises-Fischer distribution shown in Equation (55), which is a normal distribution on a sphere and has only directional data, is used as the probability density function of the mixture model. A group of unknown parameters possessed by C probability density functions is defined as C model parameters, and a group of model parameters included in the c-th distribution is called a c-th cluster. In a mixture model represented by a series sum of C probability density functions, the importance of the c-th probability density function is the mixture weight of the cluster π _c , and the variation around the average direction μ _c of the c-th cluster of all basis vectors is Define the degree of directional concentration _κc . Also, the probability that the basis vector of frequency bin f and source signal n belongs to the c-th cluster is defined as the membership degree u _cn (f). In the von Mises-Fischer distribution mixture model, three types of mixture weight π _c , directional concentration κ _c , and mean direction μ _c are defined as model parameters. Also, the probability of belonging to each cluster is expressed using the degree of belonging to the cluster of basis vectors u _cn (f). This degree of membership is defined as an auxiliary function. The objective function is designed by adding the logarithmized value of the probability density function to the value obtained by multiplying the constraint on the model parameters and the auxiliary function by the Lagrangian multiplier. The optimal solution is a combination of model parameters and auxiliary functions that maximizes this objective function. In this embodiment, instead of directly maximizing the objective function, we fix the assumption that one of the auxiliary function and the model parameters is known, maximize the objective function with respect to the other, and alternately update the solution. By iterating, the solution that maximizes the objective function is searched for.

5 and 6 are calculation flows of basis vector clustering. The cluster calculation procedure of this embodiment consists of an initial clustering section shown in FIG. 5 and a cluster determination section shown in FIG. The initial clustering unit treats the auxiliary function (degree of membership) as a continuous value, and alternately updates the auxiliary function and the model parameters to calculate the solution that maximizes the objective function. A specific calculation process will be described. Update the solution and search for the optimal solution by alternately repeating the steps of fixing one of the auxiliary function and model parameters as known, partially differentiating the objective function with respect to the other, and updating the solution. . First, the model parameters are fixed under the condition that they are known, and the objective function is maximized with respect to the auxiliary function, thereby executing an operation to update the auxiliary function. Next, the auxiliary function is fixed under the condition that it is known, and an operation is performed to update the model parameters. The optimum solution is searched for by alternately repeating these two kinds of calculations. The calculation for updating the model parameters includes the calculation for partially differentiating the objective function with respect to the mixture weight π _c and the average direction μ _c among the model parameters to search for a solution that makes the gradient zero, and the calculation for approximating the directional concentration κ _c . It consists of a formula or an operation calculated using a gradient method or a maximum likelihood estimation method. Equation (59) shows the calculation for updating the mixture weight _πc and the average direction _μc among the model parameters. Also, the procedure for calculating the directional concentration degree by the approximation formula and the gradient method is shown in formulas (60) to (63). Since the calculation method for obtaining the directional concentration degree is described in claim 8, it will be described later in detail. The degree of concentration of directions can be calculated using other approximation formulas or maximum likelihood estimation methods. is within the scope of claim 6. Expressions (58) to (63) for updating the model parameters and auxiliary variables represent the case where the basis vectors form a two-dimensional structure of frequency bins and source signals. This operation of updating the variable with the index of the two-dimensional structure of the frequency bins and the source signal is a pattern that has not existed in the past, and is a point where novelty can be claimed. Also, no attempt has been made to collect basis vectors obtained from signal feature amounts for each same source signal by a cluster analysis method based on the von Mises-Fischer distribution. In addition, each basis vector is multiplied by a weight α _n (f) to extend the method so that it can handle cases in which the importance of the basis vectors differs from each other. When cluster analysis is performed by obtaining base vectors from signal feature amounts and the like in units of time-frequency points, in the update formulas of formulas (58) to (63), the index of the frequency bin should be extended to the time frame/frequency bin. can be handled by Therefore, claim 6 of this patent mainly assumes two patterns: a case where basis vectors form a two-dimensional structure of frequency bins and source signals, and a case where basis vectors form a three-dimensional structure of time frames, frequency bins and source signals. ing.

Among the optimum solutions in the initial clustering section obtained by repeating these calculations, the model parameters are set to the initial values of the cluster determination section. In the cluster determination unit, the initial values of the model parameters are obtained, and while providing a constraint condition based on the foresight information that basis vectors with different source signal indices belonging to the same frequency bin always belong to different clusters, the objective function is determined to be the maximum. search for a solution. The cluster determination unit limits the auxiliary function (degree of membership) to discrete values of 0 or 1 when searching for a solution that maximizes the objective function. That is, instead of using equation (58) in calculating the degree of membership, the value of the degree of membership is limited to a discrete value of 0 or 1, and while setting constraints based on foresight information, the membership that maximizes the objective function Explore degree combinations. In the combination search of the degrees of membership, a brute-force search may be performed, but using the concept of set partitioning shown in Algorithm 1 makes the solution search more efficient. In this way, with the degree of membership discretized to 0 or 1, providing a constraint condition based on a priori information that basis vectors with different source signal indices belonging to the same frequency bin always belong to different clusters, An operation for alternately updating model parameters and auxiliary functions corresponds to claim 7 of this patent.

In basis vector clustering, the degree of directional concentration is also included in the modified Bessel function of the first kind, so it cannot be obtained analytically. Therefore, in claim 8, a method of calculating the directional concentration degree using the approximation formula and the gradient method according to the following procedure is shown. First, the ratio of modified Bessel functions of the first kind with different ranks by one is roughly calculated using the basis vectors and the degree of membership according to the equation (58). To this approximate value, the approximation formula of formula (61) or formula (63) obtained by expanding the ratio of the modified Bessel function of the first kind different in rank by one to continued fractions is applied, and the direction concentration degree Calculate an approximation. After that, the approximate values of the directional concentration degrees are substituted into the modified Bessel functions of the first kind having different ranks by one, and the detailed values of the ratios of the modified Bessel functions of the first kind different in rank by one are calculated. After calculating the value obtained by differentiating the detailed value of the ratio of the modified Bessel function of the first kind by the approximate value of the degree of concentration of direction by the formula (61) or (63), the detailed value and the rough ratio of the modified Bessel function of the first kind With the goal of minimizing the error in the values, we update the directional convergence index according to equation (62) applying Newton's method. The updated directional convergence degrees are again substituted into the modified Bessel functions of the first kind that differ in rank by one, and the detailed values of the ratios of the modified Bessel functions of the first kind that differ in rank by one are updated. The gradient is the value obtained by differentiating the updated detailed value of the ratio of the modified Bessel function of the first kind with a different rank by one with the updated degree of directional concentration, and the error between the detailed value and the approximate value of the ratio of the modified Bessel function of the first kind is minimized. With the goal of achieving uniformity, the degree of directional concentration is updated by Equation (62) obtained by applying Newton's method. By repeating the update calculation of the directional concentration degree by the Newton method several times, the directional concentration degree can be estimated with high accuracy. The directional concentration degree calculation procedure according to the calculation procedure described above corresponds to the content of claim 8 of this patent.

In this way, the search for solutions in the initial clustering section and the cluster determination section is a calculation procedure in which one of the auxiliary function and the model parameter is fixed on the assumption that one is known, and the step of searching for the optimum solution for the other is alternately repeated. done by In the calculation of the initial clustering part, when the number of auxiliary functions and model parameters increases due to the number of source signals, etc., the cluster determination part should prevent the solution from erroneously converging to a local solution. has the effect of being able to However, when the number of source signals is small, there is no problem even if the initial values of the model parameters are set by random numbers, and in that case, the calculation of the cluster determination section can be used alone. In these optimum solution search processes, the degree of concentration of directions is included in the modified Bessel function of the first kind, and the calculation requires some ingenuity. In the initial clustering unit, it is desirable to use the approximation formula (63) in which the directional concentration degree is uniform for all clusters. should be used.

[Source signal direction-of-arrival calculation unit]
By taking a weighted average for all frequency bins for each cluster with respect to the direction of arrival of the signal calculated from each basis vector using the correspondence between the signal number and the cluster number that were previously associated by clustering, Obtain the direction of arrival of the source signal at all frequency points.

FIG. 7 is a calculation flow of source signal direction-of-arrival calculation. The inverse matrix of the separation matrix is obtained, and the frequency transfer function for each frequency bin and source signal is obtained from the column vector using Equation (8). From the estimated value of the frequency transfer function including the information between the source signal and the observation point and the information on the arrangement position of the sensor, using the equation (74) derived based on the information on the phase difference between the sensors, the frequency bin / source The direction of arrival of the source signal for each signal is calculated. Also, a continuous value mask is calculated by equation (71) so as to reduce the influence of signals with low signal separation accuracy in each frequency bin on the signal calculation result and source signal arrival direction estimation result. By obtaining the weighted average of the direction of arrival of the source signal for each frequency bin and source signal masked by Equation (75) from the result of associating the signal number and the cluster number, the final direction of arrival of the source signal in all frequency bins can be obtained. decide. The base weights used in weighted averaging can be indexed by the power spectrum of observed values for each frequency as shown in Equation (76), but can also be all 1 if all frequencies have the same reliability. Alternatively, the sum of the eigenvalues of the signal subspace can be used as an index for the base weight, and in this case, the sum of the power of the incoming waves that directly reach the observation point without being approximately reflected is taken. become.

[Separation signal output part]
The signal information separated for each source signal is received from the signal waveform separator and displayed on the monitor in the form of time-series data of waveforms or power spectra in the frequency domain. The displayed waveform data is erased or saved according to the user's selection. In addition, by receiving the estimation result of the direction of arrival of the source signal and superimposing it on the image captured by the camera in advance, the direction of arrival of the source signal can be visualized in an easy-to-understand manner. A general-purpose computer or dedicated hardware can be used for the device that displays and stores the separated waveform data and the device that displays the estimation result of the direction of arrival of the source signal. In the former case, general-purpose computer hardware and software are configured in cooperation.

A specific calculation procedure for the direction-of-arrival source signal separation calculation in this embodiment will be described below. As described above, the computation is divided into the separation signal calculation section, the basis vector calculation section, the clustering section, and the source signal arrival direction calculation section, and each calculation process will be described.

[Separation signal calculator]
Complex number components in the frequency domain of the source signal s _n (t) in the time domain, the separated signal y _n (t), and the observed signal x _m (t) are represented by the following equations using vector representation.

ｆ＝１，…，Ｆ、τ＝１，…，Ｔは、それぞれ周波数ビン及び時間フレームのインデックスを表す。又、ｎ＝１，…，Ｎで全周波数帯域における帯域ごとの信号数の最大値、ｍ＝１，…，Ｍで観測センサの番号を表す。

f=1, . . . , F, τ=1, . Also, n=1, .

A windowed Fourier transform is applied to the observation signal x _m (τ) in the time domain while gradually moving the position of the window to generate observation data x(f, τ) for each time-frequency point. Separating x(f, τ) at time-frequency points by the separation matrix yields the separated signal y(f, τ), which is an estimate of the source signal s(f, τ).

Various methods such as independent component analysis and independent low-rank matrix analysis can be used to calculate the separation matrix. Independent component analysis is suitable when the prior information of the waveform is not known, for example, when dealing with abnormal noise in a building or when obtaining noise waveforms in different directions in a city. On the other hand, independent low-rank matrix analysis is suitable for separating sound signals with low-rank structures, such as music signals and human voices, by direction. Only the method of updating the separating matrix by independent component analysis is shown when prior information about the direction of arrival of the source signal is not obtained.

周波数ビンごとにＮ×Ｍ次元の行列構造を有する分離行列Ｗ（ｆ）を生成する。

A separation matrix W(f) having a matrix structure of N×M dimensions is generated for each frequency bin.

A separation matrix by independent component analysis is obtained using an optimization method based on statistical information of the target signal. Specifically, the separation matrix that minimizes the KL information amount of the product of the joint distribution and the marginal distribution of the output of the separation operation given by Equation (3) is searched for in the steepest descent direction using the gradient method. Calculated by

Ｅ［ｘ（ｆ，τ）］は、観測信号の時間フレームτ方向の平均値を算定する処理である。Ｅ［ｘ（ｆ，τ）ｙ^Ｈ（ｆ，τ）］は、観測信号と分離信号の時間フレーム方向の共分散行列を算定する処理を表す。ｐ_ｎ（ｙ_ｎ（ｆ，τ））は、目的信号の統計的な情報を表す確率密度関数である。確率密度関数をｙ_ｎ（ｆ，τ）で微分し、確率密度関数のＮ次元勾配ベクトルであるスコア関数φ（ｆ，τ）を算定する。式（４）に示すＫＬ情報量をＷ^Ｈ（ｆ）で偏微分することにより得られる式の勾配ベクトルを用いて、分離行列を最適化する。

E[x(f, τ)] is a process of calculating the average value of the observation signal in the time frame τ direction. E[x(f, τ)y ^H (f, τ)] represents the process of calculating the covariance matrix of the observed signal and the separated signal in the time frame direction. p _n (y _n (f, τ)) is a probability density function representing statistical information of the target signal. The probability density function is differentiated with respect to y _n (f, τ) to calculate the score function φ(f, τ), which is the N-dimensional gradient vector of the probability density function. The separation matrix is optimized using the gradient vector of the equation obtained by partially differentiating the KL information shown in equation (4) with W ^H (f).

ηは、分離行列の学習速度を制御する定数であり以降、学習率と称する。Ｉ_ＭはＭ次元単位行列である。

η is a constant that controls the learning speed of the separation matrix and is hereinafter referred to as the learning rate. I _M is the M-dimensional identity matrix.

スコア関数には、様々な種類があり、分離信号の確率密度関数がガウス分布に従うと仮定する場合は、式（５）となる。σ_ｎ（ｆ）は、分離信号の確率分布の統計的分散である。

There are various types of score functions, and when it is assumed that the probability density function of the separated signals follows a Gaussian distribution, Equation (5) is obtained. σ _n (f) is the statistical variance of the split signal probability distribution.

分離信号の確率密度関数が双曲線余弦に従うと仮定する場合は、スコア関数は、式（６）となる。この確率密度を用いると、優ガウス分布に従う源信号の高次モーメントを再現できる。

Assuming that the probability density function of the separated signals follows a hyperbolic cosine, the score function becomes Equation (6). Using this probability density, higher-order moments of the super-Gaussian source signal can be reproduced.

分離信号の確率密度関数がラプラス分布に従うと仮定する場合は、スコア関数は、式（７）となる。

Assuming that the probability density function of separated signals follows the Laplace distribution, the score function is given by Equation (7).

There is a music method as a method specialized in obtaining the direction of arrival of the source signal. If the conditions for applying the music method are met, the combination of source signal arrival directions in a specific frequency bin can be known. Also, there are times when the source signal position is physically known. In this way, we propose a method of improving signal separation accuracy by using prior information when the direction of arrival of the source signal is known at a specific frequency bin. By incorporating this technique as an option in the signal separation operation, it can be selectively used in the system.

A reference sensor number M' is selected from observation point numbers 1-M. A sensor having this sensor number is hereinafter referred to as a reference sensor. Assuming a far-field model as the signal transfer model, the difference in distance from the source signal between observation points is sufficiently small compared to the distance between each observation point and the source signal, so the arrival wave between observation points is The difference in amplitude attenuation of is very small and can be ignored. Considering this, the influence of attenuation from the source signal to the sensor is ignored, and the positional relationship between the source signal and the observation point and the wavelength λ(f) of the arriving wave at the frequency bin f are used as parameters in the equation ( 8) approximates the frequency transfer function.

ｐ_ｍはｍ番目のセンサの位置ベクトル、ｑ_ｎは基準センサから視た目的信号への方向を表す単位ベクトルである。式（８）において、基準センサの番号をＭ’とし、基準センサとその他センサとの周波数伝達関数の比をとることで、式（９）を得る。このように比をとることで、基準センサにおける波の位相を基準にして、その他のセンサにおける位相を算定していることになる。

pm is the position vector of the _m -th sensor, and _qn is a unit vector representing the direction from the reference sensor to the target signal. In equation (8), the number of the reference sensor is set to M', and the ratio of the frequency transfer functions of the reference sensor and the other sensors is obtained to obtain equation (9). By taking the ratio in this way, the phase of the wave at the reference sensor is used as a reference to calculate the phase at the other sensor.

上式においてセンサ間の位置ベクトルの差と源信号到来方向への単位ベクトルとの内積の値を波速で割った値がセンサ間の信号到達時間差（遅延差）となる。

In the above equation, the signal arrival time difference (delay difference) between sensors is obtained by dividing the value of the inner product of the position vector difference between the sensors and the unit vector in the source signal arrival direction by the wave velocity.

Consider a case in which the directions of arrival of N' source signals out of N source signals are known as a priori information. When the direction of arrival of the source signal is known, the frequency transfer function is calculated by equation (9) based on the delay difference between the sensors. Equation (10-2) defines an M-dimensional array manifold vector in which frequency transfer functions corresponding to source signals for which spatial prior information is obtained are arranged in the vertical direction with reference to sensor M'. Assuming that the mixing system is linearly time-invariant and can be represented by complex instantaneous mixing in the time-frequency domain, the array manifold vector is used to convert the M×N′-dimensional complex mixing matrix A(f) into Equation (10-1) can be defined by

ここで、到来方向が既知の源信号インデックスをｎ’＝１，…，Ｎ’としている。式（１０－２）で、ベクトル成分の分母に源信号から基準センサまでの周波数伝達関数があるが、この基準センサ番号は、分離信号を出力したいセンサ番号とすればよい。

Here, the index of the source signal whose direction of arrival is known is n'=1, . . . , N'. In equation (10-2), the denominator of the vector component is the frequency transfer function from the source signal to the reference sensor, and the reference sensor number may be the sensor number for which the separated signal is to be output.

If the frequency transfer function from the source signal to each sensor can be calculated without error, the separated signals can be calculated by using the inverse matrix of the mixing matrix created based on them as the separating matrix. However, when there is an observation error in the direction of arrival of the source signal, the mixing matrix consists of ill-conditioned equations, so even a small error has a large effect on the result. Also, since it is difficult to accurately measure the distance from the sensor to the position of the source signal in the far field, the frequency transfer function can be expressed as a relative value from the reference sensor by equation (9), ignoring the effect of signal amplitude attenuation. approximate. Therefore, if the inverse matrix of the mixing matrix is used, the separated waveforms will always include errors.

First of all, when N=N', that is, the direction of arrival of all source signals is known in advance, the statistical variance of the source signals is calculated and added to the direction of arrival information to calculate the separated signals. Summarize.

これまでと同様に、混合系の近似モデルを仮定すると、分離信号と観測信号の相関行列は、混合行列を用いて式（１１）で関連づけられる。観測信号の相関行列を固有値分解すると次の式が得られる。

As before, assuming a mixed system approximation model, the correlation matrices of the separated signal and the observed signal are related by Equation (11) using the mixing matrix. Eigenvalue decomposition of the correlation matrix of the observed signal yields the following equation.

▲Ｑ▼（ｆ）は、Ｍ次元固有ベクトルｑ_ｍ（ｆ）を並べたＭ×Ｍ次元の行列であり、Ｍ個の列ベクトルの内、Ｎ＜＝Ｍを満たす固有ベクトルが信号部分空間を張る。Ａ（ｆ）は、Ｍ個の固有値を値が大きい順に並べた対角行列である。式（１２）の固有値と固有ベクトルにおいて、源信号インデックスｎは周波数ビンごとに異なる源信号と対応していて一貫性がない。この交換不定性問題を解消するために、大きさが上位Ｎ個の固有値に対応する信号部分空間の固有ベクトルをその方向別にクラスタリングする。基底ベクトルを次の式で算定する。

▲Q▼(f) is an M×M-dimensional matrix in which M-dimensional eigenvectors q _m (f) are arranged, and among the M column vectors, eigenvectors satisfying N<=M span the signal subspace. A(f) is a diagonal matrix in which M eigenvalues are arranged in descending order. In the eigenvalues and eigenvectors of equation (12), the source signal index n is inconsistent, corresponding to different source signals for each frequency bin. To solve this exchange ambiguity problem, the eigenvectors of the signal subspace whose magnitudes correspond to the top N eigenvalues are clustered according to their directions. The basis vectors are calculated by the following formula.

Ｎ×Ｆ本の基底ベクトルｖ_ｎ（ｆ）をその方向が類視したペアごとにまとめたＮ×Ｆ個の大きさが同じクラスタに分類する。最も簡易なクラスタリング方法は、Ｋ－ｍｅａｎｓ法に基づく方法である。まず同じ源信号インデックスｎを有する基底ベクトルを一つのクラスタとみなして、各クラスタにおけるＭ次元平均方向ベクトルμ_ｃを算定する。基底ベクトルとその平均方向との差が最小すなわち内積が最大となるように、基底ベクトルを並び替え、クラスタを更新する。再度クラスタの平均方向を求める。クラスタに変化がなくなるまでこの交互に演算を繰り返す。演算過程を式（１４）に示す。

The N×F base vectors v _n (f) are classified into N×F clusters of the same size, which are groups of pairs whose directions are analogous. The simplest clustering method is based on the K-means method. First, the basis vectors having the same source signal index n are regarded as one cluster, and the M-dimensional average direction _vector μc in each cluster is calculated. The basis vectors are rearranged so that the difference between the basis vectors and their mean direction is minimized, ie, the inner product is maximized, and the clusters are updated. Find the average direction of the cluster again. This alternate operation is repeated until there is no change in clusters. The calculation process is shown in Equation (14).

ｃ＝１，…，Ｃはクラスタインデックスである。ｕ_ｃｎ（ｆ）は周波数ビンｆにおける源信号インデックスｎの基底ベクトルが第ｃクラスタに属する時に１、属さない時に０を返す帰属度変数である。基底ベクトルは必ずいずれかのクラスタに所属するため、帰属度変数のクラスタインデックスｃ方向への和は必ず１となる（式（１５－１））。又、帰属度変数の源信号インデックスｎ方向への和は、周波数ビンｆの中に源信号ｎが存在する時に１、存在しない時に０となる。今回のように信号を同数のクラスタに分類する並び替えの問題、すなわちＣ＝Ｎが成り立つ場合は、クラスタ番号と源信号番号が１対１で対応するため、源信号インデックス▲ｎ▼方向への和が１となる（式（１５－２））

c=1, . . . , C is the cluster index. u _cn (f) is a membership variable that returns 1 if the basis vector of source signal index n in frequency bin f belongs to the c-th cluster and 0 if it does not. Since a basis vector always belongs to one of the clusters, the sum of membership degree variables in the cluster index c direction is always 1 (equation (15-1)). Also, the sum of the degree-of-membership variables in the direction of the source signal index n is 1 when the source signal n exists in the frequency bin f, and 0 when it does not exist. In the rearrangement problem of classifying the signals into the same number of clusters as in this case, that is, when C=N holds, the cluster number and the source signal number correspond one-to-one. sum becomes 1 (formula (15-2))

しかし、Ｋ－ｍｅａｎｓ法によるクラスタリングでは、基底ベクトルの本数が多い時、最終的に生成されるクラスタが初期のクラスタに依存してしまい、一意に解が定まらないという問題がある。そこで、後程「クラスタリング部」に示す混合モデルとしてフォンミーゼス・フィッシャー分布を仮定した演算によりクラスタ分析を行い、帰属度変数の集合ｕ_ｃｎ（ｆ）を決定することを推奨する。クラスタリングが完了すると、その結果を用いて固有ベクトルをクラスタごとに算定し、各クラスタに

However, clustering by the K-means method has a problem that when the number of base vectors is large, the finally generated cluster depends on the initial cluster, and a unique solution cannot be determined. Therefore, it is recommended to perform cluster analysis by calculations assuming the von Mises-Fischer distribution as the mixture model described later in the “Clustering part” and determine the set u _cn (f) of membership variables. After clustering is complete, the results are used to compute the eigenvectors for each cluster, and for each cluster

式（１６）に示す平均方向ベクトルの成分を用いて基準センサＭ’から視たクラスタごとの源信号到来方向を次式で算定する。

Using the component of the average direction vector shown in Equation (16), the direction of arrival of the source signal for each cluster viewed from the reference sensor M' is calculated by the following equation.

クラスタごとの源信号到来方向が実際の到来方向ＤＯＡ_ｎと一致するように、全周波数ビンのクラスタインデックスｃを同一規則で並び替える。まず、式（１８）によりクラスタごとの源信号到来方向が実際の源信号到来方向と近い値となるように対応付けを行う。この対応付けによりクラスタ番号ｃを入力すると正しい源信号番号▲ｎ▼を返す順列集合｛Π（ｃ）｝を作成する。

The cluster indices c of all frequency bins are rearranged according to the same rule so that the direction of arrival of the source signal for each cluster matches the actual direction of arrival DOA _n . First, according to Equation (18), the source signal arrival direction for each cluster is associated with a value close to the actual source signal arrival direction. Based on this association, a permutation set {Π(c)} that returns a correct source signal number {circumflex over (n)} when a cluster number c is input is created.

この順列Π（ｃ）を用いて、固有値と固有ベクトルの並び替え（クラスタリング）を以下により行う。

Using this permutation Π(c), eigenvalues and eigenvectors are rearranged (clustered) as follows.

並び替えが完了したら、源信号の統計的分散を周波数伝達関数と固有値分解の結果を用いて求める。分離信号が独立であれば、その相関行列は、次の式に示すように源信号の統計的分散を対角に並べた行列と一致する。

Once permutation is complete, the statistical variance of the source signal is determined using the frequency transfer function and eigenvalue decomposition results. If the separated signals are independent, their correlation matrix will match the diagonal matrix of the statistical variances of the source signals as shown in the following equation.

式（１１）と式（１２）の形状を比較すると明らかなように、▲Ｑ▼（ｆ）は無相関化された分離信号の混合行列と対応する。そのため、式（２０）を式（１１）に代入し 式（１２）と係数比較すると、源信号の統計的分散は以下により算定できる。

As is clear from comparing the shapes of equations (11) and (12), Q(f) corresponds to the mixing matrix of decorrelated separated signals. Therefore, substituting equation (20) into equation (11) and comparing the coefficients with equation (12), the statistical variance of the source signal can be calculated by:

源信号到来方向の事前情報を用いて分離信号を概略的に算定する方法についてまとめる。以下に述べる統計的分散推定演算で用いられる分離信号算定方法は本特許の請求項２と対応する。

A method for roughly estimating separated signals using a priori information of the direction of arrival of the source signal is summarized. The separated signal calculation method used in the statistical variance estimation operation described below corresponds to claim 2 of this patent.

When the direction of arrival of the source signal is known, we attempt to simply calculate the separated signals using the direction of arrival as prior information. Consider the case where the source signals are mixed and received by the sensor as the observed signal. At this time, the mixture model in the frequency domain can be expressed as follows.

源信号の推定値である分離信号は、混合行列の逆行列を用いて算定できる。

Separated signals, which are estimates of the original signal, can be calculated using the inverse of the mixing matrix.

上付文字＋は、Ｍｏｏｒｅ－Ｐｅｎｒｏｓｅ疑似逆行列の表記であり、行列の行数と列数が等しい時に、通常の逆行列となる。解析対称周波数の中の低周波領域や源信号が近接した位置に存在する場合に、混合行列は、悪条件な方程式となる。そのため、式（２３）のように逆行列を用いる方法では、僅かな源信号到来方向の推定誤差が分離信号算定結果に大きな誤差をもたらす。この算定誤差をｅ（ｆ）と定義し、混合行列の悪条件性に対処するために混合モデルを等式制約とした次の最小化問題を解く。

The superscript + is the notation for the Moore-Penrose pseudo-inverse matrix, which is the normal matrix inverse when the number of rows and columns of the matrix is equal. The mixing matrix becomes an ill-conditioned equation when the low-frequency region among the analytic symmetric frequencies and the source signals exist at close positions. Therefore, in the method using an inverse matrix as in Equation (23), a slight error in estimating the direction of arrival of the source signal causes a large error in the separated signal calculation result. Defining this computational error as e(f), we solve the following minimization problem with equality constraints on the mixture model to deal with the ill-conditioned nature of the mixture matrix.

観測信号がガウス分布に従うと仮定すると、源信号の確率密度関数は、次式で表せる。

Assuming that the observed signal follows a Gaussian distribution, the probability density function of the source signal can be expressed by the following equation.

雑音が各周波数ビンで独立に分散σ_Ｅ（ｆ）のガウス分布に従って発生する時、混合行列と分離信号が既知という条件下での観測信号の確率密度関数は以下で与えられる。

When the noise occurs independently in each frequency bin according to a Gaussian distribution with variance σ _E (f), the probability density function of the observed signal under the condition that the mixing matrix and separated signals are known is given by:

ベイズの定理より、混合行列と観測信号が既知という条件下での分離信号の確率密度関数は、

From Bayes' theorem, the probability density function of the separated signals under the condition that the mixing matrix and the observed signal are known is

となり、Ｐｒ（Ａ（ｆ），ｘ（ｆ，τ））を一定とし、式（２７）に式（２２）、（２５）、（２６）を代入し、対数化すると、式（２８）が得られる。

Assuming that Pr(A(f), x(f, τ)) is constant, and substituting equations (22), (25), and (26) into equation (27) for logarithmization, equation (28) becomes can get.

上記の対数尤度を最大化する分離信号を求めるために、分離信号ｙ（ｆ，τ）で左辺を微分して極大点すなわち微分値が０となる分離信号を求めると、

In order to obtain the separated signal that maximizes the above logarithmic likelihood, if the left side is differentiated with the separated signal y(f, τ) to obtain the separated signal with the maximum point, that is, the differential value of 0,

Σ（ｆ）はＮ×Ｎ次元対角行列である。この式を解くには、源信号の分散と雑音の分散の２つを求める必要があるが、これらはいずれも未知である。そこで、雑音の分散を未知数として、式（２９）を微分した値を０と等値すると、

Σ(f) is an N×N dimensional diagonal matrix. To solve this equation, it is necessary to find two variables, the source signal variance and the noise variance, both of which are unknown. Therefore, if the noise variance is an unknown and the value obtained by differentiating equation (29) is equal to 0, then

式（２９）の形状に着目すると、雑音と源信号の分散の比率が正則化項の強さを決定することがわかる。そこで源信号の分散を用いて、雑音の分散の初期値を次式で与える。

Looking at the shape of equation (29), it can be seen that the ratio of the noise and source signal variances determines the strength of the regularization term. Therefore, using the variance of the source signal, the initial value of the variance of the noise is given by the following equation.

Ｋ_Ｆは観測ノイズの大きさを推定するパラメータであり、０．０１～０．１程度とするのがよい。但し、源信号の位置が互いに離れている場合や解析対称周波数が特定の帯域のみの場合など、混合行列が悪条件でない時は、０．０１より小さな値とすることもできる。ｐｓは、混合行列が低周波数になるほど悪条件となる性質を反映したパラメータで、０～１程度とするのがよい。計算は次の手順で行う。
まず源信号の分散として式（１２）に示す観測信号の相関行列の固有値から求めた値を用いる。源信号数が観測センサ数よりも多い過完備基底の条件では、固有値分解ができないため、源信号の分散構造に対する事前情報が得られない。そのような場合、源信号の分散は全て１とするとよい。雑音分散の初期値は、式（３０）により求める。次に、式（２９）と式（３０）を交互に計算して、分離信号の値と雑音分散の値を交互に更新する。この交互演算を全周波数ビンにおける雑音の分散の変化が小さくなるまで行う。この交互演算が収束した時の分離信号がその事後確率を最大化する最適解となる。

_KF is a parameter for estimating the magnitude of observation noise, and is preferably about 0.01 to 0.1. However, when the mixing matrix is not ill-conditioned, such as when the positions of the source signals are separated from each other or when the analysis symmetric frequency is only in a specific band, the value can be smaller than 0.01. ps is a parameter that reflects the property that the mixing matrix becomes worse as the frequency becomes lower, and is preferably set to approximately 0-1. Calculation is performed in the following procedure.
First, the value obtained from the eigenvalue of the correlation matrix of the observed signal shown in Equation (12) is used as the variance of the source signal. Under the condition of an overcomplete basis in which the number of source signals is greater than the number of observed sensors, eigenvalue decomposition cannot be performed, and prior information on the distributed structure of the source signals cannot be obtained. In such a case, the variance of the source signal should be all unity. The initial value of the noise variance is obtained by Equation (30). Equations (29) and (30) are then alternately calculated to alternately update the value of the separated signal and the value of the noise variance. This alternating operation is performed until the change in noise variance in all frequency bins becomes small. The separated signal when this alternating operation converges becomes the optimum solution that maximizes the posterior probability.

After the update is completed, amplitude correction is performed using the following equation, which corresponds to the inverse matrix of the separation matrix. The method of calculating separated signals by arithmetic operations described so far has a small amount of calculation and is fast, but the source signal waveform tends to have a sparse structure than it actually is. However, considering the fact that the characteristics of the power spectrum shape can be well captured and the simplicity of calculation, it is considered practical.

In addition to this, various methods are conceivable for roughly obtaining the separation matrix. The simplest method is to find the frequency transfer function between the reference sensor and the source signal so that the inter-sensor delay is compensated only for the source signal whose direction of arrival is known, and arrange the frequency transfer functions for each source signal. In this method, a matrix obtained by arranging the complex conjugate vectors of the vector h _n (f) arranged in column vectors is used as a separating matrix (equation (33-1)). This method should be used when the mixing matrix is ill-conditioned and the inverse matrix cannot be accurately calculated, or when the number of signal sources N′ whose arrival directions are known is smaller than the total number of signal sources N. A method of using the Moore-Penrose pseudo-inverse matrix of the mixing matrix as the separating matrix is also conceivable (equation (33-2)). This method should be used when the inverse matrix of the mixing matrix can be accurately calculated for high frequency bands. Claim 4 summarizes how to obtain the initial values of these separation matrices.

I _p is the p-dimensional identity matrix. When the source signal has a clear frequency structure such as an instrumental sound or voice, the statistical variance is updated by adding the prior information of the direction of arrival of the source signal as a regularization term, and the separation matrix and the separated signal are calculated more accurately. can do. The method is summarized below. First, separated signals are calculated using the separation matrix calculated by the rough calculation method shown so far. Assuming that the separated signals follow the signal source model of MNMF (Non-Patent Document 13), the power spectrum of the separated signals is preserved in the sense of the expected value, so the Wiener mask M _n ^(G) (f, τ) is can be designed by the formula

基準センサＭ’における観測信号のパワースペクトルをウィナーマスクによりセンサ内の源信号ごとのパワースペクトルに分解し、基準センサ内における源信号の統計的分散を次のいずれかの式で求める。

The power spectrum of the observed signal at the reference sensor M' is decomposed by the Wiener mask into the power spectrum of each source signal in the sensor, and the statistical variance of the source signal in the reference sensor is obtained by one of the following equations.

一般に上記のような分離信号のパワーの加法性は成り立たないが、観測信号の生成モデルが多変量複素ガウス分布に従う時、パワースペクトルの加法性が期待値の意味で保存される。上記に示すウィナーマスクによる分離信号の分散推定はこの性質を利用したものである。

In general, the above additivity of separated signal powers does not hold, but when the observed signal generation model follows a multivariate complex Gaussian distribution, the power spectrum additivity is preserved in the sense of the expected value. Variance estimation of separated signals by the Wiener mask shown above utilizes this property.

上記の算定式は、分離信号のパワーが源信号パワーと等しいという条件式である。源源信号到来方向が全て既知でない場合、ウィナーマスクを設計できないので、上記の算定式を用いるとよい。式（３５）又は式（３６）により算定される源信号の統計的分散と概略算定により求めた分離行列を初期値に設定する。これらを更新して源信号推定精度向上を実現する手法を示す。

The above calculation formula is a conditional formula that the power of the separated signal is equal to the power of the source signal. If the direction of arrival of the source signal is not known, the Wiener mask cannot be designed, so the above formula should be used. The statistical variance of the source signal calculated by Equation (35) or Equation (36) and the separation matrix obtained by rough calculation are set as initial values. A technique for updating these to improve the accuracy of source signal estimation is presented.

This paper summarizes a technique for highly accurate signal separation by adding a priori information on the direction of arrival of the source signal as a constraint to the objective function. MNMF uses an objective function based on a multichannel observed signal generation model (Non-Patent Documents 13 and 15) assuming that the source signal follows a multivariate complex Gaussian distribution. This objective function is equivalent to the posterior probability of the source signal under the condition that the Gaussian distributed source signal variance and the frequency transfer function between the sensor and the source signal are known. In the proposed method, the objective function of the MNMF is obtained by approximating the spatial correlation matrix with a rank-1 matrix, and the constraint condition LS (W(f)) regarding the direction of arrival of the source signal is added to the objective function L _MNMF . Define a function LL(W(f)). The objective function L _MNMF approximated by this rank-1 matrix is used in independent low-rank matrix analysis (Non-Patent Document 14).

ＭＮＭＦに基づく目的関数Ｌ_ＭＮＭＦに制約条件ＬＳ（Ｗ（ｆ））を加えた式（３７－１）の目的関数が最小となるように、分離行列を更新する。目的関数の勾配降下方向を求めるために、式（３７－１）をＷ^Ｈ（ｆ）で偏微分すると、次式が得られる。

Objective function L based on _MNMF The separation matrix is updated so that the objective function of equation (37-1) obtained by adding constraint condition LS(W(f)) to MNMF is minimized. In order to obtain the gradient descent direction of the objective function, the following equation is obtained by partially differentiating equation (37-1) with respect to W ^H (f).

右辺第１、２項がＭＮＭＦに基づく目的関数の勾配である。上式において右辺第１、２項のみＷ^Ｈ（ｆ）Ｗ（ｆ）を掛けて自然勾配に補正すると、最急降下方向は近似的に以下により与えられる。数学的には、式中の特定項のみを自然勾配に補正する演算は間違いなのだが、正則化項ＬＳ（Ｗ（ｆ））の次元が自然勾配補正後の次元と一致するように設定するため、実用上問題はない。

The first and second terms on the right side are the gradients of the objective function based on MNMF. In the above equation, if only the first and second terms on the right side are multiplied by W ^H (f) W (f) to correct for the natural gradient, the steepest descent direction is approximately given by the following. Mathematically, it is wrong to correct only a specific term in the formula to the natural gradient, but because the dimension of the regularization term LS (W (f)) is set to match the dimension after natural gradient correction , there is no practical problem.

Ｉ_Ｍは、Ｍ×Ｍ次元単位行列である。

I _M is an M×M dimensional identity matrix.

ＭＮＭＦに基づく目的関数を用いると、分離信号の確率密度関数の微分値であるスコア関数は式（３９）となる。この分布では、源信号の統計的分散が時間周波数点単位で定義されている。源信号到来方向の事前情報を基に求めた混合行列に、源信号が互いに独立であるという条件を加えて解を探索する。源信号到来方向事前情報の制約条件に基づく正則化項ＬＳ（Ｗ（ｆ））を次の式により表現し、ＭＮＭＦに基づく目的関数に加えることで、分離精度を向上させる。

Using the MNMF-based objective function, the score function, which is the differential value of the probability density function of the separated signals, is given by Equation (39). In this distribution, the statistical variance of the source signal is defined per time-frequency point. A solution is searched for by adding a condition that the source signals are independent from each other to the mixing matrix obtained based on the prior information of the direction of arrival of the source signals. A regularization term LS(W(f)) based on the constraint condition of the source signal arrival direction prior information is expressed by the following equation and added to the objective function based on the MNMF to improve the separation accuracy.

ｅ_ｎは▲ｎ▼番目の要素が１、他要素が０のＭ次元単位ベクトルである。正則化係数における上付括弧文字のｃｙｃｌｅ＝１，…，ＣＹＣＬＥは、繰り返し計算の反復回数を表すインデックスである。式（４１－１）の第１項は、同じ源信号番号の分離行列の行ベクトルと混合行列の列ベクトルの内積が１になるという制約条件であり、目的方向から到来する信号を強調する効果がある。第２項は、異なる源信号番号の分離行列の行ベクトルと混合行列の列ベクトルの内積が０になるという制約条件になっていて、目的方向以外から到来する信号を遮断する死角形成効果がある。正則化係数ρ_Ｄ ^{（ｃｙｃｌｅ）}とρ_Ｕ ^{（ｃｙｃｌｅ）}（ｆ）は、各項の制約の強さを支配する定数である。これら２種類の正則化係数は次の関係式で関連付けておくのがよい。

e _n is an M-dimensional unit vector in which the {n}th element is 1 and the other elements are 0s. The superscripted parenthesis characters cycle=1, . The first term of equation (41-1) is a constraint that the inner product of the row vector of the separation matrix and the column vector of the mixing matrix with the same source signal number is 1, and has the effect of emphasizing the signal arriving from the target direction. There is The second term is a constraint that the inner product of the row vector of the separation matrix and the column vector of the mixing matrix with different source signal numbers becomes 0, and has the effect of forming a blind spot that blocks signals coming from directions other than the target direction. . The regularization coefficients ρ _D ^(cycle) and ρ _U ^(cycle) (f) are constants that govern the strength of each term's constraint. These two types of regularization factors are preferably related by the following relational expression.

ｐｓは、０～１程度の定数である。ある特定の２つのセンサを用いて、源信号到来方向を推定したい時に、グレーティングローブが可視領域に入らないようにするためには、センサ間隔は、解析対称の最大周波数で信号伝搬速度（音の場合は音速）を割った値すなわち最小波長の半分以下になるようにすべきである。よって、Ｍ個のセンサで信号を観測し、基準センサＭ’から視た源信号到来方向を推定する場合、基準センサと基準センサから最も距離の近いセンサとの距離が最小波長の半分以下であれば、源信号到来方向が一意に定まることになる。この条件を満たす時、例えば、基準センサと基準センサから最も距離の近いセンサに着目すると、低周波数になるほど可視領域が狭くなる。一方、その他のセンサと基準センサとの間に着目すると、グレーティングローブが可視領域に入らない周波数の中の最大値よりも小さい値の周波数に対して同様の関係が成り立つ。そのため、グレーティングローブが可視領域に入らず源信号到来方向が一意に定まる時、混合行列Ａ（ｆ）は、傾向として低周波数になるほど悪条件になる。この傾向を基に低周波になるほど、目的関数の制約条件に基づく正則化項における死角形成効果を弱める役割を式（４２）の正則化係数ｐｓは担っている。式（４１）の正則化項を分離行列の複素共役で偏微分すると、次の式が得られる。

ps is a constant of about 0-1. When estimating the direction of arrival of a source signal using two specific sensors, in order to prevent grating lobes from entering the visible region, the sensor spacing should be set to the signal propagation velocity (sound It should be less than half the minimum wavelength divided by the speed of sound, if applicable. Therefore, when observing signals with M sensors and estimating the direction of arrival of the source signal viewed from the reference sensor M′, even if the distance between the reference sensor and the sensor closest to the reference sensor is less than half the minimum wavelength, Then, the direction of arrival of the source signal is uniquely determined. When this condition is satisfied, for example, focusing on the reference sensor and the sensor closest to the reference sensor, the lower the frequency, the narrower the visible region. On the other hand, focusing on the relationship between the other sensors and the reference sensor, a similar relationship holds for frequencies smaller than the maximum value among the frequencies at which the grating lobes do not fall within the visible region. Therefore, when the grating lobe does not enter the visible region and the direction of arrival of the source signal is uniquely determined, the mixing matrix A(f) tends to be ill-conditioned as the frequency becomes lower. Based on this trend, the regularization coefficient ps of Equation (42) plays a role of weakening the blind spot formation effect in the regularization term based on the constraints of the objective function as the frequency becomes lower. Partial differentiation of the regularization term in equation (41) with the complex conjugate of the separation matrix yields the following equation.

Ｏ_ｐｑはｐ×ｑ次元の零行列を表す。制約条件を分離行列で偏微分して得られるＮ×Ｍ次元の行列におけるＮ’×Ｍ次元のＮ’×Ｍ次元の非零要素Ｐ_Ｄ（ｆ）及びＰ_Ｕ（ｆ）は、Ｍ次元アレイマニフォールドベクトルａ_ｎ’（ｆ）を用いて、式（４４）で算定できる。

O _pq represents a zero matrix of dimension p×q. Non-zero elements P _D (f) and P _U (f) of N′×M dimensions in the N′×M dimensional matrix obtained by partially differentiating the constraint with the separation matrix are M-dimensional arrays It can be calculated by Equation (44) using the manifold vector a _n' (f).

Ｒ_ｎ（ｆ）は、空間相関行列をアレイマニフォールドベクトルでランク１近似したＭ×Ｍ次元エルミート行列である。このように到来方向が既知の源信号に対してのみ方向情報に基づく制約条件を加えることで、複数の源信号が混合した波形から到来方向が既知の信号波形を正確に抽出することが可能になる。

R _n (f) is an M×M-dimensional Hermitian matrix obtained by rank-1 approximating the spatial correlation matrix with the array manifold vector. In this way, by adding constraint conditions based on direction information only to source signals with known directions of arrival, it is possible to accurately extract signal waveforms with known directions of arrival from waveforms in which multiple source signals are mixed. Become.

Adding the regularization term of Equation (44) to the objective function complicates the solution space, and updating the separation matrix based on the steepest descent method by the natural gradient method may result in local solutions or unstable calculation. Diverge. Therefore, a new separation matrix update method according to claim 5 is proposed. By adding the previous update direction of the separation matrix as an inertia term to the direction of gradient descent and adjusting the amount of movement in the direction of gradient descent, abrupt changes in the direction of update are suppressed and the update of the separation matrix is stabilized. try. The separation matrix is updated for each matrix block such as row vector. The gradient descent direction is obtained from the direction of the value obtained by partially differentiating the objective function shown in Equation (37) with the complex conjugate of the separation matrix, and the separation matrix is updated. Decompose the matrix into several matrix blocks and update the solution for each matrix block. In the following, consider the case where the matrix blocks are partitioned by the same source signal number, that is, the case where the row vectors of the separation matrix are defined as matrix blocks.

In order to prevent erroneous convergence to a local solution, for each matrix block, the gradient of the objective function is multiplied by the update amount of the separation matrix at the previous time step and added as an inertia term. Due to this inertia term, for example, even in a case where there is a deep valley in the solution search area, it is possible to easily escape from the valley, making it difficult to fall into a local optimum. Furthermore, in order to stabilize the calculation, the learning rate is adjusted by taking the squared moving average of the objective function gradient and the past separation matrix update amount. Specifically, when the accumulated value of the sum of squares of the gradients of the objective function is large, the update amount is small, and when it is small, the update amount is large. When the accumulated sum of squares of is large, the update amount is large, and when it is small, the update amount is small. Adjusting the update amount of the separation matrix in this way is equivalent to automatically adjusting the learning rate in the update formula of the natural gradient method used in independent component analysis according to the shape of the solution space. lead to transformation. Specifically, the separation matrix is updated by the following formula.

分離行列の第ｎ行を式（４５）で定義する。ｅ_ｎはｎ番目の要素が１、他要素が０のＭ次元単位ベクトルである。

The n-th row of the separation matrix is defined by Equation (45). e _n is an M-dimensional unit vector in which the n-th element is 1 and the other elements are 0s.

周波数ビンｆにおける式（３９）の目的関数勾配行列の第ｎ行を抜き出したＮ次元ベクトル▽Ｌ_ｎ（ｆ）を式（４６）で定義する。これらの定義式では、行列ブロックを分離行列の行ごとに定義しているが、その他の行列ブロック形状を取ることも可能である。一番単純なのは、分離行列全体を行列ブロックとみなし、分割数を１にする方法である。これらのケースでも請求項５に記載の事項と合致する場合は、本特許の権利範囲とみなす。周波数ビンｆ、源信号ｎごとに次の更新式を用いて分離行列を更新する。

An N-dimensional vector ▽L _n (f) extracted from the n-th row of the objective function gradient matrix of Equation (39) at frequency bin f is defined by Equation (46). Although these definitions define a matrix block for each row of the separation matrix, other matrix block shapes are possible. The simplest method is to treat the entire separation matrix as a matrix block and set the number of divisions to one. Even in these cases, if the matters described in claim 5 are met, they are considered to be within the scope of this patent. The separation matrix is updated using the following update formula for each frequency bin f and source signal n.

εは、１０^－６程度の微小な定数、γは、０．９程度の定数である。μ^{（ｃｙｃｌｅ）}は、慣性項の強さを支配する定数で、交互反復計算回数ｃｙｃｌｅに応じて値を調整する。η_Ｍａｘは学習率の上限値であり、０．１程度の値とするのがよい。上付き括弧文字ｌ＝１，…，Ｌは、分離行列更新計算における繰り返しインデックスである。変数Λ_ｎ ^（ｌ）（ｆ）とΓ_ｎ ^（ｌ）（ｆ）は、それぞれタイムステップｌにおける分離行列更新量の二乗和と損失関数（目的関数の勾配）の二乗和の累積値であり、いずれも実変数である。式（４７－１）と式（４７－４）で、目的関数の勾配の二乗と分離行列更新量の二乗について移動平均をとる。式（４７－２）で周波数ビン・源信号ごとに、学習率η_ｎ（ｆ）を目的関数の勾配と１つ前のタイムステップまでの分離行列の更新量に応じて自動調整する。式（４７－２）では学習率の上限値を０．１としているが、適用ケースに応じてこの上限値は任意に定めることができる。式（４７－３）では、右辺第１項の慣性項で局所解からの脱却を図ることで、目的関数が複雑化した時の解探索を可能にしている。さらに右辺第２項で移動平均により学習率を調整することで、目的関数降下方向への分離行列更新量が急激に変化するのを抑えている。分離行列更新後は式（４７－６）で分離信号を算定する。

ε is a minute constant of about 10 ⁻⁶ and γ is a constant of about 0.9. μ ^(cycle) is a constant that governs the strength of the inertia term, and its value is adjusted according to the number of alternate iterative calculation cycles. η _Max is the upper limit of the learning rate, and is preferably set to a value of about 0.1. The superscripted parenthesis letter l=1,...,L is the iteration index in the separation matrix update computation. The variables Λ _n ^(l) (f) and Γ _n ^(l) (f) are cumulative values of the sum of squares of the separation matrix update amount and the sum of squares of the loss function (gradient of the objective function) at time step l, respectively, Both are real variables. In equations (47-1) and (47-4), a moving average is taken for the square of the gradient of the objective function and the square of the separation matrix update amount. In Equation (47-2), the learning rate η _n (f) is automatically adjusted according to the gradient of the objective function and the update amount of the separation matrix up to the previous time step for each frequency bin and source signal. Although the upper limit of the learning rate is set to 0.1 in equation (47-2), this upper limit can be arbitrarily determined according to the application case. In equation (47-3), the inertia term of the first term on the right side is used to break away from the local solution, thereby enabling solution search when the objective function becomes complicated. Furthermore, by adjusting the learning rate by moving average in the second term on the right side, a rapid change in the separation matrix update amount in the descending direction of the objective function is suppressed. After the separation matrix is updated, the separation signal is calculated by Equation (47-6).

After the separation matrix update operation has converged, the statistical variance of the source signal is then updated. The frequency spectrum of the source signal's statistical variance σ _n (f, τ) is approximated by t _b (f) and its time-varying excitation pattern is reconstructed with γ _b (τ). Assuming that the ranks of the spectra between the source signals differ greatly, we divide the B bases into the N source signals using a partitioning function β _nb to obtain a low-rank approximation of the source signal variance. The variance is updated by non-negative matrix factorization (ISNMF) (Non-Patent Document 16) using Itakura-Saito pseudorange as an error evaluation index of the statistical variance of the separated signal and the source signal. It is performed by an update formula (Non-Patent Document 17).

ｂ＝１，…，Ｂは、低ランク近似で用いる基底のインデックスである。この更新式は、源信号の統計的分散が低ランク分解できるという拘束条件を満たしつつ、式（３７）の目的関数を最小化する分散を求める最尤推定と等価である。

b=1, . . . , B are the indices of the basis used in the low-rank approximation. This update formula is equivalent to maximum likelihood estimation for finding the variance that minimizes the objective function of Equation (37) while satisfying the constraint that the statistical variance of the source signal can be decomposed into low ranks.

Using the gradient of the regularization term based on the source signal spatial information obtained by Equation (41) and the initial value of the source signal variance obtained by Equation (35) or (36), the gradient of the objective function is calculated by Equation (39). Ask. Using the separation matrix obtained from the source signal space information by the equation (33) as an initial value, the separation matrix and the separated signals are updated by the equation (47). Using the updated separation matrix and separation signal, the calculation of the objective function gradient by Equations (39) and (41) and the update of the separation matrix by Equation (47) are repeated. When the iteration index l reaches L times, the update of the separation matrix is finished, and the update of the statistical variance of the source signal is started. Update the statistical variance of the source signal using equation (48). The statistical distribution is updated repeatedly about three times at the initial update (cycle=1), and updated once thereafter (cycle≧2). Signal separation is performed by alternately repeating a two-step update of the separation matrix and the update of the statistical variance. The number of iterations of alternating computation is determined in advance, and the value of the regularization term is decreased as the iteration progresses.

例えば式（４９）のように、交互反復回数ｃｙｃｌｅが予め設定した反復回数ＣＹＣＬＥに０．５～０．８程度の定数Ｃｒを掛けた回数を超えると正則化項と慣性項がゼロになるように、各項の係数を単調減少させる。このように正則化項と慣性項を繰り返し計算が進むほど、小さくしていくことで、源信号の空間情報に基づいて、統計的分散を学習させることができる。正則化項の初期値ρ_Ｄ ^（ｌ）は、１０～３０程度の定数であり、源信号が互いに近接して存在する場合などの混合行列が悪条件の時には小さな値、源信号到来方向の信頼度が高い時に大きな値とするとよい。

For example, as shown in equation (49), the regularization term and the inertia term become zero when the number of alternating repetitions cycle exceeds the number of times obtained by multiplying the preset number of repetitions CYCLE by a constant Cr of about 0.5 to 0.8. , the coefficient of each term is monotonically decreasing. By reducing the regularization term and the inertia term as the calculation progresses in this way, the statistical variance can be learned based on the spatial information of the source signal. The initial value ρ _D ^(l) of the regularization term is a constant of about 10 to 30. When the mixing matrix is under bad conditions such as when the source signals are close to each other, the initial value ρ D (l) is small. It is better to set a large value when the degree is high.

ｂ．目的関数の勾配算定
分離信号と源信号分散及び式（４１）で求まる源信号空間情報に基づく正則化項の勾配を用いて、式（３９）により目的関数の勾配を求める。
ｃ．分離行列の更新
式（４７）で分離行列と分離信号を更新する。繰り返しインデックスｌがＬ回に達していなければ、ｌを１つ増やしてステップｂに戻る。ｌがＬ回に達したら、次の式で分離行列更新の際に用いる変数の初期値を設定し、繰り返しインデックスｌを１に初期化して源信号分散の更新に移る。

ｄ．源信号分散の更新
式（４８）で源信号分散を更新する。初回更新時（ｃｙｃｌｅ＝１）には３回程度反復更新し、それ以降の更新時（ｃｙｃｌｅ≧２）には、１回ずつ更新する。
ｅ．収束判定
繰り返しインデックスｃｙｃｌｅを１つ増やす。この時、ｃｙｃｌｅ≦ＣＹＣＬＥならばステップｂに戻る。The separation matrix calculation flow for minimizing the objective function LL(W(f)) obtained by adding the constraint condition LS(W(f)) regarding the source signal space prior information to the _LMNMF described above is summarized below.
a. Initial value setting Initialize the repetition index (superscript parenthesis character) as cycle=1, l=1. The number of updates CYCLE and L of each index is determined in advance. Initial values of the separation matrix and the source signal variance are set by the equations (32) to (36), respectively. Initialize the separation matrix update parameters with the following formula.

b. Gradient Calculation of Objective Function The gradient of the objective function is determined by Equation (39) using the gradient of the regularization term based on the separation signal and source signal variance and the source signal spatial information determined by Equation (41).
c. Updating Separation Matrix The separation matrix and the separated signal are updated by Equation (47). If the iteration index l has not reached L times, increase l by one and return to step b. When l reaches L times, the initial values of the variables used in updating the separation matrix are set according to the following equation, the iteration index l is initialized to 1, and the source signal variance is updated.

d. Update Source Signal Variance Update the source signal variance with equation (48). At the time of the first update (cycle=1), the update is repeated about three times, and at the time of subsequent updates (cycle≧2), the update is performed once.
e. Convergence judgment Increase the repetition index cycle by one. At this time, if cycle≤CYCLE, the process returns to step b.

When updating the separation matrix in step c, after the number of alternate iterations cycle exceeds half of the total number of iterations CYCLE and the regularization term becomes 0, the objective function shown in Equation (37) becomes the same as the normal MNMF. , its solution space becomes a simple shape. Therefore, after the regularization term becomes 0, it becomes possible to update the separation matrix quickly and stably by applying the iterative projection method (Non-Patent Document 18). Even in such a case, when the regularization term according to claim 3 is added to the objective function until the regularization term becomes 0, or when the separation matrix initial value according to claim 4 is used , when the separating matrix updating method described in claim 5 is applied, it is regarded as the scope of this patent.

[Base vector calculator]
If the direction of arrival of the source signal is not known, the amplitude of the complex vector y(f, τ) of the separated signal calculated by the separated signal calculator is not uniquely determined, and remains indefinite. Therefore, frequency distortion occurs when reconstructing a time domain signal from all frequency components by inverse Fourier transform. Therefore, by multiplying W(f) by the inverse matrix of the separation matrix, the separation matrix W ^(S) (f) corresponding to the M'-th sensor is obtained.

上記分離信号算定部で求めた分離行列を用いて式（５０）により、分離信号を算定する。

Separated signals are calculated according to the equation (50) using the separation matrix obtained by the separated signal calculation unit.

The isolated signals calculated in the previous section suffer from the problem of commutation ambiguity that the isolated signal with the nth source signal index does not necessarily correspond to the nth source signal. As a method to solve the ambiguity of the exchange in separated signals for each frequency bin, there is a method of constructing basis vectors using the feature values of the source signals, and a method of constructing basis vectors using information on the direction of arrival of the source signals contained in the separation matrix. can be considered as a method of configuring Of these, only the latter method is shown.

周波数伝達関数の基準センサにおける値との比を分離行列の逆行列を用いることで、式（５１）により近似する。

By using the inverse matrix of the separation matrix, the ratio of the frequency transfer function to the value of the reference sensor is approximated by Equation (51).

Ｍ’は基準センサの番号、ｄ_ｍａｘは基準センサと他のセンサの距離の最大値である。式（５２）は、フォンミーゼス・フィッシャー分布を仮定したクラスタリング手法に適用できる基底ベクトルｖ_ｎ（ｆ）の算定式であり、分離行列の逆行列から周波数伝達関数を算定して導くことができる。上記、基底ベクトルは周波数（波長）に依存せず、センサと源信号の相対位置によって値が変化する構成になっているため、周波数ごとの交換不定性解決のためのクラスタリングに用いることができる。又、ノルムが１となるように正規化を行っているため、球面上の分布となり、フォンミーゼス・フィッシャー分布を仮定したクラスタリング手法が適用可能となる。

M' is the reference sensor number, and d _max is the maximum distance between the reference sensor and the other sensors. Equation (52) is a formula for calculating basis vectors v _n (f) that can be applied to a clustering method assuming a von Mises-Fischer distribution, and can be derived by calculating a frequency transfer function from the inverse matrix of the separation matrix. Since the basis vectors do not depend on the frequency (wavelength) and the value changes depending on the relative position of the sensor and the source signal, they can be used for clustering to solve exchange ambiguity for each frequency. Also, since normalization is performed so that the norm is 1, the distribution is on a spherical surface, and a clustering method assuming a von Mises-Fischer distribution can be applied.

[Clustering part]
Consider the conditional probability Pr(v _n (f)|θ _c ) of the M-dimensional basis vectors in the group where the set of unknown parameters θ _c (hereinafter referred to as model parameters) in the c-th cluster is known. At this time, according to _Bayes ' _theorem , the model parameter The weighted log-likelihood LC of the conditional probability of is defined as follows.

ｃ＝１，…，Ｃはクラスタインデックスである。対数尤度の下限はＪｅｎｓｅｎの不等式を用いると、次のように与えられる。

c=1, . . . , C is the cluster index. Using Jensen's inequality, the lower bound of the log-likelihood is given as follows.

ｕ_ｃｎ（ｆ）は、Ｍ次元基底ベクトルのクラスタへの帰属度合を表す実変数であり、帰属度と称する。帰属度は、周波数ビンｆ、信号ｎごとの和が常に１となる。式（５４）の右辺は、重み付き対数尤度ＬＣの下限を与えるため、右辺を最大化して下限を持ち上げれば、対数尤度ＬＣを最大化する変数の組み合わせが求まる。

u _cn (f) is a real variable representing the degree of membership of the M-dimensional basis vector to the cluster, and is called the degree of membership. The degree of membership always sums to 1 for each frequency bin f and each signal n. Since the right side of equation (54) gives the lower bound of the weighted log-likelihood LC, the combination of variables that maximizes the log-likelihood LC can be obtained by maximizing the right side and raising the lower bound.

In the following, we assume the von Mises-Fischer distribution as a mixture model and expand the formula. At this time, the probability density function is given by the following equation.

この時、θ_ｃは、次の式で与えられる。

At this time, _θc is given by the following equation.

ここに、μ_ｃは各クラスタにおける基底ベクトルの平均方向を表すＭ次元ベクトルであり以降、平均方向と称する。κ_ｃは、クラスタごとの基底ベクトルの方向集中度を表す実定数であり以降、方向集中度と称する。変数Ｉ_{ｏ．ｒＭ－１}（κ_ｃ）は、（０．５Ｍ－１）階の第１種変形ベッセル関数である。

Here, μc is an M-dimensional _vector representing the average direction of basis vectors in each cluster, and is hereinafter referred to as the average direction. κ _c is a real constant representing the degree of directional convergence of basis vectors for each cluster, and is hereinafter referred to as the degree of directional convergence. Variable I _{o. rM−1} (κ _c ) is a modified Bessel function of the first kind of order (0.5M−1).

In order to maximize the right-hand side of Equation (55) while constraining the average direction to satisfy the normalization condition while considering the domain of mixture weights and degrees of membership, we define the following optimal conditions: .

Equation (58) can be derived by imposing three types of equality constraints on the objective function of Equation (57) using Lagrangian variables and then partially differentiating with the membership degree u _cn (f) to obtain the extreme value.

同様にパラメータ集合θ_ｃに含まれる変数の内、π_ｃ及びμ_ｃで偏微分して極値を求めることで、各変数を算定できる。

Similarly, among the variables included in the parameter set θc, each variable can be calculated by partially differentiating with respect to _πc and _μc and _obtaining the extreme value.

Since the degree of directional concentration of basis vectors is included in the modified Bessel function of the first kind, the value cannot be obtained as a calculation formula. Therefore, the ratio of modified Bessel functions of the first kind with different orders is defined.

この比率を連分数に展開した後、一次近似することで、方向集中度κ_ｃを算定する。

After expanding this ratio into a continued fraction, linear approximation is performed to calculate the degree of directional concentration _κc .

式（６１－１）の近似的に求めた比率と、式（６２－１）の実際に基底ベクトルの方向集中度を第１種変形ベッセル関数に代入して求めた比率との誤差が最小となるように、ニュートン法に基づく解の更新を行う。具体的には、式（６２－１）と式（６２－２）を交互に適用して、最適な基底ベクトルの方向集中度を探索する。更新は２回程度行えばよい。

The error between the approximately calculated ratio of formula (61-1) and the ratio of formula (62-1) obtained by substituting the directional convergence of the basis vectors into the modified Bessel function of the first kind is minimized. Update the solution based on Newton's method so that Specifically, by alternately applying equations (62-1) and (62-2), the optimum directional concentration of basis vectors is searched for. Updating may be performed about twice.

基底ベクトルの方向集中度は、次の式変換により、全てのクラスタで同一の値と仮定することもできる。この時、変数に対してクラスタインデックスｃを削除する置換処理を施すことにより、基底ベクトルとその重心ベクトルとの距離が近い、すなわち内積μ_ｃ ^Ｈｖ_ｎ（ｆ）が大きいほど、そのクラスタへの帰属度が大きくなるという単純な関係性が導けるため、基底ベクトルの方向集中度の推定誤差の影響を受けにくくなり、数値計算の安定性が向上する。一方で、源信号数が観測点数を超えている場合や観測ノイズが含まれる場合などの方向集中度がクラスタ間で異なる場合や、基底ベクトルの重みが互いに異なる場合は正確な解が導けない。しかし、近似的にクラスタの位置を求めたい場合には、数値安定性に優れるこの仮定は便利である。よって、源信号数とセンサ数が同数の場合において、初期クラスタ決定時には、一時的に基底ベクトルの重みを全て１とした上で、方向集中度が、全てのクラスタで同一とする仮定を用いることを推奨する。源信号数がセンサ数よりも多い場合は、クラスタの規模が互いに大きく異なる可能性が高いので、初期クラスタ決定時においてもクラスタごとに方向集中度を算定することを推奨する。方向集中度が、全てのクラスタで同一とする仮定を用いた場合、方向集中度は全クラスタの平均値に相当する次の式で近似される。

The directional convergence degree of basis vectors can also be assumed to be the same value for all clusters by the following expression conversion. At this time, by performing a replacement process for deleting the cluster index c from the variable, the closer the distance between the base vector and its centroid vector, that is, the larger the inner product μ _c ^H v _n (f), the more likely the cluster will be assigned to it. Since a simple relationship that the degree of membership is increased can be derived, the influence of the estimation error of the directional concentration degree of the basis vectors is reduced, and the stability of the numerical calculation is improved. On the other hand, when the number of source signals exceeds the number of observation points, when observation noise is included, and the degree of directional concentration is different between clusters, or when the weights of basis vectors are different, an accurate solution cannot be derived. However, this assumption, which is excellent in numerical stability, is convenient when it is desired to approximate the positions of clusters. Therefore, when the number of source signals and the number of sensors are the same, when determining the initial cluster, the weight of the basis vector is temporarily set to 1, and the directional concentration is assumed to be the same for all clusters. recommended. If the number of source signals is greater than the number of sensors, it is highly likely that the scales of the clusters will differ greatly from each other. Therefore, it is recommended to calculate the degree of directional concentration for each cluster even when determining the initial clusters. Using the assumption that the degree of directional concentration is the same for all clusters, the degree of directional concentration is approximated by the following equation, which corresponds to the average value of all clusters.

The base weight α _n (f) _is a value calculated based on the reliability of the base vector and can be obtained in various ways. It is considered to be the simplest to use it as a base weight.

When the number of source signals is unknown, the number of source signals is estimated by finding the optimal number of clusters. The optimal number of clusters is determined by the following procedure. First, several patterns of the number of clusters C are set, and cluster analysis is performed by repeated calculation of equations (58) to (63). The result is evaluated by equation (64), and the number of clusters that minimizes the evaluation index EV is assumed to be optimal.

ρは定数で２又は３程度の値を用いる。

ρ is a constant and a value of about 2 or 3 is used.

Restriction that basis vectors (separated signals) having different signal numbers in the same frequency bin do not belong to the same cluster based on the calculation results of the degree of membership and the model parameters obtained by the update formulas of formulas (58) to (63). While adding , the cluster division is performed in detail. Conventionally, in this cluster division, an evaluation value is set, and a round-robin search is performed for a combination of base vector exchanges that gives the best value. However, this method has the problem that when the number of source signals increases, the number of search combinations increases and the calculation takes a long time. Therefore, we propose a high-speed clustering method by limiting the range of optimal solution search. In order to prevent clusters from converging to a local solution when the number of source signals increases, values obtained by updating formulas (58) to (63) are used as the initial values of the degree of membership and model parameters.

In signal separation by independent component analysis or the like, it can be assumed that separation matrices having different signal numbers in the same frequency bin do not belong to the same cluster, so the number of clusters C is always greater than or equal to the number of source signals N. Also, when the degree of membership is binarized as 0 or 1, the log-likelihood lower limit value becomes 0 when the degree of membership is 0, and has a value when the degree of membership is 1. Therefore, the log-likelihood lower limit when the degree of membership is 1 is defined as the partial log-likelihood by Equation (65). The partial log-likelihood is the same as the log-likelihood before the degree of membership is given as an auxiliary variable, resulting in a three-dimensional tensor structure for each cluster/signal in each frequency bin.

確率分布がフォンミーゼス・フィッシャーズ分布に従う時、部分対数尤度は式（６６）となる。

When the probability distribution follows the von Mises-Fischer's distribution, the partial log-likelihood becomes Equation (66).

フォンミーゼス・フィッシャーズ分布において、帰属度が０か１で表され、基底ベクトルの方向集中度がクラスタ間において一定で、基底ベクトルの信頼度が全て同一のため、基底重みα_ｎ（ｆ）が全て１とできるという特定の条件が揃った時、部分対数尤度は以下により、簡易的に算定することもできる。

In the von Mises-Fischer's distribution, the degree of membership is represented by 0 or 1, the directional concentration of basis vectors is constant between clusters, and the reliability of basis vectors is all the same, so all basis weights α _n (f) are When the specific condition that it can be 1 is satisfied, the partial logarithmic likelihood can also be simply calculated by the following.

周波数ビンごとに、クラスタ番号を行成分、信号番号を列成分にもつ次の行列を考える。周波数ビンを固定して、クラスタ・信号番号の２次元テンソル構造に着目すると、クラスタ・信号番号が重複しないように、選んだ部分対数尤度のＮ個のテンソル和が最大となる組み合わせが、帰属度を０又は１と２値化した時の、対数尤度下限の最大値を与える。そのため、対数尤度を最大化する最適解を探索するには、部分対数尤度において、周波数ビンごとに、クラスタ・信号番号の２次元テンソル構造を考え、行と列が重複しないように選んだテンソルの組み合わせからその和が最大となる組み合わせをみつければよい。この時、和が最大となる部分対数尤度のテンソルの組み合わせの候補は、周波数ビンごとに値が大きい順に並び替えた時、その上位Ｎ－１位のテンソルが含まれている組み合わせから選べばよいことは明らかである。このことを念頭において、最適な組み合わせを部分対数尤度のテンソル及びその行と列番号の集合を用いて探索する手法を以下に示す。

Consider the following matrix with the cluster number as the row component and the signal number as the column component for each frequency bin. Focusing on the two-dimensional tensor structure of the cluster/signal number with the frequency bin fixed, the combination that maximizes the sum of the N tensors of the partial logarithmic likelihoods selected so that the cluster/signal number does not overlap is the attribution Gives the maximum value of the lower log-likelihood when the degree is binarized as 0 or 1. Therefore, in order to search for the optimal solution that maximizes the log-likelihood, in the partial log-likelihood, for each frequency bin, a two-dimensional tensor structure of clusters and signal numbers is considered, and rows and columns are selected so that they do not overlap. It suffices to find the combination that maximizes the sum from the combination of tensors. At this time, the candidates for the combination of partial log-likelihood tensors that maximize the sum should be selected from combinations that include the top N-1 tensors when rearranged in descending order of values for each frequency bin. The good is clear. With this in mind, a technique for searching for the optimal combination using a partial log-likelihood tensor and its set of row and column numbers is presented below.

部分対数尤度、クラスタ番号、信号番号に関する３種類の集合を定義し、それぞれ上付き文字としてＬ，Ｃ，Ｎをつける。
部分対数尤度Ｌ_ｃｎ（ｆ）を周波数ビンごとに大きい順に並べた集合をＢ^Ｌ（ｆ）、対応するクラスタ（式（６８）の行）と源信号（式（６８）の列）の番号の集合をそれぞれＢ^Ｃ（ｆ）及びＢ^Ｎ（ｆ）と定義し、それらの和集合をＢ（ｆ）とおく。又、周波数ビンごとに、部分対数尤度、クラスタ番号、源信号番号の各集合要素の内、ｉ番目の要素を集めた集合をＢ_ｉ（ｆ）と定義する。以降、使用するその他の同じ集合に対しても同様の表記規則を用いる。

We define three sets of partial logarithmic likelihoods, cluster numbers, and signal numbers, with L, C, and N as superscripts, respectively.
B ^L (f) is a set in which the partial logarithmic likelihood L _cn (f) is arranged in descending order for each frequency bin, and the number of the corresponding cluster (row of equation (68)) and source signal (column of equation (68)) are defined as B ^C (f) and B ^N (f), respectively, and their union set is B(f). Also, for each frequency bin, a set of the i-th elements among the set elements of the partial logarithmic likelihood, cluster number, and source signal number is defined as B _i (f). Similar conventions are used hereinafter for other identical sets used.

ここで、上付文字ｌ＝１～Ｌは繰り返しインデックスである。又、ｉ＝１，…，Ｎは集合要素の位置を表す。
式（６９）の集合構造を以降、Ｂ（ｆ）｛ｂ（ｆ）｝と短縮表記し、その他の集合にも同じ表記を適用する。最適解の条件を満たす要素すなわち部分対数尤度とそのクラスタ、信号番号を格納する集合をΩ^（ｌ）｛ω^（ｌ）｝、条件を満たさない要素を格納する集合をΨ｛ψ｝を定義する。アルゴリズム内の記号φは空

Here, the superscripts l=1 to L are repetition indices. Also, i=1, . . . , N represent the positions of set elements.
The set structure of equation (69) is hereafter abbreviated as B(f){b(f)}, and the same notation applies to the other sets. Define Ω ^(l) {ω ^(l) } as the set that stores the elements that satisfy the conditions of the optimal solution, that is, the partial logarithmic likelihood and its cluster, and the signal number, and Ψ {ψ} as the set that stores the elements that do not satisfy the conditions. do. The symbol φ in the algorithm is empty

１．繰り返しインデックスと最適解候補格納用集合を初期化し、探索開始解を登録する。

２．最適解でない要素の格納用集合Ψを初期化する。

３．集合Ｂ（ｆ）の要素に対して、最適解条件を満たすものと満たさないものに分ける。

４．最適解条件を満たさない集合Ψ内の要素が過去の繰り返し計算で最適解集合Ωに含まれているか調べ、含まれている時は集合Ψから該当要素を削除する。

５．集合Ψが空集合ならば、ステップ７へ、そうでなければ、ステップ６へ進む。
６．空集合でない集合Ψの要素の内、要素番号が最小の要素が最適解条件を満たす集合の最初の要素Ω_１に保存して、ステップ２へ戻る。

ステップ２へ戻る。
７．最適解条件を満たす集合Ωの内、評価値の最も高い集合を最適解として帰属度に反映する。

８．ｆ＜Ｆならステップ１へ戻る。ｆ＝Ｆなら繰り返し計算を終了する。・Algorithm 1
0. Initialize the frequency index

1. Initialize the iteration index and the set for storing optimal solution candidates, and register the search start solution.

2. Initialize the storage set Ψ of non-optimal elements.

3. The elements of the set B(f) are divided into those that satisfy the optimal solution condition and those that do not.

4. It is checked whether an element in the set Ψ that does not satisfy the optimum solution condition is included in the optimum solution set Ω in past repeated calculations, and if it is included, the corresponding element is deleted from the set Ψ.

5. If the set Ψ is an empty set, go to step 7, else go to step 6.
6. Among the elements of the non-empty set Ψ, the element with the smallest element number is stored in the _first element Ω1 of the set that satisfies the optimum solution condition, and the process returns to step 2.

Return to step 2.
7. Among the sets Ω that satisfy the optimal solution conditions, the set with the highest evaluation value is reflected in the degree of membership as the optimal solution.

8. If f<F, return to step 1. If f=F, terminate the iterative calculation.

The calculation flow for generating the optimal cluster that maximizes the logarithmic likelihood when the degree of membership is binarized as 0 or 1 is summarized below. First, by solving equations (59) to (63) using the calculated value of the degree of membership of equation (58), the model parameter θ _c that satisfies the posterior probability maximization condition is obtained. Using the model parameters, the partial log-likelihoods are calculated, and Algorithm 1 searches for the combination of partial log-likelihoods that maximizes the log-likelihood. This optimal combination is reflected in the degree of membership. By repeating this calculation, the optimal combination of clusters that maximizes the log-likelihood can be identified. Specifically, the following iterative calculation is performed.

i) Calculation of variable set θ _C The variable set θ _C that increases the logarithmic likelihood under the condition that the degree of membership is given is updated and calculated using equations (58) to (63).
ii) Calculation of Partial Log-Likelihood From the variable set _θC , the partial log-likelihood L _ci (f) for each frequency bin/signal in each cluster is calculated by Equation (66) or Equation (67).
iii) Determining the degree of membership Based on the calculation results of the partial log likelihood, search for a combination that allocates the basis vectors to the clusters so as to maximize the log likelihood by steps 2 to 6 of Algorithm 1, and select the optimal allocation combination. , the degree of membership is determined as 0 or 1 by step 7 of Algorithm 1.
iv) In the iterative calculation process, the calculations i) to iii) are repeated until there is no change in the degree of membership defined by 0 or 1 or the mean direction.

Separation matrices calculated by optimization may contain various estimation errors, such as estimation errors of source signal waveforms due to noise, etc., and errors due to unsuccessful learning itself. Therefore, a continuous value mask M _n (f) is defined to evaluate the estimation accuracy of the separation matrix and reduce the influence of the separation signal calculated from the low-reliability separation matrix. The calculation method of the continuous value mask differs depending on the signal separation environment and its purpose, so it is not specified, but an effective calculation procedure for obtaining main signal data for each direction is shown as an example. The correctly calculated inverse matrix of the separation matrix corresponds to the mixing matrix having the frequency transfer function for each source signal as a column vector. The column vectors of the mixing matrix have the property that they are orthogonal to each other in the whitened space. This property is used to determine whether the separation matrix has been calculated correctly. First, the demultiplexing matrix is inverted by the whitening matrix V(f) to compute the N×N-dimensional matrix G(f) transferred to the whitening space at each frequency bin. The whitening matrix is calculated from the eigenvalues and eigenvectors of the correlation matrix E[x(f, τ)x ^H (f, τ)] of the observed signal. , is used when decorrelating between signals.

行列Ｇ（ｆ）の列ベクトルｇ_ｎ（ｆ）が互いになす角度をその内積の大きさにより評価する。評価は周波数ビンごとに行い、同一周波数ビンの他列ベクトルとの内積の２乗和の平均値が一定値を超えた源信号インデックスを異常とみなし、異常度合に応じてマスクを０に近づける。それ以外の正常値については、連続値マスクを１に近づける。この判定をシグモイド関数により連続値として行うと、次の式により連続値マスクが求まる。

The angle between the column vectors g _n (f) of the matrix G(f) is evaluated by the size of the inner product. The evaluation is performed for each frequency bin, and the source signal index where the average value of the sum of squares of the inner products of the same frequency bin with other column vectors exceeds a certain value is regarded as abnormal, and the mask is brought closer to 0 according to the degree of abnormality. For other normal values, the continuous value mask is brought closer to 1. When this determination is performed as a continuous value using a sigmoid function, a continuous value mask is obtained from the following equation.

Ｋｉｊｕｎは、列ベクトルの内積の和が異常と判定する基準値であり、０～１の値をとる。Ｋｏｕｂａｉは、正の実数で、シグモイド関数が１から０に切り替わる速さを表し、値が大きいほど０か１のどちらかの値のみを出力する傾向が強まる。

Kijun is a reference value for determining that the sum of inner products of column vectors is abnormal, and takes a value of 0-1. Koubai is a positive real number and represents the speed at which the sigmoid function switches from 1 to 0, and the larger the value, the stronger the tendency to output either 0 or 1 only.

連続値マスクの値が閾値Ｔｈｒｅｓｈｏｌｄを下回った分離信号を除外するために、バイナリマスクＭ_ｎ ^（Ｂ）（ｆ）を定義する。

A binary mask M _n ^(B) (f) is defined to exclude isolated signals whose continuous value mask values are below a threshold Threshold.

The separated signals are multiplied by a degree of membership defined as 0 or 1 and summed clusterwise to compute the separated signal y _c ^(ICA) (f, τ) with one source signal in one row vector.

以上のクラスタリング方法により、効率的に交換の不定性を解決できる。

The above clustering method can efficiently solve the ambiguity of exchange.

The separated signal y _c ^(ICA) (f, τ) contains a separated signal with low calculation accuracy, and if it is output as it is, the error in the source signal estimation result will increase. Therefore, after summarizing the source signal direction-of-arrival calculation method, a method for reducing the influence of separated signals with low calculation accuracy is summarized as an arithmetic expression in the separated signal correction unit.

[Source signal direction-of-arrival calculation unit]
First, the source signal arrival direction vector q̂ _n (f) is estimated using the Moore-Penrose pseudo-inverse matrix P ⁺ of the observation point coordinate vector (Non-Patent Document 12). In equation (74), the position of the reference sensor is the position of the origin of the direction-of-arrival vector of the source signal, that is, the viewpoint.

行列Ｐ及びベクトルｒ＾_ｎ（ｆ）はそれぞれ、観測点番号として基準センサ以外の（Ｍ－１）通りの位置ベクトルに対して、基準センサを起点とする観測点座標ベクトル及び周波数伝達関数の偏角ａｒｇを羅列したものである。

Matrix P and vector r̂ _n (f) are observation point coordinate vectors and deviations of frequency transfer functions from the reference sensor with respect to (M−1) position vectors other than the reference sensor as observation point numbers, respectively. It is a list of angles arg.

源信号及び周波数ビン番号に対して源信号到来方向推定値を加算することでクラスタごとの信号方向ベクトルＱ＾_ｃを式（７５）により算定する。この源信号到来方向推定値に対してカメラにより事前に撮影した画像を重ね合わせることで、信号到来方向の推定精度を向上させることもできる。
Ｍ_ｎ（ｆ）は、連続値マスクであり、周波数ビン・源信号ごとの分離行列の信頼度に基づいて設定する。上記算定部で掛けている連続値マスクは、背景ノイズのほとんどない場合など分離行列の信頼度が高い場合などは、省略することもできる。Ｗｈｅｉｇｈｔ ｆ）は周波数ビン・源信号ごとの源信号到来方向ベクトルの重みであり、どの周波数も同じ信頼度なら全て１となる。重みとして周波数ごとの観測値のパワースペクトルを指標とする時は式（７６）で算定される。他にも重みとしては、信号部分空間の固有値の和を指標とすることもでき、近似的に直接、観測点まで届いた到来波のパワーの総和をとっていることになる。

The signal direction _vector Q̂c for each cluster is calculated by Equation (75) by adding the source signal direction-of-arrival estimate to the source signal and frequency bin number. By superimposing an image captured in advance by a camera on the source signal arrival direction estimation value, it is possible to improve the estimation accuracy of the signal arrival direction.
M _n (f) is a continuous value mask and is set based on the reliability of the separating matrix for each frequency bin and source signal. The continuous value mask multiplied by the calculation unit can be omitted when the separation matrix is highly reliable, such as when there is almost no background noise. Wheight f) is the weight of the direction-of-arrival vector of the source signal for each frequency bin and source signal, and is all 1 if all frequencies have the same reliability. When the power spectrum of the observed value for each frequency is used as the index as the weight, it is calculated by Equation (76). Alternatively, the sum of the eigenvalues of the signal subspace can be used as an index for the weight, which is approximately the sum of the powers of the incoming waves that reach the observation point directly.

ＰＳＤ（ｆ）は、全センサのパワースペクトルを時間フレーム方向に加算した値である。
〔分離信号補正部〕

PSD(f) is a value obtained by adding the power spectra of all sensors in the time frame direction.
[Separation signal corrector]

クラスタ番号ｋを混合行列の列番号すなわち源信号番号へ変換する順列Π_ｆ（ｋ）を式（７７）で定義する。この順列は、周波数ビンｆにクラスタ番号ｋと対応する源信号が存在しない時に０を返すため、源信号の存在を調べる指標として用いることもできる。周波数ビンｆに存在する源信号の数Ｎ^{（Ａｃｔ）}（ｆ）は、式（７８）で帰属度を源信号番号及びクラスタ番号の方向へ加算することで求まる。

A permutation Π _f (k) that transforms the cluster number k into the column number of the mixing matrix, that is, the source signal number, is defined by equation (77). Since this permutation returns 0 when the source signal corresponding to the cluster number k does not exist in the frequency bin f, it can also be used as an index for checking the existence of the source signal. The number of source signals N ^(Act) (f) existing in the frequency bin f is obtained by adding the degree of belonging in the direction of the source signal number and cluster number in Equation (78).

式（７５）で得られる源信号到来方向推定結果を用いて、Ｍ×Ｎ^{（Ａｃｔ）}（ｆ）次元の混合行列Ａ＾（ｆ）を算定する。

Using the direction-of-arrival estimation result of the source signal obtained by equation (75), an M×N ^(Act) (f)-dimensional mixing matrix Â(f) is calculated.

この混合行列の逆行列を求めることで、簡易的に分離行列と分離信号が算定できる。

By obtaining the inverse matrix of this mixing matrix, the separation matrix and the separated signal can be easily calculated.

ｅ_ｎ（ｆ）はｎ番目の要素が１、他要素が０のＮ^{（Ａｃｔ）}（ｆ）次元単位ベクトルである。式（８０－３）の上付括弧文字＋は、ムーアペンローズの擬似逆行列を表す。源信号到来方向より得られる分離信号ｙ_ｃ ^{（ＤＯＡ）}（ｆ，τ）をベクトル形式でまとめると、

e _n (f) is an N ^(Act) (f)-dimensional unit vector in which the n-th element is 1 and the other elements are 0s. The superscript character + in equation (80-3) represents the Moore-Penrose pseudo-inverse. Putting together the separated signals y _c ^(DOA) (f, τ) obtained from the direction of arrival of the source signal in vector form,

When the number of clusters {C} and the number of source signals N are the same, the computation for calculating the separated signals from the direction of arrival of the source signals shown in Equations (77) to (81) can be summarized simply by the following equations. .

Using the separated signal calculation result by ICA shown in Equation (73) and the separated signal calculation result based on the direction of arrival of the source signal, the separated signal y ^(MOD) (f , τ).

この式ではバイナリマスクを用いて、ＩＣＡにより算定される分離信号ｙ^{（ＩＣＡ）}（ｆ，τ）の内、算定精度の低い分離信号を取り除き、代わりに源信号到来方向に基づく分離信号ｙ^{（ＤＯＡ）}（ｆ，τ）を加算して、取除分を埋め合わせている。

In this equation, a binary mask is used to remove poorly calculated separated signals from the separated signals y ^(ICA) (f, τ) calculated by ICA, and instead replace the separated signals y ^{(DOA )} (f, τ) is added to make up for the subtraction.

If the source signal spatial prior information is obtained, the separated signal of formula (50) is obtained, and if not obtained, the separated signal of formula (83) is subjected to inverse Fourier transform, multiplied by a window function, and overlapped to obtain separation A time waveform of the signal is obtained.

The present invention is not limited to the embodiments of the invention described above, and various modifications are possible within the scope of the invention described in the claims, and these are also included in the scope of the invention. be.

The text of this patent describes the case of using independent component analysis and independent low-rank matrix analysis as the separating matrix, but even when using other separating matrices, the clustering method or source signal arrival direction estimation in the present invention Any use of a method is considered to be within the scope of the present invention.

As a method for source signal waveforms, there is independent deep learning matrix analysis as a method for making artificial intelligence pre-learn the statistical variance information of the source signal and performing signal separation with high accuracy. In this method, the separation matrix is updated using an iterative projection method based on the auxiliary function method. However, the independent deep learning matrix analysis has a problem that the statistical variance information of the source signal cannot be correctly learned unless the information of the source signal is obtained at a certain level in advance. By using the separation matrix update method based on prior information constraint addition of Equation (37) in which the information on the direction of arrival of the source signal is added as a constraint condition shown in Equation (41) instead of the iterative projection method in independent deep learning matrix analysis , it becomes possible to estimate the source signal waveform with high accuracy even when there is little prior information of the separated signal. Therefore, application in independent deep learning matrix analysis (Non-Patent Document 19) is also conceivable. Even in such a case, when adding a constraint with the direction of arrival of the source signal as a priori information to the objective function, the constraint is defined as "a row vector of the separation matrix and a column vector of the mixing matrix having the same source signal number. The inner product is 1, and the inner product of the row vector of the separation matrix and the column vector of the mixing matrix with different source signal numbers is 0." considered within the scope of the present invention. Also, when the value obtained by multiplying the update amount of the separation matrix at the previous time step by a constant is used as the inertia term, or the update amount in the gradient descent direction for each matrix block, that is, the learning rate, is used as the previous separation matrix When the items described in claim 5 are met, such as when adjusting according to the update amount of or the gradient of the objective function, it is considered to belong to the scope of the present invention.

Even when applying the separation matrix update method with prior information such as the direction of arrival of the source signal as a constraint to the independent deep learning matrix analysis, the operation flow is the same as the independent low-rank matrix analysis, and the statistical variance of the signal The only difference is the update method. The basic operation flow is shown below.

分離信号のパワースペクトル構造に関する事前情報が得られている場合、ニューラルネットワークにより信号波形の特徴量から統計的分散を上式により推定することができる。ニューラルネットワークとして独立深層学習行列分析を用いた場合の具体的な分散推定の演算方法は、非特許文献１９を参照されたい。ニューラルネットワークによる分散推定の方法には、様々な手法があるが、いずれの手法を用いた場合でも制約条件の設定方法や分離行列の更新方法は変わらない。まず式（３７）に示す源信号到来方向に関する事前情報制約を付加した目的関数を用い、式（４７）のの分離行列更新ステップを数回、繰り返して分離行列を更新する。更新が完了した分離行列を用いて分離信号を算定し、ニューラルネットワークに入力する。ニューラルネットワークで分離信号波形の入力データと過去に学習した分離信号のパワースペクトル構造を比較し、最も適合する構造を選択してその分散を時間周波数点ごとに求める。得られた分散を用いて、再び上記計算ステップで分離行列を更新し、分離信号を算定する。この演算フローを繰り返すことで、分離信号を正確に推定することができる。この場合、本特許における請求項３に記載の分離行列更新方法と請求項５に制約条件の記載を用いていることになる。これらのいずれかを用いている場合は、本特許の請求範囲とする。

When prior information about the power spectrum structure of the separated signals is obtained, the neural network can estimate the statistical variance from the feature values of the signal waveform using the above equation. See Non-Patent Document 19 for a specific calculation method for variance estimation when independent deep learning matrix analysis is used as a neural network. There are various methods for estimating the variance using a neural network, but the method of setting constraints and the method of updating the separation matrix are the same regardless of which method is used. First, the separation matrix is updated by repeating the separation matrix update step of Expression (47) several times using the objective function to which the prior information constraint on the direction of arrival of the source signal shown in Expression (37) is added. Using the updated separation matrix, separated signals are calculated and input to the neural network. A neural network compares the input data of separated signal waveforms with the power spectrum structure of previously learned separated signals, selects the most suitable structure, and obtains its variance for each time-frequency point. Using the obtained variance, the separation matrix is updated again in the above calculation step, and the separation signal is calculated. By repeating this operation flow, the separated signals can be accurately estimated. In this case, the separating matrix updating method described in claim 3 and the description of the constraint condition in claim 5 of this patent are used. Use of any of these is covered by this patent.

In the clustering calculation according to the calculation flow of Figures 5 and 6, basis vectors are generated for each frequency bin f and source signal n. can also be replaced by time-frequency points to perform cluster analysis. Even in such a case, when cluster analysis is performed with an objective function based on the von Mises-Fischer distribution and the parameters are updated based on the optimality condition shown in Equation (57), claims 6 and 7 Or if it matches the clustering method of claim 8, it is considered to belong to the scope of the present invention.

The cluster analysis calculation proposed in the "clustering unit" according to claims 6 to 8 is performed by estimating the direction of arrival of the source signal for each frequency bin by the music method (Non-Patent Document 7) over the entire analysis symmetric frequency range, It can also be used for clustering. The application method is shown below.

水平角θと仰角Φを変えながら、式（８５）により源信号到来方向と対応する周波数伝達関数を算定する。式（８６）で算定されるミュージックスペクトルが局所的にピークになる時の水平角θと仰角Φを源信号到来方向として、源信号到来方向に対応する周波数ビンごとの基底ベクトルを式（１３）により、源信号ごとに算定する。

While changing the horizontal angle .theta. and the elevation angle .PHI. Assuming that the horizontal angle θ and elevation angle Φ at which the music spectrum calculated by equation (86) peaks locally is the direction of arrival of the source signal, the basis vector for each frequency bin corresponding to the direction of arrival of the source signal is expressed by equation (13). is calculated for each source signal by

Ｑ_Ｎ（ｆ）は、観測信号の相関行列の最小の固有値に対応した固有ベクトルである。周波数ビンごとに、源信号到来方向別に算定された基底ベクトルに対して、本特許のクラスタリング部に記載の演算（図－５及び図－６）を適用すると、解析対称周波数全域にわたって、源信号インデックスｎを対応付けることが可能となる。このような場合でも、クラスタ分析に、フォンミーゼス・フィッシャー分布に基づく目的関数を設定し、式（５８）～式（６３）の最適性条件の基でパラメータを更新している場合など、請求項６、請求項７又は請求項８に記載のクラスタリング方法と合致する場合は、本発明の範囲に属するとみなす。

Q _N (f) is an eigenvector corresponding to the minimum eigenvalue of the correlation matrix of the observed signal. Applying the calculations described in the clustering section of this patent (Figs. 5 and 6) to the basis vectors calculated for each frequency bin and for each source signal arrival direction yields the source signal index n can be associated. Even in such a case, cluster analysis, set the objective function based on the von Mises-Fischer distribution, and update the parameters under the optimality conditions of formulas (58) to (63), such as when the claim If the clustering method of claim 6, claim 7 or claim 8 is met, it is considered to belong to the scope of the present invention.

The cluster analysis calculation proposed in the "clustering unit" according to claims 6 to 8 is based on the KR product-compressed music method (Non-Patent Document 20), which is an extension of the music method. It can also be used when clustering the estimation results over the entire analytically symmetrical frequency range. The application method is shown below.

Using the M×M-dimensional correlation matrix for each time frame of the observed signal, the N ² ×T-dimensional data matrix Y ^˜ (f) is defined for each frequency bin by the following equation.

Ｖｅｃ（・）は行列のベクトル化を表す。データ行列の特異値分解を行い、特異値０に対応する左特異ベクトルから成る行列をＱ^～（ｆ）と定義する。この時、一般化ミュージックスペクトルＳｐｅｃ（ｆ，θ，Φ）は次の式で求まる。

Vec(·) represents matrix vectorization. Singular value decomposition of the data matrix is performed, and the matrix consisting of left singular vectors corresponding to singular value 0 is defined as Q ^~ (f). At this time, the generalized music spectrum Spec(f, .theta., .PHI.) is obtained by the following formula.

番号の列ベクトルにクロネッカー積を適用する演算であり、一般にＫＲ積と呼ばれる。推定源信号到来方向集合Ｃは、データ行列に対して、事前に直交マッチング追跡（非特許文献２１）などの圧縮センシングのアルゴリズムを適用することで、得られるＮ－Ｔ個の源信号到来方向推定値の集合である

The operation of applying the Kronecker product to a column vector of numbers, commonly referred to as the KR product. The estimated source signal direction-of-arrival set C is obtained by applying a compression sensing algorithm such as orthogonal matching tracking (Non-Patent Document 21) in advance to the data matrix, and NT source signal direction-of-arrival estimates are obtained. is a set of values

While changing the horizontal angle .theta. and the elevation angle .PHI. Assuming that the horizontal angle θ and the elevation angle Φ at which the music spectrum calculated by Equation (88) peaks locally is the source signal arrival direction, the basis vector for each frequency bin corresponding to the source signal arrival direction is given by Equation (13). is calculated for each source signal by Applying the calculations (Eqs. (58) to (63)) described in the clustering section of this patent to the basis vectors calculated for each frequency bin and for each direction of arrival of the source signal, the source It becomes possible to associate the signal index n. Even in such a case, in the cluster analysis, an objective function based on the von Mises-Fischer distribution is set, and the parameters are updated under the optimality condition shown in Equation (57). If it is consistent with the clustering method described in claim 7 or claim 8, it is considered to be within the scope of the present invention.

In addition, when using a compressed sensing algorithm such as orthogonal matching tracking alone, or when using the KR music method (Non-Patent Document 22) or the compressed music method (Non-Patent Document 23), each frequency bin Then, the direction of arrival of the source signal is obtained, and using the clustering operation invented in this patent, it is possible to associate the source signal index n over the entire analytical symmetric frequency range. Even in such a case, in the cluster analysis, an objective function based on the von Mises-Fischer distribution is set, and the parameters are updated under the optimality condition shown in Equation (57). If it is consistent with the clustering method described in claim 7 or claim 8, it is considered to be within the scope of the present invention.

The clustering calculations of equations (58)-(63) can also be applied to the separation of speech signals assuming a harmonic structure. The application method is shown below.

式（８９）の観測信号の時間フレームτと周波数ビンｆを時間周波数点のインデックスｋ＝１，…，Ｋを定義することで、同一次元にしてベクトル化する。この時、時間周波数点の各点ｋにおける観測信号ベクトルは、センサ数をＭとすると、次のＭ次元ベクトルで定義される。

The time frame τ and the frequency bin f of the observed signal in equation (89) are vectorized with the same dimension by defining the index k=1, . . . , K of the time-frequency point. At this time, the observed signal vector at each time-frequency point k is defined by the following M-dimensional vector, where M is the number of sensors.

時間周波数点の各点ｋごとに、ノルムが１のＭ次元埋め込みベクトルｖ^～ _ｎ（ｋ）を考える。

Consider an M-dimensional embedding vector v ^∼ _n (k) with a norm of 1 for each point k of time-frequency points.

ｎ＝１，…，Ｎは、源信号のインデックスを表し、各時間周波数点において、多くてもＮ個の源信号がアクティブと仮定している。埋め込みベクトルを用いると、時間周波数点の各点ｋにおいて、支配的な源信号波形を観測信号波形から算定することが可能となる。埋め込みベクトルは時間周波数点ごとに独立して生成されるため、異なる時間周波数点において、源信号インデックスｎの一貫性がなくなる。例えば、源信号数が２つの場合を考えると、時間周波数点１と時間周波数点２で、源信号の組み合わせは２パターン考えられる。そこで、全時間周波数点において、埋め込みベクトルをクラスタリングすることで、対応関係を明確化する。クラスタリングは、全源信号数をクラスタ数Ｃと定義し、ｆ→ｋとインデックスを置き換えて、「クラスタリング部」に示す演算方法すなわち式（５８）～式（６３）及びアルゴリズム１を適用することで、実施することができる。よって調波構造を仮定した音声信号の分離を行う場合においても、クラスタ分析に、混合モデルとしてフォンミーゼス・フィッシャー分布に基づく目的関数を設定し、式（５７）に示す最適性条件の基でパラメータを更新している場合など、請求項６又は請求項７又は請求項８に記載のクラスタリング演算手法と合致する場合は、本発明の範囲に属するとみなす。

n=1, . . . , N represents the index of the source signals, assuming that at each time-frequency point at most N source signals are active. Using the embedding vector, at each time-frequency point k, the dominant source signal waveform can be calculated from the observed signal waveform. Since the embedding vector is generated independently for each time-frequency point, the source signal index n will be inconsistent at different time-frequency points. For example, when the number of source signals is two, two combinations of source signals are conceivable at time-frequency point 1 and time-frequency point 2 . Therefore, the correspondence relationship is clarified by clustering the embedding vectors at all time-frequency points. Clustering is performed by defining the total number of source signals as the number of clusters C, replacing the index with f → k, and applying the calculation method shown in the "clustering section", that is, equations (58) to (63) and algorithm 1. , can be implemented. Therefore, even when separating speech signals assuming a harmonic structure, an objective function based on the von Mises-Fischer distribution is set as a mixture model in the cluster analysis, and the parameter If it matches with the clustering operation method described in claim 6, claim 7, or claim 8, such as updating the , it is considered to belong to the scope of the present invention.

1 signal acquisition unit 2 signal waveform separation unit 3 calculation result display unit

Claims (8)

When defining the vibration of air or an object as a source signal, and defining a signal obtained by a microphone or a sensor in a state in which the source signals arriving from multiple directions are mixed as an observation signal,
A direction-specific signal waveform separation method for separating the observed signal into the source signal waveforms by direction of arrival and estimating the direction of arrival of the source signal,
In the case where a matrix for separating the observed signal into the source signals for each direction of arrival is defined as a separation matrix, the observed signal is multiplied by the separation matrix to calculate a separated signal that is an estimated value of the source signal. A separated signal calculation unit characterized by a feature, a clustering unit that classifies the feature amount of the direction of arrival or the waveform shape of the source signal into categories by a cluster analysis method based on the von Mises-Fischer distribution, and the source signal based on the signal separation result. characterized by having a source signal arrival direction calculation unit for estimating the arrival direction as a calculation function,
When several pieces of spatial a priori information (a priori information on directions of arrival) are obtained from among the plurality of source signals, the separated signal calculation unit uses the observed signals and the a priori spatial information to performing an operation of spatial prior information-based separated signal calculation to separate signals by direction of arrival to calculate separated signals;
When the spatial a priori information of the source signal is not obtained, after the separated signal calculation unit calculates the separated signals based on the a priori information on the variance of the observed signal and the source signal, the clustering unit and characterized by performing the calculation of the direction-of-arrival calculation unit for the source signal,
Signal waveform separation method. A computing method for split signal estimation based on the spatial a priori information of claim 1, comprising:
The separated signal calculation method is a method of calculating separated signals y (f, τ) for each time frame τ and frequency bin f using the direction of arrival of the source signal as prior information,
calculating a frequency transfer function from the sensor to the source signal from the sensor arrangement position and the source signal arrival direction;
When defining a matrix in which the frequency transfer function is arranged so that the sensor number is a row vector and the source signal number is a column vector as a mixing matrix A(f),
The separation signal is calculated with the goal of maximizing the probability density (posterior probability) of the separation signal under the condition that the mixing matrix and the observed signal x (f, τ) at the sensor are known. year,
Regularization calculated by the ratio of the variance σ _E (f) of the noise and the variance σ _n (f) of the separated signal by assuming that the probability density of noise other than the target signal and the separated signal follows a Gaussian distribution calculating the inverse matrix of the matrix ^AH (f)A(f)+Σ(f) obtained by adding the term Σ(f) to the square of the mixing matrix;
The separation matrix W (f) is calculated from the inverse matrix by formula (92-2), and the separated signal is calculated by multiplying the observed signal by the separation matrix by formula (92-1). to be
The separation signal calculation method.

A computing method for split signal estimation based on the spatial a priori information of claim 1, comprising:
The separation matrix has the source signal number as a row vector and the sensor number as a column vector,
calculating a frequency transfer function from the sensor to the source signal from the sensor arrangement position and the source signal arrival direction;
When defining a matrix in which the frequency transfer function is arranged so that the sensor number is a row vector and the source signal number is a column vector as a mixing matrix,
characterized by setting information that the mixing matrix is an inverse matrix of the separation matrix as a constraint,
The constraints are such that the inner product of the row vector of the separation matrix and the column vector of the mixing matrix having the same source signal number is 1, and the row vector of the separation matrix having the source signal number different from each other. characterized by including information that the inner product of the column vector of the mixing matrix is 0,
Calculating the separation matrix by optimizing the objective function including the constraint,
Calculating the separated signal by multiplying the observed signal by the separation matrix,
A computational method for calculating the separation signal. A computing method for split signal estimation based on the spatial a priori information of claim 1, comprising:
In the process of calculating the separation matrix that optimizes the objective function for obtaining the accurate source signal waveform,
calculating a gradient descent direction by partially differentiating the objective function with the separation matrix;
As the separation matrix roughly calculated from the prior information of the direction of arrival of the source signal, a matrix having, as tensor components, an inverse matrix of the mixing matrix or a matrix obtained by Hermitian transposing the mixing matrix, or the formula according to claim 2 ( 92-2) using the matrix as the initial value of the separation matrix,
calculating the separation matrix that optimizes the objective function by iteratively updating the separation matrix for the gradient descent direction;
Calculating the separated signal by multiplying the observed signal by the separation matrix,
A computational method for calculating the separation signal. A method of searching for the optimal separation matrix using the objective function including the constraint according to claim 3,
When updating the separation matrix based on the objective function for obtaining the accurate source signal waveform, when the objective function has a complicated solution space, the separation matrix is stably updated without falling into a local solution. A method for updating, comprising:
The separation matrix has the source signal number as a row vector and the sensor number as a column vector,
When a set of tensors obtained by dividing the separation matrix into row vectors or a set of tensors obtained by dividing the separation matrix by an arbitrary range is defined as a matrix block, the update amount is adjusted for each matrix block to divide the separation matrix. characterized by updating
defining the sum of squares of the tensor components for each block of the matrix as the block-wise sum of squares of the matrix;
A gradient obtained by partially differentiating the objective function with the separation matrix is defined as a loss function,
When defining a coefficient to be multiplied by the loss function as a learning rate in order to adjust the update amount in the gradient descent direction of the objective function,
Obtaining the block-wise sum of squares of the loss function and adding it to the cumulative value of the block-wise sum of squares of the loss function up to the previous time step, thereby obtaining the loss function at the current time step. After calculating the cumulative value of the block-wise sum of squares,
Using the accumulated value of the block-wise sum of squares of the loss function at the current time step and the accumulated value of the block-wise sum of squares of the update amount of the separation matrix at the previous time step, the matrix block Calculate the learning rate for each
Updating the separation matrix at the current time step by adding the value obtained by multiplying the loss function by the learning rate to the value obtained by multiplying the update amount of the separation matrix at the previous time step by a constant. After calculating the amount for each matrix block,
Obtaining the block-wise sum of squares of the update amount of the separation matrix and adding it to the accumulated value of the block-wise sum of squares of the update amount of the separation matrix at the previous time step, thereby obtaining Based on the four-stage variable update step of calculating the cumulative value of the block-wise sum of squares of the update amount of the separation matrix in
Searching for the separation matrix that optimizes the objective function by repeating the operation of updating the separation matrix many times,
A search method for the separating matrix. The computing method in the clustering unit according to claim 1,
When defining a cluster as a collection of unknown parameters in a mixture model represented by a linear sum of multiple probability density functions,
calculating a signal basis vector as a feature quantity of the direction of arrival or waveform shape of the source signal for each frequency bin/source signal or for each time frame/frequency bin/source signal;
using the von Mises-Fischer distribution as the probability density function of the mixture model,
Define the membership probability of the signal basis vector to the cluster in the mixture model as the degree of membership,
The degree of importance for each cluster is defined as a mixture weight, the average value of the signal basis vectors for each cluster is defined as an average direction, and the degree of variation of the signal basis vectors around the average direction is defined as a direction concentration degree. case,
With the goal of maximizing the probability density function of the von Mises-Fischer distribution,
An operation for fixing the mixture weight, the average direction, and the direction concentration degree to update the degree of membership, fixing the degree of membership and updating the mixture weight and the average direction, and Using a likelihood estimation method, by alternately repeating the calculation of updating the direction concentration degree, calculating the optimum degree of membership, the mixture weight, the average direction, and the direction concentration degree,
Clustering calculation method. 2. The calculation method in the clustering unit according to claim 1, which can be used when the value of the degree of membership is limited to a discrete value of 0 or 1,
With the goal of maximizing the probability density function of the von Mises-Fischer distribution,
While imposing a constraint that the signal basis vectors belonging to the same frequency bin or time frame/frequency bin do not belong to the same cluster, the attribution is performed by fixing the mixture weight, the average direction, and the direction concentration degree. an operation that updates the degree to either a value of 0 or 1;
Fixing the degree of membership, updating the mixture weight and the average direction, and using an approximation formula, a gradient method, or a maximum likelihood estimation method to alternately repeat the operation of updating the degree of concentration in the direction, Calculating the degree of membership, the mixture weight, the average direction and the degree of concentration in the direction,
Clustering calculation method. 8. The clustering operation according to claim 6 or claim 7, wherein an operation method for calculating the degree of directional concentration _κc or κ,
The degree of directional concentration is also included in the modified Bessel function of the first kind I _0.5M (κ _c ) included in the von Mises-Fischer distribution formula,
Based on the signal basis vector v _n (f) and the signal basis vector, the basis weight α _n (f) set based on the degree of importance or reliability and the degree of membership u _cn (f), the rank 1 Calculate an approximate ratio, which is an approximate value of the ratio of the two different modified Bessel functions of the first kind, using formula (93) or formula (94),

Using the approximation shown in Equation (95) or Equation (96), which is obtained by expanding the ratio of the modified Bessel functions of the first kind that differ in order by one to continued fractions, the directional concentration degree calculate an approximation,

By substituting the approximation values of the degree of directional concentration into the modified Bessel functions of the first kind that differ by one rank, the detailed ratios, which are the detailed values of the ratios of the modified Bessel functions of the first kind that differ by one rank, are expressed by the formula ( 97) calculated by

With the goal of minimizing the error between the detailed ratio and the general ratio, a Newton method update formula (equation (98)) is applied, in which the gradient is the value obtained by differentiating the detailed ratio with the approximate value of the degree of directional concentration. and updating the degree of directional concentration,

The directional concentration degree updated in the previous update step is substituted into the modified Bessel function of the first kind having a rank different by one, and the detailed ratio is updated by the formula (97), and updated with the updated detailed ratio. Substituting the directional concentration degree into the Newton method update formula (equation (98)) and updating the directional concentration degree by alternately repeating the calculation process using the equations (97) and (98) several times. By
characterized by calculating the degree of directional concentration,
The method for calculating the degree of concentration of directions in a clustering operation.

Images