WO2004079388A1 - 位置情報推定装置、その方法、及びプログラム - Google Patents

位置情報推定装置、その方法、及びプログラム Download PDF

Info

Publication number
WO2004079388A1
WO2004079388A1 PCT/JP2004/002610 JP2004002610W WO2004079388A1 WO 2004079388 A1 WO2004079388 A1 WO 2004079388A1 JP 2004002610 W JP2004002610 W JP 2004002610W WO 2004079388 A1 WO2004079388 A1 WO 2004079388A1
Authority
WO
WIPO (PCT)
Prior art keywords
matrix
position information
signal
signal source
frequency
Prior art date
Application number
PCT/JP2004/002610
Other languages
English (en)
French (fr)
Inventor
Hiroshi Sawada
Ryo Mukai
Shoko Araki
Shoji Makino
Original Assignee
Nippon Telegraph And Telephone Corporation
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Nippon Telegraph And Telephone Corporation filed Critical Nippon Telegraph And Telephone Corporation
Priority to US10/508,708 priority Critical patent/US7039546B2/en
Priority to EP04716703A priority patent/EP1600789B1/en
Priority to DE602004029867T priority patent/DE602004029867D1/de
Priority to JP2005503053A priority patent/JP3881367B2/ja
Publication of WO2004079388A1 publication Critical patent/WO2004079388A1/ja

Links

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S3/00—Direction-finders for determining the direction from which infrasonic, sonic, ultrasonic, or electromagnetic waves, or particle emission, not having a directional significance, are being received
    • G01S3/02—Direction-finders for determining the direction from which infrasonic, sonic, ultrasonic, or electromagnetic waves, or particle emission, not having a directional significance, are being received using radio waves
    • G01S3/14—Systems for determining direction or deviation from predetermined direction
    • G01S3/46—Systems for determining direction or deviation from predetermined direction using antennas spaced apart and measuring phase or time difference between signals therefrom, i.e. path-difference systems
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S3/00—Direction-finders for determining the direction from which infrasonic, sonic, ultrasonic, or electromagnetic waves, or particle emission, not having a directional significance, are being received
    • G01S3/02—Direction-finders for determining the direction from which infrasonic, sonic, ultrasonic, or electromagnetic waves, or particle emission, not having a directional significance, are being received using radio waves
    • G01S3/74—Multi-channel systems specially adapted for direction-finding, i.e. having a single antenna system capable of giving simultaneous indications of the directions of different signals
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S3/00—Direction-finders for determining the direction from which infrasonic, sonic, ultrasonic, or electromagnetic waves, or particle emission, not having a directional significance, are being received
    • G01S3/80—Direction-finders for determining the direction from which infrasonic, sonic, ultrasonic, or electromagnetic waves, or particle emission, not having a directional significance, are being received using ultrasonic, sonic or infrasonic waves
    • G01S3/8006—Multi-channel systems specially adapted for direction-finding, i.e. having a single aerial system capable of giving simultaneous indications of the directions of different signals
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S3/00—Direction-finders for determining the direction from which infrasonic, sonic, ultrasonic, or electromagnetic waves, or particle emission, not having a directional significance, are being received
    • G01S3/80—Direction-finders for determining the direction from which infrasonic, sonic, ultrasonic, or electromagnetic waves, or particle emission, not having a directional significance, are being received using ultrasonic, sonic or infrasonic waves
    • G01S3/802—Systems for determining direction or deviation from predetermined direction
    • G01S3/808—Systems for determining direction or deviation from predetermined direction using transducers spaced apart and measuring phase or time difference between signals therefrom, i.e. path-difference systems

Definitions

  • the present invention observes signals radiated from a plurality of signal sources such as sound sources and radio sources and mixed with each other with a plurality of sensors, and obtains positional information of each signal source, at least a parameter representing the position.
  • the present invention relates to an apparatus, a method, and a program for estimating information including one, detecting the direction of arrival of a signal, and applying the separation and restoration to signals for each signal source. Background art
  • the direction of arrival of a signal is often estimated in the frequency domain.
  • This includes observation signals Signal) ”(t) is Fourier transformed in a short time to obtain a time-series signal ⁇ ” ( ⁇ , m) in the frequency domain.
  • Direction of arrival vector with direction as a variable
  • ⁇ ( ⁇ , m) W ( ⁇ ) X ( ⁇ , m)
  • ⁇ ( ⁇ , m) ⁇ , ⁇ , m), ⁇ , ⁇ ( ⁇ , m)] ⁇ ⁇ ⁇ ⁇ (6)
  • Wi (o>) produces the separated signal Yi (, m). Therefore, Wi (o>) emphasizes and extracts one of the source signals S! W, m),..., S, (ft) ( m) and suppresses the others.
  • Wi ((W)) By analyzing the directional characteristics formed by Wi ((W)), it is possible to extract the signal arriving from which direction and to suppress in which direction the signal arriving is suppressed. The direction of arrival of the signal st) can be estimated.
  • This ( ⁇ , ⁇ ) is considered to be the frequency response from the source signal in direction S to the separated signal ⁇ ; ( ⁇ , m).
  • Figure 2 shows the gain I ⁇ ; ( ⁇ ,) I of the directional characteristic after independent component analysis at a frequency of 3156 Hz.
  • the horizontal axis represents ⁇
  • the vertical axis represents gain.
  • the solid line is the directional characteristic IB ft), ⁇ ) I given by the first row of the separation matrix, and the dashed line is the directional characteristic IB 2 ( ⁇ , I) given by the second row.
  • FIG. 3 shows the gain IBw, ⁇ ) I of the directional characteristic at the frequency of 2734Hz after ICA in this case.
  • the solid line in Fig. 3 is the directivity given by the first row of the separation matrix, the broken line is the second row, and the dashed line is the directivity given by the third row is there.
  • each source signal is emphasized by one row of the separation matrix and suppressed by the other two rows.
  • the two lines to be suppressed are not necessarily suppressed in the same direction.
  • i B 2 I and IB 3 I are both minimum around 45 °, w! W) takes out the source signal around 45 °, and w 2 (w) and w 3 (cw) are This source signal is suppressed.
  • I ⁇ , I and IB 3 I are minimized around 90 °, the source signal around 90 ° is extracted at w 2 (o)), and this signal is suppressed at w! Cw) and w 3 (a>). You can read it.
  • the blind signal separation is a technique for estimating the source signal from the observed mixed signal.
  • a mixed signal in which I source signals are mixed is observed by J sensors will be described as an example.
  • a jj is an impulse response from the signal source i to the sensor j
  • ⁇ is a convolution operator.
  • the convolutional mixture in the time domain can be transformed into multiple instantaneous mixtures in the frequency domain. That is, the above equations (7) and (8) are expressed by the above equations (1) and (5), respectively.
  • W ( ⁇ ) is a separation matrix, and This is the solution that Y k ( ⁇ , m) and Y k ′ ( ⁇ , m) are calculated to be independent from each other, that is, the solution of ICA.
  • the problems in performing the above-described blind signal separation in the frequency domain are the problem of permutation and the problem of scaling.
  • the separation matrix W ( ⁇ ) is a solution of the independent component analysis even if the rows are exchanged.
  • W ( ⁇ ) is a solution of IC ⁇
  • any arbitrary mutation matrix this matrix is any matrix ⁇ ( ⁇ ) D ( ⁇ ) W ( ⁇ ) is also the solution of ICA when ⁇ ⁇ ( ⁇ ) is replaced by ⁇ ( ⁇ ).
  • This is based on the fact that ICA separates source signals only under the condition of statistical independence between source signals, and the degree of freedom of the solution given by D ( ⁇ ) is called scaling degree of freedom.
  • the degree of freedom of the solution given by ⁇ ( ⁇ ) is called the permutation degree of freedom.
  • an appropriate solution W ( ⁇ ) must be specified as a separation matrix from the solution of I C ⁇ .
  • the appropriate solution W ( ⁇ ) is specified by multiplying the arbitrarily determined solution of IC ⁇ by appropriate D ( ⁇ ) and ⁇ ( ⁇ ), and defining the result as the appropriate solution W ( ⁇ ). This is done by:
  • the problem of appropriately determining D ( ⁇ ) is called the scaling problem, and the problem of properly determining ⁇ ( ⁇ ) is called the permutation problem.
  • Permutation is a bijective function from ⁇ 1,2, ⁇ , 1 ⁇ to ⁇ 1,2, ⁇ , 1 ⁇ ⁇ : ⁇ 1,2,-, 1 ⁇ ⁇ ⁇ 1,2,-, 1 ⁇ , which has a one-to-one correspondence with the permutation matrix.
  • the degree of freedom of scaling corresponds to the degree of freedom of the filter that changes the frequency characteristics in the time domain. Therefore, in order to generate a distortion-free separated signal in the time domain, D must be determined appropriately for all ⁇ .
  • the separated signal ( ⁇ , The separated signal ( ⁇ 2 , m) may be output as an estimated value of the source signal S 2 ( ⁇ 2 , m) for another angular frequency ⁇ 2 .
  • the time domain source signal s, (t) and the component of s 2 (t) are mixed in the time domain output signal y, (t), and the separated signal is correctly extracted. Cannot be generated.
  • the estimation accuracy of the direction of arrival of a signal is estimated by searching for and estimating the direction in which the gain of the directional characteristic is low Differs depending on the position of the signal source, and the error increases especially when the signal arrival direction is close to the straight line connecting the pair of sensors 1 "and 1".
  • the microphones are placed at a distance of 2.83 cm as the two sensors 101 and 102, and a fixed distance (about 1 cm) from the center (original) of sensors 101 and 102 is obtained.
  • Two signal sources 1 1 1 and 1 1 2 apart from each other at an angle interval of 20 ° are provided with sound sources, and the origin is determined with reference to the direction from sensor 101 to sensor 102 (0 °).
  • the sound sources 111 and 112 were moved while maintaining the constant distance and angular interval until the angle S of the signal source 111 viewed from 10 ° to 150 °.
  • FIG. 4B shows an example of the result of the blind signal separation performed while moving the sound sources 111, 112 as described above.
  • the vertical axis in Fig. 4B shows the signal to interference ratio (Signal to Interference Ratio).
  • the separation matrix is obtained by using the independent component analysis, the directional pattern is obtained from each row, and the direction of the signal source (the direction of the incoming signal) is obtained by searching for the direction with the low gain.
  • the direction of the signal source the direction of the incoming signal
  • An object of the present invention is to provide a position information estimating apparatus, a method and a program thereof, which require a short calculation time for estimating position information of a signal source by obtaining a separation matrix by independent component analysis. Disclosure of the invention
  • ⁇ ( ⁇ )) of the mixing matrix A ( ⁇ ⁇ ), '', ⁇ ( ⁇ ) including the degrees of freedom, and ⁇ ⁇ ( ⁇ ⁇ ) ( ⁇ 1,-, ⁇ )
  • i is a parameter representing the signal source
  • One parameter of the position information of the signal source i for example, a conical surface or a curved surface where the signal source exists is calculated.
  • FIG. 1 is a diagram for explaining the relationship between the sensor array and the arrival time difference between the arrival signal and each sensor.
  • Fig. 2 is a diagram showing the gain directional characteristics of each row of the separation matrix calculated by ICA for the two sound source mixed signals.
  • Fig. 3 is a diagram showing the gain directivity of each row of the separation matrix calculated by ICA for the three sound source mixed signals.
  • FIG. 4A is a diagram for explaining the relationship between the sensor and the signal source used in the preliminary experiment.
  • FIG. 4B shows the results of the preliminary experiment.
  • FIG. 4C is a diagram showing the relationship between the estimated angle and its error and declination error.
  • FIG. 5 is a block diagram showing an example of a functional configuration of the first embodiment in which the present invention is applied to estimation of a signal arrival direction.
  • FIG. 6 is a flowchart illustrating an example of a processing procedure according to the first embodiment.
  • FIG. 7 is a block diagram showing an example of a specific functional configuration of the angle calculation unit in FIG.
  • FIG. 8 is a flowchart showing an example of a specific processing procedure of step S4 in FIG.
  • FIG. 9 is a diagram showing an experimental result of the direction estimation according to the first embodiment.
  • Fig. 1 OA is a diagram showing experimental results of direction estimation according to the first embodiment when there are two sound sources and three microphones.
  • FIG. 10B is a diagram showing an experimental result obtained by estimating the direction by the MU SIC method under the same conditions as the experiment for FIG. 1OA.
  • FIG. 11A is a diagram showing experimental results obtained by performing direction estimation according to the first embodiment when there are two sound sources and three microphones.
  • FIG. 11B is a diagram showing an experimental result obtained by estimating the direction by the MUSIC method under the same conditions as the experiment for FIG. 11A.
  • FIG. 12 is a block diagram showing a functional configuration example of a second embodiment in which the present invention is applied to blind signal separation.
  • FIG. 13 is a flowchart illustrating an example of a processing procedure according to the second embodiment.
  • FIG. 14 is a block diagram showing an example of a specific processing procedure of step S14 in FIG.
  • FIG. 15 is a diagram for explaining a common linear direction of a plurality of conical surfaces.
  • FIG. 16 is a block diagram showing another example of a specific functional configuration of the arrival direction determining unit 16 in FIG.
  • FIG. 17 is a block diagram showing a functional configuration example of a third embodiment in which the present invention is applied to blind signal separation.
  • FIG. 18 is a flowchart illustrating an example of a processing procedure according to the third embodiment.
  • FIG. 19 is a flowchart showing an example of a specific processing procedure of step S35 in FIG.
  • FIG. 20 is a diagram for explaining the relationship between the sensor arrangement, the signal source position, and the estimated curved surface.
  • FIG. 21 is a diagram illustrating an example of an estimated spherical surface.
  • FIG. 22 is a diagram showing the relationship between the room, microphone arrangement, and sound source used in the experiment.
  • FIG. 23 is a histogram of the estimation direction using a small-interval microphone pair.
  • FIG. 24 is a diagram showing the distribution of the estimation direction with respect to the frequency.
  • Figure 2 5 shows a distribution to the radius of the frequency of the estimated spherical for the sound source 2 4 and 2 5.
  • FIG. 26 is a diagram showing the experimental results of each method.
  • FIG. 27 is a block diagram showing a functional configuration example of a main part of the fourth embodiment of the present invention.
  • FIG. 28 is a flowchart showing an example of a processing procedure according to the fifth embodiment in which the present invention is applied to blind signal separation.
  • FIG. 29 is a view showing an experimental result for the fifth embodiment.
  • Figure 3 OA is a flowchart showing the main part of the processing procedure for finding the common linear direction of a plurality of conical surfaces and solving the permutation of the separation matrix.
  • FIG. 30B is a flowchart showing a main part of a processing procedure for solving the permutation of the separation matrix by using the conical surface estimation and the spherical surface estimation.
  • FIG. 30C is a flowchart showing a main part of another processing procedure for solving the permutation of the separation matrix by using the conical surface estimation and the spherical surface estimation.
  • the direction of a signal source that is, the arrival direction of a source signal radiated from the signal source is obtained.
  • FIG. 5 shows a functional configuration diagram of the first embodiment
  • FIG. 6 shows a flowchart of a part of the processing procedure.
  • ⁇ ( ⁇ ) ( ⁇ ( ⁇ ( ), ⁇ ( ⁇ 2 ), ..., ⁇ ( ⁇ ⁇ ))
  • this inverse matrix is calculated as a pseudo-inverse when J> I.
  • a Moore-Penrose generalized inverse matrix is used as the pseudo inverse matrix.
  • Two elements of each row i of the H (omega) at least inverse of one frequency in ⁇ in this embodiment ( ⁇ ⁇ ) ( ⁇ ⁇ ) and ⁇ ".
  • Polarization of the ratio of the i (omega eta) is calculated by the angle calculator 14 (FIG. 6, step S4).
  • One unselected column i in the inverse matrix ⁇ ( ⁇ ⁇ ) of the angular frequency ⁇ ⁇ is selected by the selecting unit 14a (FIG. 8, step S4a), and two elements ⁇ ⁇ ”from the i-th column; ( ⁇ ⁇ ) and ⁇ ; ( ⁇ ⁇ ) are selected (FIG. 8, step S 4 b).
  • the argument calculating unit 14 b calculates the argument arg [Hji ( ⁇ ⁇ ) / H r ; ( ⁇ ⁇ )] of the ratio of the selected element, and the interval calculating unit 14 stores sensor information.
  • the phase rotation calculator 14 d calculates the product of the distance d j-dj. And the amount of phase rotation of the signal per unit distance at the angular frequency ⁇ ⁇ o) n / c (c is the speed of the source signal).
  • phase rotation amount of the signal at the interval dj—dj. Is calculated, and the argument arg [ ⁇ ( ⁇ ⁇ ) ⁇ ”.i ( ⁇ ⁇ )] is divided by the division unit 14 e using the phase rotation amount. ( Figure 8, step S4d).
  • ⁇ i ( ( y n ) cos " 1 G i ( ⁇ ): For I Gi ((W n ) II, ⁇ ' ⁇ (9) I Gi ( ⁇ jo)
  • ⁇ , ( ⁇ ⁇ ), ⁇ ⁇ ( ⁇ ⁇ ), ⁇ , ⁇ ! ( ⁇ ⁇ ) are the source signals S l (t), s 2 (t), ⁇ , s J (t) Corresponding to one of the directions of arrival (directions of signal sources 1, 2,..., I).
  • D ( ⁇ ) is a diagonal matrix indicating the degree of freedom of scaling
  • ⁇ ( ⁇ ) is a permutation matrix indicating the degree of freedom of permutation
  • I is a unit matrix.
  • the mixed row including the degrees of freedom of scaling and permutation A column estimate is obtained.
  • the inverse matrix ⁇ ( ⁇ ) is obtained by rearranging the columns of the mixing matrix ⁇ ( ⁇ ) by ⁇ ( ⁇ ), and further multiplying each column by the diagonal element of D ( ⁇ ).
  • ⁇ is a permutation corresponding to multiplying the permutation matrix ⁇ ( ⁇ ) from the right.
  • the estimated value of the mixing matrix ⁇ ( ⁇ ) including the degrees of freedom of scaling and permutation is calculated using the inverse matrix ⁇ ( ⁇ ) of the separation matrix W ( ⁇ ).
  • the angle Si ( ⁇ ⁇ ) may be obtained for each column, and the direction of arrival may be determined from the total of these angles. That is, the arrival direction of each frequency estimated by the angle calculation unit 14 is divided for each angle by the sorting unit 32 in FIG. 5 (FIG. 6, step S5). For example, they are arranged in descending order of direction (angle).
  • Angular frequency In the order each component in the ⁇ ( ⁇ ) is larger for (( ⁇ ⁇ ), ⁇ 2 ((Wj), ⁇ , ⁇ ! ( ⁇ ,)) and side by side, ⁇ for the ⁇ 2 ( ⁇ 2) in ( ⁇ 2 ), ⁇ ⁇ ( ⁇ 2 ), ⁇ , ⁇ !
  • the classification angle ( ⁇ ,), ⁇ ( ⁇ 2 ),..., ⁇ ⁇ ( ⁇ ⁇ ) is integrated into one angle i and one angle (direction of arrival); 6, Step S 6).
  • the average value of each separation angle (( ⁇ ⁇ ), ⁇ i ( ⁇ 2 ), ⁇ , ⁇ ⁇ ( ⁇ ⁇ )) is taken as the integration angle, or the classification angle (Si ( ⁇ ⁇ ), ⁇ ⁇ ( ⁇ 2 ), ⁇ , ⁇ ⁇ ( ⁇ )), the mode, the median, etc. may be set as the integrated angle 0 ;.
  • the direction of arrival can be estimated more correctly than when the direction of arrival is estimated only from the inverse matrix H ( ⁇ ) of one certain frequency.
  • FIG. 1OA shows the results obtained by the method of this embodiment
  • FIG. 10B shows the results obtained by the MUSIC method. From these results, it can be seen that both methods estimate the direction fairly correctly.
  • the results of the fractional averaging of the results according to this embodiment in two directions are 45 ° and 123 °, and the MUSIC method is 45 ° and 122 °, respectively.
  • the same experiment was performed for the case where two sound sources were placed with the sound source directions of 105 ° and 119 ° relatively close to each other.
  • FIGS. 11A and 11B The results obtained by the method of this embodiment and the results obtained by the MU SIC method are shown in FIGS. 11A and 11B, respectively.
  • the MUS IC method cannot estimate the sound source direction at most frequencies, in this embodiment, it may be possible to perform angle calculation at most frequencies. ° and 124 °, and the MUS IC method is 94 ° and 128 °, indicating that the direction estimation according to this embodiment is fairly correct.
  • the estimation direction can be obtained only by substituting the value into Equation (9), compared with the case of searching in the direction where the gain of the directivity pattern is low.
  • the separation matrix W ( ⁇ ) obtained by ICA includes the degrees of freedom of scaling and the degrees of freedom of permutation.Therefore, the inverse matrix of the separation matrix W '( ⁇ ) that solves these is calculated. In other words, it may be possible to calculate the true mixing matrix ⁇ ( ⁇ ) and estimate the direction of arrival using the ratio of the two elements for each column of the mixing matrix ⁇ ( ⁇ ).
  • the true mixture matrix ⁇ ( ⁇ ) itself is, for example, the average power of the signal S i (t) of the signal source. It cannot be done unless constraints such as one are set to one. In the wireless communication field, such a restriction may be imposed on the signal source.For example, when the signal S i (t) of the signal source is a voice signal directly uttered by a human, such a restriction may be imposed. Is virtually impossible.
  • the freedom of scaling is obtained by taking the ratio of two elements for each column of the inverse matrix ⁇ ( ⁇ ) of the separation matrix W ( ⁇ ) including the degrees of freedom of scaling and permutation.
  • This method solves the above problem and can be applied to any signal source, and furthermore, it is not necessary to calculate a separation matrix that solves the above two problems, and the calculation time is correspondingly short. Furthermore, as described above, the permutation problem can be easily solved by separating the estimated directions obtained for each frequency in a predetermined order. Even if the number is the same as the number J of sensors, the direction of each signal source can be estimated. Even if the sound source directions are relatively close, it is possible to estimate fairly accurately.
  • direction information is obtained as one piece of position information of a signal source. At least three sensors arranged two-dimensionally are used to determine the direction of the signal source in any direction.
  • the direction information is determined by estimating the conical surface based on the direction information and estimating the common straight line of the plurality of conical surfaces.
  • FIG. 12 shows a functional configuration in which the second embodiment is applied to a blind signal separation device
  • FIG. 13 shows a processing procedure thereof.
  • four sensors 1 lt 1 2, 1 3, 1 4 are arranged at equal intervals on the same circle, the interval of any two sensors, is less than half the shortest wavelength of the original signal.
  • the number of sensors is J, and it is assumed that J ⁇ 3.
  • the inverse matrix calculator 13 calculates the inverse matrix ⁇ of each separation matrix W ( ⁇ ) for each frequency (step S13).
  • the conical surface estimating unit 14 uses two sensor pairs corresponding to these elements from the ratio of each element for a plurality of element pairs different for each column in each inverse matrix H ( ⁇ ) for each frequency.
  • the sensor pair axis as the central axis, a conical surface on which one of the signal sources exists, that is, a plurality of conical surfaces is estimated for each column of one mixing matrix ⁇ ( ⁇ ) (step S14) .
  • the functional configuration of the conical surface estimation unit 14 is almost the same as that of the angle calculation unit 14 in FIG. 15, and the processing procedure is similar to that of FIG.
  • a specific example of the processing of step S14 for the inverse matrix ⁇ ( ⁇ ) of a certain frequency performed in the conical surface estimation unit 14 will be described with reference to FIG.
  • control parameters i and ⁇ stored in a register in the conical surface estimation unit 14 are initialized to 0 (step S 20).
  • i corresponds to the number of each signal source
  • p is the number of estimated conical surfaces for each i.
  • this selection is based on the fact that the sensor pair axis (the straight line passing through the sensor j and the sensor j ') specified by the selected j, j' is determined by the j, j 'selected earlier in this routine processing. It is desirable that the conical surface estimating section 14 estimates a plurality of conical surfaces whose central axes do not overlap within a certain error range within a certain error range.
  • the determination as to whether or not each sensor pair axis is on the same straight line may be made, for example, by storing a vector indicating the position of each sensor in the sensor information storage unit 15 and starting from each sensor. This is done by extracting the vector information that indicates the position of.
  • the vector dj. Indicating the position is extracted (step S24).
  • j-th row and i-column element Hji (a>) and j'-row and i-column element H ⁇ ⁇ ) are selected and extracted from the mixing matrix H ( ⁇ ) (step S25).
  • the above processing is performed by the selection unit 14a in FIG. Therefore, a register for storing i, p, j, j 'and the selected j, j', etc. is also provided in the selection unit 14a.
  • Equation (9) is for a case where the sensors are arranged in one dimension and the angle of arrival of the signal is two-dimensional. This is extended when the sensors can be arranged in two or three dimensions and the angle of the signal arrival direction can be in three-dimensional space. Therefore, equation (9 ') includes equation (9).
  • the angle 0 ⁇ u ( ⁇ ) estimated by the equation (9 ′) and its i, j, j ′ are temporarily recorded as cone surface information in a register (memory) in the cone surface estimation unit 14 (step S 27). As shown by the broken line in FIG.
  • step S26 is the argument calculation unit 14b, interval calculation unit 14c, phase rotation calculation unit 14d, division unit 14e, determination unit 14f, and inverse cosine calculation unit in FIG. According to 14 g, steps S 4 c, S 4 d, S 4 e and
  • This P is the number of conical surfaces to be estimated for each i.
  • Sensor 1 3 is arranged in the arrangement direction perpendicular to the direction of the sensor 1 i and 1 2, between sensor 1 i and 1 2, the distance between the sensor 1 2 and 1 3 are the same.
  • i 2 of the signal source 2 2
  • the conical surface 4 12 around axis sensor pair axis 3 12 with sensor 1, and 1 2 ( ⁇ ⁇ 2, 12 ( ⁇ )) is, the sensor 1 2 and 1 3 conical surface 4 23 around axis sensor pair axis 3 23 ( ⁇ 2.
  • 23 ( ⁇ )) is a conical surface 4 13 around axis sensor pair shaft 313 by the sensor 11 and 1 3 ( ⁇ ⁇ 2, 13 ( ⁇ )) are estimated respectively.
  • the conical surface 4 12 and 4 13 are considered to be in line contact with each other at a common straight line 5 2.
  • This common linear 5 2 towards direction u 2 (omega) determines the direction of the signal source 2 arrival direction u 2 of 2 radiated signals from clogging the signal source 2 2.
  • substantially linear contact Suruga with conical surface 4 12, 4 13 if also translate conical surface 4 23, conical surface 4 23 a fourth embodiment described later will be discarded.
  • the normalized axis vector calculator 16a normalizes the axis vector (dj (p) —d (p)) connecting the positions of the pair of sensors to a length of 1, respectively. I mean
  • this system of equations has no solution or is not uniquely determined. Therefore, u that minimizes II Vu ⁇ c ⁇ ( ⁇ ) II is obtained and is set as the solution of the simultaneous equations, that is, the direction of the common line 5 i ⁇ ( ⁇ ). The calculation for u to minimize this error is performed in the arithmetic unit 16b. Since this direction Ui ( ⁇ ) is a three-dimensional direction, the azimuth ( ⁇ ) and elevation ⁇ i> i ( ⁇ ) can be obtained by polar coordinates.
  • sorting is performed on the arrival azimuth angle ⁇ as follows, and the unresolved columns are similarly sorted by the arrival angle 0i.
  • the columns of the inverse matrix ⁇ ( ⁇ ) are exchanged so that the calculated arrival azimuth ( ⁇ ) is in a predetermined order at any frequency, for example, 0 or 0 2 ,..., ⁇ 1 in ascending order.
  • the permutation matrix that replaces the columns of the inverse matrix ⁇ ( ⁇ ) is generated by the permutation matrix generation unit 17a so that the elevation angle becomes ascending order at any frequency for the columns for which the permutation was not performed.
  • the matrix P ( ⁇ ) is generated by the inverse matrix generator 17, and the rearranger 17 c multiplies the inverse matrix ⁇ ( ⁇ ) and the separation matrix W ( ⁇ ) from the left. Is calculated.
  • the permutation matrix generation unit and the inverse matrix generation unit 17b constitute a permutation matrix generation unit.
  • the separation matrix W ′ ( ⁇ ) is converted to a separation filter coefficient group in the time domain by, for example, an inverse Fourier transform in the time domain transformation unit 18.
  • the separation filter coefficient group is set in the signal separation unit 19.
  • the signal separation unit 19 calculates the equation (8) using the observation signals (t),-, X j (t) from each sensor and the separation filter coefficient group, and performs the separation signal (t),-,, ⁇ (t) is output.
  • Time domain signal (t), ..., y! (t) may be generated.
  • the sorting unit 32 in FIG. 5 sorts out each direction range.
  • the integration may be performed by the integration unit 33.
  • the conical surface information (. ( ⁇ ), j ( ⁇ )
  • the amount of calculation is small.
  • a plurality of conical surfaces are estimated and the direction of arrival of the signal is calculated from the common straight line. From 360. No matter in which direction the signal source exists, the direction of signal source 3 can be uniquely estimated.
  • the prediction matrix P ( ⁇ ) since the estimated direction is used to determine the prediction matrix P ( ⁇ ), the prediction problem can be properly solved regardless of the position of the signal source.
  • the third embodiment uses a curved surface on which a signal source exists based on a ratio of distances between a pair of sensors and one signal source as position information.
  • Each signal source was located far from the sensor, and the signal from the signal source was processed assuming that it arrived at the sensor as a plane wave. However, when the distance between the signal source and the sensor is relatively short, the signal arrives at the sensor as a spherical wave. From these points, the ratio of the elements of the mixing matrix A ( ⁇ ) ⁇ ⁇ ) ⁇ ”. ⁇ ) is interpreted by a spherical wave (near-field) model, and it is possible to estimate information other than the direction of the signal source i. it can.
  • the frequency response ⁇ ( ⁇ ) can be expressed as follows.
  • qi is a vector indicating the position of the signal source i.
  • the permutation problem with respect to the separation matrix is solved by using the position information indicating the curved surface on which the signal source exists and the directional information.
  • FIG. 17 shows a functional configuration of the third embodiment
  • FIG. 18 shows a processing procedure thereof.
  • Sensor spacing of sensors 1 2 and 1 3 to distance sensor 1 i and 1 2 as shown in FIG. 2 0 In example embodiment
  • three or more are arranged in a two-dimensional or three-dimensional 1 It is 0 to 20 times, preferably about 15 times.
  • each observation signal is 0 to 20 times, preferably about 15 times.
  • step S 11 A separation matrix W ( ⁇ ) is calculated for each (step S12), and a matrix ⁇ ( ⁇ ) is calculated as an inverse matrix of each separation matrix W ( ⁇ ) (step S13).
  • step S14 one or preferably a plurality of conical surfaces are estimated using element pairs selected from each column of the inverse matrix H ( ⁇ ) for each frequency (step S14).
  • the distance ratio calculator 31 uses an element pair selected for each column from the inverse matrix ⁇ ( ⁇ ) for each frequency, and uses the element pair selected for each column to determine the distance between the corresponding sensor and one signal source i.
  • the following equation (10 ′) is calculated from the ratio, that is, equation (10) and equation (21) (step S35).
  • step S35 A specific example of the processing in step S35 performed in the distance ratio calculation unit 31 will be described with reference to FIG.
  • the parameter i is initialized to 0 (step S20), then i is incremented by 1 (step S21), and natural numbers less than J are randomly selected as parameters j and j ', for example, j is selected. ⁇ , and the pair selected once is not selected (step S23).
  • the position vector dj of the sensor j and the position vector dj of the sensor j ' are extracted (step S24), and from the i-th column of the inverse matrix H ( ⁇ ), the elements Hji w) and H r And are selected (step S25).
  • the ratio DRi, j ( ⁇ ) of these two selected elements is calculated (Step S41).
  • the distance ratio information DRi, jr ( ⁇ ) calculated by the distance ratio calculator 31 is sent to the permutation resolver 17, and the permutation resolver 17 is estimated by the arrival direction determiner 16.
  • the separation matrix calculated by the separation matrix calculation unit 12 is performed. Solve the simulation problem.
  • the permutation problem is solved by replacing the line of W ( ⁇ ).
  • the distance estimator 17d calculates the distance II qi ( ⁇ ) II of the signal source 2i from the direction information and the distance ratio information (step S36).
  • the permutation resolution unit 17 determines the separation matrix W ( ⁇ ) such that the signal source distance II qi ( ⁇ ) II obtained for each frequency is given in a predetermined order, for example, in ascending order. The line of is replaced. That is, the permutation matrix ⁇ ( ⁇ ) is calculated to perform this replacement (step S37).
  • the calculation of the permutation matrix ⁇ ( ⁇ ) may be performed in the same manner as the calculation of the permutation matrix performed by the permutation solving unit in the second embodiment.
  • the output separation matrix W '( ⁇ ) is sent to the time domain transform section 18 and used for signal separation in the time domain transform section 18.
  • the permutation solving unit 17 uses the distance ratio DRi, ( ⁇ ) calculated as shown in parentheses in step S37 in FIG. 18 to calculate the inverse matrix ⁇ ( ⁇ ) at all frequencies.
  • a permutation matrix is generated such that the distance ratio DR jT ( ⁇ ) to be obtained is in ascending order in the order of the columns, the second column, the third system U,..., the I column, and the permutation matrix ⁇ ( ⁇ ) May be generated.
  • the distance estimating unit 17 d in FIG. 17 is omitted.
  • the direction information Ui ( ⁇ ) cannot be calculated by the direction-of-arrival determination unit 16 or only the signal source i and the sensor ij ′ for which the same Ui ( ⁇ ) is calculated for a plurality of columns of the inverse matrix ⁇ ( ⁇ ),
  • the distance ratio DRi,, ( ⁇ ) may be calculated, and the rows of the separation matrix W ( ⁇ ) may be exchanged by themselves or by the obtained II q t ( ⁇ ) II. That is, the rows of the separation matrix W ( ⁇ ) are rearranged using the direction information Ui ( ⁇ ), and then, for rows that cannot be rearranged by Uj ( ⁇ ), the distance ratio DRi, ( ⁇ ) or II ( ⁇ )
  • a permutation matrix ⁇ ( ⁇ ) that performs both exchanges simultaneously is generated.
  • the distance ratio DRi, j ( ⁇ ) is obtained, the rows of the separation matrix W ( ⁇ ) are replaced, and the direction information Ui ( ⁇ ) is used for those for which the distance ratio DRi,, ( ⁇ ) could not be obtained. Further, the separation matrix W ( ⁇ ) may be replaced. In this case as well, a permutation matrix ⁇ ( ⁇ ) that performs both exchanges simultaneously is generated.
  • Ui ( ⁇ ) can be obtained with higher precision than DRi, j ( ⁇ ). Therefore, it is better to replace mainly by Ui ( ⁇ ) and to replace those that cannot be performed by DRi, jr ( ⁇ ).
  • Equation (10 ') is expressed as the distance ratio DRi, j ( ⁇ ) on the right-hand side, and solving for cii, the center, ( ⁇ ) and the radius Ri, j ( ⁇ ) are The spherical surface is given by (1 1) and (1 2).
  • the state of the spherical surface determined by Eq. (9) is as shown in Fig. 21. That is, the signal source i exists on the spherical surface given by the equations (11) and (12) as position information. Therefore, the distance estimating unit 1 ⁇ d in the permutation solving unit 17 in Fig. 17 is a curved surface estimating unit as shown in parentheses, and each of i and jj for which the direction information Ui ( ⁇ ) could not be obtained.
  • the surface estimator 17 d calculates the radius Ri, j ( ⁇ ) and the center Oi, ( ⁇ ) of the equation (11) as shown in parentheses in step S36 in FIG. each calculated, these Ri, j (omega) and the center Oi, and the j (omega) the inverse matrix of which frequencies Eta (omega) even be the same order, be determined by Pamyute one sucrose emissions matrix [rho (omega) Good. .
  • the correlation calculation unit 34 calculates the correlation again for the fank component in (t), and repeats the same until the calculated correlation becomes maximum. If the maximum value of the correlation does not reach the predetermined value, furthermore , the component of the fank and the fundamental or harmonic (generally, the harmonic ) for the fank in the frequency component for which the permutation matrix is obtained. The calculation of the sum of the correlations with those having a relationship is performed by the correlation calculation unit 34, and the correction unit 17e performs the permutation of the separation matrix with respect to fank so that the correlation sum becomes large. Repeat until the correlation sum is maximum return. As described in the second embodiment, when the separation matrix W ( ⁇ ) is applied to the frequency domain signal X ( ⁇ , m) to obtain the frequency domain separated signal ⁇ ( ⁇ , m), ( ⁇ , m) may be used for the calculation of the correlation calculator 34.
  • an omnidirectional microphone li ⁇ as a sensor is located at a height of about 135 cm at the center in a horizontal plane with a reverberation time of 130 ms, 355 cm x 445 cm x 250 cm (height).
  • the center of the microphone port Hong disposed as an origin, a 1 20 cm away from one 30 ° direction to the sound source 2i from the origin, + 30 ° direction to the sound source 2 2 a + a 90 ° direction to the sound source 2 3 one 1 respectively arranged sound source 2 6 50 ° direction, the sound source 2 4 to the distance 60 cm from the origin at 1 50 ° direction, distance 1 80 cm to the tone generator 2 5 were arranged respectively.
  • the sampling frequency was kHz
  • the data length was 6 seconds
  • the frame length was 2048 samples (256 milliseconds)
  • the frame shift was 512 samples (64 milliseconds).
  • FIG. 22 shows the estimated histogram.
  • the horizontal axis in Fig. 23 is the estimation direction, and the vertical axis is the number of signals. From this, there are five groups (clusters) of the estimated direction, of which 150. Nearby clusters are twice as large as other clusters. This is presumed that two of the six sound sources came from the same direction (around 150 °). For the four sound sources, permutations were resolved based on the estimated directions.
  • Figure 24 shows the results. In FIG. 24, the horizontal axis represents frequency, and the vertical axis represents direction (degrees).
  • the incoming signal is determined by the radius of the spherical surface where the sound source is located using the widely spaced microphone pairs 17 and 15 , 17 and 18 , 16 and 15 , and 16 and 18. Were distinguished.
  • Figure The horizontal axis of 25 represents frequency, and the vertical axis represents radius (m).
  • D + C solves permutation by direction (conical surface) estimation, and applies correlation method to those that cannot be solved D + S + C
  • the permutation is solved by the direction (conical surface) estimation and the sphere estimation, and the unresolved frequencies are calculated using the correlation method. According to this latter method, all six sound sources could be separated, and the SIR improvement averaged 17.1 dB from the input SIR.
  • the sensors are two-dimensionally arranged.
  • the spherical surface estimated by the pair of sensors appears symmetrically with the center bisector of these sensors, so that the signal source exists three-dimensionally.
  • the sensor cannot be identified by the two-dimensional arrangement sensor, so the sensor must also be three-dimensionally arranged.
  • the amount of calculation is small because the conical surface information is estimated by Expression (9 ′) and the curved surface information is estimated by Expression (10 ′). Moreover, since the permutation is solved by the conical surface and the distance ratio DRi, ". ( ⁇ ) or the distance II qi ( ⁇ )
  • the reliability of the estimated conical surface is verified, and the conical surface determined to be highly reliable is used for solving permutation of the separation matrix.
  • the conical surface information estimated by the conical surface estimating unit 14 ⁇ (. ( ⁇ ) is determined by the conical surface verifying unit 41 based on the allowable angle information from the allowable angle information storage unit 42. Is verified.
  • the angle S ⁇ ( ⁇ ) obtained by equation (9 ′) is a relative angle with respect to the array direction of the sensors 1 ′′ and 1 ′′ used to obtain this, and as described with reference to FIG.
  • the angle ⁇ i, ( ⁇ ) is around 0 degree or around 180 degrees, it is estimated that correct permutation is not performed.
  • the minimum value ⁇ min and the maximum value ⁇ mx of the angle estimated to be able to perform permutation correctly are stored in the allowable angle information storage unit 42 as allowable angle information, and the conical surface verification unit 4 1 If the estimated conical surface information S ⁇ j ( ⁇ ) is between ⁇ 9 nin and nax, it is determined to be a reliable conical surface and output, that is, used for permutation resolution, but ⁇ min and S nax If and, j ( ⁇ ) is discarded as unreliable-not used for permutation resolution. For example a circular conical surface 4 13 of FIG. 1 5 is discarded.
  • the conical surface information verified to be reliable by the conical surface verification unit 41. ''''. ( ⁇ ) is sent to the arrival direction determination unit 16 in FIG. 12 or as described in FIG. Is sent directly to the permutation resolution unit 17.
  • estimated conical surface Step S 1 00 as indicated by the broken line in the next step S 1 4 ⁇ ⁇ ;, jr ( ⁇ ) verification of whether there is a reliable Will be carried out, and only those that are reliable will be moved to the next step.
  • the unreliable information in the estimated conical surface information is discarded, so that the direction of arrival can be accurately estimated without being adversely affected by this, so that the correct permutation can be obtained.
  • the matrix P ( ⁇ ) can be generated, which means that the SIR (performance) of the separated signal is improved.
  • a distance ratio or spherical information estimated from the distance ratio is used as position information.
  • FIG. FIG. 28 schematically shows the processing procedure of this embodiment.
  • the sensors are arranged at a relatively wide interval, for example, 30 cm when the signal source shown in FIG. 22 is a sound source, and is arranged at least two-dimensionally.
  • the observation signal in the time domain is converted into a signal in the frequency domain as shown in FIG. 17, and the separation matrix W ( ⁇ ) is generated by the separation matrix generation unit 12 (step S 5).
  • the inverse matrix ⁇ (») is generated by the inverse matrix generation unit 13 (step S52).
  • the spherical information is estimated for each column of the inverse matrix H ( ⁇ ) at each frequency (step S53).
  • This spherical surface information is calculated in the same manner as in the method described in the third embodiment.
  • the spherical information is the distance ratio DRi. ( ⁇ ) calculated by the distance calculator 31 or II q calculated by the distance estimator or the curved surface estimator 17 d ; IL or the radius Ri, ( ⁇ ) and the center. Oi. ( ⁇ ).
  • a permutation matrix for the separation matrix W ( ⁇ ) is generated so that the arrangement order is a predetermined order, and the rows of the separation matrix W ( ⁇ ) are replaced (step S54).
  • This processing in the third embodiment is performed by the permutation solution unit 17, however, in this fifth embodiment, only spherical information is used. If there is a frequency for which the permutation cannot be resolved only with this spherical information, the frequency is resolved by the above-described correlation method (step S55).
  • the results show that the method using only the sound source position has poor overall performance, and the method using only the correlation results in variable separation performance. When both are used, stable high performance can be obtained.
  • the threshold value Ath is preferably a relatively large value, and it can be seen that a frequency of about 1/5 to 1/10 of 8 to 16 can be determined by the position of the sound source.
  • the calculation amount is small because the spherical surface information is obtained by the equation (10 ′).
  • the distance ratio DRi. "" ( ⁇ ) is preferable.
  • the sixth embodiment solves the problem of permutation based on, for example, one estimated conical surface.
  • the conical surface ⁇ ⁇ ;, '( ⁇ ⁇ ),... ⁇ ⁇ i, jj. ( ⁇ ⁇ ) Go to Resolution 1 7
  • Figure 30 ⁇ , Figure 30 0, and Figure 30C show the method of generating a separation matrix that resolves the permutation in blind signal separation described above, and the processing after generation of the inverse matrix.
  • Figure 3 OA estimates the conical surface from the element ratio of each column of the inverse matrix (step S61), discards the unreliable conical surface if necessary (step S62), The direction of the common line is determined from the conical surface (step S63), the permutation matrix P ( ⁇ ) is generated using the common line direction, and the separation matrix is replaced (step S64). Rows of the separation matrix are replaced by the correlation method for the frequencies that have not been found (step S65). As shown by the broken line in FIG. 3OA, the process may immediately proceed from step S6 to step S64, and the permutation matrix P ( ⁇ ) may be generated directly using the conical surface.
  • the conical surface is estimated in step S61, the conical surface is discarded in step S62 if necessary, and the common straight line direction is estimated in step S63.
  • the conical surface estimated in step S61 or the conical surface processed in step S62 may be directly used for the processing in step S66.
  • the conical surface could not be estimated or the common linear direction could not be determined, or if the conical surface or the common linear direction corresponded to the uncertain one or had the same value, the element of the column of the inverse matrix was used.
  • the spherical information is estimated from the ratio (step S67), and the inverse is calculated using the estimated spherical information.
  • the permutation matrix P ( ⁇ ) is generated by exchanging the matrices, and the rows of the separation matrix are exchanged (step S64). Furthermore, a method based on correlation is applied to the frequencies for which the permutation matrix ⁇ ( ⁇ ) could not be created (step S65).
  • the spherical information is first estimated from the element ratio of each column of the inverse matrix (step S68), and the columns of the inverse matrix are replaced using this spherical information (step S69), and the replacement of this column is performed.
  • a conical surface was estimated from the element ratio of the column of the inverse matrix (step S70), and the common straight line direction of these multiple conical surfaces was determined.
  • Step S71 the permutation matrix ⁇ ( ⁇ ) is generated by directly using the determined direction or the conical surface estimated in Step S70, and the permutation matrix ⁇ ( ⁇ ) is generated.
  • Exchange is performed (step S64).
  • the correlation method is applied to the frequencies for which a permutation matrix could not be created (step S65).
  • a signal is obtained by using a frequency-domain separation matrix W ( ⁇ ) and an observation signal X ( ⁇ ). After the separation, the separated frequency domain signal ⁇ ( ⁇ ) may be converted to a time domain signal y (t).
  • the direction of arrival of a signal in two dimensions is estimated.
  • the present invention can be applied to the estimation of the direction of arrival of a signal in three dimensions.
  • the second to sixth embodiments can also be applied to two-dimensional signal separation.
  • the estimated conical surface ⁇ 9 is two directions symmetrical with respect to sensor pair axis of the sensor pair used in the estimation, estimation spherical Ri, jj. ( ⁇ ) or DRi. J ( ⁇ ) is the radius of the circle Or it will be almost corresponding to this.
  • the devices shown in FIG. 5, FIG. 12, and FIG. 17, respectively, and the device according to the fifth embodiment can also be caused to function by a computer.
  • the corresponding processing procedure that is, a magnetic disk, a semiconductor memory, a CD ROM, etc., which stores a program for causing a computer to execute the steps shown in FIGS. 6, 13, and 18 etc.
  • the program may be downloaded to a memory in the computer, and the computer may execute the program.

Description

明 細 書 位置情報推定装置、 その方法、 及びプログラム 技術分野
この発明は音源や電波源などの信号源の複数から放射され、 互いに混合された 信号を複数のセンサで観測して、 各信号源の位置情報、 正しくは位置を表わすパ ラメ一タの少くとも 1つを含む情報を推定し、 信号の到来方向の検出や、 信号源 ごとの信号に分離復元に適用される装置、 その方法、 及びプログラムに関する。 背景技術
複数の信号源からの各信号が空間内で混合されて複数のセンサに到来し、 これ らセンサで観測された到来信号から、 各源信号の到来方向の推定や各源信号を分 離することを、その源信号の混合系の情報を知らずに、独立成分分析（Independent Component Analysis、 以後、 I C Aと記述する場合もある） により行うことが提 案されている。 空間内での前記混合は信号源からセンサまでの到達遅延及び減衰 度が直接波と障害物などによる複数の反射波で異なるため ある信号は複数の時 間遅れを持って混合された畳み込み混合となる。 時間領域で直接分離フィルタを 求める I C Aは最終的な解への収束が非常に遅いため、 周波数領域で周波数毎に I C Aを適用する方法が現実的である。
[到来方向推定]
周波数領域での I C Aを用いて信号源の方向を位置情報として推定する従来技 術を図 1を参照して簡単に説明する。 J個のセンサ 1ぃ 12, ···, が直線状に 配列されている。 センサ 1」（ j = 1 , 2, ···, J)の位置を dj、 そのセンサ 1 j での観測信号を Xj(t )とする。 センサ 1！ , ···, 1」·の配列方向と垂直な方向を 90° とし、 源信号 s 1:)の到来方向を0° ≤6i≤ 180° とする。 I個の源 信号 S l(t)， ···， S l(t)の混合信号を J個のセンサ l ,〜ljにより、 観測信号 x^t), x t)として検出しているものとする。
信号の到来方向の推定は周波数領域で行われることが多い。 これには、 観測信 号 χ」 （t)を短時間フ一リェ変換して周波数領域における時間系列信号 Χ』（ω, m)を求める。 ωは角周波数 （周波数を f とすると ω = 2 π f ) mは時刻を表す番 号である。 源信号 S i(t) ( i = 1 , -, I ) も、 同様に、 周波数領域における時 間系列 S , m)に変換したものとすると、 観測信号 X』（ω, m)は、 信号 s iの 信号源からセンサ 1』 への周波数応答 Aji )を用い Xj ， m)=∑i=1 IAji (ω)
Figure imgf000004_0001
と表現することができる。 ベクトルと行列を用いると、
X (ω, m) =Α (ω) S (ω, m) · · · ( 1 ) となる。 ここで、
X (ω, m) = [Χ,(ω, m)， ···, Χ}(ω, m)] τ · · · (2)
S (ω, m) = CS^a), m), ···, S! , m)] T · · · (3) はそれぞれ、 J個のセンサによる観測信号および I個の源信号をべクトル表現し たものである。 A (ω) は周波数応答 A (ω) を要素とする J χ I行列であり、 信号混合系の周波数応答を表現するので混合行列と呼ばれる。 [a] τはベク トル 又は行列 aの転置を表わす。
図 1において、 方向 ；から到来する源信号は、 di = 0の位置のセンサ 1 ,に 対し センサ 1 j にもて ji== c— cos<9iだけ速く到達する。 cは源信号 S iの速 度である。 従って、 直接波のみを考慮すると角周波数 における周波数応答は Ajj ) =exp ( j ω c d j cos Θ；) · · · (4;
とモデル化することができる。 方向 を変数とする到来方向ベク トルを
a (ω, θ ) = [exp 、 j ω c— 1 d !cos (9 )、 exp ( j ω c _1 d 2cos Θ ), ···, exp ( j ω
Figure imgf000004_0002
Tと表現すると、 観測信号は
X (ω, m) =∑ i=1 1 a (ω , θ；) S , (ω, m)
とも近似表現することができる。
独立成分分析を用いて信号源の方向を推定する方法については例えば S. Kurita, H. Saruwatari, S Kaj ita, K. Takeda, and F. Itakura, "Evaluation of blind signal separation method using directivity pattern under reverberant conditions, "in Proc. ICASSP2000, 2000, pp.3140-3143 (文献 1 という） に示され ている。 この方法を以下に簡単に述べる。 観測信号 X (ω, m) は X (ω, m) =Α (ω) S (ω, m) であり、 源信号 S (ω, m) が混合されたものであるので、 互いに独立ではない。 X (ω, m) に独立成分分析
Υ (ω, m) =W (ω) X (ω, m) · · · (5) を適用すると、 互に独立な分離信号
Υ (ω, m) = ίΥ,ίω, m), ···, Υ (ω, m)] τ · · · (6) が得られる。
W (ω) は要素が Wi」. (ω) である I χ Jの行列であり、 分離行列と呼ばれる。 例えば、 1 =3 = の場合、 独立成分分析は、 Υ,(ω, m)と Υ2(ω， m)が互い に独立になる様に
Υ^ω,ιτι) Ψη ω)Ψ12_ω)
Υοΐω,ΐΉ) W21(ro) W22(o) Χ2(ω,ηι) を満足する分離行列 W (ω) を求める。 源信号 Sid m). ···, S,((w, m)が 互いに独立であれば、 分離信号 Υ^ω, m), ···, Υ,(ω, m)はそれぞれ、 源信 号の何れかに対応する。 但し、 独立成分分析は信号の独立性のみに基づいている ところから、 分離信号の順序と大きさに関して任意性がある。 即ち 分離行列 W (ω) の行を入れ換えても、 W (ω) の行を定数倍しても、 独立成分分析の解で ある。 なお後で述べるが順序の任意性はパーミュテーションの問題となり、 大き さの任意性はスケーリングの問題となる。
分離行列 W (ω) の i行目は、
Figure imgf000005_0001
である。 この Wi(o>)は分離信号 Yi( , m)を作り出していることが分かる。 従 つて、 Wi(o>)は、 源信号 S! w, m)， ···， S,(ft)( m)の何れか 1つを強調して 取り出し、 それ以外を抑圧している。 Wi((W)が形成する指向特性を解析すること により、 何れの方向から到来する信号を取り出し、 何れの方向から到来する信号 を抑圧しているかを解析することができ、 この解析により、 源信号 s t)の到来 方向を推定することができる。 これを全ての Wi(oj)， i = 1 , ■··, Iに対して行 うと、 分離行列の W (ω)各行 Wi(w)が取り出している源信号の到来方向 Θ (ω) = ίθχ{ω), ···, （ω)] τを推定することができる。
Wi(w)が形成する指向特性は、到来方向べク トル a (ω, Θ)を用いて Bi (ω, Θ) =Wi(o)) a (ω, Θ) と表現することができる。 この (ω, Θ) は方向 Sにある源信号から、 分離信号 Υ;(ω, m)への周波数応答と考えられる。 周波数 3 1 56Hzにおいて独立成分分析後の指向特性のゲイン I Β;(ω， ） I特性を 図 2に示す。 図 2の横軸は ^を、 縦軸はゲインを表わす。 実線は分離行列の 1行 目が与える指向特性 I B ft), θ) Iであり、 破線は 2行目が与える指向特性 I B 2(ω, Iである。 実線は 55° でゲインが最小となっており、 破線は 1 2 1 ° でゲインが最小となっている。 このことから、 分離行列の 1行目は 55° から到 来する信号を抑圧して 1 2 1 ° から到来する信号を取り出し、 逆に分離行列の 2 行目は 1 2 1 ° から到来する信号を抑圧して 55。 から到来する信号を取り出し ている。 従って、 Θ (3 1 56H z) = [ 1 2 1 ° ， 55° ] τと推定すること ができる。
複数の信号源方向を複数センサを用いて、 センサの観測信号を周波数領域に変 換して推定する方法として、 MUS I C (Multiple Signal Classification) 法 かある (例 ば S. Unnikrishna Pi 1 lai , "Array Signal Processing", Springer -Verlag, 1989, ISBN 0-387-96951-9. ISBN 3-540- 9695ト 9 参照）。 この方法はセン ザの数 Jより 1個少ない ( J - 1 ) 個までしか信号源の方向を推定できない。 し かし前記独立成分分析 ( I CA) を用いる方法 （以下単に I C A法ともいう） に よれば 2つのセンサで 2信号の混合に対応することができ、 この点で、 MUS I C法より優れている。 しかし、 上記の I C Aを用いる方法では 3信号以上の混合 に対処するには以下に示すように困難を伴う。 また、 指向特性のゲインの最小値 を求める操作に多くの計算を要することも欠点である。
3信号の混合に対して 3つのセンサを用いて I C A法を適用した状況を説明す る。 この場合の I C Aは分離行列が 3x3になるだけで同様に行えるが、 指向特 性のゲインを解析するのに困難を伴う。 この場合の I C A後の周波数 2734H zにおいて指向特性のゲイン I B w, θ) Iを図 3に示す。 図 3の実線は分離行 列の 1行目が、 破線は 2行目が、 1点鎖線は 3行目がそれぞれ与える指向特性で ある。 この場合、 各源信号は分離行列のある行によって強調され、 他の 2つの行 によって抑圧される箦である。 しかし、 抑圧する 2つの行が同じ方向で抑圧して いるとは限らない。 例えば、 図 3においては 45° 付近で i B2 I と I B3 I が共 に極小となり、 w! w)が 45°付近の源信号を取り出し、 w2(w)及び w3(cw)はこ の源信号を抑圧している。 同様に 90° 付近で I Β, I 及び I B3 I が極小となり、 w2(o))で 90° 付近の源信号を取り出し、 w!Cw)及び w3(a>)でこの信号を抑圧 していると読み取れる。 しかし λ¥3(ω)が 1 20° 付近の源信号を取り出し、 w ω)および w2(cw)がこの源信号を抑圧していると考えられるが、 1 20° 付近 での I B, I 及び I B,2 I の極小がかなりずれている。 この様なずれが大きくなる と、 何れの方向の抑圧が何れの源信号に対応するか判然としなくなる。 従って、 I C A法による従来技術を 3信号以上の状況に対して適用するのは困難であると 考えられる。
〔ブラインド信号分離〕
次に I CAを用いたプラインド信号分離の従来技術を説明する。 プラインド信 号分離とは、 観測された混合信号から源信号を推定する技術である。 以下では、 I個の源信号が混在した混合信号が J個のセンサで観測される場合を例にとって 説明する。
信号源 iから発生した源信号を s i (t) ( i = 1 , -, I . tは時刻）、 センサ jで観測された混合信号を χ』 (t) ( j = 1，'·'， J) とすると、 混合信号 χ」· (t) は次式で表わせる。
xj(t) =∑i=1 I (ajj* S i) (t) · · · (7)
ここで、 a jjは信号源 iからセンサ jへのィンパルス応答、 氺は畳み込み演算 子である。 ブラインド信号分離の目的は、 観測信号 χ』 (t) のみを用いて、 分離 のためのフィルタ wkj及び分離信号 yk ( t ) (k = 1 , ···, I ) を次式により求め ることである。
yk(t) =∑j=1 J (wkj*x ]·) (t) · · · (8)
時間領域での畳み込み混合は、 周波数領域における複数の瞬時混合に変換する ことができる。 すなわち、 前記の式 （7) 及び式 （8) は、 それぞれ前記式 （1 ) 及び式 （5) で表わされ、 前述したように、 W (ω) は分離行列であり、 I CA を用いて Yk (ω, m) と Yk' (ω, m) が互いに独立となるように計算された もの、 つまり I C Aの解である。
以上の周波数領域でのブラインド信号分離を行うに際し問題となるのは、 パー ミュテーシヨンの問題とスケーリングの問題である。
前述したように分離行列 W (ω) は行の入れ替えを行っても独立成分分析の解 となる。 つまりある角周波数 ω において、 W (ω) が I C Αの解だとし、 複素数 を要素とする任意の対角行列を D (ω) とし、任意のパ一ミュテーシヨン行列（こ の行列を任意の行列の左から掛けた結果は、 この任意の行列の行を入れ替えた行 列となる) を Ρ (ω) とした場合における Ρ (ω) D (ω) W (ω) もまた I C Aの解となる。 これは、 I C Aが源信号間の統計的独立性のみを条件として源信 号の分離を行うことに基づくものであり、 この D (ω) によって与えられる解の 自由度をスケーリングの自由度と呼び、 Ρ (ω) によって与えられる解の自由度 をパーミュテーシヨンの自由度と呼ぶ。
そのため、 適切なブラインド信号分離を行うためには、 全ての ω について、 I C Αの解から分離行列として適切な解 W (ω) を特定しなければならない。 一般 に、 この適切な解 W (ω) の特定は、 任意に求めた I C Αの解に、 適切な D (ω) や Ρ (ω) を乗じ その結果を適切な解 W (ω) とすることによって行われる。 そして、 全ての について、 この D (ω) を適切に決める問題をスケーリングの 問題と呼び、 Ρ (ω) を適切に決める問題をパーミュテーシヨンの問題と呼ぶ。 また、 パーミュテーシヨンとは、 {1,2, ···, 1} から {1,2, ···, 1} への全単射な関数 Ζ ： {1,2,-, 1} → {1,2, -, 1} であり、 パーミュテーション行列と 1対 1対応す る。
スケーリングの自由度は、 時間領域において周波数特性を変化させるフィルタ の自由度に相当する。 それゆえ、 時間領域において歪のない分離信号を生成する ためには、 全ての ω について、 D を適切に決める必要がある。 このスケー リングの問題は、 例えば、 D (ω) = d i a g (W"1 (ω)) とすることで容易 に解決することができる。 d i a g (a) は行列 αの対角化 （対角要素の他は全 ての要素を 0とする） を表わす。 つまり任意に求めた I C Αの解 W。 （ω) の逆行 列を求め、 それを対角化した行列を D (ω) とし、 D (ω) W。 （ω) を、 適切な 分離行列 W (ω) として特定する。 このことは既に知られている。 例えば参考文 m： K. Matsuoka and S. Nakashima, "Minimal Distortion Principle for Blind Source Separation" ,Proc. ICA 2001. pp.722- 727.に示されている。
一方、パーミュテーションの自由度のため、前記の式（5) の演算結果として、 例えば、 ある角周波数 については分離信号 (ω,, m) が源信号 Si {ω , m) の推定値として出力され、別の角周波数 ω2については分離信号 (ω2, m) が源信号 S2 (ω2, m) の推定値として出力されることもありうる。 このような 場合、 時間領域の出力信号 y, (t) の中に、 時間領域の源信号 s, (t) の成分 と s2 (t) の成分とが混在してしまい、 分離信号を正しく生成することができな い。 そのため、 時間領域における出力信号 y〗 (t) が正しく源信号 (t) の 推定値となるためには、 全ての ω について、 Υ1 (ω, m) が Si (ω， m) の推 定値となるように P (ω) を適切に決める必要がある。
従来のパーミュテーション問題の代表的な解法として、 前記文献 1に示す信号 の到来方向を推定する方法がある。 つまり、 図 2を参照して説明したように、 各 周波数における分離行列 W (ω) の各行と対応する指向特性を求める （図 2では f = 3 1 56 H zの指向特性のみを示している）。これら各周波数の指向特性にお いて例えば分離行列 W (ω) の 1行目が与える指向特性のゲイン最小値が 55° で 2行目が与える指向特性のゲイン最小値が 1 2 1 ° になるようにする。 つま りある角周波数 ωηではその分離行列 W (ωπ) の 1行目が与える指向特性のゲイ ン最小値が 1 2 1 ° で、 2行目が与える指向特性のゲイン最小値が 55° であれ ば、 その分離行列 W (ωη) の 1行目と 2行目とを入れ替える。 つまりこの入れ替 えを行うパーミュテーシヨン行列 Ρ (ωη) を W (ωη) にその左から掛算する。 このパ一ミュテーション問題の解法は先にも述べたように指向特性のゲイン最 小値を求める操作に多くの計算を要し、 しかも信号源の数 Iが 3以上では全角周 波数の W (ω) について適当に並べ替えを行ってみるという試行錯誤が必要とな る。 更に図 3を参照して説明したように W (ω) のある行がある方向の信号 Si (ω, m) を取り出している場合は、 その W (ω) の他の行は全てその方向の信 号 Si (ω, m) を抑圧する状態に必ずしもならない。
また指向特性のゲインが低い方向を探索して推定する信号到来方向の推定精度 は信号源の位置によって異なり、特に一対のセンサ 1」と 1』.を結んだ直線（以下、 センサペア軸という） に信号到来方向が近い場合に誤差が大きくなる。 このこと を実験により示す。 図 4 Aに示すように、 2つのセンサ 1 01, 102としてマ イク口ホンを距離 2. 83 cm離して配置し、 センサ 101と 1 02の中点 （原 点） から一定の距離 （約 1 5 O cm) 離れ、 かつ角度間隔が 20° 離れた 2つの 信号源 1 1 1, 1 1 2として音源を設け、 センサ 1 01からセンサ 1 02を見た 方向を基準 （0° ) として、 原点からみた信号源 1 1 1の角度 Sが 10° から 1 50° となるまで音源 1 1 1と 1 12を前記一定距離及び角度間隔を保持した状 態で移動させた。
このように音源 1 1 1, 1 1 2を移動させながら行ったブラインド信号分離結 果を図 4 Bに例示する。 図 4 Bの縦軸は、 信号対妨害音比 （Signal to Interfer ence Ratio)を示しており、分離音中に含まれる目的信号及び妨害信号を用いて、
S I R l Ologw (目的信号のパワー/妨害信号のパワー） （dB)
のように計算したものである。 また、 図 4 Bの横軸は、 原点からみた音源 1 1 1 の角度 0を表し、 実線はこの実験結果を示しており、 点線はパーミュテーシヨン に正解を与えたときの S I Rを示している。
この図より信号源 1 1 1が、 センサペア軸に近い方向 （0° 或いは 180° 付 近） に近づくと、 実験結果の S I Rはパーミュテ一シヨンに正解を与えたときの S I Rに比べ大きく低下していることがわかる。 このことは信号源 1 1 1の方向 がセンサペア軸方向と近いとパーミュテ一ションが誤っていることに基づくと考 えられる。
以上述べたように独立成分分析を用いて分離行列を求め、 その各行から指向特 性パターンを求めて、 その利得が低い方向探索により、 信号源の方向 (到来信号 方向） を求め、 更にこれを利用してブラインド信号分離を行う場合は、 前記指向 特性パターンを求めて、 利得が低い方向を探索するために多くの計算時間を必要 とした。
この発明の目的は独立成分分析により分離行列を求めて信号源の位置情報を推 定するための計算時間が短かい位置情報推定装置、 その方法及びプログラムを提 供することにある。 発明の開示
この発明によれば周波数領域の分離行列 W (ωι)， ···， W (ωΝ) の逆行列 （ I < Jの場合は擬似逆行列） を計算してスケーリングとパ一ミュテ一シヨンの自由 度を含む混合行列 A (ωχ), ·'·, Α (ω ) の推定値 Η ( ···， Η (ωΝ) を生 成し、 これらの周波数ごとの Η (ωη) (η= 1, -, Ν) について各列ごとにそ の 2つの要素 Η」Ί (ωη) と H i (ωη) (j、 j ' はセンサを表すパラメータ、 i は信号源を表すパラメータ）の比から、信号源 iの位置情報の一つのパラメータ、 例えば信号源が存在する円錐面あるいは曲面を計算する。
このように要素比に基づき表された計算式を演算すればよく、 分離された信号 の指向性パターンを求め、 更にその極小値角度の探索を行う場合より少ない計算 量で済む。 また要素比をとるため、 前述したスケーリングの自由度による影響が なくなる。 図面の簡単な説明
図 1はセンサァレーと到来信号の各センサへの到来時間差の関係を説明するた めの図である。
図 2は 2音源混合信号に対し I C Aにより計算した分離行列の各行のゲィン 指向特性を示す図である。
図 3は 3音源混合信号に対し、 I C Aにより計算した分離行列の各行のゲイン 指向特性を示す図である。
図 4 Aは予備実験に用いたセンサと信号源の関係を説明するための図である。 図 4 Bは前記予備実験の結果を示す図である。
図 4 Cは推定角度とその誤差及び偏角誤差との関係を示す図である。
図 5はこの発明を信号到来方向の推定に適用した第 1実施形態の機能構成例を 示すブロック図である。
図 6は第 1実施形態の処理手順の例を示す流れ図である。
図 7は図 5中の角度計算部の具体的機能構成例を示すプロック図である。
図 8は図 6中のステップ S 4の具体的処理手順の例を示す流れ図である。
図 9はこの第 1実施形態により方向推定した実験結果を示す図である。 図 1 O Aは音源が 2個、 マイクロホンが 3個の場合に第 1実施形態により方向 推定した実験結果を示す図である。
図 1 0 Bは図 1 O Aについての実験と同一条件で MU S I C法により方向推定 した実験結果を示す図である。
図 1 1 Aは音源が 2個、 マイクが 3個の場合に第 1実施形態により方向推定し た実験結果を示す図である。
図 1 1 Bは図 1 1 Aについての実験と同一条件で M U S I C法により方向推定 した実験結果を示す図である。
図 1 2はこの発明をブラインド信号分離に適用した第 2実施形態の機能構成例 を示すブロック図である。
図 1 3は第 2実施形態の処理手順の例を示す流れ図である。
図 1 4は図 1 3中のステップ S 1 4の具体的処理手順の例を示すプロック図で める。
図 1 5は複数円錐面の共通直線方向を説明するための図である。
図 1 6は図 1 2中の到来方向決定部 1 6の具体的機能構成の他の例を示すプロ ック図である。
図 1 7はこの発明をプラインド信号分離に適用した第 3実施形態の機能構成例 を示すプロック図である。
図 1 8は第 3実施形態の処理手順の例を示す流れ図である。
図 1 9は図 1 8中のステップ S 3 5の具体的処理手順の例を示す流れ図である。 図 2 0はセンサ配置と信号源位置と推定曲面との関係を説明するための図であ る。
図 2 1は推定球面の例を示す図である。
図 2 2は実験に用いた室とマイクロホン配置と音源との関係を示す図である。 図 2 3は小間隔マイクロホン対を用いた推定方向のヒストグラムである。
図 2 4はその周波数に対する推定方向の分布を示す図である。
図 2 5は音源 24及び 25に対する推定球面の半径の周波数に対する分布を示す 図である。
図 2 6は各方法での実験結果を示す図である。 図 27はこの発明の第 4実施形態の要部の機能構成例を示すプロック図である。 図 28はこの発明をブラインド信号分離に適用した第 5実施形態の処理手順の 例を示す流れ図である。
図 29は第 5実施形態に対する実験結果を示す図である。
図 3 OAは複数の円錐面の共通直線方向を求めて分離行列のパーミュテーショ ンを解決する処理手順の要部を示す流れ図である。
図 30Bは円錐面の推定と球面の推定を利用して分離行列のパーミュテーショ ンを解決する処理手順の要部を示す流れ図である。
図 30 Cは円錐面の推定と球面の推定を利用して分離行列のパーミュテ一ショ ンを解決する他の処理手順の要部を示す流れ図である。 発明を実施するための最良の形態
まず信号源の位置情報として方向情報を推定する場合にこの発明を適用した実 施形態を説明する。 なお以下の説明において同一又は対応するものについて図面 中に同一の参照番号を付けて、 重複説明を省略する。
[第 1実施形態]
この第 1実施形態は信号源の方向 つまりその信号源から放射された源信号の 到来方向を求める。
図 5に第 1実施形態の機能構成図を、 図 6にその処理手順の一部の流れ図をそ れぞれ示す。
信号源の数 I以上の個数である J個のセンサ 1い 12, …， が、 例えば図 1 に示したように配列されている。 隣接センサの間隔は通常は、 源信号の最も短か い波長の 1 /2以下である。 これらセンサ 1」 （ j = 1， 2， ···， J) で観測され た信号 Xj (t) はそれぞれ周波数領域変換部 1 1」 において、 例えば短時間フー リェ変換により周波数領域信号 Χ」·（ω， m)に変換される（図 6、ステップ S 1 )。 これら周波数領域信号 Xj (ω, m) に対する各角周波数 ωηごとの分離行列 W (ω η) (η= 1, 2, -, Ν) が分離行列算出部 1 2において独立成分分析により算 出される （図 6、 ステップ S 2)。
W (ω) = (W ίωχ), W ( 2)， ···， W N)) の各周波数ごとの分離行列 W (ωη) の逆行列が逆行列算出部 1 3で算出され 逆行列 Η (ωη) が求まる （図 6、 ステップ S 3)。
Η (ω) = (Η (ω{), Η (ω2), ···, Η (ωΝ))
なおこの逆行列の計算は J > Iの場合は疑似逆行列を計算することになる。 疑 似逆行列としては例えば Mo o r e -P e n r o s e型一般逆行列が用いられる。 この実施形態においては H (ω) 中の少くとも 1つの周波数の逆行列 Η (ωη) の各列 iの 2つの要素 (ωη) と Η」. i (ωη) との比の偏角から方向情報、 つま り源信号の到来方向が角度計算部 1 4で計算される （図 6、 ステップ S 4)。 この 角度計算部 1 4の具体的機能構成例と処理手順例を図 7及び図 8を参照して説明 する。 選択部 14 aにより角周波数 ωηの逆行列 Η (ωη) 中の未選択の一つの列 iが選択され （図 8、 ステップ S 4 a)、 その i列目から 2つの要素 Η』； (ωη) と Η ; (ωη) が選択される (図 8、 ステップ S 4 b)。
偏角計算部 1 4 bで、選択された要素の比の偏角 arg[Hji (ωη) /Hr； (ωη)] が計算され、 また間隔計算部 1 4じで、 センサ情報格納部 1 5内のセンサ ljと 1 j. の位置情報 dj と dj. が取り出され、 これらの差からセンサ 1 j と 1 の間隔 dj. - d が計算される （図 8、 ステップ S 4 c)。 位相回転計算部 1 4 dで、 間隔 d j-dj. と 角周波数 ωηでの単位距離当りの信号の位相回転量 o)n/c ( cは源信 号の速度） との積が計算されて間隔 dj— dj. での信号の位相回転量が求められ、 この位相回転量により偏角 arg [Η (ωη) ΖΗ」·. i (ωη)] が除算部 1 4 eで割 算される (図 8、 ステップ S 4 d)。 なお単位距離当たりの位相回転量 は例 えば f = 680 H zの音波 (速度 c = 34 Om/s) では、 1 m当たり 2 ττ · 2 の位相が回転することになる。
この除算部 1 4 eの割算結果 Gi (ωη) の絶対値が 1以下か否かかが判定部 1 4 f で判定され （図 8、 ステップ S 4 e)、 1以下であれば Gi (ωη) の逆余弦 S i n) =cos-1Gi (ωη) が逆余弦計算部 1 4 gで計算される （図 8、 ステップ S 4 f )0 つまり角度計算部では次の計算が行われる。
Figure imgf000014_0001
[Η (ωη) ノ Η s (ωη)] / (ω (ά}- άν) / c))
^ i((yn)=cos"1Gi (ω) ： I Gi((Wn) I Iに対して、 ·'· （9) ステップ S 4 eで I Gi (ω„) | が 1以下でなければ角度 ^ i (ωη) が虚数にな るので、 別の組み合せを選択するため、 判定部 1 4 f で選択した i列目のすべて の要素の組み合せを選択したか判定され （図 8、 ステップ S 4 g)、 選択していな い組み合せがあればステップ S 4 bに戻り、 全ての組み合せを選択していれば、 判定部 1 4 f ですベての列を選択したか判定され （図 8、 ステップ S 4 h)、 選択 していない列があればステップ S 4 aに戻り、 全ての列を選択していれば角度計 算処理を終了する。 なお図 7中の間隔計算部 1 4 cは各周波数の角度計算部 1 4 に対して共通に用いられる。
角度計算部 1 4から H (ω) 中の選択した角周波数 ωηの逆行列 Η (ωη) にお ける各列と対応して、 その選択順に、 もし第 1列目から順次選択すればその順に 式 （7) の計算結果、 つまり I個の信号源の方向 （信号到来方向） Θ (ωη) = (Θ ！ (ωη), θ2 (ωη), ···, Θ！ (ωη)) が出力される。 これら θ、 (ωη)， θζ (ωη)， ···, θ ! (ωη) は各源信号 S l (t)， s2 (t), ···, s J (t) の到来方向 (信号 源 1， 2, ···, Iの方向) のいずれかの 1つと対応する。
この実施形態において信号到来方向が推定できるしくみを以下に説明する。 独 立成分分析 （ I CA) により分離が達成されていれば、 I CAにより算出した分 離行列 W (ω) と真の混合行列 Α (ω) とは Ρ (ω) D (ω) W (ω) Α (ω) = Iの関係にある。 ここで D (ω) はスケーリングの自由度を示す対角行列、 Ρ (ω) はパーミュテーシヨンの自由度を示すパーミュテーシヨン行列、 Iは単 位行列である。 I C Αを用いても、 混合行列 A (ω) そのものは一般には算出で きない。 しかし、 W (ω) の逆行列 Η (ω) =W— 1 (ω) =Α (ω) Ρ (ω) D (ω) を算出すれば、 スケーリングとパーミュテーションの自由度を含む混合行 列の推定値が得られる。 すなわち、 逆行列 Η (ω) は、 混合行列 Α (ω) の列を Ρ (ω) により並べ替え、 さらに各列にそれぞれ D (ω) の対角要素を掛け合わ せたものになる。
この実施形態では、 逆行列 Η (ω) の同じ列 iから 2つの要素 H」i (ω) と Hj' i (ω) を取り出してそれらの比 H」i (ω) /Η- (ω) を求めることで、 算出で きない D (ω) によるスケーリングの自由度を取り除く。 すなわち、 ¾(ω) [Α(ω)ρ(ω)θ(ω)1ί Ajz(i)((o
Η』 ） [Α(ω)Ρ(ω)0(ω¾, Α』,ζ(0(ω)
となる。 ここで、 Ζはパーミュテ一シヨン行列 Ρ (ω) を右から掛けることに対 応ずる順列である。 逆行列 Η (ω) の全ての列 i に対して式 （2 1) に関わる操 作を行うことで、 P (ω) によるパーミュテ一シヨン Ζにかかわらず、 全ての信 号の到来方向が推定できる。
背景技術の説明では、 混合行列 Α (ω) の要素を Aw (ω) =exp (ja) c-'djC os とモデル化した。 しかし、 この実施形態では、 分離行列 W (ω) の逆行列 Η (ω) を用いて、 スケーリングとパーミュテーションの自由度を含む混合行列 Α (ω) の推定値を算出しているので、 この単純なモデルでは不十分である。 そ こで、 振幅の減衰度 a (実数） と原点における位相差 exp を用いて、 A ji {ω) = a jjex j ；) exp j ω c 1 d jcos Θ；) とモテル化する。 このモテル を用いて Α」ζω (ω) /Ayz(i) (ω) を計算すると、 式 (2 1) の関係から、 ¾W =½ ) =¾ ωο -! (d d )cos6 ) (22)
Hj-ilro) Α/Ζ(ί)(ω) /ζ(ί)
が得られる。 その結果、 上記の Gi (ω) は
G; (ω) =arg [Η (ω) /Hj>i (ω)] / (ω c' (dj—d )) =cos θ z(i) となり、 I Gj (ω) I <.lであれば Szw cos-'Gi (ω) が実数となり、 到来方 向を推定することができる。 先に述べたように、 パーミュテ一シヨン Ζの自由度 により、 I個の方向 Θ (ω) = ίθζω ( )， ···， θζω (ω)] 全体を適当に並べ 替えたものが、 信号 s!, ···, の方向に対応する。
H (ω) 中の複数の周波数又は全ての周波数の各逆行列 H (ωη) (η= 1 , …，
Ν) について先に述べたように各列ごとに角度 Si (ωη) を求め、 これらの全体 により各到来方向を決定してもよい。 つまり角度計算部 1 4で推定された各周波 数の到来方向は図 5中のソート部 32においてその角度ごとに分けられる（図 6、 ステップ S 5)。 例えばそれぞれ方向 （角度） の大きい順に並べられる。 角周波数 についての Θ (ω ) 中の各成分が大きい順に （ （α^), θ2 ((Wj), ···, θ ! (ω,)) と並べ、 ω2についての Θ (ω2) 中の各成分が大きい順に （ （ω2), θζ (ω2)， ···， Θ！ (ω2)) と並べ、 ···， ωΝについての Θ (ωΝ) 中の各成分が大 きい順に (θ γ (ωΝ), θ2 (ωΝ), ■··, θ！ (ωΝ)) と並べられる。 このように大 きい角度順に分別された角度は、 大きい順が同一であっても周波数間でばらつき がある、 つまり角度 i ( ,), θ{ (ωζ), ···, θ； (ωΝ) ( i = 1 , ···, I ) はば らついている。
よって統合部 33で、 分別角度 (ω,), θ (ω2), …， θ{ (ωΝ) ごとに 1 つの角度 iに統合されて 1つの角度 （到来方向） ；とされる （図 6、 ステップ S 6)。 この統合の方法としては例えば分別角度 （ {ωλ), Θ i (ω2), ···, θ { (ωΝ)) ごと平均値が統合角度 とされ、 あるいは分別角度 （Si (ωχ), θχ (ω 2), ···, θ { (ω )) ごとの最頻値や中央値などを統合角度 0;としてもよい。 この ように分別角度ごとに 1つの角度に統合することにより、 1つのある周波数の逆 行列 H (ωη) からのみ到来方向を推定する場合より、 正しく到来方向を推定する ことができる。
なお図 8に示した角度 （方位） 計算推定処理において、 その 1つでも選択した 列について 角度 0i (α)η) を求めることができなかった時は そこでその逆行 列 H (α)η) に対する処理を終了し.， 他の周波数の逆行列 H (ωπ) に対する処理 を行い、 最初に全ての列について角度 (方向) を計算することができたら、 その 時の計算結果を推定 （方向） 0い -, , としてもよい。 あるいは各周波数の逆 行列 H (ωη) 中で全ての列について角度計算ができた結果を角度ごとに分別し、 統合してもよい。 あるいは 1つの列について角度 0 i (ωη) を計算することがで きないものがあれば、 その状態が最初に発生した時に、 その後の処理を全て終了 し、 あらためて観測信号を求める処理からやりなおすことにより、 推定結果の信 頼性を高めるようにしてもよい。
この第 1実施形態の実験例を述べる。 残響時間が 1 9 Omsの部屋に 3つのマ イク口ホンを 56. 6mm間隔で 1列に配列し、 その配列方向を基準として角度 48° 、 73° 、 1 1 9° の方向に 3つの音源を配置した。 これら音源からの音 響信号の混合時間を 6秒、 観測信号の標本化周波数を 8 kH z、 空間エイリアシ ングが生じない最大周波数を 3 kHz, 短時間フーリエ変換のフレームを 1 02 4サンプルとした。 前述したようにして、 各周波数ごとにそれぞれ計算した各角 度を図 9に、 横軸を周波数、 縦軸を方向として示す。 図 9中の◊， +, □はそれ それ 3つの音源方向の推定計算値を示す。 この結果を 3つの角度範囲に分別し、 各分別された角度平均値は 45° , 74° , 1 23° となった。 3音源の方向を 3つのマイクロホンで推定することは MU S I C法ではできないが、 この実施形 態によれば可成り正しく各音源方向が推定されていることがわかる。
この実施形態の方法と MUS I C法とを比較するため、 音源方向を 48° と 1 19° と比較的大きく角差がある状態で 2つの音源を配置した場合について同様 に実験を行った。 図 1 OAにこの実施形態の方法により得られた結果を、 図 1 0 Bに MUS I C法により得られた結果をそれぞれ示す。 これらよりいずれの方法 も、 かなり正しく方向を推定していることがわかる。 この実施形態による結果を 2つの方向範囲で分別平均した結果は 45° , 1 23° であり、 MUS I C法は それぞれ 45° , 1 22° である。 次に、 音源方向を 1 05° と 1 1 9° と方向 が比較的接近している状態で 2つの音源を配置した場合について同様に実験を行 つた。 この実施形態方法による結果及び MU S I C法による結果を図 1 1 A及び 図 1 1 Bにそれぞれ示す。 MUS I C法では大部分の周波数で音源方向の推定が できないが この実施形態では大部分の周波数で角度計算を行うことができ し かも、 これらを角度範囲で分別して、 それぞれ平均した値は 1 05° , 1 24° であり、 MUS I C法は 94° , 1 28° でありこの実施形態による方向推定が かなり正しいことがわかる。
以上述べたようにこの実施形態によれば従来方法において、 指向性パターンの 利得が低い方向探索する場合と比較して、 式 (9) に値を代入するだけで推定方 向が求まるため計算時間がかなり短かい。 なお前述したように I C Aにより求め た分離行列 W (ω) はスケーリ ングの自由度及びパーミュテーシヨンの自由度を 含んでいるためこれらを解決した分離行列 W' (ω) の逆行列を計算し、 つまり真 の混合行列 Α (ω) を計算し、 その混合行列 Α (ω) の各列について二つの要素の 比を用いて到来方向を推定することが考えられるかもしれない。 しかし、 真の混 合行列 Α (ω) そのものを求めることは、 例えば信号源の信号 S i(t)の平均パヮ 一を 1にするなどの制約を与えない限りできない。 無線通信分野においては信号 源に対しそのような制約を与えることができるかも知れないが、 例えば信号源の 信号 S i( t )が人間が直接発声した音声信号である場合は、そのような制約は実質 的にはできない。 一方この第 1実施形態によれば、 スケーリング及びパーミュテ ーシヨンの自由度を含んでいる分離行列 W (ω) の逆行列 Η (ω) の各列について 二つの要素の比をとることによりスケーリングの自由度の問題を解決しておりど のような信号源に対しても適用でき、 しかも前記両問題を解決した分離行列を計 算する必要がなくそれだけ計算時間が短い。 更に前述したように、 各周波数につ いて求めた推定方向を予め決めた順に分別すればパーミュテーションの問題も簡 単に解決できる。 センサの数 Jと同数であっても各信号源の方向を推定すること ができる。 また音源方向が比較的接近していてもかなり正しく推定することがで ぎる。
[第 2実施形態]
この第 2実施形態は信号源の一つの位置情報として方向情報を求めるものであ り、 かつ少なく とも 2次元に配置された少なくとも 3個のセンサを用いて、 信号 源がいずれの方向に位置していても方向を推定でき、 従ってブラインド信号分離 におけるパーミュテーションの問題も比較的箇単に解決するものである。 つまり 方向情報に基づく円錐面を推定し 複数の円錐面の共通直線を推定して、 方向情 報を決定する。
この第 2実施形態をプラインド信号分離装置に適用した機能構成を図 1 2に、 その処理手順を図 1 3に示す。 例えば 4個のセンサ 1 lt 12, 13, 14が同一円上 で等間隔に配置され、 そのいずれの 2つのセンサの間隔も、 源信号の最短波長の 1/2以下とされる。以下の説明はセンサの数は J個であり、 J ≥ 3として行う。 第 1実施形態と同様に各センサ j ( j = l , -, J ) で観測された観測信号 χ』 ( t ) は、 それぞれ周波数領域変換部 1 1で例えば、 短時間フーリエ変換によつ て、 周波数領域の信号 Χ』 (ω, τη) にそれぞれ変換される （ステップ S 1 1 )。 これら周波数領域の信号 Χ』 (ω, τη) から、 分離行列算出部 1 2で、 独立成分 分析により、 各周波数ごとの分離行列
Figure imgf000020_0001
が算出される （ステップ S 1 2)。
逆行列算出部 1 3でこれら周波数ごとの各分離行列 W (ω) の逆行列 Η がそれぞれ算出される （ステップ S 1 3)。
Figure imgf000020_0002
この第 2実施形態では、 円錐面推定部 1 4において、 周波数ごとの各逆行列 H (ω) における各列ごとに異なる複数の要素対について各要素の比からこれら要 素と対応する二つのセンサのセンサペア軸を中心軸とし、 何れかの信号源が存在 する円錐面、 つまり一つの混合行列 Η (ω) の各列ごとに複数の円錐面がそれぞ れ推定される (ステップ S 1 4)。
円錐面推定部 1 4の機能構成は図 1 5中の角度計算部 1 4とほぼ同様であり、 その処理手順も図 8の手順と似ている。 円錐面推定部 1 4において行われるある 一つの周波数の逆行列 Η (ω) に対するステップ S 1 4の処理の具体例を図 1 4 を参照して説明する。
まず、 例えば、 円錐面推定部 1 4内のレジスタ内に格納されている制御パラメ 一夕 i と Ρを 0に初期化する (ステップ S 2 0)。 ここで、 iは各信号源の番号と 対応し、 pは i ごとの推定済み円錐面の個数である。
iを + 1 し （ステップ S 2 1)、 gを + 1 し （ステップ S 22)、 制御パラメ一 タ j， j ' はそれぞれ互いに異なる j以下の自然数を、例えばランダムに選択する (ステップ S 23)。 この制御パラメータ j , j ' の組は、 同一の iに対しては一 度選択された の組み合わせの選択は行わないようにする。 すなわち、 j ≠ j ' であり、 また、 例えば、 i = 1について、 一度 （ j , j ' ) = ( 1 , 2) が選択 されると、 i = lについては、 この処理が終了するまで、 再びステップ S 23に おいて （ j， j ' ) = ( 1 , 2) という選択は行われないものとする。 （また、 この 選択は、選択する j , j ' によって特定されるセンサペア軸（センサ j とセンサ j ' とを通る直線） が、 このルーチン処理においてそれ以前に選択された j , j ' によ つて特定されるセンサペア軸と同一直線上に配置されないように行われることが 望ましい。 これにより、 円錐面推定部 1 4は、 一定の誤差範囲内で中心軸が重な らない複数の円錐面を推定することになる。 なお、 各センサペア軸が同一直線上 となるか否かの判断は、 例えば、 センサ情報格納部 1 5に各センサの位置を示す べクトルを格納しておき、 ここから各センサの位置を示すべクトルの情報を抽出 して行う。）
次に、 センサ情報格納部 1 5から、 ステップ S 23で選択したパラメータ j に 対応する j番目のセンサ jの位置を示すベク トル djと、 変数 Γ に対応する j ' 番目のセンサ； j ' の位置を示すベク トル dj. とを抽出する （ステップ S 24)。 ま た、 混合行列 H (ω) から、 j行 i列要素 Hji(a>)と、 j ' 行 i列要素 H ^ω) とを選択して抽出する （ステップ S 2 5)。以上の 理は図 7中の選択部 1 4 aに より行われる。 従って、 i , p, j , j ' や選択した j , j ' などを格納するレ ジスタも選択部 1 4 aに設けられることになる。
抽出したこれらの情報を用い、 次式 （9' ) を演算する （ステップ 2 6)。
-1 ^[Η^(ω)/Η^(ω)]
θ;』, (ω) = cos -1 (9')
d j.— d.
J ここで ll dj— d II はセンサ 1」 と 1」. との間隔（距離）である。 この第 2実施形 態においては複数のセンサが 2次元又は 3次元に配置されているため各センサの 位置情報は例えば図 1 2中のセンサ 1 !〜 14が配置されている円の中心を原点と する 2又は 3要素の座標ベク トルで表わされる。 つまり式 （9) はセンサが 1次 元に配置され、 信号到来方向の角度が 2次元の場合であるが、 この式 （9' ) はセ ンサが 2次元あるいは 3次元に配置されていてもよく信号到来方向の角度が 3次 元空間でもよい場合に拡張されたものである。 従って式 （9') は式 （9) を包含 する。 この式 （9') により推定された角度 0 Λ u (ω) 及びその i , j , j ' を円錐面情報として円錐面推定部 1 4内のレジスタ （メモリ）に一時記録する（ス テツプ S 27)。なお図 1 2中に破線で示すように各周波数に対する間隔計算部 1
4 cは共通に用いられる。 また式 （9) で演算した角度 0i (ω) は三次元空間に おいてセンサペア軸 （センサ 1』 と 1』. とを結ぶ直線） に対する角度が Si (ω) となる直線の無数の集合、 つまり円錐面に信号源 iが存在することを推定してい ることになる。 この拡張された式 （9') の演算結果を単なる ； (ω) ではなく θ ^；, j (ω) と表わした。 ステップ S 26で行う演算は図 7中の偏角計算部 1 4 b、 間隔計算部 1 4 c、 位相回転計算部 1 4 d、 除算部 1 4 e、 判定部 1 4 f 及び逆余弦計算部 1 4 gにより、 図 8中のステップ S 4 c， S 4 d, S 4 e及び
54 f の処理により行う。
次に、 p = Pであるか否かを判断する (ステップ S 28)。 この Pは、 i ごとに 推定すべき円錐面の個数であり、 このステップでは、 その i について円錐面を P 個推定したか否かを判断する。 p = Pでなければステップ S 22に戻り、 p = P であれば i = Iであるか否かを判断する （ステップ S 29)。 すなわち すべて の iについて円錐面の推定が終了したか否かを判断し、 i = Iでなければステツ プ S 2 1に戻り、 i = 1であれば処理を終了する (ステップ S 1 4の具体例の説 明終了。）。 第 1実施形態では選択した i について角度 Si (ω) を 1つ推定すれば 次の iについての角度推定に移るがこの第 2実施形態では各 i ごとに複数 P個の 組の角度 (円錐面) Θ ^ jr (ω) を推定する。 ステップ S 28及び S 29の処 理は図 7中の判定部 1 4 f で行われる。
図 1 2中の到来方向決定部 1 6で、 円錐面推定部 1 4で推定された複数の円錐 面の情報（この例では、 i , j， j '、 Θ ^ i, jj. (ω))から源信号の到来方向 Ui (ω) = (方位角 (ω), 仰角 (ω)) ( i =l ,-, I ) が決定される （図 1 2、 ス テツプ S 1 5)。具体的には、例えば、 iについて推定された複数の円錐面のうち、 互いに線接触しているとみなされる共通直線の方向を、 信号源 i に対応する信号 の到来方向 (ω) として算出する。 この信号の到来方向 Ui (ω) の推定方法を図 1 5を参照して説明する。
センサ 1 i と 12の配列方向と直角方向にセンサ 13が配され、 センサ 1 i と 12 間、 センサ 12と 13間の各間隔は同一とされている。 i = 2の信号源 22に対し、 センサ 1 ,及び 12によりセンサペア軸 312を中心軸とする円錐面 412 ( ^ ^ 2, 12 (ω))が、センサ 12と 13によりセンサペア軸 323を中心軸とする円錐面 423 ( Λ 2. 23 (ω)) が、 センサ 11と 13によりセンサペア軸 313を中心軸とする円錐面 413 (θ Λ 2, 13 (ω)) がそれぞれ推定されている。 これらのうち、 円錐面 412と 4 13は、 共通直線 52で互いに線接触しているとみなされる。 この共通直線 52の方 向 u2 (ω) が、 信号源 22から放射された信号の到来方向 u2つまり信号源 22の 方向と決定する。 なお円錐面 423も平行移動すれば円錐面 412, 413とほぼ線接触 するが、 後述する第 4実施形態により円錐面 423は廃棄される。
図 1 2中の到来方向決定部 1 6におけるこの共通直線 5 iの方向 Si (ω) を決 定するための具体的計算方法の例を説明する。角周波数 ωに対して信号源 2；が存 在すると推定された複数の円錐面を 4j (1 ), ···, 4jy (P) とし、 これら円錐 面 4』」. (p) (p = 1 , ···, Ρ) の推定に用いられた 1対のセンサの位置情報を d j (p), dj. (p) とし、 角周波数 ωについて推定された円錐面 4」』. （p) と対 応する角度を ^』' （ω， ρ) とし、 円錐面 4jr (p) を表わすベクトルを uと する。
正規化軸べク トル計算部 1 6 aで前記各一対のセンサの位置間を結ぶ軸べク ト ル (dj (p) — d (p)) をそれぞれ長さ 1に正規化する。 つまり
vp= (dj (p) — dj. (p)) /II dj (p) -dr (p) II
をそれぞれ演算する。この vpと円錐面べク トルの内積はこれらべク トルのなす角 度の余弦とする。 つまり
νρ τ · uノ il u II = cos (9 jj. (ω, ρ )
関係が成立する。 知りたいのは共通直線 5 iの方向のみであるから円錐面べク ト ル uは単位べクトル、 つまり II u II = 1 とする。 全ての円錐面と共通する直線の 方向を求めるには
V= (vf νΡ)τ, c Λ (o)) = (cos^ ^jj. (ω, 1 )'"cos^ Λ ' （ , Ρ))τ
として次式の連立方程式を uについて解けば良い。 Vu= c " (ω) ( II u II = 1 )
一般にはこの連立方程式は解が存在しないか、 もしくは一意には決まらない。 そこで II Vu— c - (ω) IIを最小化する uを求めて前記連立方程式の解、つまり 共通直線 5 iの方向 ^ (ω) とする。 この誤差を最小化する uを求める計算が演 算部 1 6 bで行われる。 この方向 Ui (ω) は 3次元での方向であるから、 極座標 表示により方位角 (ω) 及び仰角 <i>i (ω) が求まることになる。
もっと計算を簡便にするには次のようにしてもよい。 図 1 6に示すように正規 化軸べクトル計算部 1 6 aで円錐面推定に用いた各センサ対についての正規化軸 ベク トル vp ( p = 1 , ·■·, P) を求め、 逆行列計算部 1 6 cで V= (Vい ···, vP) τの Mo o r e -P e n r o s e型一般逆行列 V +を計算し、 この V+と推定 円錐面と対応する余弦べクトル c " (ω)= (cos θ " ' (ω, 1 )--cos θ " jr (ω, P)) Tとを用いて最小ノルムな最小 2乗誤差解を求め、 大きさを正規化して近似 解としてもよい。 つまり演算部 1 6 dで
Ui (ω) =V+ c " (ω)/ II V+ c " (ω) ||
を計算する。
このようにして複数の推定円錐面の共通直線とみなされる直線の方向が各周波 数ごとに、 かつ各信号源ごとに求められる。
図 1 2中のパーミュテーション解決部 1 7は、 到来方向決定部 1 6で決定され た到来方向 u i= (θ ,) を用いて、 分離行列算出部 1 2で算出された分離行 列 W (ω) の行の並ぴ替えを行ってパーミュテーション問題が解決された分離行 列を生成する。
パーミュテーション解決部 1 7で行われる具体例としては、到来方位角 Θ につ いて以下のように並べ替えを行い、 解決できなかった列について同様に到来仰角 0 iによる並べ替えを行う。 つまり算出決定された到来方位角 (ω) がどの周 波数においても既定の順序、 例えば、 0い 02,…， θ 1が昇順となるように逆行列 Η (ω) の列を入れ替え、入れ替えができなかった列について仰角 がどの周波 数でも昇順になるように逆行列 Η (ω) の列を入れ替える並替行列が並替行列生 成部 1 7 aで生成され、 その並替行列の逆行列 P (ω) が逆行列生成部 1 7わで 生成され、 並替部 1 7 cでその逆行列 Ρ (ω) 、 分離行列 W (ω) に左から掛 算される。 並替行列生成部及び逆行列生成部 1 7 bはパーミュテーシヨン行列生 成部を構成している。
並替行列生成部 1 7 aでの処理を具体的に述べる。 この例では逆行列 H (ω) の第 1列目で （ (ω), , (ω)) が、 第 2列目で （ (ω), Φ2 (ω)) が， ···、 第 I列目で (θ ^ (ω), J (ω)) がそれぞれ式 （9') に基づき計算され たとする。 いま到来方向決定部 1 6 1から入力された到来方向を、 前記昇順に並 ベたものと区別するために Sおよび (^の各添字にダッシュ に」 を付けて、
(ω), φν (ω))， (θ2· (ω), 2· (ω)), ···, (θ Γ (ω), φ Γ (ω)) と表 わし、 これらを例えばまず r (ω) について昇順に並べた時に、 例えば （ 3. (ω)' φΆ. (ω)) > (Θ！. (ω), φ ν (ω)) > (θ2· (ω), φζ· (ω)) >··· であれば、 逆行列 Η (ω) の 3列目を 1列目に、 I列目を 2列目に、 2列目を 3 列目になるように移動し、 残りの列も同様に移動し、 θ (ω) が同一のものに ついては^, (ω) が昇順になるように列を移動する。 そのような列の移動並べ 替えを行う並替行列を作る。 このような並べ替えを行う行列を作ることは従来か らよく行われていることである。 このようにして全ての周波数について得られた 到来方向 (Θ i (ω)， φ！ (ω))， ···, (θ 1 (ω), φ2 (ω)) を用いて、 並替行 列が作られ 更にその逆行列つまりパーミュテーシヨン行列 Ρ (ω) が算出され る （図 1 3 ステップ S 1 6)。
この算出したパ一ミ ュテーシヨン行列 Ρ (ω) が並替部 （1 7 c) で分離行列 W (ω) の左から掛けられ、 その結果 W' (ω) =Ρ (ω) W (ω) がパーミ ュテ —ション問題が解決された分離行列として出力される (ステップ S 1 7)。つまり この分離行列 W' (ω) はいずれの周波数のものでも、 1行目は信号源 2!の信号 を分離する要素、 2行目は信号源 22の信号を分離する要素、以下同様に同一行の 要素は同一信号源からの信号を分離する要素となる。
分離行列 W' (ω) は、 時間領域変換部 1 8で、 例えば、 逆フーリエ変換によつ て、 時間領域の分離フィルタ係数群 Wn … wn
Wn … wn
に変換され、 これら分離フィルタ係数群は信号分離部 1 9に設定される。
信号分離部 1 9は、 各センサからの観測信号 (t) ,-, X j (t) と分離フ ィルタ係数群とにより式 （8) の演算を行い分離信号 (t) , -, Υ Ϊ (t) を 出力する。
図 1 2中に破線で示すように、 パーミュテーション解決部 1 7よりの並替えさ れた分離行列 W' (ω)と周波数領域変換部 1 1よりの周波数領域観測信号 Χ(ω, m) とにより式 (5) の演算を周波数領域分離信号生成部 1 9' で行い、 その結 果の周波数領域分離信号 Υ (ω, m) =W (ω) X (ω, m) を時間領域信号変 換部 1 8' で時間領域信号 (t), ···, y！ (t) を生成してもよい。 また到来 方向決定部 1 6で决定された各周波数ごとの到来方向 （01 (ω), φι (ω)), ···,
(Θ！ (ω), ！ (ω)) を第 1実施形態と同様に、 つまり図 5中のソート部 32 で方向範囲ごとに分別し、 各分別されたものごとの全周波数における到来方向を 統合部 33でそれぞれ統合してもよい。
以上説明したように、 この実施形態においても指向特性パターンの低ゲイン方 向を探索することなく、 式 （9') の演算により円錐面情報 ( . (ω),
Figure imgf000026_0001
j (ω)) を推定しているため計算量が少ない。 しかも同一信号源に対し、 円錐面 を複数推定してその共通直線から信号の到来方向を算出しているため、 0。 から 360。 までのどの方向に信号源が存在しても、信号源 3；の方向を一意に推定す ることができる。 またその推定方向をパ一ミュテーシヨン行列 P (ω) の決定に 利用しているため、 信号源の位置に関わらず、 パ一ミュテーシヨン問題を適切に 解決することができる。
〔第 3実施形態〕
第 3実施形態は、 位置情報として一対のセンサと一つの信号源との各距離の比 に基づくその信号源が存在する曲面を用いる。 第 1及び第 2実施形態においては 各信号源はセンサから遠方にあって、 信号源からの信号はセンサには平面波とし て到来することを想定して処理した。 しかし信号源とセンサとの距離が比較的短 かい場合はセンサには信号が球面波として到来する。 このような点から混合行列 A (ω) の要素の比 Α ω)ΖΑ」. ω)を球面波 （近距離場） モデルによって解 釈すると、 信号源 iの方向以外の情報を推定することができる。
すなわち、 近距離場モデルを用いると、 周波数応答 Αη(ω)は以下のように表 わせる。
A (ω) = (1/11 q; - dj il ) exp ( j o e—1 ( II q - d, I ))
ここで、 qiは信号源 iの位置を示すベクトルである。
このように表現された周波数応答について混合行列の同一列の 2つの要素比 A ji (ω) /Ay i (ω) をとり、 その絶対値を計算すると次式となる。
II Qi- dy 11 /11 qj-dj ll = I Aji (ω) /Ar ; (ω) I - ( 1 0) なお、 I 0 I は ι8 の絶対値を表わす。
この式（ 1 0)を満す q iの点の無数の集合は信号源 iが存在する曲面を与え、 遠距離場 （平面波） モデルを用いて推定された方向 （又は円錐面） とあわせて用 いることにより、 センサから信号源 iまでの距離を推定することができる。 これ により.， 2つ以上の信号源が同一方向もしくは互いに近い方向にある場合であつ ても、 距離が違っていれば それらを区別でき パーミュテーシヨン問題を適切 に解決することが可能となる。
この第 3実施形態の例では、 この信号源が存在する曲面を表わす位置情報と前 記方向情報とを用いて分離行列に対するパーミュテーシヨン問題を解決する。 第 3実施形態の機能構成を図 1 7に、 その処理手順を図 1 8にそれぞれ示す。 センサは 3個以上が 2次元又は 3次元に配置されるがこの実施形態においては例 えば図 2 0に示すようにセンサ 1 iと 12の間隔に対しセンサ 12と 13の間隔は 1 0〜20倍好ましくは 1 5倍程度とされる。 先の実施形態と同様、 各観測信号
( t ), …， j ( t ) は、 それぞれ周波数領域の信号 (ω, m) . --. Xj (ω, m) に変換され （ステップ S 1 1)、 更に独立成分分析により、 各周波数ごとに分 離行列 W (ω) が算出される （ステップ S 1 2)、 その各分離行列 W (ω) の逆行 列として行列 Η (ω) が算出される （ステップ S 1 3)。 またこの例では第 2実施 形態と同様に周波数ごとの逆行列 H (ω) の各列から選択した要素対を用い、 円 錐面が一つ好ましくは複数推定される （ステップ S 1 4)。 この第 3実施形態にお いては距離比算出部 3 1において、 周波数ごとの逆行列 Η (ω) から列ごとに選 択した要素対を用いその対応するセンサと 1つの信号源 i との距離比、 つまり式 (1 0) と式 （2 1 ) から次式 （1 0') を算出する （ステップ S 35)。
II qz(i)- dj- II / II qz(i)- dj 1 = I Ajz(i) (ω) / y z(i) (ω) 1 = 1 H (ω)
Figure imgf000028_0001
距離比算出部 3 1において行われるステップ S 35の処理の具体例を図 1 9を 参照して説明する。 この処理は図 1 4の処理の流れとほぼ同様である。 パラメ一 タ iを 0に初期化し （ステップ S 20)、 次に、 i を + 1 し （ステップ S 2 1 )、 パラメ一夕 j , j ' として J以下の自然数を、 例えばランダムに選択し j ≠ 、 かつ 1度選択した組は選択しない（ステップ S 23)。センサ jの位置べク トル d jと、センサ j ' の位置べク トル dj.とを抽出し （ステップ S 24)、逆行列 H (ω) の i列目から、 要素 Hji w)と、 Hr とを、 選択する (ステップ S 25)。 この実施形態ではこれら選択した 2つの要素の比 DRi, j (ω) を計算する （ス 'テツプ S 41)。 次に、 i = Iであるか否かを判断し （ステップ S 29)。 i = I でなければステップ S 2 1に戻り、 i =2であれば処理を終了する。
距離比算出部 3 1で算出された距離比情報 DRi, jr (ω) は パーミュテ一シ ョン解決部 1 7に送られ、 パーミュテ一ション解決部 1 7は、 到来方向决定部 1 6で推定された方向情報 U i ( ω ) と 距離比算出部 3 1で算出された距離比情報 DRi, jj. (ω) とを用いて、 分離行列算出部 1 2で算出された分離行列を行って パ一ミュテ一ション問題を解決する。
W (ω) の行替えを行ってパーミュテーシヨンの問題を解決する。 例えば方向 情報と距離比情報とから信号源 2iの距離 II q i (ω) IIを距離推定部 1 7 dで計 算する （ステップ S 36)。
この距離 II qs (ω) IIの計算方法を図 20を参照して説明する。 信号源 2,及 ぴ 22がセンサ 1ぃ 12からみて同一方向 Βにある。 この場合センサ 1 !， 12及び 遠距離場モデルで推定される信号源 2,及び 22の方向 U2は同一の直線を定 める。 一方、 間隔が大きいセンサ 12， 13及び近距離場モデルを用いて信号源 2 i が存在すると曲面 6 ,は距離比
DR,.23 (ω) = I Η21 (ω) /Η31 (ω) I = II q「 d3|IZII q1-d2ll から推定でき、 つまり II q, (ω) IIが推定できる。 また信号源 22が存在する曲 面 62は距離比
DR2, 23 (ω) = I Η22 (ω) /Η32 (ω) I = 11 q2-d3II II q2-d2ll から推定でき、 つまり II q2 (ω) IIが推定できる。
信号源 2!の位置は、 直線 u i^と、 曲面 6, との共通部分上にあると推定で き、 信号源 22の位置は、 直線 U l=u2と、 曲面 62との共通部分上にあると推定 できる。 例えば直線 1^= 1^を表す式と曲面 6i及び 62をそれぞれ表わす式との 連立方程式を解くことにより II (ω) II及び ll q2 (W) ii を求めればよい。 このようにして信号源の方向が同一や近似する場合であっても、 信号源位置を区 別することができる。
パーミュテーション解決部 1 7は、 このようにして得られた各周波数ごとの信 号源の距離 II qi (ω) IIが、 所定の順序が例えば昇順になるように、 分離行列 W (ω) の行の入れ替えを行う。 つまりこの入れ替えを行うためパーミュテーショ ン行列 Ρ (ω)を算出する (ステップ S 37)。そのパーミュテ一ション行列 Ρ (ω) の算出は実施形態 2中のパ一ミュテーシヨン解決部で行ったパーミュテーシヨン 行列の算出と同様に行えばよい。 算出したパーミ ュテーション行列 Ρ (ω) を分 離行列 W (ω) の左から掛け、 その結果 W' (ω) =Ρ (ω) W (ω) を分離行列 として出力する (ステップ S 38)。
出力された分離行列 W' (ω) は、 時間領域変換部 1 8に送られ、 時間領域変換 部 1 8での信号分離に用いられる。
図 20から理解されるように、 センサ 12と 13についてみると、 信号源 2 ,及び センサ 12間の距離と信号源 2 i及びセンサ 13間の距離との差が大きく、一方信号 源 22及びセンサ 12間の距離と信号源 22及びセンサ 13間の距離との差が小さい。 従って D RL 23= II q!-dgll/ll Q!-dJIと数値 1 との差の絶対値は比較的大 きく、 DR2, 23= I q2-d3IIXII q2-d2|| と数値 1 との差の絶対値は小さい。 センサ 12と 13の間隔が大きい程、 距離比 DRL 23 (ω) と D R2. 23 (ω) との差 が大きくなる。 このため間隔が大きいセンサ対、 この例では 12と 13が設けられ る。
パーミュテーション解決部 1 7において、 図 1 8中のステップ S 3 7中に括弧 で示すように計算した距離比 DRi, （ω) を用いて、 すべての周波数において 逆行列 Η (ω) の 1列目， 2列目， 3歹 U目， …， I列目の順に、 求まる距離比 D R jT (ω) が昇順になるように並替行列を生成し、 更に前記パーミュテーショ ン行列 Ρ (ω) を生成してもよい。 この場合、 図 1 7中の距離推定部 1 7 dは省 略される。 また到来方向決定部 1 6で方向情報 Ui (ω) を計算できなかったり、 逆行列 Η (ω) の複数列について同一の Ui (ω) が計算された信号源 iとセンサ i j 'についてのみ、 距離比 DRi, ， （ω) を計算し、 これら自体により、 あるい は更に求めた II qt (ω) IIにより分離行列 W (ω) の行の入れ替えを行ってもよ レ、。つまり、方向情報 Ui (ω) を用いて分離行列 W (ω) の行の入れ替えを行い、 その後、 Uj (ω) によって入れ替えができなかった行については距離比 DRi, (ω) 又は II (ω) ||を用いて W (ω) の行の入れ替えを更に行ってもよい。 実際には両入れ替えを同時に行うパーミュテーシヨン行列 Ρ (ω) を生成する。 距離比 DRi, j (ω) をまず求め、 分離行列 W (ω) の行の入れ替えを行い、 距離比 DRi, ， (ω)を求めることができなかったものについて方向情報 Ui (ω) を用いて分離行列 W (ω) の入れ替えを更に行ってもよい。 この場合も両入れ替 えを同時に行うパーミュテーシヨン行列 Ρ (ω) を生成する。 一般には DRi, j (ω) より Ui (ω) の方が高い精度で得られるため、 Ui (ω) による入れ替え を主とし、できないものについて DRi, jr (ω)による入れ替えを行う方がよい。 また式 （1 0') をその右辺の値を距離比 DRi, j (ω) と表わして、 ciiにつ いて解くと、 中心の ， (ω) 及び半径 Ri, j (ω) がそれぞれ下記式 ( 1 1 ) 及び ( 1 2) で与えられる球面となる。
O j (ω)
Figure imgf000030_0001
i}. (ω) — 1) … （11)
R .， （ω) = II DRi, ' （ ω ) · ( d — d」） Z (D R2 j (ω) 一 1 ) ||
… (12) 例えばセンサ 1 1 rの各位置が d j= ( 0， 0. 1 5, 0 ) , d = ( 0， ― 0. 1 5 , 0) (単位はメートル） である場合、 半径 Ri, をパラメータとした とき、 式 （9) によって決まる球面の様子は図 2 1に示すようになる。 つまり信号源 iは位置情報としての式 （1 1 ) 及び （1 2) で与えられる球面 上に存在することになる。 従って図 1 7中のパーミュテーション解決部 1 7中の 距離推定部 1 Ί dを括弧書きで示すように曲面推定部とし、 方向情報 Ui (ω) を 求めることができなかった各 i , j j ' について曲面推定部 1 7 dで、 式 （1 1 ) の半径 Ri, j (ω) 及び中心 Oi, （ω) を、 図 1 8中のステップ S 3 6中に括 弧書きで示すようにそれぞれ計算し、 これら Ri, j (ω) 及び中心 Oi, j (ω) がどの周波数の逆行列 Η (ω) でも同一順になるようにして、 パーミュテ一ショ ン行列 Ρ (ω) で求めてもよい。 .
この円錐面情報 Θ ，』」·. (ω) 又は方向情報 Ui (ω) と球面情報 DRi,」Ί. (ω) 又は Ri. jj. (ω)を用いてもパーミュテーシヨン行列 P (ω) を求めることができ ない場合はその周波数について、従来の相関法（例えば、 H. Sawada他 "A robust and precise method for solving the permutation problem of frequency-domain bl ind source separation. ", in Proc. Intl. Symp. on Independent Component Analysis and Blind Signal Separation (ICA 2003) , 2003, pp.505- 510参照) を適用すると よい。 図 1 7中に破線で示すように、 信号分離部 1 9で分離された出力信号 y i
(t) y T (t ) を周波数領域変換部 33で例えば周波数領域変換部 1 1 と 同様の方法で周波数領域の信号 Υ, (ω' m), ···, Y j (ω, m) に変換し 相関 算出部 34において、 これら周波数領域信号 (ω, m)， ■··, Y , (ω, m) 中 のパーミ ュテ一ション解決部 1 7でパーミュテ一ション行列を作ることができな かった周波数成分 f ank成分と、 パーミュテ一シヨン行列が得られた周波数中の f Mkと隣接した周波数成分との相関を計算し、 パーミュテーション解決部 1 7では この相関が大きくなるように周波数 f ankの分離行列の行の入れ替えを行い、 この 行入れ替えした分離行列に基づく分離出力信号 (t) y！ (t ) 中の f ank 成分について再び相関計算部 34で相関を計算し、 同様のことを計算した相関が 最大になるまで繰り返す。 この相関の最大値が所定値に達しなかった場合は、 更 にその f ankの成分と、 パーミュテーション行列が得られている周波数成分中の f ankに対する基本波又は高調波（一般には高調波）関係にあるものとの'各相関の和を 相関計算部 34で計算し、 この相関和が大きくなるように、 f ankに対する分離行 列の行入れ替えを修正部 1 7 eで行うことを、 前記相関和が最大になるまで繰り 返す。 なお第 2実施形態で説明したように周波数領域の信号 X (ω, m) に対し 分離行列 W (ω) を適用し、 周波数領域の分離信号 Υ (ω, m) を求める場合は、 この Υ (ω, m) を相関計算部 34の計算に用いればよい。
以下にこの第 3実施形態の実験例を信号源として、 室内で実測したィンパルス 応答を畳み込んで作成した混合音声を用いて分離実験を行った場合について説明 する。 図 22に示すように残響時間が 1 30ミリ秒、 355 cmx445 cmx2 50 cm (高さ） の室内の水平面内でほぼ中心部で高さ 1 35 cmの所にセンサ としての無指向性マイクロホン l i〜 l4を直径 4 cmの円上に等間隔で配し （こ れらマイクロホン 1 i〜 14は図中の左上に拡大して示した）、更に無指向性マイク 口ホン 15〜 18を直径 30 c mの円上に等間隔かつ、 マイクロホン 1い 13， 15 及び 17が、 またマイクロホン 12， 14, 16及び 18がそれぞれ一直線上に配列す るように配置した。 マイクロホン 1 iから 15への方向を基準 (0°) とし、 マイク 口ホン配置の中心を原点とし、 原点から 1 20 cm離れ 一 30°方向に音源 2i を、 + 30°方向に音源22を、 + 90°方向に音源 23を、 一 1 50°方向に音源 2 6をそれぞれ配置し、 1 50°方向で原点からの距離 60 cmに音源 24を、 距離 1 80 cmに音源 25をそれぞれ配置した。サンプリング周波数 kH z、 データ長 6 秒 フ レーム長 2048サンプル （ 256ミ リ秒） フ レームシフ ト 5 1 2サンプ ル （64ミ リ秒） とした。
周波数領域での分離行列 W (ω) の各行のうち、 マイクロホン対 1 ,と 13、 12 と 14、 1！と 12、 12と 13に相当する行を用いて音源方向 （円錐面） 推定を行つ た。 推定したヒストグラムを図 22に示す。 図 23の横軸は推定方向を、 縦軸は 信号数である。 これより推定方向は 5つの群 (クラスタ） があり、 そのうち 1 5 0。付近のクラスタは他のクラスタより 2倍の大きさである。これは 6個の音源中 の 2個が同一方向 （1 50°付近） から到来していると推測される。 4つの音源に ついては推定した方向をもとにパーミュテ一シヨンを解決した。 その結果を図 2 4に示す。 図 24の横軸は周波数、 縦軸は方向 （度） を表わす。
残りの 2個の音源については、広い間隔のマイクロホン対 17と 15、 17と 18、 16と 15、 16と 18を用いて音源が存在する球面の半径によって到来信号を区別 した。 マイクロホン対 17と 15を用いて推定した球面の半径を図 25に示す。 図 25の横軸は周波数、 縦軸は半径 （m) を表わす。
残響や推定誤差の影響により、 位置情報だけではパーミュテーションを完全に は解決することはできない。 従って推定した位置情報に基づいて信号が矛盾なく 分類できた周波数はその情報に基づいてパ一ミュテーション行列を生成し、 残り の周波数については相関に基づく方法を用いてパーミュテーションを解決した。 最後に時間領域の分離フィルタ係数を求める際にスぺクトルの平滑化を行った。 スぺク トル平滑化については例えば H. Sawada 他、 "Spectral smoothing for frequency-domain blind source separation" , in Proc. IWAENC 2003, 2003, PP.311- 314参照されたい。 分離性能 （S I R) の評価結果を図 26に示す。 図 2 6中の数値の単位は dBである。 Cは相関法のみによってパーミュテーシヨンを 解決したもの、 D + Cは方向 （円錐面） 推定でパーミュテーシヨンを解決し、 解 決できないものに相関法を適用したもの D + S + Cは方向 （円錐面） 推定と球面 推定とによりパーミュテーション解決を行った後、 解決できなかつた周波数につ いて、 相関法を用いたものである。 この後者の方法によれば、 6音源すベての分 離ができ、 S I Rの改善は入力 S I Rから平均 1 7. l dBとなった。
前記した第 2実施形態の例ではセンサは 2次元配置したが、 一対のセンサによ り推定する球面はこれらセンサの中心 2等分線と対称に現われるため 信号源が 3次元に存在している場合は 2次元配置センサでは判別できなくなるため セン サも 3次元配置とする必要がある。
以上説明したように、 この第 3実施形態においても式 （9') により円錐面情報 を、また式（1 0 ')により曲面情報をそれぞれ推定しているため計算量が少ない。 しかも円錐面と、 距離比 DRi,」』. (ω)、 又は距離 II qi (ω) || あるいは球面の 半径 Ri. (ω) とによりパーミュテーションを解決しているために 2つ以上の 信号源が同一方向もしくは互いに近い方向にある場合であっても、 それらを区別 できる。 更にこれに相関法を加えれば一層確実に分離が可能である。 なお球面情 報としては計算が簡単であるから DRi. j (ω) が好ましい。
〔第 4実施形態〕
この第 4実施形態では、 推定した円錐面の信頼性を検証し、 信頼性が高いと判 断した円錐面を分離行列のパーミュテーションの解決に用いる。 例えば図 27に 示すように円錐面推定部 14で推定された円錐面情報 」Ί. (ω) は円錐面検 証部 4 1で許容角情報格納部 42よりの許容角情報に基づき、 信頼性があるか否 かの検証がなされる。 つまり式 （9') により求まる角度 S Λ (ω) はこれを 求めるために用いたセンサ 1』及び 1」.の配列方向に対する相対角度であり、図 4 Bを参照して説明したように、 角度 Λ i, （ω) が 0度付近及び 1 80度付近 の場合、 正しいパーミュテーションが行われないと推定される。
従って、 パーミュテーシヨンを正しく行うことができると推定される角度の最 小値 Θ minと最大値 Θ mxが許容角情報として許容角情報格納部 42に格納され、 円 錐面検証部 4 1で推定円錐面情報 S Λ j (ω) が <9ninと naxとの間にあれば信 頼できる円錐面と判定されて出力され、 つまりパーミ ュテーシヨン解決に用いら れるが、 ^minと Snaxとの間になければ、 その ， j (ω) は信頼性がないもの として破棄され-. パーミュテーシヨン解決に用いられない。 例えば図 1 5中の円 錐面 413は破棄される。
円錐面検証部 4 1で信頼性があると検証された円錐面情報 .」」. (ω) は図 1 2中の到来方向決定部 1 6へ送られ、 又は図 1 7で説明したように、 パーミュ テーション解決部 1 7へ直接送られる。 図 1 3及び図 1 8の処理手順において、 ステップ S 1 4の次に破線で示すようにステップ S 1 00として推定円錐面 θ Λ ；, jr (ω) が信頼性があるか否かの検証が行われ、 信頼性があるものについての み次のステップに移るようにされる。 式 （9') において計算された偏角 arg (H ノ H i) に含まれる誤差を AargH 、 推定された角度 ；に含まれる誤差を Δ Θ とすると、 これらの比 I Δ Θ ソ AargH Iは式 (9') を偏微分して次 式のように近似できる。
I Δ Θ ソ AargH Λ 1 = 1 1/ (ω c"1 I d「 d I sin^ ） I (13) 式 （1 3) をいくつかの周波数について計算した結果を図 4 Cに示す。 この図 4 Cから推定された角度 < の方向がセンサ lj, 1」. のセンサペア軸に近い場 合は arg (Hji/Hj. ；) に含まれる誤差 Δ argH Λは推定角度 Λ i に対し、 大き な誤差を生じさせることがわかる。 つまりセンサペア軸に近い方向の推定角度 Θ を用いて分離行列 W ( ω )に対するパーミュテーシヨンの問題を解決した場合、 これは正しい解にならない可能性が高い。 図 4 Β及び図 4 Cから例えば^ ninは 2 0° 程度、 Smaxは 1 60° 程度とされる。 図 4 Cからわかるように I Θ / argH Λ Iは周波数により可成り異なり、 低い周波数では AargH Λが到来方向の 誤差に大きく影響する。 従って、 低周波数については、 信号源が存在する球面の 推定に基づく情報 DRi, j 、 II q; Ik R;. j を用いて、 あるいは相関法によりパ ーミュテーシヨン問題を解決するとよい。
この第 4実施形態によれば、 推定円錐面情報中の信頼性のないものは破棄され るため、 これにより悪影響を受けることなく、 それだけ正確に到来方向を推定で き、 従って正しいパーミュテーシヨン行列 P (ω) を生成することができ、 つま り分離信号の S I R (性能） が向上する。
〔第 5実施形態〕
第 5実施形態では位置情報として距離比、 又はこれより推定した球面情報を用 いる。 この機能構成は図 1 7中に示されている。 図 28にこの実施形態の処理手 順を簡略に示す。 この場合はセンサは比較的広い間隔、 例えば図 22に示した信 号源が音源の場合、 30 cmとして少くとも 2次元に配置される。
先の実施形態と同様に時間領域の観測信号は図 1 7中に示すように周波数領域 の信号に変換され、 分離行列生成部 1 2により分離行列 W (ω) が生成され （ス テツプ S 5 1) これより逆行列生成部 1 3で逆行列 Η ( ») が生成される （ステ ップ S 52)。 各周波数における逆行列 H (ω) の各列ごとに球面情報が推定され る （ステップ S 53)。 この球面情報は第 3実施形態で説明した方法と同様に算出 される。つまり球面情報は距離算出部 3 1により算出された距離比 DRi. (ω)、 あるいは距離推定部又は曲面推定部 1 7 dで算出された II q; IL又は半径 Ri, (ω) 及ぴ中心 Oi. (ω) である。
これら球面情報を用いて、 その配列順が予め決めた順になるよう分離行列 W (ω) に対するパーミュテ一シヨン行列が生成され分離行列 W (ω) の行の入れ 替えが行われる （ステップ S 54)。第 3実施形態におけるこの処理はパーミュテ ーシヨン解決部 1 7で行われるが、 ただしこの第 5実施形態では球面情報のみが 用いられる。 この球面情報のみでパ一ミュテーションを解決できない周波数があ れば、 その周波数に対しては前述した相関法による解決が行われる （ステップ S 55)。 図 22に示した室内に配したマイクロホン 16及び18を用ぃ、 1 20° 方向に おいて原点から 60 cm及ぴ 1 50 c mの所に音源 24及び 25を配し、 実験室内 で測定したインパルス応答に話者 4名の音声を畳み込んで作成した 1 2通りの組 合せの混合音声と対する分離実験を行った。 パ一ミュテーション解決は推定され た半径 R4, 68 (ω) と R5. 68 (ω) を比較して R4, 68 (ω) ≤R5, 68 (ω) となるよ うにして行った。 R4, 68及び R5, 68中の最大値が、 R4, 68及び R5, 68中の最小値に しきい値 Athを乗算した値以上、 つまり ma X (R4, 68, R5, 68) ≥Ath · m i n (R4, 68- R5, 68) を満す周波数は音源位置 （R4, eg (ω), R5.68 (ft))) によるパ ーミュテ一ションの解決が信頼できると判断して、音源位置による方法を適用し、 それ以外の周波数では相関による方法を適用した。 従って Ath= 1. 0のときは 全て音源位置による方法を適用し A th が無限大のときは全て相関による方法を 適用することになる。 しきい値 Athを変化させながら分離性能として S I R (信 号干渉比) を音声の組合せごとにプロッ ト した結果を図 29に示す。 この結果よ り音源位置のみによる方法では全体的に性能が悪く、 相関のみによる方法では分 離性能がばらつき、両者を併用すると安定して高い性能が得られることがわかる。 またしきい値 Athとしては比較的大きい値が好ましく、 8〜 1 6で全体の 1 /5 〜 1 /1 0程度の周波数を音源位置によって決定できることがわかる。
また図 22に示した実験条件で音源位置による方法と相関による方法を併用し た実験結果を図 26中の D + Cの欄に示す。 これよりこの併用方法によれば可成 り性能よく分離が行われることがわかる。 この結果は式 （9 ') を用いる到来方向 (円錐面情報） による方法と相関による方法とを併用した場合と同様の傾向があ る。なお、到来方向による方法と相関による方法とを併用した場合の実験結果は、 相関による方法の説明で引用した文献に示されている。
この第 5実施形態によれば式（ 1 0 ') により球面情報を求めているため計算量 が少ない。 なお球面情報としては距離比 DRi.』」. (ω) が好ましい。
[第 6実施形態]
第 6実施形態は例えば一つの推定円錐面に基づき、 パーミュテ一ションの問題 を解決する。 図 1 2中に破線で示すように円錐面推定部 1 4で推定された円錐面 θ Λ ;, ' (ω α), · · · θ Λ i, jj. (ωΝ)、 はパーミュテ一ション解決部 1 7へ 直接入力され、 図 1 3中に破線で示すようにステップ s 1 4で設定された円錐面 θ ^ (, jj. ( ω η) ( η = 1、 · · ·、 Ν) がステップ S 1 5の処理をすることなくス テツプ S 1 6でどの周波数においても例えば Θ ^ j ( ω η) が昇順になるような パーミュテーシヨン行列 Ρ ( ω ) が算出される。
この場合、 ステップ S 1 4で信号源 i ごとに 1つの円錐面を推定するのみでも よい。 また図に示していないが、 行入れ替えができなかったものについては、 第 3実施形態で説明したように相関法による行入れ替えを得るようにしてもよい。 この第 6実施形態によれば逆行列 H ( ω ) の各列ごとの二つの要素比をとるこ とによりスケーリングの問題が簡単に除去され、 式 （9 ' ) の演算をすればよく計 算時間が短くて済む。
[信号分離のまとめ]
以上述べたブラインド信号分離におけるパーミュテーションを解決した分離行 列の生成方法を、 逆行列の生成以後の処理を図 3 0 Α、 図 3 0 Β、 図 3 0 Cに示 す。 図 3 O Aでは逆行列の各列の要素比から円錐面を推定し (ステップ S 6 1 )、 必要に応じて信頼性のない円錐面を破棄し （ステップ S 6 2 )、 これらの複数の円 錐面から共通直線の方向を決定し （ステップ S 6 3 )、共通直線方向を用いてパー ミュテーシヨン行列 P ( ω ) を生成して分離行列の入れ替えを行い （ステップ S 6 4 ) ーミュテーションを解決できなかった周波数に対しては相関法による分 離行列の行入れ替えを行う (ステップ S 6 5 )。 図 3 O A中に破線で示すように、 ステップ S 6からステップ S 6 4に直ちに移り、 円錐面を直接用いてパーミュテ ーション行列 P ( ω ) を生成してもよい。
図 3 0 Βでは前記ステップ S 6 1の円錐面推定、 必要に応じて前記ステップ S 6 2の円錐面破棄を行い、 ステップ S 6 3の共通直線方向推定を行い、 これら共 通直線方向を用いて逆行列の列の入れ替えを行う （ステップ S 6 6 )。 この際、 図 中破線で示すようにステップ S 6 1で推定された円錐面又はステップ S 6 2で処 理後の円錐面を直接、 ステップ S 6 6の処理に用いてもよい。 次に円錐面の推定 又は共通直線方向の決定ができなかった、 又は円錐面或いは共通直線方向が不確 かなものと対応する、 或いは同一値となったものについて、 その逆行列の列の要 素比から球面情報を推定し（ステップ S 6 7 )、 この推定した球面情報を用いて逆 行列の入れ替えを行ってパーミュテーシヨン行列 P (ω) を生成して分離行列の 行入れ替えを行う （ステップ S 64)。 更にこのパーミュテーシヨン行列 Ρ (ω) を作ることができなかった周波数に対しては相関による方法を適用する （ステツ プ S 65)。
図 30 Cはまず逆行列の各列の要素比から球面情報を推定し（ステップ S 6 8)、 この球面情報を用いて逆行列の列の入れ替えを行い（ステップ S 69)、 この列の 入れ替えを行うことができなかった、 又は球面情報が不確かなものについて、 そ の逆行列の列の要素比から円錐面を推定し （ステップ S 70)、 これらの複数の円 錐面の共通直線方向を決定し （ステップ S 7 1)、 この決定した方向又はステップ S 70で推定した円錐面を直接用いて、 逆行列の列の入れ替えを行ってパーミュ テーション行列 Ρ (ω)を生成し、分離行列の入れ替えを行う (ステップ S 64)。 パーミュテーシ aン行列を作ることができなかった周波数に対して相関法を適用 する (ステップ S 65)。
これらに対し更に図 28に示した方法がある。
第 3実施形態、 第 5実施形態及び第 6実施形態においても、 第 2実施形態で説 明したように、 周波数領域の分離行列 W (ω) と観測信号 X (ω) とを用いて信 号分離を行い その後、分離された周波数領域信号 Υ (ω)を時間領域信号 y ( t ) に変換してもよい。
第 1実施形態の説明では 2次元における信号の到来方向の推定を行ったが、 第 2実施形態中でも説明したように、 3次元における信号の到来方向の推定にも適 用できる。 また第 2乃至第 6実施形態を 2次元における信号分離に適用すること もできる。 この場合は推定円錐面 <9 (ω) はその推定に用いたセンサ対の センサペア軸に対し対称の 2方向となり、推定球面 Ri, jj. (ω)又は DRi. j (ω) は円の半径又はこれとほぼ対応するものになる。
図 5、 図 1 2、 図 1 7にそれぞれ示した装置、 また第 5実施形態による装置は それぞれコンピュータに機能させることもできる。 この場合、 それぞれと対応し た処理手順、 つまり図 6、 図 1 3、 図 1 8などに示した各過程をコンピュータに 実行させるためのプログラムを記録した磁気ディスク、 半導体メモリ、 CD— R OMなどから、 コンピュータ内のメモリにインストールし、 又は通信回線を介し てコンピュータ内のメモリに前記プログラムをダウンロードして、 コンピュータ にそのプログラムを実行させればよい。

Claims

請求の範囲
1 . I箇所の信号源から放射された信号を、 J個のセンサで検出して上記信号源 の位置情報を求める信号源位置情報推定装置であって、 Iは 2以上の整数、 Jは I以上の整数であり、
上記各センサの観測信号を周波数領域の信号に変換する周波数領域変換手段と、 上記周波数領域の信号から、 独立成分分析により周波数ごとに、 信号源信号を 分離する第 1分離行列をそれぞれ算出する分離行列算出手段と、
上記各第 1分離行列の逆行列をあるいは擬似逆行列 「以下、 両者を単に逆行列と いう」 をそれぞれ算出する逆行列算出手段と、
上記周波数ごとの逆行列の少なくとも 1つに対し、 その各列ごとにその二つの 要素比から上記信号源の一つの位置情報を計算する位置情報計算手段とを具備す る。
2 . 請求の範囲第 1項記載の装置において、
上記位置情報計算手段は、 上記周波数ごとの逆行列の複数について、 上記各列 ごとの要素比からの位置情報計算により上記各信号源の位置情報を求める手段で あり
上記周波数ごとの信号源ごとの位置情報から、 上記複数の逆行列と対応する周 波数の分離行列における位置情報と対応する行が予め决めた順序になるように行 の入れ替えを行うパーミュテーション行列を生成するパーミュテーション行列生 成手段と、
上記パーミュテーション行列と上記第 1の分離行列とを乗算して行を入れ替え た第 2の分離行列を求める並替手段とを備える。
3 . 請求の範囲第 2項記載の装置において、 '
上記 Jは 3以上であり、 かつ J個のセンサは少なくとも 2次元に配置され、 上記位置情報は上記センサから上記信号源への方向を含み、 その信号源が存在 する円錐面であり、
上記位置情報計算手段は、 各列ごとの上記要素比から上記円錐面の計算を異 なる二つの要素組の複数について行う手段と、 周波数ごとの上記複数の円錐面の 共通直線の方向を、 上記位置情報として推定する到来方向決定手段を含む。
4. 請求の範囲第 2項記載の装置において、
上記 Jは 3以上であり、 上記位置情報は一対のセンサから上記信号源への方向 を含み、 その信号源が存在する円錐面及びその信号源が存在する曲面であり、 上記位置情報計算手段は、 上記二つの要素の比から上記円錐面を計算する手段 と、 上記円錐面の計算に用いた二つの要素と対応するセンサ間隔よりも大きいセ ンサ間隔のセンサと対応する二つの要素の比からこれらセンサと上記信号線との 距離比を計算する手段と、 上記距離比から上記曲面を計算する手段を備え、 上記パーミュテーション行列生成手段は周波数ごとに上記円錐面及び上記曲面 に基づいて上記パーミュテ一ション行列を生成する手段である。
5 . 請求の範囲第 2項、 第 3項又は第 4項に記載の装置において、
上記 Jは 3以上であり-. 上記 J個のセンサは少なくとも 2次元に配置され-, 上記位置情報は上記センサから上記信号源への方向を含み、 その信号源が存在 する円錐面であり、
上記円錐面を表わす角度が予め決めた第 1角度と第 2角度との間であるか否 かを判定し、 これら第 1角度及び第 2角度の間の円錐面を有効とする判定手段を 備える。
6 . 請求の範囲第 2項記載の装置において
上記 Jは 3以上であり、 上記 J個のセンサは少なくとも 2次元に配置され、 上記位置情報は信号源が存在する球面の半径であり、
上記位置情報計算手段は上記二つの要素比から距離比を計算する手段である。
7 . 請求の範囲第 2項に記載の装置において、
上記位置情報は上記センサから上記信号源への方向を含みその信号源が存在 する円錐面である。
8 . 請求の範囲第 2項、 第 3項、 第 4項、 第 6項又は第 7項に記載の装置におい て、
上記観測信号に対し、 上記第 2分離行列を用いて分離した信号の周波数領域の 信号中の、 上記パーミュテーシヨン行列生成手段で、 パ一ミュテーシヨン行列が 得られない周波数成分と、 パーミュテーション行列が得られた周波数成分との相 関を計算する相関計算手段と、
上記相関が大きくなるように上記パーミュテーション行列が得られない周波数 ' の分離行列の行を入れ替えるパーミュテーション行列を生成する修正手段とを備 んる。
9 . 請求の範囲第 1項記載の装置において、
上記位置情報は上記センサから上記信号源への方向情報であり、
上記位置情報計算手段は、 上記比の偏角を計算する手段と、 単位距離当たりの 位相回転量と上記二つの要素とそれぞれ対応する上記センサ間の距離との積を計 算する手段と、 上記偏角を上記積の結果で除算する手段と、 その除算結果の逆コ サインを計算して上記方向情報を出力する手段とを備える。
1 0 . 請求の範囲第 9項記載の装置において、
上記位置情報計算手段は、 一つの信号源の上記方向情報を、 周波数ごとに計算 する手段であり、 上記各一つの信号源ごとに計算された周波数ごとの方向情報を 統合して方向情報を確定する統合手段を備える。
1 1 . I箇所の信号源から放射された信号を、 J個のセンサで検出して上記信号 源の位置情報を求める信号源位置情報推定方法であって、 Iは 2以上の整数、 J は I以上の整数であり
上記各センサの観測信号を周波数領域の信号に変換する周波数領域変換過程と 上記周波数領域の信号から、 独立成分分析により周波数ごとに、 信号源信号を 分離する第 1分離行列をそれぞれ算出する分離行列算出過程と、
上記各第 1分離行列の逆行列あるいは擬似逆行列（以下、両者を単に逆行逆列と いう）をそれぞれ算出する逆行列算出過程と、
上記周波数ごとの逆行列の少なくとも 1つにおける列ごとの異なる二つの要素 比から上記信号源の一つの位置情報を計算する位置情報計算過程とを有する。
1 2 . 請求の範囲第 1 1記載の方法において、
上記位置情報計算過程は、 上記周波数ごとの逆行列の複数について、 各列ごと の上記要素比からの位置情報計算により上記各信号源の位置情報を求める過程で あり、
上記周波数ごとの信号源ごとの位置情報から、 上記複数の逆行列と対応する分 離行列における位置情報と対応する行が予め決めた順序になるように行の入れ替 えを行うパ一ミュテーション行列を生成するパ一ミュテ一ション行列生成過程と、 上記パーミュテーシヨン行列と上記第 1の分離行列とを乗算して行を入れ替え た第 2分離行列を求める分離行列置換過程とを有する。
1 3 . 請求の範囲第 1 2項記載の方法において、
J≥ 3であり、 かつ J個のセンサは少なくとも 2次元に配置され、
上記位置情報は上記センサから上記信号源への方向を含み、
その信号源が存在する円錐面であり、
上記位置情報計算過程は、 各列ごとの上記要素比から上記円錐面の計算を異な る二つの要素組の複数について行う過程であり、 周波数ごとの上記複数の円錐面 の共通直線の方向を、 上記位置情報として推定する共通線推定過程を含む。
1 4 . 請求の範囲第 1 2項記載の方法において、
J≥ 3であり、 上記位置情報は一対のセンサから上記信号源への方向を含み その信号源が存在する円錐面及びその信号源が存在する曲面であり、
上記位置情報計算過程は、 上記 2つの要素の比から上記円錐面を計算する過程 と、 上記円錐面の計算に用いた 2つの要素と対応する一対のセンサの間隔より大 きい間隔の一対のセンサと対応する 2つの要素の比からこれら一対のセンサと一 つの信号源との距離比を計算する過程と 上記距離比から上記曲面を計算する過 程を有し、
上記パーミュテーション行列生成過程は周波数ごとに上記円錐面及び上記球面 に基づいて上記パーミュテーシ 3ン行列を生成する過程である。
1 5 . 請求の範囲第 1 4項記載の方法において、
上記パ―ミュテーション行列生成過程は、 円錐面及び上記球面の一方について 上記位置情報計算過程を実行させ、 得られた上記一方の位置情報に基づいて、 各 周波数での逆行列の各列と対応する上記一方の位置情報が予め決めた順になるよ うに列の入れ替えを行う第 1の入替行列を生成する過程と、
この過程により、 列の入れ替えを行うことができない列について上記位置情報 の他方について上記位置情報計算過程を実行させ、 得られた上記他方の位置情報 に基づいて上記逆行列の列の入れ替えを行うように上記第 1の入替行列を修正し て第 2の入替行列を生成する過程と、 上記第 2の入替行列の逆行列を計算して上 記パーミュテ一ションとを有する。
1 6 . 請求の範囲第 1 2項乃至第 1 5項のいずれかに記載の方法において、
J 3であり、 上記 J個のセンサは少なくとも 2次元に配置され、
上記位置情報は上記センサから上記信号源への方向を含み、 その信号源が存在 する円錐面であり、
上記円錐面を表わす角度が予め決めた第 1角度と第 2角度との間であるか否か を判定し、 これら角度の間にない円錐面を破棄する判定過程を有する。
1 7 . 請求の範囲第 1 2項に記載の方法において、
J≥ 3であり、 上記 J個のセンサは少なくとも 2次元に配置され、
上記位置情報は信号源が存在する球面の半径であり、
上記位置情報計算過程は上記二つの要素比から距離比を計算する過程である。
1 8. 請求の範囲第 1 2項に記載の方法において、
上記位置情報は上記センサから上記信号源への方向を含み、 その信号源が存在 する円錐面である。
1 9 . 請求の範囲第 1 2項〜第 1 5項のいずれか又は第 1 7項又は第 1 8項に記 載の方法において
上記パーミュテーション行列生成過程でパーミ ュテーション行列を生成できな い周波数があれば、
上記観測信号を上記第 2分離行列で分離した信号の周波数領域の信号中の上記 パーミュテ一ション行列を生成できた周波数成分と上記パーミ ュテーション行列 を生成できない周波数成分との相関を計算する過程と、
この計算した相関が大きくなるように、 上記生成できなかった周波数の分離行 列に対するパーミュテーション行列を生成する過程とを有する。
2 0 . 請求の範囲第 1 1項記載の方法において、
上記位置情報は上記センサから上記信号源への方向情報であり、
上記位置情報計算過程は、 上記比の偏角を、 単位距離当たりの位相回転量と上 記二つの要素とそれぞれ対応する上記センサ間の距離との積で除算し、 その除算 結果の逆コサインを計算して上記方向情報を出力する過程である。
2 1 . 請求の範囲第 2 0項記載の方法において、
上記位置情報計算過程は、 一つの信号源の上記方向情報を、 周波数ごとに計算 する過程であり、 上記各一つの信号源ごとに計算された周波数ごとの方向情報を 統合して方向情報を確定する統合過程を有する。
2 2 . 請求の範囲第 1 1項乃至第 1 9項のいずれかに記載した信号源位置情報推 定方法の各過程をコンピュータに実行させるためのプログラム。
PCT/JP2004/002610 2003-03-04 2004-03-03 位置情報推定装置、その方法、及びプログラム WO2004079388A1 (ja)

Priority Applications (4)

Application Number Priority Date Filing Date Title
US10/508,708 US7039546B2 (en) 2003-03-04 2004-03-03 Position information estimation device, method thereof, and program
EP04716703A EP1600789B1 (en) 2003-03-04 2004-03-03 Position information estimation device, method thereof, and program
DE602004029867T DE602004029867D1 (de) 2003-03-04 2004-03-03 Positionsinformationsschätzeinrichtung, verfahren dafür und programm
JP2005503053A JP3881367B2 (ja) 2003-03-04 2004-03-03 位置情報推定装置、その方法、及びプログラム

Applications Claiming Priority (4)

Application Number Priority Date Filing Date Title
JP2003-057070 2003-03-04
JP2003057070 2003-03-04
JP2003-297580 2003-08-21
JP2003297580 2003-08-21

Publications (1)

Publication Number Publication Date
WO2004079388A1 true WO2004079388A1 (ja) 2004-09-16

Family

ID=32964882

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/JP2004/002610 WO2004079388A1 (ja) 2003-03-04 2004-03-03 位置情報推定装置、その方法、及びプログラム

Country Status (5)

Country Link
US (1) US7039546B2 (ja)
EP (1) EP1600789B1 (ja)
JP (1) JP3881367B2 (ja)
DE (1) DE602004029867D1 (ja)
WO (1) WO2004079388A1 (ja)

Cited By (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
EP1792366A2 (en) * 2004-09-23 2007-06-06 Interdigital Technology Corporation Blind signal separation using array deflection
JP2007226036A (ja) * 2006-02-24 2007-09-06 Nippon Telegr & Teleph Corp <Ntt> 信号分離装置、信号分離方法、信号分離プログラム及び記録媒体、並びに、信号到来方向推定装置、信号到来方向推定方法、信号到来方向推定プログラム及び記録媒体
JP2007240209A (ja) * 2006-03-06 2007-09-20 Kddi Corp 信号到来方向推定装置及び方法、並びに信号分離装置及び方法、コンピュータプログラム
JP2008039693A (ja) * 2006-08-09 2008-02-21 Toshiba Corp 方向探知システム及び信号抽出方法
JP2009074990A (ja) * 2007-09-21 2009-04-09 Toshiba Corp 方向測定装置
US7647209B2 (en) 2005-02-08 2010-01-12 Nippon Telegraph And Telephone Corporation Signal separating apparatus, signal separating method, signal separating program and recording medium
US8577054B2 (en) 2009-03-30 2013-11-05 Sony Corporation Signal processing apparatus, signal processing method, and program

Families Citing this family (26)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US7265715B2 (en) * 2004-09-23 2007-09-04 Interdigital Technology Corporation Communications device with adaptive decoding and associated methods
JP4278610B2 (ja) * 2004-12-28 2009-06-17 富士通株式会社 数値解析支援装置，数値解析支援方法，数値解析支援プログラムおよび同プログラムを記録したコンピュータ読取可能な記録媒体
US7812718B1 (en) 2005-11-21 2010-10-12 The Hong Kong University Of Science And Technology Distributed position estimation for wireless sensor networks
US7769118B2 (en) * 2006-02-10 2010-08-03 Interdigital Technolgy Corporation Method and apparatus for equalizing received signals based on independent component analysis
JP2009529699A (ja) * 2006-03-01 2009-08-20 ソフトマックス，インコーポレイテッド 分離信号を生成するシステムおよび方法
JP4952979B2 (ja) * 2006-04-27 2012-06-13 独立行政法人理化学研究所 信号分離装置、信号分離方法、ならびに、プログラム
KR100780210B1 (ko) * 2006-07-20 2007-11-27 삼성전기주식회사 휴대형 보안 송신 장치 및 보안 인증 시스템
KR101434200B1 (ko) * 2007-10-01 2014-08-26 삼성전자주식회사 혼합 사운드로부터의 음원 판별 방법 및 장치
GB0720473D0 (en) * 2007-10-19 2007-11-28 Univ Surrey Accoustic source separation
KR101394338B1 (ko) * 2007-10-31 2014-05-30 삼성전자주식회사 무선 센서 네트워크의 토폴로지 정보 표시 방법 및 장치 및이를 위한 시스템
US8321214B2 (en) * 2008-06-02 2012-11-27 Qualcomm Incorporated Systems, methods, and apparatus for multichannel signal amplitude balancing
JP5195652B2 (ja) * 2008-06-11 2013-05-08 ソニー株式会社 信号処理装置、および信号処理方法、並びにプログラム
EP2159593B1 (en) * 2008-08-26 2012-05-02 Nuance Communications, Inc. Method and device for locating a sound source
US20100265800A1 (en) * 2009-04-16 2010-10-21 Graham Paul Eatwell Array shape estimation using directional sensors
JP5299233B2 (ja) * 2009-11-20 2013-09-25 ソニー株式会社 信号処理装置、および信号処理方法、並びにプログラム
US8521477B2 (en) * 2009-12-18 2013-08-27 Electronics And Telecommunications Research Institute Method for separating blind signal and apparatus for performing the same
KR101419377B1 (ko) * 2009-12-18 2014-07-15 배재대학교 산학협력단 암묵신호 분리 방법 및 이를 수행하는 장치
JP2012018066A (ja) * 2010-07-07 2012-01-26 Panasonic Electric Works Sunx Co Ltd 異常検査装置
WO2012105385A1 (ja) * 2011-02-01 2012-08-09 日本電気株式会社 有音区間分類装置、有音区間分類方法、及び有音区間分類プログラム
KR101241842B1 (ko) * 2011-10-12 2013-03-11 광운대학교 산학협력단 3차원 측위 장치 및 방법
US20140303929A1 (en) * 2013-04-03 2014-10-09 Umm Al-Qura University Method to obtain accurate vertical component estimates in 3d positioning
US10657958B2 (en) * 2015-03-18 2020-05-19 Sogang University Research Foundation Online target-speech extraction method for robust automatic speech recognition
US10991362B2 (en) * 2015-03-18 2021-04-27 Industry-University Cooperation Foundation Sogang University Online target-speech extraction method based on auxiliary function for robust automatic speech recognition
US11694707B2 (en) 2015-03-18 2023-07-04 Industry-University Cooperation Foundation Sogang University Online target-speech extraction method based on auxiliary function for robust automatic speech recognition
US10334390B2 (en) * 2015-05-06 2019-06-25 Idan BAKISH Method and system for acoustic source enhancement using acoustic sensor array
DE102020213286A1 (de) * 2020-10-21 2022-04-21 Robert Bosch Gesellschaft mit beschränkter Haftung Verfahren zur Bestimmung einer Phasenlage eines Drehratensignals oder eines Quadratursignals, Verfahren zur Anpassung einer Demodulationsphase und Drehratensensor

Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO1998058450A1 (en) * 1997-06-18 1998-12-23 Clarity, L.L.C. Methods and apparatus for blind signal separation
JP2000242624A (ja) * 1999-02-18 2000-09-08 Retsu Yamakawa 信号分離装置
WO2000054404A1 (en) * 1999-03-08 2000-09-14 Telefonaktiebolaget Lm Ericsson (Publ) Method and device for separating a mixture of source signals
JP2002084593A (ja) * 2000-09-08 2002-03-22 Inst Of Physical & Chemical Res 信号分離システムおよび信号分離方法
JP2003090871A (ja) * 2001-09-19 2003-03-28 Toshiba Corp 受信装置
JP2003099093A (ja) * 2001-09-25 2003-04-04 Nippon Telegr & Teleph Corp <Ntt> 雑音抑圧方法及びその装置、雑音抑圧プログラム並びにそのプログラム記録媒体

Family Cites Families (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US5737431A (en) * 1995-03-07 1998-04-07 Brown University Research Foundation Methods and apparatus for source location estimation from microphone-array time-delay estimates
JP3975153B2 (ja) 2002-10-28 2007-09-12 日本電信電話株式会社 ブラインド信号分離方法及び装置、ブラインド信号分離プログラム並びにそのプログラムを記録した記録媒体

Patent Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO1998058450A1 (en) * 1997-06-18 1998-12-23 Clarity, L.L.C. Methods and apparatus for blind signal separation
JP2000242624A (ja) * 1999-02-18 2000-09-08 Retsu Yamakawa 信号分離装置
WO2000054404A1 (en) * 1999-03-08 2000-09-14 Telefonaktiebolaget Lm Ericsson (Publ) Method and device for separating a mixture of source signals
JP2002084593A (ja) * 2000-09-08 2002-03-22 Inst Of Physical & Chemical Res 信号分離システムおよび信号分離方法
JP2003090871A (ja) * 2001-09-19 2003-03-28 Toshiba Corp 受信装置
JP2003099093A (ja) * 2001-09-25 2003-04-04 Nippon Telegr & Teleph Corp <Ntt> 雑音抑圧方法及びその装置、雑音抑圧プログラム並びにそのプログラム記録媒体

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
See also references of EP1600789A4 *

Cited By (10)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
EP1792366A2 (en) * 2004-09-23 2007-06-06 Interdigital Technology Corporation Blind signal separation using array deflection
EP1792366A4 (en) * 2004-09-23 2007-09-12 Interdigital Tech Corp BLIND SEPARATION OF SIGNALS WITH THE HELP OF ARRAY SLOWING
KR100848846B1 (ko) * 2004-09-23 2008-07-28 인터디지탈 테크날러지 코포레이션 어레이 편향을 사용하는 암묵 신호 분리
US7647209B2 (en) 2005-02-08 2010-01-12 Nippon Telegraph And Telephone Corporation Signal separating apparatus, signal separating method, signal separating program and recording medium
JP2007226036A (ja) * 2006-02-24 2007-09-06 Nippon Telegr & Teleph Corp <Ntt> 信号分離装置、信号分離方法、信号分離プログラム及び記録媒体、並びに、信号到来方向推定装置、信号到来方向推定方法、信号到来方向推定プログラム及び記録媒体
JP4630203B2 (ja) * 2006-02-24 2011-02-09 日本電信電話株式会社 信号分離装置、信号分離方法、信号分離プログラム及び記録媒体、並びに、信号到来方向推定装置、信号到来方向推定方法、信号到来方向推定プログラム及び記録媒体
JP2007240209A (ja) * 2006-03-06 2007-09-20 Kddi Corp 信号到来方向推定装置及び方法、並びに信号分離装置及び方法、コンピュータプログラム
JP2008039693A (ja) * 2006-08-09 2008-02-21 Toshiba Corp 方向探知システム及び信号抽出方法
JP2009074990A (ja) * 2007-09-21 2009-04-09 Toshiba Corp 方向測定装置
US8577054B2 (en) 2009-03-30 2013-11-05 Sony Corporation Signal processing apparatus, signal processing method, and program

Also Published As

Publication number Publication date
EP1600789A1 (en) 2005-11-30
US20050203981A1 (en) 2005-09-15
JPWO2004079388A1 (ja) 2006-06-08
US7039546B2 (en) 2006-05-02
JP3881367B2 (ja) 2007-02-14
EP1600789B1 (en) 2010-11-03
EP1600789A4 (en) 2006-09-20
DE602004029867D1 (de) 2010-12-16

Similar Documents

Publication Publication Date Title
WO2004079388A1 (ja) 位置情報推定装置、その方法、及びプログラム
JP6335985B2 (ja) マルチセンサ音源定位
Chen et al. Maximum-likelihood source localization and unknown sensor location estimation for wideband signals in the near-field
Bialer et al. Performance advantages of deep neural networks for angle of arrival estimation
Herzog et al. Eigenbeam-ESPRIT for DOA-vector estimation
US20070253561A1 (en) Systems and methods for audio enhancement
CN113126028B (zh) 一种基于多个麦克风阵列的噪声源定位方法
Talmon et al. Supervised source localization using diffusion kernels
Jo et al. Direction of arrival estimation using nonsingular spherical ESPRIT
JP5654980B2 (ja) 音源位置推定装置、音源位置推定方法、及び音源位置推定プログラム
JP2018063200A (ja) 音源位置推定装置、音源位置推定方法、及びプログラム
Traa et al. Robust source localization and enhancement with a probabilistic steered response power model
Chen et al. Broadband sound source localisation via non-synchronous measurements for service robots: A tensor completion approach
Lin et al. Beamforming pointing error of a triaxial velocity sensor under gain uncertainties
Krause et al. Data diversity for improving DNN-based localization of concurrent sound events
JP4738284B2 (ja) ブラインド信号抽出装置、その方法、そのプログラム、及びそのプログラムを記録した記録媒体
Yao et al. A fast multi-source sound DOA estimator considering colored noise in circular array
Hu et al. Analytical geometry calibration for acoustic transceiver arrays
Jarrett et al. Eigenbeam-based acoustic source tracking in noisy reverberant environments
Matsumoto et al. Multiple signal classification by aggregated microphones
Kijima et al. Tracking of multiple moving sound sources using particle filter for arbitrary microphone array configurations
Chen et al. A multiplanes geometric approach for sound source localization with TDOA
Hioka et al. Multiple-speech-source localization using advanced histogram mapping method
Markovich-Golan et al. A probability distribution model for the relative transfer function in a reverberant environment
Wei et al. Estimating Angle of Arrival (AoA) of multiple Echoes in a Steering Vector Space

Legal Events

Date Code Title Description
WWE Wipo information: entry into national phase

Ref document number: 2005503053

Country of ref document: JP

AK Designated states

Kind code of ref document: A1

Designated state(s): AE AG AL AM AT AU AZ BA BB BG BR BW BY BZ CA CH CN CO CR CU CZ DE DK DM DZ EC EE EG ES FI GB GD GE GH GM HR HU ID IL IN IS JP KE KG KP KR KZ LC LK LR LS LT LU LV MA MD MG MK MN MW MX MZ NA NI NO NZ OM PG PH PL PT RO RU SC SD SE SG SK SL SY TJ TM TN TR TT TZ UA UG US UZ VC VN YU ZA ZM ZW

AL Designated countries for regional patents

Kind code of ref document: A1

Designated state(s): BW GH GM KE LS MW MZ SD SL SZ TZ UG ZM ZW AM AZ BY KG KZ MD RU TJ TM AT BE BG CH CY CZ DE DK EE ES FI FR GB GR HU IE IT LU MC NL PL PT RO SE SI SK TR BF BJ CF CG CI CM GA GN GQ GW ML MR NE SN TD TG

WWE Wipo information: entry into national phase

Ref document number: 2004716703

Country of ref document: EP

Ref document number: 10508708

Country of ref document: US

Ref document number: 20048001009

Country of ref document: CN

121 Ep: the epo has been informed by wipo that ep was designated in this application
WWP Wipo information: published in national office

Ref document number: 2004716703

Country of ref document: EP