WO2020003342A1 - 波源方向推定装置、波源方向推定方法、およびプログラム記録媒体 - Google Patents

波源方向推定装置、波源方向推定方法、およびプログラム記録媒体 Download PDF

Info

Publication number
WO2020003342A1
WO2020003342A1 PCT/JP2018/023970 JP2018023970W WO2020003342A1 WO 2020003342 A1 WO2020003342 A1 WO 2020003342A1 JP 2018023970 W JP2018023970 W JP 2018023970W WO 2020003342 A1 WO2020003342 A1 WO 2020003342A1
Authority
WO
WIPO (PCT)
Prior art keywords
frequency
wave
cross spectrum
direction information
wave source
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Ceased
Application number
PCT/JP2018/023970
Other languages
English (en)
French (fr)
Inventor
友督 荒井
裕三 仙田
玲史 近藤
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
NEC Corp
Original Assignee
NEC Corp
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 NEC Corp filed Critical NEC Corp
Priority to PCT/JP2018/023970 priority Critical patent/WO2020003342A1/ja
Priority to JP2020526724A priority patent/JP7056739B2/ja
Priority to US17/252,391 priority patent/US20210263125A1/en
Publication of WO2020003342A1 publication Critical patent/WO2020003342A1/ja
Anticipated expiration legal-status Critical
Ceased legal-status Critical Current

Links

Images

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/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
    • H—ELECTRICITY
    • H04—ELECTRIC COMMUNICATION TECHNIQUE
    • H04R—LOUDSPEAKERS, MICROPHONES, GRAMOPHONE PICK-UPS OR LIKE ACOUSTIC ELECTROMECHANICAL TRANSDUCERS; ELECTRIC HEARING AIDS; PUBLIC ADDRESS SYSTEMS
    • H04R3/00—Circuits for transducers
    • H04R3/005—Circuits for transducers for combining the signals of two or more microphones
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01H—MEASUREMENT OF MECHANICAL VIBRATIONS OR ULTRASONIC, SONIC OR INFRASONIC WAVES
    • G01H3/00—Measuring characteristics of vibrations by using a detector in a fluid
    • 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/801—Details
    • 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
    • G01S3/8083—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 determining direction of source
    • G—PHYSICS
    • G10—MUSICAL INSTRUMENTS; ACOUSTICS
    • G10L—SPEECH ANALYSIS TECHNIQUES OR SPEECH SYNTHESIS; SPEECH RECOGNITION; SPEECH OR VOICE PROCESSING TECHNIQUES; SPEECH OR AUDIO CODING OR DECODING
    • G10L25/00—Speech or voice analysis techniques not restricted to a single one of groups G10L15/00 - G10L21/00
    • G10L25/03—Speech or voice analysis techniques not restricted to a single one of groups G10L15/00 - G10L21/00 characterised by the type of extracted parameters
    • G10L25/06—Speech or voice analysis techniques not restricted to a single one of groups G10L15/00 - G10L21/00 characterised by the type of extracted parameters the extracted parameters being correlation coefficients
    • G—PHYSICS
    • G10—MUSICAL INSTRUMENTS; ACOUSTICS
    • G10L—SPEECH ANALYSIS TECHNIQUES OR SPEECH SYNTHESIS; SPEECH RECOGNITION; SPEECH OR VOICE PROCESSING TECHNIQUES; SPEECH OR AUDIO CODING OR DECODING
    • G10L25/00—Speech or voice analysis techniques not restricted to a single one of groups G10L15/00 - G10L21/00
    • G10L25/03—Speech or voice analysis techniques not restricted to a single one of groups G10L15/00 - G10L21/00 characterised by the type of extracted parameters
    • G10L25/18—Speech or voice analysis techniques not restricted to a single one of groups G10L15/00 - G10L21/00 characterised by the type of extracted parameters the extracted parameters being spectral information of each sub-band
    • H—ELECTRICITY
    • H04—ELECTRIC COMMUNICATION TECHNIQUE
    • H04R—LOUDSPEAKERS, MICROPHONES, GRAMOPHONE PICK-UPS OR LIKE ACOUSTIC ELECTROMECHANICAL TRANSDUCERS; ELECTRIC HEARING AIDS; PUBLIC ADDRESS SYSTEMS
    • H04R1/00—Details of transducers, loudspeakers or microphones
    • H04R1/20—Arrangements for obtaining desired frequency or directional characteristics
    • H04R1/32—Arrangements for obtaining desired frequency or directional characteristics for obtaining desired directional characteristic only
    • H04R1/40—Arrangements for obtaining desired frequency or directional characteristics for obtaining desired directional characteristic only by combining a number of identical transducers
    • H04R1/406—Arrangements for obtaining desired frequency or directional characteristics for obtaining desired directional characteristic only by combining a number of identical transducers microphones
    • H—ELECTRICITY
    • H04—ELECTRIC COMMUNICATION TECHNIQUE
    • H04R—LOUDSPEAKERS, MICROPHONES, GRAMOPHONE PICK-UPS OR LIKE ACOUSTIC ELECTROMECHANICAL TRANSDUCERS; ELECTRIC HEARING AIDS; PUBLIC ADDRESS SYSTEMS
    • H04R2430/00—Signal processing covered by H04R, not provided for in its groups
    • H04R2430/20—Processing of the output signals of the acoustic transducers of an array for obtaining a desired directivity characteristic
    • H04R2430/21—Direction finding using differential microphone array [DMA]

Definitions

  • the present invention relates to a wave direction estimation device, a wave direction estimation method, and a program.
  • the present invention relates to a source direction estimating device, a source direction estimating method, and a program that estimate a source direction based on signals acquired by a plurality of sensors.
  • Patent Literature 1 and Non-Patent Literature 1 disclose a method of estimating the direction of a sound source from the arrival time difference between sound reception signals of two microphones.
  • the probability density functions of the arrival time difference of sound waves are obtained for each frequency, the arrival time difference is calculated from the probability density function obtained by superimposing them, and the sound source direction is estimated.
  • Patent Document 2 discloses a search method in which sound and vibration transmitted to a predetermined observation point are sampled to search whether or not the sound source of the sound is a sound generated from a vibration source.
  • the sound transmitted from the sound source and the vibration of the surface wave transmitted from the vibration source are simultaneously measured.
  • the direction of the sound source obtained from the data of the sound pressure level of the sound is compared with the direction of the vibration source obtained from the data of the vibration level of the vibration, and the sound from the sound source is transmitted from the vibration source accompanied by sound generation. Determine if it is a sound.
  • Patent Literature 1 and Non-Patent Literature 1 in a frequency band having a high signal-to-noise ratio (SNR), the probability density function of the arrival time difference forms a sharp peak. Even if the band is small, the arrival time difference can be accurately estimated.
  • SNR signal-to-noise ratio
  • Patent Literature 1 and Non-Patent Literature 1 when the probability density functions of the arrival time differences for each frequency are superimposed, even if no sound source exists, the probability density functions superimposed by accidentally aligning the phases are superposed. A peak is generated. Therefore, the methods of Patent Literature 1 and Non-Patent Literature 1 have a problem that a virtual image sound source is erroneously estimated.
  • the method of Patent Document 2 it is determined whether the sound from a sound source is a sound from a vibration source accompanied by sound generation or a sound from a sound source without vibration, or the vibration source does not accompany a sound. It is possible to accurately determine whether or not it is a vibration source.
  • the method of Patent Literature 2 has a problem that the arrival time difference of a virtual image sound source in a direction different from that of a sound source may be calculated due to the accidental alignment of phases between different microphones, and the virtual image sound source may be erroneously estimated. Was.
  • An object of the present invention is to provide a wave source direction estimating apparatus that solves the above-described problems, reduces occurrence of erroneous estimation of a virtual image sound source, and can estimate the direction of a sound source with high accuracy.
  • the wave source direction estimating apparatus combines a plurality of input units that acquire an electric signal converted from a wave acquired by a plurality of sensors as an input signal, and at least two input signals from the plurality of input signals.
  • a signal selection unit that selects at least two pairs, and a relative delay time calculation unit that calculates, as a relative delay time, an arrival time difference of the wave for each wave source search direction between at least two input signals that form the input signal pair,
  • At least one frequency-dependent estimated direction information generating unit that generates the estimated direction information of the wave source of the wave for each frequency using the input signal pair and the relative delay time, and for each frequency generated by the frequency-based estimated direction information generating unit
  • an integrating unit for integrating the estimated direction information.
  • the information processing device obtains, as an input signal, an electric signal converted from a wave obtained by the plurality of sensors, and obtains at least two input signals from the plurality of input signals. Select at least two of the combined pairs, calculate the arrival time difference of the wave for each wave source search direction as a relative delay time between at least two input signals constituting the input signal pair, and input signal pair and relative delay time. Is used to generate at least one estimated direction information of a wave source for each frequency, and integrates the estimated direction information for each frequency.
  • a program includes a process of acquiring an electric signal converted from a wave acquired by a plurality of sensors as an input signal, and at least two pairs of a combination of at least two input signals from the plurality of input signals.
  • the process of selecting, the process of calculating the arrival time difference of the wave for each wave source search direction as a relative delay time between at least two input signals constituting a pair of input signals, and using the input signal pair and the relative delay time A process for generating at least one estimated direction information of the wave source of the wave for each frequency, a process for integrating the estimated direction information for each frequency, and a computer are executed.
  • the present invention it is possible to provide a wave source direction estimating apparatus capable of reducing the occurrence of erroneous estimation of a virtual image sound source and estimating the direction of the sound source with high accuracy.
  • FIG. 6 is a conceptual diagram for describing an example of a process of a relative delay time calculation unit in the wave source direction estimation device according to the first embodiment of the present invention. It is a conceptual diagram for demonstrating another example of the process of the relative delay time calculation part in the wave source direction estimation apparatus which concerns on 1st Embodiment of this invention.
  • FIG. 2 is a block diagram illustrating an example of a configuration of a frequency-specific estimated direction information generation unit included in the wave source direction estimation device according to the first embodiment of the present invention.
  • FIG. 3 is a block diagram illustrating an example of a configuration of a cross-spectrum generating unit for each frequency included in the wave source direction estimating apparatus according to the first embodiment of the present invention. It is a block diagram showing an example of composition which added at least one sensor to a wave source direction estimating device concerning a 1st embodiment of the present invention.
  • 5 is a flowchart for explaining an outline of an operation of the wave source direction estimation device according to the first embodiment of the present invention.
  • 6 is a flowchart for explaining an operation of a frequency-specific estimated direction information generation unit of the wave source direction estimation device according to the first embodiment of the present invention.
  • FIG. 5 is a flowchart for explaining an operation of a frequency-specific cross spectrum generating unit of the frequency-specific estimated direction information generating unit of the wave source direction estimating device according to the first embodiment of the present invention. It is a block diagram showing an example of composition of a wave source direction estimating device concerning a 2nd embodiment of the present invention.
  • FIG. 3 is a block diagram illustrating an example of a hardware configuration for realizing a wave source direction estimation device according to each embodiment of the present invention.
  • the wave source direction estimating apparatus of the present embodiment estimates the source of a sound wave that is a vibration wave of air or water. Therefore, the wave source direction estimating device of the present embodiment verifies the vibration wave converted into the electric signal by the microphone.
  • the estimation target of the wave source direction estimating apparatus of the present embodiment is not limited to a sound wave generation source, and can be used to estimate an arbitrary wave generation source (also referred to as a wave source) such as a vibration wave or an electromagnetic wave.
  • FIG. 1 is a block diagram illustrating a configuration of a wave source direction estimating apparatus 10 according to the present embodiment.
  • the wave source direction estimating device 10 includes an input terminal 11, a signal selecting unit 12, a relative delay time calculating unit 13, a frequency-specific estimated direction information generating unit 15, and an integrating unit 17.
  • the wave source direction estimating apparatus 10 includes p input terminals 11 (p is an integer of 2 or more).
  • the wave source direction estimating apparatus 10 includes R frequency-specific estimated direction information generating units 15 (R is an integer of 1 or more).
  • R is an integer of 1 or more.
  • numbers from 1 to p are attached with a hyphen between them.
  • numbers 1 to R are attached after a symbol with a hyphen between them.
  • Each of the input terminals 11-1 to 11-p (also referred to as an input unit) is connected to a microphone not shown (hereinafter also referred to as a microphone). Electrical signals converted from sound waves (also referred to as sound signals) collected by microphones arranged at different positions are input to each of the input terminals 11-1 to 11-p as input signals.
  • the input signal input to the m-th input terminal 11-m at the time t is referred to as x m (t) (t: real number, m: an integer of 1 to p).
  • a microphone collects sound waves in which sound generated from a target sound source and various noises generated around the microphone are mixed, and converts the collected sound waves into a digital signal (also referred to as a sample value sequence).
  • the microphones are arranged at different positions corresponding to the input terminals 11-1 to 11-p in order to collect sounds from the target sound source.
  • the input signal converted from the sound wave collected by the m-th microphone is supplied to the m-th input terminal 11-m.
  • the input signal supplied to the m-th input terminal 11-m is also referred to as an “m-th microphone input signal”.
  • the signal selector 12 selects two input signals from among the P input signals supplied to the input terminals 11-1 to 11-p.
  • the signal selection unit 12 outputs the selected two input signals to the frequency-dependent estimated direction information generation units 15-1 to 15-R, and outputs the position information of the microphone that supplies those input signals (hereinafter, also referred to as microphone position information).
  • the number R of the frequency-specific estimated direction information generators 15 corresponds to the number R of combinations of input signals.
  • the signal selection unit 12 may select all combinations or some combinations. When all combinations are selected, R is represented by the following equation 1.
  • the wave source direction estimating apparatus 10 estimates the direction of the sound source using the time difference between the sound from the target sound source and the two microphones. If the distance between the microphones (hereinafter, also referred to as the microphone distance) is too large, the sound from the target sound source is not observed as the same sound due to the influence of a medium such as air or water, so that the direction estimation accuracy decreases. On the other hand, if the microphone interval is too small, the difference in arrival time of the sound wave between the two microphones becomes too small, and the direction estimation accuracy decreases. Therefore, as shown in Expression 2, it is preferable that the signal selection unit 12 selects an input signal of a microphone in which the microphone interval d falls within a certain range (d min , d max : real numbers).
  • the signal selection unit 12 may select two input signals having the maximum microphone interval d.
  • the signal selecting unit 12 arranges the microphone intervals d in descending order and selects a combination of input signals up to the R-order (r ⁇ C (p, 2)). Is also good.
  • the selection of some combinations by the signal selection unit 12 leads to a reduction in the amount of calculation in addition to preventing a decrease in the direction estimation accuracy.
  • Microphone position information is also important when determining the difference in arrival time of sound from a target sound source to two microphones. Therefore, the signal selection unit 12 outputs the microphone position information to the relative delay time calculation unit 13 in addition to the input signal.
  • the relative delay time calculation unit 13 calculates the relative delay time between all the microphone pairs selected by the signal selection unit 12 using the microphone position information and the sound source search target direction.
  • the relative delay time is a difference in arrival time of a sound wave uniquely determined based on a microphone interval and a sound source direction.
  • the sound source search target direction is set at predetermined angle intervals. That is, the relative delay time is calculated for the direction of the sound source search target.
  • the relative delay time calculation unit 13 sets the calculated sound source search target direction and the calculated relative delay time and outputs the set to the frequency-specific estimated direction information generation unit 15.
  • the method of calculating the relative delay time differs depending on the positional relationship between the microphone pairs.
  • two positional relationships of the microphone pairs will be described, and a calculation method of the relative delay time will be described for each of the positional relationships of the microphone pairs.
  • FIG. 2 shows an example in which all microphones are arranged on the same straight line.
  • the sound speed is c
  • the microphone interval is d r
  • the sound source search target direction also referred to as a sound source direction
  • the sound source direction ⁇ is at least one angle set for estimating the direction of the sound source 100.
  • the relative delay time ⁇ r ( ⁇ ) with respect to the sound source direction ⁇ can be calculated using the following Expression 3.
  • the microphone interval d differs depending on the combination of input signals selected by the signal selection unit 12. Therefore, the relative delay time ⁇ r ( ⁇ ) differs for each combination number r. For example, assuming that the distance between the microphone pairs AB in FIG. 2 is d 1 , the relative delay time ⁇ 1 ( ⁇ ) can be calculated using the following Expression 4.
  • the relative delay time ⁇ 2 ( ⁇ ) can be calculated using the following Expression 5.
  • the relative delay time ⁇ r ( ⁇ ) for a certain sound source is proportional to the microphone interval d, but the sound source direction ⁇ can be regarded as the same from any microphone.
  • FIG. 3 shows an example in which two microphone pairs are arranged on a straight line orthogonal to each other.
  • the sound source direction ⁇ differs depending on the microphone pair.
  • the relative delay time ⁇ 1 ( ⁇ ) between the microphone pairs AB in FIG. 3 can be calculated using the following Expression 6.
  • the relative delay time ⁇ 2 ( ⁇ ) between the microphones CD in FIG. 3 can be calculated by the following Expression 7.
  • the relative delay time ⁇ r ( ⁇ ) of another microphone pair can be generalized as a function of the sound source direction ⁇ as viewed from the reference microphone pair, as in the following Expression 8. Note that any reference microphone pair may be selected.
  • the relative delay time calculator 13 calculates relative delay times for all sound source search target directions. For example, if the search range in the sound source direction is 0 to 90 degrees in increments of 10 degrees, that is, 0, 10, 20,. Calculate the relative delay time. Then, the relative delay time calculator 13 outputs the sound source search target direction and the relative delay time to the frequency-specific estimated direction information generator 15.
  • the frequency-dependent estimated direction information generation units 15-1 to 15-R receive the input signal of one of the microphone pairs selected by the signal selection unit 12 and the relative delay time supplied from the relative delay time calculation unit 13. Is entered.
  • the frequency-specific estimated direction information generation units 15-1 to 15-R generate frequency-specific estimated direction information between input signals of the microphone pair using the input signal of the microphone pair and the relative delay time.
  • FIG. 4 is a block diagram of the frequency-specific estimated direction information generator 15.
  • the frequency-specific estimated direction information generation unit 15 includes a conversion unit 151, a cross spectrum calculation unit 152, an average calculation unit 153, a variance calculation unit 154, a frequency-specific cross spectrum generation unit 155, an inverse conversion unit 156, and a frequency-specific estimated direction information calculation.
  • a unit 157 is provided.
  • Conversion section 151 receives two input signals (input signal A and input signal B) from signal selection section 12.
  • the converter 151 converts the two input signals supplied from the signal selector 12 into a converted signal (also called a frequency domain signal).
  • the conversion unit 151 performs conversion for decomposing an input signal into a plurality of frequency components. For example, the conversion unit 151 decomposes the input signal into a plurality of frequency components using Fourier transform.
  • Conversion section 151 outputs a conversion signal to cross spectrum calculation section 152.
  • the conversion unit 151 receives two types of input signals x m (t).
  • m is the number of the input terminal 11.
  • the converter 151 cuts out a waveform of an appropriate length from the input signal supplied from the input terminal 11 while shifting the waveform at a constant cycle.
  • the signal section thus cut out is called a frame
  • the length of the cut out waveform is called a frame length
  • the cycle of shifting the frame is called a frame cycle.
  • conversion section 151 converts the cut-out signal into a frequency domain signal using Fourier transform.
  • n is a frame number
  • the Fourier transform X m (k, n) of the input signal x m (t, n) can be calculated using the following equation 9.
  • K represents a frequency bin number and is an integer of 0 or more and K-1 or less.
  • k is also referred to simply as a frequency, not as a frequency bin number.
  • the converted signal is input from the converter 151 to the cross spectrum calculator 152.
  • the cross spectrum calculator 152 calculates a cross spectrum using the converted signal supplied from the converter 151.
  • the cross spectrum calculation section 152 outputs the calculated cross spectrum to the average calculation section 153.
  • the cross spectrum calculator 152 calculates the product of the complex conjugate of the converted signal X 2 (k, n) and the converted signal X 1 (k, n) to calculate the cross spectrum.
  • the cross spectrum of the converted signal is S 12 (k, n).
  • the cross spectrum calculation unit 152 calculates the cross spectrum using the following Expression 10.
  • conj (X 2 (k, n)) represents the complex conjugate of X 2 (k, n).
  • a cross spectrum normalized by an amplitude component may be used.
  • the cross spectrum calculation unit 152 calculates the cross spectrum using the following Expression 11.
  • the average calculator 153 receives the cross spectrum from the cross spectrum calculator 152.
  • the average calculation unit 153 calculates an average of the cross spectra supplied from the cross spectrum calculation unit 152 (also referred to as an average cross spectrum).
  • the average calculation section 153 outputs the calculated average cross spectrum to the variance calculation section 154 and the cross spectrum generation section 155 for each frequency.
  • the average calculation unit 153 calculates the average cross spectrum for each frequency bin from the cross spectrum input in the past.
  • the average calculation unit 153 may calculate the average cross spectrum in units of subbands in which a plurality of frequency bins are bundled, instead of in units of frequency bins.
  • the cross spectrum in the frequency bin k of the n-th frame is defined as S 12 (k, n).
  • the average calculation unit 153 uses the following equation 12 to calculate the average cross spectrum obtained from the past L frames. Calculate SS 12 (k, n).
  • the average calculation unit 153 may calculate the average cross spectrum SS 12 (k, n) using the following leak integration.
  • ⁇ is a real number larger than 0 and smaller than 1.
  • the variance calculator 154 receives the average cross spectrum from the average calculator 153.
  • the variance calculation unit 154 calculates variance using the average cross spectrum supplied from the average calculation unit 153.
  • the variance calculation unit 154 outputs the calculated variance to the frequency-specific cross spectrum generation unit 155.
  • the variance calculation unit 154 calculates the variance V 12 (k, n) using the following Expression 14.
  • the variance calculation unit 154 may calculate the variance V 12 (k, n) using the following Expression 15.
  • the variance calculation unit 154 calculates the variance V 12 (k, n) using Expression 16 below.
  • FIG. 5 is a block diagram illustrating an example of the configuration of the frequency-specific cross spectrum generation unit 155.
  • the frequency-specific cross spectrum generation unit 155 includes a frequency-specific basic cross spectrum calculation unit 551, a kernel function spectrum generation unit 552, and a multiplication unit 553.
  • the average cross spectrum from the average calculator 153 is input to the basic cross spectrum calculator 551 for each frequency.
  • the frequency-specific basic cross spectrum calculator 551 calculates a frequency-specific basic cross spectrum using the average cross spectrum supplied from the average calculator 153.
  • the frequency-specific basic cross spectrum calculation unit 551 outputs the calculated frequency-specific basic cross spectrum to the multiplication unit 553.
  • the frequency-specific basic cross spectrum calculation unit 551 uses the average cross spectrum SS 12 (k, n) supplied from the average calculation unit 153 to calculate a cross corresponding to each frequency k of the average cross spectrum SS 12 (k, n). Calculate the spectrum (also called the basic cross spectrum for each frequency).
  • the frequency-specific basic cross spectrum calculation unit 551 outputs the calculated frequency-specific basic cross spectrum to the multiplication unit 553.
  • the basic cross spectrum for each frequency is calculated to calculate a correlation function for each frequency component.
  • the frequency-specific basic cross spectrum calculation unit 551 calculates a frequency-specific basic cross spectrum for calculating a correlation function (also referred to as a frequency-specific correlation function) corresponding to a certain frequency k in a subsequent stage.
  • the frequency-specific basic cross spectrum calculation unit 551 calculates the frequency-specific basic cross spectrum of the frequency k.
  • the basic cross spectrum calculating unit 551 obtains the phase component and the amplitude component separately beforehand. Integrate. Assuming that the basic cross spectrum U k (w, n) for each frequency of the frequency k, its amplitude component is
  • w represents a frequency, and is an integer of 0 or more and W-1 or less.
  • the frequency-based basic cross spectrum calculation unit 551 uses the average cross spectrum SS 12 (k, n) of the frequency k to calculate the amplitude component
  • phase component arg (U k (w, n)) of a frequency obtained by multiplying k by an integer a value obtained by multiplying the average cross spectrum SS 12 (k, n) of the frequency k by a constant is used.
  • the frequency k, 2k, 3k, and each of phase components of 4k is, arg (SS 12 (k, n)), 2arg (SS 12 (k, n)), 3arg (SS 12 (k, n) ), 4 arg (SS 12 (k, n)).
  • the phase component of a frequency that is a non-integer multiple of k is set to zero. Therefore, the phase component arg (U k (w, n)) of the basic cross spectrum for each frequency corresponding to the frequency k is calculated using the following equation (20).
  • p is an integer of 1 or more and P or less (P> 1).
  • the frequency-specific basic cross spectrum calculation unit 551 integrates the amplitude component calculated using Expressions 18 and 19 and the phase component calculated using Expression 20 using Expression 17, and calculates the frequency of the frequency k. Another basic cross spectrum U k (w, n) is obtained.
  • the amplitude component and the phase component are separately obtained, and then the fundamental cross spectrum for each frequency is calculated.
  • the power of the cross spectrum shown in the following Expression 21 is used, the basic cross spectrum U k (w, n) for each frequency can be obtained without obtaining the amplitude component and the phase component.
  • kernel function spectrum generator The variance is input from the variance calculator 154 to the kernel function spectrum generator 552.
  • the kernel function spectrum generator 552 calculates the kernel function spectrum using the variance supplied from the variance calculator 154.
  • the kernel function spectrum is obtained by Fourier-transforming a kernel function and taking its absolute value.
  • the kernel function spectrum may be squared instead of taking the absolute value of the Fourier transform of the kernel function.
  • the kernel function spectrum may be obtained by performing a Fourier transform on the kernel function and obtaining the square of its absolute value.
  • the kernel function spectrum generator 552 outputs the calculated kernel function spectrum to the multiplier 553.
  • the kernel function spectrum is G (w), and the kernel function is g ( ⁇ ).
  • a Gaussian function is used as a kernel function. At this time, the Gaussian function is given by the following Expression 22.
  • g 1, g 2, and g 3 are positive real numbers.
  • g 1 controls the magnitude of the Gaussian function
  • g 2 controls the position of the peak of the Gaussian function
  • g 3 controls the spread of the Gaussian function.
  • g 3 for adjusting the spread of the Gaussian function is important because a significant effect on the sharpness of the peak of the frequency-correlation function. That is, Equation 22 shows that the spread of the Gaussian function increases as g 3 increases.
  • Equation 23 the probability density function of the logistic distribution of Expression 23 below may be used as a kernel function.
  • Equation 23 g 4 and g 5 are positive real numbers.
  • the probability density function of the logistic distribution has the same shape as the Gaussian function, but has a longer tail than the Gaussian function.
  • g 5 for adjusting the spread of the probability density function of the logistic distribution similarly to the case of g 3 in the Gaussian function of Equation 22, is a parameter which greatly affects the sharpness of the peak of the frequency-correlation function.
  • a cosine function or a uniform function may be used for the kernel function.
  • g 3 and g 5 that affect the spread of the kernel function are determined using the variance input from the variance calculation unit 154.
  • these parameters are called spread control parameters and are expressed as q (k, n). Therefore, if the kernel function is a Gaussian function, g 3 is q (k, n). If the variance is small, the parameters are changed so that the peak of the frequency-dependent correlation function is sharp and the tail is narrow. Therefore, the spread control parameter is reduced.
  • the spread control parameter can be calculated by converting a variance value using a predetermined mapping function. For example, if the variance exceeds a certain threshold, the spread control parameter is set to a large value (for example, 10), and if the variance is below a certain threshold, the spread control parameter is set to a small value (for example, 0.01).
  • the variance is V 12 (k, n)
  • the threshold is p th .
  • the spread control parameter q (k, n) in the frequency bin k of the n-th frame can be calculated using the following Expression 24. However, in Expression 24, q 1 and q 2 are positive real numbers that satisfy q 1 > q 2 .
  • the spread control parameter q (k, n) may be calculated using a linear function as in the following Expression 25.
  • q 3 is a real number larger than 0, and q 4 is a real number.
  • Equation 27 the values shown in Equation 27 may be used.
  • L represents the number of frames averaged when the average calculation unit 153 calculates the average cross spectrum. Since the error of the average cross spectrum is inversely proportional to the number L of averaged frames, the spread control parameter can be obtained by using Equations 26 and 27 in consideration of the error (reliability) of the average cross spectrum.
  • a variance function represented by a linear mapping function, a higher-order polynomial function, a non-linear function, or the like can be used for the variance calculation.
  • the variance may be used as it is as the spread control parameter.
  • the function for obtaining the spread control parameter may be a function of the frequency k as well as the variance. For example, a function that decreases as the frequency k increases can be used. As a typical example, an example using the reciprocal of k is given. In this case, the spread control parameter q (k, n) can be calculated using the function of the following equation 28 instead of the equation 25.
  • the spread control parameter q (k, n) can be calculated using the function of Expression 29 below.
  • the multiplier 553 receives the basic cross spectrum for each frequency from the basic cross spectrum calculator for frequency 551, and receives the kernel function spectrum from the kernel function spectrum generator 552.
  • the multiplication unit 553 calculates a product of the basic cross spectrum for each frequency supplied from the basic cross spectrum for each frequency calculation unit 551 and the kernel function spectrum supplied from the kernel function spectrum generation unit 552 to calculate the cross spectrum for each frequency. I do.
  • the multiplying unit 553 outputs the calculated cross spectrum for each frequency to the inverse transform unit 156.
  • the basic cross spectrum for each frequency supplied from the basic cross spectrum calculation unit for frequency 551 is U k (w, n), and the kernel function spectrum supplied from the kernel function spectrum generation unit 552 is G (w).
  • the multiplication unit 553 calculates the frequency-specific cross spectrum UM k (w, n) using Expression 30 below.
  • the inverse transform unit 156 receives the frequency-specific cross spectrum from the multiplication unit 553 of the frequency-specific cross spectrum generation unit 155. For example, when the transform unit 151 uses Fourier transform, the inverse transform unit 156 performs inverse transform using inverse Fourier transform. The inverse transform unit 156 calculates an inverse transform of the frequency-specific cross spectrum supplied from the frequency-specific cross spectrum generation unit 155.
  • the frequency-specific cross spectrum supplied from the frequency-specific cross spectrum generation unit 155 is UM k (w, n).
  • the inverse transform unit 156 inversely transforms UM k (w, n) and calculates a frequency-dependent cross-correlation function u k ( ⁇ , n) using Expression 31 below.
  • the frequency-dependent estimated direction information calculator 157 receives the frequency-dependent cross-correlation function from the inverse transformer 156 and the relative delay time from the relative delay time calculator 13.
  • the frequency-specific estimated direction information calculation unit 157 uses the frequency-dependent cross-correlation function supplied from the inverse transform unit 156 and the relative delay time supplied from the relative delay time calculation unit 13 to determine the correspondence between the direction and the correlation value. The relationship is obtained as frequency-wise estimated direction information.
  • the frequency-specific estimated direction information calculation unit 157 outputs the obtained frequency-specific estimated direction information to the integration unit 17.
  • the frequency-specific estimated direction information calculation unit 157 calculates the frequency-specific estimated direction information H k, r ( ⁇ , n) using Expression 32 below.
  • the integrating unit 17 receives estimated frequency direction information from the frequency estimated direction information generating units 15-1 to 15-R.
  • the integration unit 17 integrates the frequency-specific estimated direction information supplied from the frequency-specific estimated direction information generation units 15-1 to 15-R to calculate integrated estimated direction information.
  • the integrating unit 17 obtains one piece of estimated direction information by mixing or superimposing a plurality of individually obtained estimated direction information for each frequency.
  • the integrating unit 17 outputs the calculated integrated estimated direction information. For example, the integration unit 17 outputs integrated estimation direction information to a higher-level system (not shown).
  • the integrating unit 17 first integrates the frequency-specific estimated direction information H k, r ( ⁇ , n) for the number of combinations (R) of the input signals, thereby integrating the frequency-specific integrated estimated direction information H k ( ⁇ , n). n) is calculated. Then, the integrating unit 17 calculates integrated estimated direction information H ( ⁇ , n) by integrating the calculated frequency-based integrated estimated direction information for all frequencies.
  • the integrating unit 17 calculates the frequency-based integrated estimated direction information H k ( ⁇ , n) by calculating the total power of the frequency-based estimated direction information H k, r ( ⁇ , n). At this time, the integrating unit 17 calculates the frequency-based integrated estimation direction information H k ( ⁇ , n) using Expression 33 below.
  • the integrating unit 17 may calculate the frequency-specific integrated estimation direction information H k ( ⁇ , n) by calculating the sum of the frequency-specific estimated direction information H k, r ( ⁇ , n). At this time, the integration unit 17 calculates the integrated estimation direction information H k ( ⁇ , n) for each frequency using the following Expression 34.
  • the integration unit 17 calculates the sum or the square of the frequency-based integrated estimated direction information H k ( ⁇ , n) for each frequency.
  • the integration unit 17 calculates the sum of the frequency-based integrated estimation direction information H k ( ⁇ , n) for the frequency k as the integrated estimation direction information H ( ⁇ , n) using Expression 35 below.
  • the integrating unit 17 calculates the total power of the frequency-based integrated estimation direction information H k ( ⁇ , n) for the frequency k as the integrated estimation direction information H ( ⁇ , n) using Expression 36 below. I do.
  • the integrating unit 17 obtains the integrated estimated direction information using only the frequency-based integrated estimated direction information corresponding to the frequency. May be.
  • the integrating unit 17 may control the degree of influence of the frequency-based integrated estimation direction information on the integration in the form of weighting. For example, assuming that a set of frequencies at which the target sound exists is ⁇ , the integration unit 17 can obtain integrated estimation direction information H ( ⁇ , n) by selecting a frequency using Expression 37 below.
  • the integrating unit 17 can calculate integrated estimated direction information H ( ⁇ , n) using Expression 38 below.
  • a and b are real numbers satisfying a>b> 0.
  • the integrating unit 17 may calculate the integrated estimated direction information H ( ⁇ , n) using another calculation method. For example, the integrating unit 17 first calculates integrated estimation direction information H r ( ⁇ , n) for each input signal in which the estimation direction information H k, r ( ⁇ , n) for each frequency is integrated for all frequencies. Then, the integrating unit 17 may calculate integrated estimated direction information H ( ⁇ , n) obtained by integrating integrated estimated direction information for each input signal combination for all combinations of input signals.
  • each sensor 110 is connected to one of the input terminals 11 of the wave source direction estimation device 10 via a network such as the Internet or an intranet or a cable.
  • the sensor 110 when detecting the sound wave, is realized by a microphone.
  • the sensor 110 is realized by a vibration sensor when detecting a vibration wave.
  • the sensor 110 is realized by an antenna when detecting an electromagnetic wave. Note that the form of the sensor 110 is not limited as long as it can convert a wave to be detected into an electric signal.
  • the wave source direction estimating apparatus 10 inputs electric signals (also referred to as input signals) from a plurality of microphones (step S111).
  • the wave source direction estimating apparatus 10 selects two input signals from input signals corresponding to a plurality of microphones (Step S112).
  • the wave source direction estimating apparatus 10 calculates a relative delay time based on an interval (also referred to as a microphone interval) between two microphones that are two selected input signal supply sources and a set sound source search target direction. (Step S113).
  • the wave source direction estimating apparatus 10 generates estimated direction information (also referred to as frequency-based estimated direction information) for each frequency using the two selected input signals and the relative delay time (step S114).
  • estimated direction information also referred to as frequency-based estimated direction information
  • the wave source direction estimation device 10 integrates the estimated direction information for each frequency to calculate integrated estimated direction information (step S115).
  • the wave source direction estimating apparatus 10 outputs integrated estimated direction information (step S116).
  • the frequency-specific estimated direction information generation unit 15 inputs the two input signals selected by the signal selection unit 12 and the relative delay times of those input signals (step S121).
  • the frequency-specific estimated direction information generating unit 15 converts the two input signals into a frequency domain signal (also referred to as a converted signal) (Step S122).
  • the frequency-specific estimated direction information generator 15 calculates a cross spectrum using the converted signal (step S123).
  • the frequency-specific estimated direction information generation unit 15 calculates an average cross spectrum using the cross spectrum (step S124).
  • the frequency-specific estimated direction information generator 15 calculates the variance using the average cross spectrum (step S125).
  • the frequency-specific estimated direction information generation unit 15 calculates a frequency-specific cross spectrum using the average cross spectrum and the variance (step S126).
  • the frequency-specific estimated direction information generating unit 15 calculates a frequency-specific cross-correlation function using the frequency-specific cross spectrum (step S127).
  • the frequency-specific estimated direction information generation unit 15 calculates frequency-specific estimated direction information using the frequency-dependent cross-correlation function and the relative delay time (step S128).
  • the frequency-specific estimated direction information generation unit 15 outputs the frequency-specific estimated direction information to the integration unit 17 (Step S129).
  • the frequency-specific cross spectrum generation unit 155 inputs the average cross spectrum from the average calculation unit 153 and the variance from the variance calculation unit 154 (step S131).
  • the frequency-specific cross spectrum generation unit 155 calculates a frequency-specific basic cross spectrum using the average cross spectrum (step S132).
  • Frequency-based cross spectrum generating section 155 calculates a kernel function spectrum using the variance (step S133). Note that the processing of step S132 and the processing of step S133 may be performed in parallel or sequentially.
  • the frequency-specific cross spectrum generation unit 155 calculates the product of the frequency-specific basic cross spectrum and the kernel function spectrum, and calculates the frequency-specific cross spectrum (step S134).
  • the frequency-specific cross spectrum generation unit 155 outputs the calculated frequency-specific cross spectrum to the inverse transform unit 156 (step S135).
  • the wave source direction estimating apparatus of the present embodiment includes a plurality of input units, a signal selecting unit, a relative delay time calculating unit, at least one frequency-based estimated direction information generating unit, and an integrating unit.
  • the plurality of input units obtain, as input signals, electric signals converted from the waves obtained by the plurality of sensors.
  • the signal selection unit selects at least two pairs obtained by combining at least two input signals from the plurality of input signals.
  • the relative delay time calculation unit calculates, as a relative delay time, the arrival time difference of the wave in each of the wave source search directions between at least two input signals forming a pair of input signals.
  • At least one frequency-dependent estimated direction information generation unit generates, for each frequency, estimated direction information of a wave source of a wave using a pair of input signals and a relative delay time.
  • the integrating unit integrates the estimated direction information for each frequency generated by the frequency-specific estimated direction information generating unit.
  • the signal selection unit selects a pair in which at least two input signals are combined from a plurality of input signals based on the sensor interval.
  • the relative delay time calculation unit based on the source search direction for a pair of sensors that are the source of one input signal pair, the relative delay time of all input signal pairs selected by the signal selection unit, It is calculated as a function of the reference wave source search direction.
  • the frequency-specific estimated direction information generation unit includes a conversion unit, a cross spectrum calculation unit, an average calculation unit, a variance calculation unit, a frequency-specific cross spectrum generation unit, an inverse conversion unit, and an estimated direction information calculation unit.
  • the conversion unit converts at least two pairs of input signals into a frequency-domain conversion signal.
  • the cross spectrum calculator calculates a cross spectrum using the converted signal converted by the conversion unit.
  • the average calculator calculates an average cross spectrum using the cross spectrum calculated by the cross spectrum calculator.
  • the variance calculator calculates the variance using the average cross spectrum calculated by the average calculator.
  • the frequency-specific cross spectrum generation unit calculates a frequency-specific cross spectrum using the average cross spectrum calculated by the average calculation unit and the variance calculated by the variance calculation unit.
  • the inverse transform unit performs an inverse transform of the frequency-specific cross spectrum calculated by the frequency-specific cross spectrum generation unit to calculate a frequency-specific cross-correlation function.
  • the estimated direction information calculation unit calculates estimated direction information for each frequency-based estimated frequency using the frequency-dependent cross-correlation function calculated by the inverse transform unit and the relative delay time.
  • the frequency-specific cross spectrum generation unit includes a frequency-specific basic cross spectrum calculation unit, a kernel function spectrum generation unit, and a multiplication unit.
  • the frequency-specific basic cross spectrum calculation unit obtains the average cross spectrum from the average calculation unit, and calculates a frequency-specific basic cross spectrum using the obtained average cross spectrum.
  • the kernel function spectrum generator acquires a variance from the variance calculator, and calculates a kernel function spectrum using the acquired variance.
  • the multiplication unit calculates the product of the frequency-specific basic cross spectrum calculated by the frequency-specific basic cross spectrum calculation unit and the kernel function spectrum calculated by the kernel function spectrum generation unit to calculate the frequency-specific cross spectrum.
  • the integration unit calculates integrated estimation direction information for each frequency by integrating estimation direction information for each of a plurality of frequencies for a plurality of pairs of input signals. Then, the integration unit integrates the calculated integrated estimation direction information for each frequency for all frequencies to calculate integrated estimation direction information.
  • the integration unit calculates integrated estimation direction information for each input signal combination in which the estimation direction information for each of a plurality of frequencies is integrated for all frequencies.
  • the integrating unit integrates the calculated integrated estimation direction information for each input signal combination for all combinations of input signals to calculate integrated estimation direction information.
  • the wave source direction estimating device includes a sensor arranged corresponding to each of the plurality of input units.
  • the wave source direction estimating device of the present embodiment obtains estimated direction information from a cross-correlation function between microphone pairs, and integrates the estimated direction information among a plurality of microphone pairs.
  • the wave source direction estimating device of the present embodiment the false peak of the estimated direction information in the direction other than the sound source direction, which is caused by the accidental alignment of the phases between the microphone pairs, is reduced, and the erroneous detection of the virtual image sound source is reduced.
  • Generation of estimation can be reduced, and the direction of the sound source can be estimated with high accuracy.
  • the estimation target of the wave source direction estimating device of the present embodiment is not limited to a sound source that is a vibration wave of air or water.
  • the wave source direction estimating apparatus according to the present embodiment can also be applied to estimating the direction of a generation source of a vibration wave using a solid as a medium such as an earthquake or a landslide.
  • a vibration sensor can be used instead of a microphone for a device that converts a vibration wave into an electric signal.
  • the wave source direction estimating apparatus of the present embodiment can be applied to the case of estimating the direction using not only the vibration waves of gas, liquid, and solid but also radio waves.
  • an antenna can be used as a device that converts radio waves into electric signals.
  • the integrated estimation direction information estimated by the wave source direction estimation device of the present embodiment can be used in various forms. For example, when the integrated estimated direction information has a plurality of peaks, it is estimated that there are a plurality of sound sources having each peak as an arrival direction. Therefore, if the integrated estimated direction information is used, not only can the direction of each sound source be estimated at the same time, but also the number of sound sources can be estimated.
  • the wave source direction estimating device of the present embodiment has a configuration in which a wave source direction calculating unit is added to the wave source direction estimating device of the first embodiment.
  • FIG. 10 is a block diagram showing the configuration of the wave source direction estimating device 20 of the present embodiment.
  • the wave source direction estimating device 20 includes an input terminal 21, a signal selecting unit 22, a relative delay time calculating unit 23, a frequency-specific estimated direction information generating unit 25, an integrating unit 27, and a wave source direction calculating unit 28.
  • the input terminal 21, the signal selection unit 22, the relative delay time calculation unit 23, the frequency-specific estimated direction information generation unit 25, and the integration unit 27 are the same as the corresponding components of the wave source direction estimation device 10 of the first embodiment. Therefore, detailed description is omitted.
  • the integrated estimated direction information is input to the wave source direction calculation unit 28 from the integration unit 27.
  • the wave source direction calculation unit 28 calculates a wave source direction using the integrated estimated direction information.
  • the wave source direction calculator 28 outputs the calculated wave source direction.
  • the calculation method of the wave source direction in the wave source direction calculation unit 28 will be described in detail below.
  • the integrated estimation direction information input from the integrating unit 27 has higher reliability (the possibility of the existence of a sound source) as the peak is larger. Therefore, for example, when the number of sound sources can be assumed to be one in advance, the wave source direction calculation unit 28 outputs the direction in which the integrated estimated direction information is the maximum as the estimated direction.
  • the integrated estimated direction information input from the integration unit 27 is set to H ( ⁇ , n).
  • the wave source direction calculation unit 28 uses the following Expression 39 to calculate a set having an argument of the integrated estimated direction information H ( ⁇ , n) such that the integrated estimated direction information H ( ⁇ , n) takes the maximum value. It can be calculated as the wave source direction ⁇ .
  • ⁇ represents all the wave source directions or wave source direction candidates.
  • the wave source direction calculation unit 28 may regard the direction having the peak exceeding the threshold value as the sound source and output the direction exceeding the threshold value as the estimated direction.
  • the wave source direction estimating device of the present embodiment can also estimate, as the sound source direction, the direction corresponding to the time at which the integrated estimated direction information becomes maximum at regular time intervals T. However, it is assumed that the direction of the sound source does not change during the fixed time T, or that the magnitude of the change is negligibly small. By making such an assumption, the estimation accuracy of the wave source direction can be improved.
  • the wave source direction estimating device of the present embodiment includes the wave source direction calculating means for calculating the wave source direction of the wave based on the integrated estimated direction information calculated by the integrating means.
  • the wave source direction calculating means calculates a direction corresponding to a time at which the integrated estimated direction information becomes maximum at regular time intervals as a wave source direction.
  • the direction of the sound source can be estimated with high accuracy without erroneously estimating the virtual image sound source.
  • the information processing device 90 includes a processor 91, a main storage device 92, an auxiliary storage device 93, an input / output interface 95, and a communication interface 96.
  • the interface is abbreviated as I / F (Interface).
  • the processor 91, the main storage device 92, the auxiliary storage device 93, the input / output interface 95, and the communication interface 96 are connected to each other via a bus 99 so as to be able to perform data communication.
  • the processor 91, the main storage device 92, the auxiliary storage device 93, and the input / output interface 95 are connected to a network such as the Internet or an intranet via a communication interface 96.
  • the processor 91 loads the program stored in the auxiliary storage device 93 or the like into the main storage device 92 and executes the loaded program.
  • a configuration using a software program installed in the information processing device 90 may be used.
  • the processor 91 executes processing by the wave source direction estimating device according to the present embodiment.
  • the main storage device 92 has an area where the program is expanded.
  • the main storage device 92 may be a volatile memory such as a DRAM (Dynamic Random Access Memory). Further, a nonvolatile memory such as an MRAM (Magnetoresistive Random Access Memory) may be configured and added as the main storage device 92.
  • DRAM Dynamic Random Access Memory
  • MRAM Magnetic Random Access Memory
  • the auxiliary storage device 93 stores various data.
  • the auxiliary storage device 93 is configured by a local disk such as a hard disk or a flash memory. Note that it is also possible to store various data in the main storage device 92 and omit the auxiliary storage device 93.
  • the input / output interface 95 is an interface for connecting the information processing device 90 and peripheral devices.
  • the communication interface 96 is an interface for connecting to an external system or device via a network such as the Internet or an intranet based on standards and specifications.
  • the input / output interface 95 and the communication interface 96 may be shared as an interface for connecting to an external device.
  • the information processing device 90 may be configured to connect input devices such as a keyboard, a mouse, and a touch panel as needed. These input devices are used for inputting information and settings.
  • input devices such as a keyboard, a mouse, and a touch panel as needed.
  • the display screen of the display device may be configured to also serve as an interface of the input device.
  • Data communication between the processor 91 and the input device may be mediated by the input / output interface 95.
  • the information processing device 90 may be provided with a display device for displaying information.
  • the information processing device 90 includes a display control device (not shown) for controlling display on the display device.
  • the display device may be connected to the information processing device 90 via the input / output interface 95.
  • the information processing device 90 may be provided with a disk drive as needed.
  • the disk drive is connected to the bus 99.
  • the disk drive mediates between a processor 91 and a recording medium (program recording medium) not shown, such as reading of a data program from the recording medium and writing of a processing result of the information processing device 90 to the recording medium.
  • the recording medium can be realized by an optical recording medium such as a CD (Compact Disc) and a DVD (Digital Versatile Disc).
  • the recording medium may be realized by a semiconductor recording medium such as a USB (Universal Serial Bus) memory or an SD (Secure Digital) card, a magnetic recording medium such as a flexible disk, or another recording medium.
  • USB Universal Serial Bus
  • SD Secure Digital
  • the above is an example of a hardware configuration for enabling the wave source direction estimation device according to each embodiment.
  • the hardware configuration in FIG. 11 is an example of a hardware configuration for executing the arithmetic processing of the wave source direction estimation device according to each embodiment, and does not limit the scope of the present invention.
  • a program that causes a computer to execute the process related to the wave source direction estimation device according to each embodiment is also included in the scope of the present invention.
  • a program recording medium on which a program according to each embodiment is recorded is also included in the scope of the present invention.
  • the components of the wave source direction estimating device of each embodiment can be arbitrarily combined. Also, the components of the wave source direction estimation device of each embodiment may be realized by software or may be realized by a circuit.
  • a plurality of input means for obtaining as an input signal an electric signal converted from the wave obtained by the plurality of sensors, Signal selecting means for selecting at least two pairs of a combination of at least two of the input signals from a plurality of the input signals, Relative delay time calculating means for calculating as the relative delay time the arrival time difference of the wave for each wave source search direction between at least two of the input signals constituting the pair of input signals, At least one frequency-based estimated direction information generating unit that generates estimated direction information of the wave source of the wave for each frequency using the pair of the input signals and the relative delay time,
  • a wave source direction estimating device comprising: integrating means for integrating the estimated direction information for each frequency generated by the frequency-specific estimated direction information generating means.
  • the frequency-specific estimated direction information generating means Conversion means for converting at least two of the paired input signals into a frequency-domain conversion signal; Cross spectrum calculation means for calculating a cross spectrum using the converted signal converted by the conversion means, Average calculation means for calculating an average cross spectrum using the cross spectrum calculated by the cross spectrum calculation means, A variance calculation unit that calculates variance using the average cross spectrum calculated by the average calculation unit, The average cross spectrum calculated by the average calculation means, and a frequency-specific cross spectrum generation means for calculating a frequency-specific cross spectrum using the variance calculated by the variance calculation means, Inverse transforming means for inversely transforming the frequency-specific cross spectrum calculated by the frequency-specific cross spectrum generating means to calculate a frequency-specific cross-correlation function, Any one of Supplementary notes 1 to 3, further comprising: frequency-specific estimated direction information calculating means for calculating the estimated direction information for each frequency using the frequency-dependent cross-correlation function calculated by the inverse transform means and the relative delay time.
  • a wave source direction estimating apparatus according to claim 1.
  • the frequency-specific cross spectrum generating means Obtain the average cross spectrum from the average calculation means, frequency-based basic cross spectrum calculation means to calculate a frequency-based basic cross spectrum using the obtained average cross spectrum, Kernel function spectrum generating means for obtaining the variance from the variance calculating means and calculating a kernel function spectrum using the obtained variance, Multiplication that calculates the product of the frequency-specific basic cross spectrum calculated by the frequency-specific basic cross spectrum calculation means and the kernel function spectrum calculated by the kernel function spectrum generation means to calculate the frequency-specific cross spectrum
  • a wave source direction estimating apparatus according to supplementary note 4 including means.
  • the integrating means includes: Calculating frequency-based integrated estimation direction information in which the estimation direction information for each of a plurality of frequencies is integrated for a plurality of input signal pairs, and integrating the calculated frequency-based integrated estimation direction information for all frequencies to obtain an integrated estimation direction
  • the wave source direction estimating apparatus according to any one of supplementary notes 1 to 5, which calculates information.
  • the integration means includes: The integrated estimation direction information for each input signal combination in which the estimation direction information for each of the plurality of frequencies is integrated for all frequencies is calculated, and the calculated integrated estimation direction information for each input signal combination is integrated for all combinations of input signals. 6.
  • the wave source direction estimating apparatus according to any one of supplementary notes 1 to 5, which calculates integrated estimated direction information.
  • the information processing device is Obtain an electric signal converted from the wave obtained by the plurality of sensors as an input signal, Selecting at least two pairs that combine at least two of the input signals from a plurality of the input signals, Calculating as a relative delay time the arrival time difference of the wave for each wave source search direction between at least two of the input signals constituting the pair of the input signals, Using the input signal pair and the relative delay time, generate at least one estimated direction information of the wave source of the wave for each frequency, A wave source direction estimating method for integrating the estimated direction information for each frequency.
  • (Appendix 12) A process of obtaining an electric signal converted from the wave obtained by the plurality of sensors as an input signal, A process of selecting at least two pairs in which at least two of the input signals are combined from a plurality of the input signals, A process of calculating the arrival time difference of the wave for each wave source search direction between at least two of the input signals constituting the pair of the input signals as a relative delay time, A process of generating at least one estimated direction information of the wave source of the wave for each frequency using the pair of the input signals and the relative delay time, A program recording medium storing a process for integrating the estimated direction information for each frequency and a program to be executed by a computer.

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Health & Medical Sciences (AREA)
  • Otolaryngology (AREA)
  • Signal Processing (AREA)
  • Acoustics & Sound (AREA)
  • General Physics & Mathematics (AREA)
  • Remote Sensing (AREA)
  • Radar, Positioning & Navigation (AREA)
  • General Health & Medical Sciences (AREA)
  • Spectroscopy & Molecular Physics (AREA)
  • Computational Linguistics (AREA)
  • Audiology, Speech & Language Pathology (AREA)
  • Human Computer Interaction (AREA)
  • Multimedia (AREA)
  • Measurement Of Velocity Or Position Using Acoustic Or Ultrasonic Waves (AREA)
  • Measurement Of Mechanical Vibrations Or Ultrasonic Waves (AREA)
  • Circuit For Audible Band Transducer (AREA)

Abstract

虚像音源の誤推定の発生を低減し、高精度に音源の方向を推定するために、複数のセンサによって取得された波動から変換される電気信号を入力信号として取得する複数の入力部と、複数の入力信号から少なくとも二つの入力信号を組み合わせたペアを少なくとも二つ選択する信号選択部と、入力信号のペアを構成する少なくとも二つの入力信号間において波源探索方向ごとの波動の到達時間差を相対遅延時間として算出する相対遅延時間計算部と、入力信号のペアと相対遅延時間とを用いて波動の波源の推定方向情報を周波数ごとに生成する少なくとも一つの周波数別推定方向情報生成部と、周波数別推定方向情報生成部によって生成される周波数ごとの推定方向情報を統合する統合部とを備える波源方向推定装置とする。

Description

波源方向推定装置、波源方向推定方法、およびプログラム記録媒体
　本発明は、波源方向推定装置、波源方向推定方法、およびプログラムに関する。特に、本発明は、複数のセンサが取得した信号に基づいて波源方向を推定する波源方向推定装置、波源方向推定方法、およびプログラムに関する。
　特許文献１および非特許文献１には、２つのマイクロフォンの受音信号の到達時間差から音源の方向を推定する方法が開示されている。特許文献１および非特許文献１の方法では、音波の到達時間差の確率密度関数を周波数別に求め、それらの重ね合わせにより得られた確率密度関数から到達時間差を算出し、音源方向を推定する。
　特許文献２には、所定の観測点に伝搬される音と振動とを採取して、音の音源が振動源から発生した音か否かを探査する探査方法について開示されている。特許文献２の方法では、音源から伝搬する音と振動源から伝搬する表面波の振動とを同時に測定する。そして、音の音圧レベルのデータから得られた音源の方向と振動の振動レベルのデータから得られた振動源の方向とを比較し、音源からの音が音の発生を伴う振動源からの音かどうかを判定する。
国際公開第２０１８／００３１５８号 特開２０１０－２３６９４４号公報
M. Kato, Y. Senda, R. Kondo, "TDOA estimation based on phase-voting cross correlation and circular standard deviation," 25th European Signal Processing Conference (EUSIPCO), EURASIP, August 2017, p.1230-1234
　特許文献１および非特許文献１の方法によれば、信号対雑音比（ＳＮＲ：Signal-to-Noise Ratio）が高い周波数帯域では、到達時間差の確率密度関数が鋭いピークを形成するため、高ＳＮＲ帯域が少なくても精度よく到達時間差を推定することができる。しかしながら、特許文献１および非特許文献１の方法においては、周波数別の到達時間差の確率密度関数を重ね合わせる際に、音源が存在しなくても、偶然位相が揃うことで重ね合わせた確率密度関数にピークが生成される。そのため、特許文献１および非特許文献１の方法には、虚像音源を誤推定してしまうという問題点があった。
　特許文献２の方法によれば、音源からの音が音の発生を伴う振動源からの音か、あるいは、振動を伴わない音源からの音か否かの判定や、振動源が音を伴わない振動源であるか否かの判定を的確に行うことができる。しかしながら、特許文献２の方法では、異なるマイクロフォンの間で位相が偶然揃うことによって音源とは異なる方向の虚像音源の到達時間差を算出し、虚像音源を誤推定する可能性があるという問題点があった。
　本発明の目的は、上述した課題を解決し、虚像音源の誤推定の発生を低減し、高精度に音源の方向を推定することができる波源方向推定装置を提供することにある。
　本発明の一態様の波源方向推定装置は、複数のセンサによって取得された波動から変換される電気信号を入力信号として取得する複数の入力部と、複数の入力信号から少なくとも二つの入力信号を組み合わせたペアを少なくとも二つ選択する信号選択部と、入力信号のペアを構成する少なくとも二つの入力信号間において波源探索方向ごとの波動の到達時間差を相対遅延時間として算出する相対遅延時間計算部と、入力信号のペアと相対遅延時間とを用いて波動の波源の推定方向情報を周波数ごとに生成する少なくとも一つの周波数別推定方向情報生成部と、周波数別推定方向情報生成部によって生成される周波数ごとの推定方向情報を統合する統合部とを備える。
　本発明の一態様の波源方向推定方法においては、情報処理装置が、複数のセンサによって取得された波動から変換される電気信号を入力信号として取得し、複数の入力信号から少なくとも二つの入力信号を組み合わせたペアを少なくとも二つ選択し、入力信号のペアを構成する少なくとも二つの入力信号間において波源探索方向ごとの波動の到達時間差を相対遅延時間として算出し、入力信号のペアと相対遅延時間とを用いて波動の波源の推定方向情報を周波数ごとに少なくとも一つ生成し、周波数ごとの推定方向情報を統合する。
　本発明の一態様のプログラムは、複数のセンサによって取得された波動から変換される電気信号を入力信号として取得する処理と、複数の入力信号から少なくとも二つの入力信号を組み合わせたペアを少なくとも二つ選択する処理と、入力信号のペアを構成する少なくとも二つの入力信号間において波源探索方向ごとの波動の到達時間差を相対遅延時間として算出する処理と、入力信号のペアと相対遅延時間とを用いて波動の波源の推定方向情報を周波数ごとに少なくとも一つ生成する処理と、周波数ごとの推定方向情報を統合する処理とコンピュータに実行させる。
　本発明によれば、虚像音源の誤推定の発生を低減し、高精度に音源の方向を推定することができる波源方向推定装置を提供することが可能になる。
本発明の第１の実施形態に係る波源方向推定装置の構成の一例を示すブロック図である。 本発明の第１の実施形態に係る波源方向推定装置における相対遅延時間計算部の処理の一例について説明するための概念図である。 本発明の第１の実施形態に係る波源方向推定装置における相対遅延時間計算部の処理の別の一例について説明するための概念図である。 本発明の第１の実施形態に係る波源方向推定装置が有する周波数別推定方向情報生成部の構成の一例を示すブロック図である。 本発明の第１の実施形態に係る波源方向推定装置が有する周波数別クロスススペクトル生成部の構成の一例を示すブロック図である。 本発明の第１の実施形態に係る波源方向推定装置に少なくとも一つのセンサを追加した構成の一例を示すブロック図である。 本発明の第１の実施形態に係る波源方向推定装置の動作の概略について説明するためのフローチャートである。 本発明の第１の実施形態に係る波源方向推定装置の周波数別推定方向情報生成部の動作について説明するためのフローチャートである。 本発明の第１の実施形態に係る波源方向推定装置の周波数別推定方向情報生成部の周波数別クロススペクトル生成部の動作について説明するためのフローチャートである。 本発明の第２の実施形態に係る波源方向推定装置の構成の一例を示すブロック図である。 本発明の各実施形態に係る波源方向推定装置を実現するハードウェア構成の一例を示すブロック図である。
　以下に、本発明を実施するための形態について図面を用いて説明する。ただし、以下に述べる実施形態には、本発明を実施するために技術的に好ましい限定がされているが、発明の範囲を以下に限定するものではない。なお、以下の実施形態の説明に用いる全図においては、特に理由がない限り、同様箇所には同一符号を付す。また、以下の実施形態において、同様の構成・動作に関しては繰り返しの説明を省略する場合がある。また、図面中の矢印の向きは、一例を示すものであり、ブロック間の信号の向きを限定するものではない。
　（第１の実施形態）
　まず、本発明の第１の実施形態に係る波源方向推定装置について図面を参照しながら説明する。以下においては、本実施形態の波源方向推定装置が、空気や水の振動波である音波の発生源を推定する例について説明する。そのため、本実施形態の波源方向推定装置は、マイクロフォンによって電気信号に変換された振動波を検証する。なお、本実施形態の波源方向推定装置の推定対象は、音波の発生源に限定されず、振動波や電磁波などの任意の波動の発生源（波源とも呼ぶ）を推定することにも使用できる。
　（構成）
　図１は、本実施形態の波源方向推定装置１０の構成を表すブロック図である。波源方向推定装置１０は、入力端子１１、信号選択部１２、相対遅延時間計算部１３、周波数別推定方向情報生成部１５、および統合部１７を備える。
　波源方向推定装置１０は、ｐ個の入力端子１１を備える（ｐは２以上の整数）。また、波源方向推定装置１０は、Ｒ個の周波数別推定方向情報生成部１５を備える（Ｒは１以上の整数）。図１においては、個々の入力端子１１を区別するために、符号の後にハイフンを挟んで１～ｐの番号を付す。同様に、図１においては、個々の周波数別推定方向情報生成部１５を区別するために、符号の後にハイフンを挟んで１～Ｒの番号を付す。
　〔入力端子〕
　入力端子１１－１～ｐ（入力部とも呼ぶ）のそれぞれは、図示しないマイクロフォン（以下、マイクとも呼ぶ）に接続される。入力端子１１－１～ｐのそれぞれには、異なる位置に配置されたマイクによって集音された音波（音信号とも呼ぶ）から変換された電気信号が入力信号として入力される。以下において、時刻ｔにおいてｍ番目の入力端子１１－ｍに入力した入力信号をｘm（ｔ）と記載する（ｔ：実数、ｍ：１以上ｐ以下の整数）。
　マイクは、目標音源で発生した音と、そのマイクの周囲で発生した様々な雑音とが混在した音波を集音し、集音した音波をデジタル信号（サンプル値系列とも呼ぶ）に変換する集音装置である。マイクは、目標音源からの音を集音するために、入力端子１１－１～ｐのそれぞれに対応させて異なる位置に配置される。以下においては、ｍ番目のマイクによって集音された音波から変換された入力信号は、ｍ番目の入力端子１１－ｍに供給されるものとする。また、以下において、ｍ番目の入力端子１１－ｍに供給された入力信号のことを「ｍ番目のマイクの入力信号」とも呼ぶ。
　〔信号選択部〕
　信号選択部１２は、入力端子１１－１～ｐに供給されるＰ個の入力信号のうちの２個の入力信号を選択する。信号選択部１２は、選択した２個の入力信号を周波数別推定方向情報生成部１５－１～Ｒに出力し、それらの入力信号の供給元であるマイクの位置情報（以下、マイク位置情報とも呼ぶ）を相対遅延時間計算部１３に出力する。ここで、周波数別推定方向情報生成部１５の数Ｒは、入力信号の組み合わせの数Ｒに相当する。信号選択部１２は、２個の入力信号を選択する際に、全ての組み合わせを選択してもよいし、一部の組み合わせを選択してもよい。全ての組み合わせを選択する場合、Ｒは以下の式１で表される。
Figure JPOXMLDOC01-appb-I000001
　波源方向推定装置１０は、目標音源からの音が２つのマイクに到達する時間差を利用して音源の方向を推定する。マイクの間隔（以下、マイク間隔とも呼ぶ）が大きすぎると、目標音源からの音が空気や水などの媒質の影響で同一の音として観測されないために方向推定精度が低下する。また、マイク間隔が小さすぎると、２つのマイクの間の音波の到達時間差が小さくなりすぎるために方向推定精度が低下する。そのため、信号選択部１２は、式２に示すように、マイク間隔ｄが一定の範囲内に収まるマイクの入力信号を選択するのがよい（ｄmin、ｄmax：実数）。
Figure JPOXMLDOC01-appb-I000002
　例えば、信号選択部１２は、マイク間隔ｄが十分に小さい場合、マイク間隔ｄが最大の２個の入力信号を選択してもよい。また、信号選択部１２は、マイク間隔ｄが十分に小さい場合、マイク間隔ｄが大きい方から順に並べ、上位Ｒ位（ｒ＜Ｃ（ｐ，２））までの入力信号の組み合わせを選択してもよい。このように、信号選択部１２が一部の組み合わせを選択することは、方向推定精度の低下を防ぐこと以外に、計算量を低減することに繋がる。
　目標音源からの音の２つのマイクへの到達時間差を求める際には、マイク位置情報も重要である。そのため、信号選択部１２は、入力信号に加えて、マイク位置情報を相対遅延時間計算部１３に出力する。
　〔相対遅延時間計算部〕
　相対遅延時間計算部１３には、信号選択部１２からマイク位置情報が入力される。相対遅延時間計算部１３は、マイク位置情報と音源探索対象方向とを用いて、信号選択部１２が選択した全てのマイクペア間の相対遅延時間を算出する。相対遅延時間とは、マイク間隔と音源方向とに基づいて一意に定まる音波の到達時間差のことである。例えば、音源探索対象方向は、所定の角度刻みで設定される。すなわち、相対遅延時間は、音源探索対象方向の分だけ算出される。相対遅延時間計算部１３は、算出した音源探索対象方向と相対遅延時間とをセットにして周波数別推定方向情報生成部１５に出力する。
　相対遅延時間は、マイクペアの位置関係によって算出方法が異なる。以下において、マイクペアの位置関係を二つ挙げ、それらのマイクペアの位置関係ごとに相対遅延時間の算出方法を示す。
　図２は、全てのマイクが同一直線上に配置される例である。図２の例では、マイクが３つの場合について説明する。ここで、音速をｃ、マイク間隔をｄr、音源探索対象方向（音源方向とも呼ぶ）をθとする。音源方向θは、音源１００の方向を推定するために設定される少なくとも一つの角である。このとき、音源方向θに対する相対遅延時間τr（θ）は、以下の式３を用いて計算できる。
Figure JPOXMLDOC01-appb-I000003
　マイク間隔ｄは、信号選択部１２によって選択される入力信号の組み合わせにより異なる。そのため、相対遅延時間τr(θ)は、組み合わせ番号ｒごとに異なる。例えば、図２のマイクペアＡＢ間の距離をｄ1とすると、相対遅延時間τ1（θ）は、以下の式４を用いて計算できる。
Figure JPOXMLDOC01-appb-I000004
　また、図２のマイクペアＡＣ間の距離をｄ2とすると、相対遅延時間τ2（θ）は、以下の式５を用いて計算できる。
Figure JPOXMLDOC01-appb-I000005
　このように、全てのマイクが同一直線上に位置する場合、ある音源に関する相対遅延時間τr(θ)はマイク間隔ｄに比例するが、音源方向θはどのマイクから見ても同一とみなせる。
　図３は、二組のマイクペアが互いに直交する直線上に配置される例である。図３の例では、音源方向θはマイクペアによって異なる。図３のマイクペアＡＢ間の相対遅延時間τ1（θ）は、以下の式６を用いて計算できる。
Figure JPOXMLDOC01-appb-I000006
　一方、図３のマイクＣＤ間の相対遅延時間τ2（θ）は、以下の式７で計算できる。
Figure JPOXMLDOC01-appb-I000007
　このように、あるマイクペアを基準として、別のマイクペアの相対遅延時間τｒ（θ）は、基準のマイクペアから見た音源方向θの関数として、以下の式８のように一般化できる。なお、基準とするマイクペアはどれを選んでもかまわない。
Figure JPOXMLDOC01-appb-I000008
　相対遅延時間計算部１３は、全ての音源探索対象方向に対して相対遅延時間を計算する。例えば、相対遅延時間計算部１３は、音源方向の探索範囲が１０度刻みで０度から９０度まで、つまり０度、１０度、２０度、・・・、９０度である場合、１０種類の相対遅延時間を計算する。そして、相対遅延時間計算部１３は、音源探索対象方向と相対遅延時間とを周波数別推定方向情報生成部１５に出力する。
　〔周波数別推定方向情報生成部〕
　周波数別推定方向情報生成部１５－１～Ｒには、信号選択部１２によって選択された全マイクペアのうち１組のマイクペアの入力信号と、相対遅延時間計算部１３から供給された相対遅延時間とが入力される。周波数別推定方向情報生成部１５－１～Ｒは、入力されたマイクペアの入力信号と相対遅延時間とを用いて、そのマイクペアの入力信号間の周波数別推定方向情報を生成する。
　ここで、図４を参照しながら、周波数別推定方向情報生成部１５の詳細な構成について説明する。図４は、周波数別推定方向情報生成部１５のブロック図である。周波数別推定方向情報生成部１５は、変換部１５１、クロススペクトル計算部１５２、平均計算部１５３、分散計算部１５４、周波数別クロススペクトル生成部１５５、逆変換部１５６、および周波数別推定方向情報計算部１５７を備える。
　〔変換部〕
　変換部１５１には、信号選択部１２から２つの入力信号（入力信号Ａおよび入力信号Ｂ）が入力される。変換部１５１は、信号選択部１２から供給された２つの入力信号を変換信号（周波数領域信号とも呼ぶ）に変換する。変換部１５１は、入力信号を複数の周波数成分に分解するための変換を行う。例えば、変換部１５１は、フーリエ変換を用いて、入力信号を複数の周波数成分に分解する。変換部１５１は、クロススペクトル計算部１５２に変換信号を出力する。
　変換部１５１には２種類の入力信号ｘm（ｔ）が入力される。ここで、ｍは、入力端子１１の番号である。変換部１５１は、入力端子１１から供給された入力信号から、適当な長さの波形を一定の周期でずらしながら切り出す。このように切り出された信号区間をフレーム、切り出された波形の長さをフレーム長、フレームをずらす周期をフレーム周期と呼ぶ。そして、変換部１５１は、フーリエ変換を用いて、切り出された信号を周波数領域信号に変換する。ここで、ｎをフレーム番号とし、切り出される入力信号をｘm（ｔ，ｎ）とする（ｔ＝０、１、・・・、Ｋ―１）。このとき、入力信号ｘｍ（ｔ,ｎ）のフーリエ変換Ｘm（ｋ,ｎ）は、以下の式９を用いて計算できる。
Figure JPOXMLDOC01-appb-I000009
　なお、上記の式９において、ｊは虚数単位、ｅｘｐは指数関数を表す。また、ｋは、周波数ビン番号を表し、０以上Ｋ－１以下の整数である。以下、簡略化のため、ｋのことを周波数ビン番号ではなく、単に周波数とも呼ぶ。
　〔クロススペクトル計算部〕
　クロススペクトル計算部１５２には、変換部１５１から変換信号が入力される。クロススペクトル計算部１５２は、変換部１５１から供給される変換信号を用いてクロススペクトルを計算する。クロススペクトル計算部１５２は、算出したクロススペクトルを平均計算部１５３に出力する。
　クロススペクトル計算部１５２は、変換信号Ｘ2（ｋ，ｎ）の複素共役と変換信号Ｘ1（ｋ,ｎ）との積を計算してクロススペクトルを計算する。ここで、変換信号のクロススペクトルをＳ12（ｋ,ｎ）とする。このとき、クロススペクトル計算部１５２は、以下の式１０を用いてクロススペクトルを計算する。
Figure JPOXMLDOC01-appb-I000010
　なお、ｃｏｎｊ（Ｘ2（ｋ，ｎ））は、Ｘ2（ｋ,ｎ）の複素共役を表す。また、式１０ではなく、振幅成分で正規化したクロススペクトルを用いてもよい。クロススペクトル計算部１５２は、振幅成分で正規化する場合、以下の式１１を用いてクロススペクトルを計算する。
Figure JPOXMLDOC01-appb-I000011
　〔平均計算部〕
　平均計算部１５３には、クロススペクトル計算部１５２からクロススペクトルが入力される。平均計算部１５３は、クロススペクトル計算部１５２から供給されたクロススペクトルの平均（平均クロススペクトルとも呼ぶ）を計算する。平均計算部１５３は、算出した平均クロススペクトルを分散計算部１５４および周波数別クロススペクトル生成部１５５に出力する。
　ここで、平均計算部１５３が過去に入力されたクロススペクトルから周波数ビンごとに平均クロススペクトルを計算する例について説明する。なお、平均計算部１５３は、周波数ビン単位ではなく、複数の周波数ビンを束ねたサブバンド単位で平均クロススペクトルを計算してもよい。ここで、第ｎフレームの周波数ビンｋにおけるクロススペクトルをＳ12（ｋ，ｎ）とする、このとき、平均計算部１５３は、以下の式１２を用いて、過去Ｌフレームから求めた平均クロススペクトルＳＳ12（ｋ，ｎ）を計算する。
Figure JPOXMLDOC01-appb-I000012
　また、平均計算部１５３は、次のリーク積分を用いて平均クロススペクトルＳＳ12（ｋ，ｎ）を計算してもよい。ただし、αは、０よりも大きく、１よりも小さい実数である。
Figure JPOXMLDOC01-appb-I000013
　〔分散計算部〕
　分散計算部１５４には、平均計算部１５３から平均クロススペクトルが入力される。分散計算部１５４は、平均計算部１５３から供給された平均クロススペクトルを用いて分散を計算する。分散計算部１５４は、算出した分散を周波数別クロススペクトル生成部１５５に出力する。
　ここで、平均クロススペクトルをＳＳ12（ｋ，ｎ）とする。このとき、分散計算部１５４は、クロススペクトルの位相の分散の計算において円周分散を用いた場合、以下の式１４を用いて、分散Ｖ12（ｋ，ｎ）を計算する。
Figure JPOXMLDOC01-appb-I000014
　なお、分散計算部１５４は、以下の式１５を用いて分散Ｖ12（ｋ，ｎ）を計算してもよい。
Figure JPOXMLDOC01-appb-I000015
　また、分散計算部１５４は、円周標準偏差を用いる場合、以下の式１６を用いて、分散Ｖ12（ｋ，ｎ）を計算する。
Figure JPOXMLDOC01-appb-I000016
　〔周波数別クロススペクトル生成部〕
　ここで、周波数別クロススペクトル生成部１５５の構成について図面を参照しながら説明する。図５は、周波数別クロススペクトル生成部１５５の構成の一例を示すブロック図である。図５のように、周波数別クロススペクトル生成部１５５は、周波数別基本クロススペクトル計算部５５１、カーネル関数スペクトル生成部５５２、および乗算部５５３を有する。
　〔周波数別基本クロススペクトル計算部〕
　周波数別基本クロススペクトル計算部５５１には、平均計算部１５３から平均クロススペクトルが入力される。周波数別基本クロススペクトル計算部５５１は、平均計算部１５３から供給された平均クロススペクトルを用いて周波数別基本クロススペクトルを計算する。周波数別基本クロススペクトル計算部５５１は、算出した周波数別基本クロススペクトルを乗算部５５３に出力する。
　周波数別基本クロススペクトル計算部５５１は、平均計算部１５３から供給される平均クロススペクトルＳＳ12（ｋ，ｎ）を用いて、平均クロススペクトルＳＳ12（ｋ，ｎ）の各周波数ｋに対応するクロススペクトル（周波数別基本クロススペクトルとも呼ぶ）を計算する。周波数別基本クロススペクトル計算部５５１は、算出する周波数別基本クロススペクトルを乗算部５５３に出力する。周波数別基本クロススペクトルは、周波数成分ごとに相関関数を計算するために算出される。周波数別基本クロススペクトル計算部５５１は、ある周波数ｋに対応する相関関数（周波数別相関関数とも呼ぶ）を後段で求めるための周波数別基本クロススペクトルを計算する。
　ここで、周波数別基本クロススペクトル計算部５５１が周波数ｋの周波数別基本クロススペクトルを計算する例について詳細に説明する。周波数別基本クロススペクトル計算部５５１は、周波数ｋの平均クロススペクトルＳＳ12（ｋ，ｎ）を用いて周波数別基本クロススペクトルを計算する際、位相成分と振幅成分とを予め別々に求めたのちに統合する。周波数ｋの周波数別基本クロススペクトルＵk（ｗ，ｎ）、その振幅成分を｜Ｕk（ｗ，ｎ）｜、位相成分をａｒｇ（Ｕk（ｗ，ｎ））とすると、以下の式１７の関係が成り立つ。ただし、式１７において、ｗは、周波数を表し、０以上Ｗ－１以下の整数である。
Figure JPOXMLDOC01-appb-I000017
　以下において、周波数別基本クロススペクトル計算部５５１が、周波数ｋの平均クロススペクトルＳＳ12（ｋ，ｎ）を用いて、周波数別基本クロススペクトルの振幅成分｜Ｕk（ｗ，ｎ）｜と位相成分ａｒｇ（Ｕk（ｗ,ｎ））とを求める方法について説明する。
　ｋの整数倍の周波数の振幅成分｜Ｕk（ｗ，ｎ）｜には１．０が用いられる。一方、ｋの非整数倍の周波数の位相成分はゼロにする。これを数式で表現すると、振幅成分｜Ｕk（ｗ，ｎ）｜は以下の式１８で与えられる。ただし、式１８において、pは1以上Ｐ以下の整数である。
Figure JPOXMLDOC01-appb-I000018
　波源方向推定を行うときに重要な情報は位相成分であるため、振幅成分には、式１８のように適当な定数を用いる。なお、ｋの整数倍の周波数の振幅成分｜Ｕk（ｗ，ｎ）｜として、１．０の代わりに｜ＳＳ12（ｋ，ｎ）｜を用いてもよい。つまり、振幅成分｜Ｕk（ｗ，ｎ）｜を以下の式１９を用いて求めてもよい。
Figure JPOXMLDOC01-appb-I000019
　ｋを整数倍した周波数の位相成分ａｒｇ（Ｕk（ｗ，ｎ））には、周波数ｋの平均クロススペクトルＳＳ12（ｋ，ｎ）を定数倍したものを用いる。例えば、周波数ｋ、２ｋ、３ｋ、および４ｋの位相成分には、周波数ｋの位相成分ａｒｇ（ＳＳ12（ｋ，ｎ））のそれぞれを同一の倍率で整数倍したものが用いられる。すなわち、周波数ｋ、２ｋ、３ｋ、および４ｋの位相成分のそれぞれには、ａｒｇ（ＳＳ12（ｋ，ｎ））、２ａｒｇ（ＳＳ1２（ｋ，ｎ））、３ａｒｇ（ＳＳ12（ｋ，ｎ））、４ａｒｇ（ＳＳ12（ｋ，ｎ））が用いられる。一方、ｋの非整数倍の周波数の位相成分はゼロにする。したがって、周波数ｋに対応する周波数別基本クロススペクトルの位相成分ａｒｇ（Ｕk（ｗ，ｎ））は、以下の式２０を用いて計算される。なお、pは１以上Ｐ以下の整数である（Ｐ＞１）。
Figure JPOXMLDOC01-appb-I000020
　周波数別基本クロススペクトル計算部５５１は、式１８や式１９を用いて算出される振幅成分と、式２０を用いて算出される位相成分とを、式１７を用いて統合し、周波数ｋの周波数別基本クロススペクトルＵk（ｗ，ｎ）を得る。
　これまで説明してきた方法では、振幅成分と位相成分を別々に求めてから周波数別基本クロススペクトルを算出した。しかし、以下の式２１に示すクロススペクトルのべき乗を用いれば、振幅成分と位相成分を求めることなく、周波数別基本クロススペクトルＵk（ｗ，ｎ）を求めることができる。
Figure JPOXMLDOC01-appb-I000021
　〔カーネル関数スペクトル生成部〕
　カーネル関数スペクトル生成部５５２には、分散計算部１５４から分散が入力される。カーネル関数スペクトル生成部５５２は、分散計算部１５４から供給された分散を用いてカーネル関数スペクトルを算出する。カーネル関数スペクトルとは、カーネル関数をフーリエ変換し、その絶対値を取ったものである。なお、カーネル関数スペクトルには、カーネル関数をフーリエ変換し、その絶対値を取る代わりに二乗してもよい。また、カーネル関数スペクトルには、カーネル関数をフーリエ変換し、その絶対値の二乗としてもよい。カーネル関数スペクトル生成部５５２は、算出したカーネル関数スペクトルを乗算部５５３に出力する。
　ここで、カーネル関数スペクトルをＧ（ｗ）とし、カーネル関数をｇ（τ）とする。また、カーネル関数としてガウス関数を用いる。このとき、ガウス関数は、以下の式２２で与えられる。
Figure JPOXMLDOC01-appb-I000022
　ただし、式２２において、ｇ1、ｇ2、およびｇ3は正の実数である。ｇ1はガウス関数の大きさを制御し、ｇ2はガウス関数のピークの位置を制御し、ｇ3はガウス関数の広がりを制御する。特に、ガウス関数の広がりを調整するｇ3は、周波数別相関関数のピークの鋭さに大きな影響を与えるので重要である。すなわち、式２２は、ｇ3が大きくなればガウス関数の広がりが大きくなることを示す。
　また、以下の式２３のロジスティック分布の確率密度関数をカーネル関数として用いてもよい。ただし、式２３において、ｇ4とｇ5は正の実数である。
Figure JPOXMLDOC01-appb-I000023
　ロジスティック分布の確率密度関数は、ガウス関数と同様の形状をしているが、ガウス関数よりも裾が長い。特に、ロジスティック分布の確率密度関数の広がりを調整するｇ5は、式２２のガウス関数におけるｇ3の場合と同様に、周波数別相関関数のピークの鋭さに大きな影響を与えるパラメータである。また、コサイン関数や一様関数をカーネル関数に用いてもよい。
　カーネル関数のパラメータのうち、カーネル関数の広がりに影響を与えるｇ3やｇ5は、分散計算部１５４から入力される分散を用いて決定される。ここで、これらのパラメータのことを広がり制御パラメータと呼び、ｑ（ｋ，ｎ）と表現する。したがって、カーネル関数がガウス関数の場合には、ｇ3はｑ（ｋ，ｎ）である。分散が小さければ、周波数別相関関数のピークが鋭く、裾が狭くなるようにパラメータを変化させる。したがって、広がり制御パラメータを小さくする。
　広がり制御パラメータは、予め設定した写像関数を用いて分散の値を変換することによって算出できる。例えば、分散がある閾値を上回れば広がり制御パラメータを大きな値（例えば１０）にし、分散がある閾値を下回れば小さな値（例えば０．０１）に設定する。ここで、分散をＶ12（ｋ，ｎ）とし、閾値をｐthとする。このとき、第ｎフレームの周波数ビンｋにおける広がり制御パラメータｑ（ｋ，ｎ）は、以下の式２４を用いて計算できる。ただし、式２４において、ｑ1およびｑ2は、ｑ1＞ｑ2を満たす正の実数である。
Figure JPOXMLDOC01-appb-I000024
　また、広がり制御パラメータｑ（ｋ，ｎ）は、以下の式２５のように一次関数を用いて計算してもよい。ただし、式２５において、ｑ3は０よりも大きな実数であり、ｑ4は実数である。
Figure JPOXMLDOC01-appb-I000025
　ｑ3、ｑ4としては、例えば、および式２７で示す値を用いてもよい。
Figure JPOXMLDOC01-appb-I000026

Figure JPOXMLDOC01-appb-I000027
Ｌは平均計算部１５３が平均クロススペクトルを求める際に平均化したフレーム数を表す。平均クロススペクトルの誤差は平均化フレーム数Ｌに反比例するため、式２６および式２７を用いることで、平均クロススペクトルの誤差（信頼性）を考慮して広がり制御パラメータを求めることができる。
　また、線形写像関数や高次の多項式関数、非線形関数などで表される分散の関数を分散の計算に用いることも可能である。また、分散をそのまま広がり制御パラメータにしてもよい。
　広がり制御パラメータを求める関数は、分散だけでなく、周波数ｋの関数としてもよい。例えば、周波数ｋが大きくなるにつれて小さくなる関数を用いることができる。このような代表例としては、ｋの逆数を用いる例が挙げられる。この場合、式２５の代わりに、以下の式２８の関数を用いて広がり制御パラメータｑ（ｋ，ｎ）を計算できる。
Figure JPOXMLDOC01-appb-I000028
　また、式２６の代わりとしては、以下の式２９の関数を用いて広がり制御パラメータｑ（ｋ，ｎ）を計算できる。
Figure JPOXMLDOC01-appb-I000029
　〔乗算部〕
　乗算部５５３には、周波数別基本クロススペクトル計算部５５１から周波数別基本クロススペクトルが入力され、カーネル関数スペクトル生成部５５２からカーネル関数スペクトルが入力される。乗算部５５３は、周波数別基本クロススペクトル計算部５５１から供給される周波数別基本クロススペクトルと、カーネル関数スペクトル生成部５５２から供給されるカーネル関数スペクトルとの積を計算し、周波数別クロススペクトルを算出する。乗算部５５３は、算出した周波数別クロススペクトルを逆変換部１５６に出力する。
　ここで、周波数別基本クロススペクトル計算部５５１から供給される周波数別基本クロススペクトルをＵk（ｗ，ｎ）、カーネル関数スペクトル生成部５５２から供給されるカーネル関数スペクトルをＧ（ｗ）とする。このとき、乗算部５５３は、以下の式３０を用いて、周波数別クロススペクトルＵＭk（ｗ，ｎ）を計算する。
Figure JPOXMLDOC01-appb-I000030
　〔逆変換部〕
　逆変換部１５６には、周波数別クロススペクトル生成部１５５の乗算部５５３から周波数別クロススペクトルが入力される。例えば、変換部１５１がフーリエ変換を用いる場合、逆変換部１５６は逆フーリエ変換を用いて逆変換を行う。逆変換部１５６は、周波数別クロススペクトル生成部１５５から供給された周波数別クロススペクトルの逆変換を求める。
　ここで、周波数別クロススペクトル生成部１５５から供給される周波数別クロススペクトルをＵＭk（ｗ，ｎ）とする。このとき、逆変換部１５６は、以下の式３１を用いて、ＵＭk（ｗ，ｎ）を逆変換して周波数別相互相関関数ｕk（τ，ｎ）を計算する。
Figure JPOXMLDOC01-appb-I000031
　〔周波数別推定方向情報計算部〕
　周波数別推定方向情報計算部１５７には、逆変換部１５６から周波数別相互相関関数が入力され、相対遅延時間計算部１３から相対遅延時間が入力される。周波数別推定方向情報計算部１５７は、逆変換部１５６から供給された周波数別相互相関関数と、相対遅延時間計算部１３から供給された相対遅延時間とを用いて、方向と相関値との対応関係を周波数別推定方向情報として求める。周波数別推定方向情報計算部１５７は、求めた周波数別推定方向情報を統合部１７に出力する。
　ここで、周波数別相互相関関数をｕk（τ，ｎ）とし、相対遅延時間をτr（θ）とする。このとき、周波数別推定方向情報計算部１５７は、以下の式３２を用いて、周波数別推定方向情報Ｈk,r（θ，ｎ）を算出する。
Figure JPOXMLDOC01-appb-I000032
　式３２を用いれば、各方向θに対して相関値が定まるので、相関値が高い方向に音源が存在する可能性が高いと判断できる。
　〔統合部〕
　統合部１７には、周波数別推定方向情報生成部１５－１～Ｒから周波数別推定方向情報が入力される。統合部１７は、周波数別推定方向情報生成部１５－１～Ｒから供給される周波数別推定方向情報を統合して統合推定方向情報を算出する。統合部１７は、個別に求めた複数の周波数別推定方向情報を混合したり、重ね合わせたりすることによって一つの推定方向情報を求める。統合部１７は、算出した統合推定方向情報を出力する。例えば、統合部１７は、上位システム（図示しない）に統合推定方向情報を出力する。
　例えば、統合部１７は、まず、入力信号の組み合わせ数（Ｒ個）分の周波数別推定方向情報Ｈk,r（θ，ｎ）を統合することによって周波数別統合推定方向情報Ｈk（θ，ｎ）を算出する。そして、統合部１７は、算出した周波数別統合推定方向情報を全周波数について統合することによって統合推定方向情報Ｈ（θ，ｎ）を算出する。
　例えば、統合部１７は、周波数別推定方向情報Ｈk,r（θ，ｎ）の総乗を計算することによって周波数別統合推定方向情報Ｈk（θ，ｎ）を算出する。このとき、統合部１７は、以下の式３３を用いて周波数別統合推定方向情報Ｈk（θ，ｎ）を算出する。
Figure JPOXMLDOC01-appb-I000033
　また、例えば、統合部１７は、周波数別推定方向情報Ｈk,r（θ，ｎ）の総和を計算することによって周波数別統合推定方向情報Ｈk（θ，ｎ）を算出してもよい。このとき、統合部１７は、以下の式３４を用いて周波数別統合推定方向情報Ｈk（θ，ｎ）を算出する。
Figure JPOXMLDOC01-appb-I000034
　統合部１７は、統合推定方向情報Ｈ（θ，ｎ）の算出においては、周波数別統合推定方向情報Ｈk（θ，ｎ）の周波数ｋについての総和または総乗を計算する。
　例えば、統合部１７は、以下の式３５を用いて、周波数別統合推定方向情報Ｈk（θ，ｎ）の周波数ｋについての総和を統合推定方向情報Ｈ（θ，ｎ）として計算する。
Figure JPOXMLDOC01-appb-I000035
　また、例えば、統合部１７は、以下の式３６を用いて、周波数別統合推定方向情報Ｈk（θ，ｎ）の周波数ｋについての総乗を統合推定方向情報Ｈ（θ，ｎ）として計算する。
Figure JPOXMLDOC01-appb-I000036
　目標音の存在する周波数や、目的音のパワーが大きい周波数が予め判明している場合、統合部１７は、その周波数に対応する周波数別統合推定方向情報だけを用いて、統合推定方向情報を求めてもよい。また、統合部１７は、重み付けという形で、統合における周波数別統合推定方向情報の影響度を制御してもよい。例えば、目標音の存在する周波数の集合をΩとすると、統合部１７は、以下の式３７を用いて、周波数の選択により統合推定方向情報Ｈ（θ，ｎ）を求めることができる。
Figure JPOXMLDOC01-appb-I000037
　また、重み付けを用いる場合には、統合部１７は、以下の式３８を用いて統合推定方向情報Ｈ（θ，ｎ）を算出できる。ただし、式３８において、ａおよびｂは、ａ＞ｂ＞０を満たす実数である。
Figure JPOXMLDOC01-appb-I000038
　このように、目的音が存在する周波数の周波数別統合推定方向情報を重点的に用いて統合すると、雑音などの非目的音の影響が小さい相関関数を生成できるため、方向推定精度が向上する。
　また、統合部１７は、別の算出方法を用いて統合推定方向情報Ｈ（θ，ｎ）を算出してもよい。例えば、統合部１７は、まず、全周波数について周波数別推定方向情報Ｈk,r（θ，ｎ）を統合した入力信号組み合わせ別統合推定方向情報Ｈr（θ，ｎ）を算出する。そして、統合部１７は、入力信号の全組み合わせについて入力信号組み合わせ別統合推定方向情報を統合した統合推定方向情報Ｈ（θ，ｎ）を算出してもよい。
　以上が、本実施形態の波源方向推定装置１０の構成についての説明である。
　なお、図６のように、波源方向推定装置１０にマイクなどのセンサ１１０を少なくとも一つ追加した構成も本実施形態の範囲に含まれる。それぞれのセンサ１１０は、インターネットやイントラネットなどネットワークやケーブルを介して、波源方向推定装置１０のいずれかの入力端子１１に接続される。
　例えば、センサ１１０は、音波を検出する場合、マイクロフォンによって実現される。例えば、センサ１１０は、振動波を検出する場合、振動センサによって実現される。例えば、センサ１１０は、電磁波を検出する場合、アンテナによって実現される。なお、センサ１１０は、検知対象の波動を電気信号に変換できるものであれば、その形態に限定を加えない。
　（動作）
　次に、本実施形態の波源方向推定装置１０の動作について図面を参照しながら説明する。
　〔波源方向推定〕
　まず、図７のフローチャートを用いて、波源方向推定装置１０の動作の概略について説明する。なお、図７のフローチャートに沿った説明においては、波源方向推定装置１０を動作の主体として説明する。
　図７において、まず、波源方向推定装置１０は、複数のマイクから電気信号（入力信号とも呼ぶ）を入力する（ステップＳ１１１）。
　次に、波源方向推定装置１０は、複数のマイクに対応する入力信号から２つの入力信号を選択する（ステップＳ１１２）。
　次に、波源方向推定装置１０は、選択した２つの入力信号の供給元である２つのマイクの間隔（マイク間隔とも呼ぶ）と、設定された音源探索対象方向とに基づいて相対遅延時間を算出する（ステップＳ１１３）。
　次に、波源方向推定装置１０は、選択された２つの入力信号と相対遅延時間とを用いて、周波数ごとに推定方向情報（周波数別推定方向情報とも呼ぶ）を生成する（ステップＳ１１４）。
　次に、波源方向推定装置１０は、周波数ごとの推定方向情報を統合して統合推定方向情報を算出する（ステップＳ１１５）。
　そして、波源方向推定装置１０は、統合推定方向情報を出力する（ステップＳ１１６）。
　以上が、波源方向推定装置１０の動作の概略についての説明である。
　〔周波数別推定方向情報生成〕
　次に、図８のフローチャートを用いて、波源方向推定装置１０の周波数別推定方向情報生成部１５の動作について説明する。図８のフローチャートの処理は、図７のフローチャートのステップＳ１１４を細分化したものである。なお、図８のフローチャートに沿った説明においては、周波数別推定方向情報生成部１５を動作の主体として説明する。
　図８において、まず、周波数別推定方向情報生成部１５は、信号選択部１２によって選択された２つの入力信号と、それらの入力信号の相対遅延時間とを入力する（ステップＳ１２１）。
　次に、周波数別推定方向情報生成部１５は、２つの入力信号を周波数領域信号（変換信号とも呼ぶ）に変換する（ステップＳ１２２）。
　次に、周波数別推定方向情報生成部１５は、変換信号を用いてクロススペクトルを算出する（ステップＳ１２３）。
　次に、周波数別推定方向情報生成部１５は、クロススペクトルを用いて平均クロススペクトルを算出する（ステップＳ１２４）。
　次に、周波数別推定方向情報生成部１５は、平均クロススペクトルを用いて分散を算出する（ステップＳ１２５）。
　次に、周波数別推定方向情報生成部１５は、平均クロススペクトルおよび分散を用いて周波数別クロススペクトルを算出する（ステップＳ１２６）。
　次に、周波数別推定方向情報生成部１５は、周波数別クロススペクトルを用いて周波数別相互相関関数を算出する（ステップＳ１２７）。
　次に、周波数別推定方向情報生成部１５は、周波数別相互相関関数と相対遅延時間とを用いて周波数別推定方向情報を算出する（ステップＳ１２８）。
　そして、周波数別推定方向情報生成部１５は、周波数別推定方向情報を統合部１７に出力する（ステップＳ１２９）。
　以上が、周波数別推定方向情報生成部１５の動作についての説明である。
　〔周波数別クロススペクトル生成〕
　次に、図９のフローチャートを用いて、波源方向推定装置１０の周波数別推定方向情報生成部１５に含まれる周波数別クロススペクトル生成部１５５の動作について説明する。図９のフローチャートの処理は、図８のフローチャートのステップＳ１２５を細分化したものである。なお、図９のフローチャートに沿った説明においては、周波数別クロススペクトル生成部１５５を動作の主体として説明する。
　図９において、まず、周波数別クロススペクトル生成部１５５は、平均計算部１５３から平均クロススペクトルを入力し、分散計算部１５４から分散を入力する（ステップＳ１３１）。
　次に、周波数別クロススペクトル生成部１５５は、平均クロススペクトルを用いて周波数別基本クロススペクトルを算出する（ステップＳ１３２）。
　また、周波数別クロススペクトル生成部１５５は、分散を用いてカーネル関数スペクトルを算出する（ステップＳ１３３）。なお、ステップＳ１３２の処理と、ステップＳ１３３の処理とは、並行して行ってもよいし、逐次的に行ってもよい。
　次に、周波数別クロススペクトル生成部１５５は、周波数別基本クロススペクトルとカーネル関数スペクトルとの積を計算し、周波数別クロススペクトルを算出する（ステップＳ１３４）。
　そして、周波数別クロススペクトル生成部１５５は、算出した周波数別クロススペクトルを逆変換部１５６に出力する（ステップＳ１３５）。
　以上が、周波数別クロススペクトル生成部１５５の動作についての説明である。
　以上のように、本実施形態の波源方向推定装置は、複数の入力部、信号選択部、相対遅延時間計算部、少なくとも一つの周波数別推定方向情報生成部、および統合部を備える。複数の入力部は、複数のセンサによって取得された波動から変換される電気信号を入力信号として取得する。信号選択部は、複数の入力信号から少なくとも二つの入力信号を組み合わせたペアを少なくとも二つ選択する。相対遅延時間計算部は、入力信号のペアを構成する少なくとも二つの入力信号間において波源探索方向ごとの波動の到達時間差を相対遅延時間として算出する。少なくとも一つの周波数別推定方向情報生成部は、入力信号のペアと相対遅延時間とを用いて波動の波源の推定方向情報を周波数ごとに生成する。統合部は、周波数別推定方向情報生成部によって生成される周波数ごとの推定方向情報を統合する。
　例えば、信号選択部は、センサの間隔に基づいて、複数の入力信号から少なくとも二つの入力信号を組み合わせたペアを選択する。
　例えば、相対遅延時間計算部は、一つの入力信号のペアの供給元であるセンサのペアに対する波源探索方向を基準として、信号選択部によって選択された全ての入力信号のペアの相対遅延時間を、基準である波源探索方向の関数として計算する。
　例えば、周波数別推定方向情報生成部は、変換部、クロススペクトル計算部、平均計算部、分散計算部、周波数別クロススペクトル生成部、逆変換部、および推定方向情報計算部を有する。変換部は、ペアをなす少なくとも２つの入力信号を周波数領域の変換信号に変換する。クロススペクトル計算部は、変換手段によって変換された変換信号を用いてクロススペクトルを計算する。平均計算部は、クロススペクトル計算部によって算出されるクロススペクトルを用いて平均クロススペクトルを計算する。分散計算部は、平均計算部によって算出される平均クロススペクトルを用いて分散を計算する。周波数別クロススペクトル生成部は、平均計算部によって算出される平均クロススペクトルと、分散計算部によって算出される分散とを用いて周波数別クロススペクトルを計算する。逆変換部は、周波数別クロススペクトル生成部によって算出される周波数別クロススペクトルを逆変換して周波数別相互相関関数を計算する。推定方向情報計算部は、逆変換部によって算出される周波数別相互相関関数と相対遅延時間とを用いて周波数別推定周波数ごとに推定方向情報を計算する。
　例えば、周波数別クロススペクトル生成部は、周波数別基本クロススペクトル計算部、カーネル関数スペクトル生成部、および乗算部を含む。周波数別基本クロススペクトル計算部は、平均計算部から平均クロススペクトルを取得し、取得した平均クロススペクトルを用いて周波数別基本クロススペクトルを計算する。カーネル関数スペクトル生成部は、分散計算部から分散を取得し、取得した分散を用いてカーネル関数スペクトルを計算する。乗算部は、周波数別基本クロススペクトル計算部によって算出される周波数別基本クロススペクトルと、カーネル関数スペクトル生成部によって算出されるカーネル関数スペクトルとの積を計算して周波数別クロススペクトルを計算する。
　例えば、統合部は、複数の周波数ごとの推定方向情報を複数の入力信号のペアについて統合した周波数別統合推定方向情報を算出する。そして、統合部は、算出した周波数別統合推定方向情報を全ての周波数について統合して統合推定方向情報を算出する。
　例えば、統合部は、複数の周波数ごとの推定方向情報を全ての周波数について統合した入力信号組み合わせ別統合推定方向情報を算出する。統合部は、算出した入力信号組み合わせ別統合推定方向情報を全ての入力信号の組み合わせについて統合して統合推定方向情報を算出する。
　例えば、波源方向推定装置は、複数の入力部のそれぞれに対応して配置されるセンサを備える。
　本実施形態の波源方向推定装置は、マイクペア間の相互相関関数から推定方向情報を求め、複数のマイクペア間で推定方向情報を統合する。その結果、本実施形態の波源方向推定装置によれば、マイクペア間で位相が偶然揃うことに起因して発生する、音源方向以外の方向における推定方向情報の偽ピークが小さくなり、虚像音源の誤推定の発生を低減し、高精度に音源の方向を推定することができる。
　本実施形態の波源方向推定装置の推定対象は、空気や水の振動波である音波の発生源に限定されない。本実施形態の波源方向推定装置は、地震や地滑りなどのように固体を媒質とする振動波の発生源の方向推定にも適用できる。その場合、振動波を電気信号に変換する装置には、マイクロフォンではなく振動センサを用いることができる。さらに、本実施形態の波源方向推定装置は、気体・液体・固体の振動波だけでなく、電波を用いて方向を推定する場合にも適用できる。電波を用いて方向推定する場合、電波を電気信号に変換する装置としてアンテナを用いることができる。
　本実施形態の波源方向推定装置によって推定される統合推定方向情報は、様々な形で利用できる。例えば、統合推定方向情報が複数のピークを有する場合には、各ピークを到来方向とする複数の音源が存在すると推定される。したがって、統合推定方向情報を用いれば、各音源の方向を同時に推定できるだけでなく、音源の数を推定できる。
　（第２の実施形態）
　次に、本発明の第２の実施形態に係る波源方向推定装置について図面を参照しながら説明する。本実施形態の波源方向推定装置は、第１の実施形態の波源方向推定装置に波源方向算出部を追加した構成を有する。
　図１０は、本実施形態の波源方向推定装置２０の構成を表すブロック図である。波源方向推定装置２０は、入力端子２１、信号選択部２２、相対遅延時間計算部２３、周波数別推定方向情報生成部２５、統合部２７、および波源方向算出部２８を備える。なお、入力端子２１、信号選択部２２、相対遅延時間計算部２３、周波数別推定方向情報生成部２５、および統合部２７は、第１の実施形態の波源方向推定装置１０の対応する構成と同様であるため、詳細な説明を省略する。
　〔波源方向算出部〕
　波源方向算出部２８には、統合部２７から統合推定方向情報が入力される。波源方向算出部２８は、統合推定方向情報を用いて、波源方向を計算する。波源方向算出部２８は、計算した波源方向を出力する。
　波源方向算出部２８における波源方向の計算方法を以下で詳細に説明する。統合部２７から入力された統合推定方向情報は、ピークが大きいほど信頼性（音源の存在可能性）が高い。そのため、例えば音源数が１つと事前に仮定できる場合には、波源方向算出部２８は、統合推定方向情報が最大となる方向を推定方向として出力する。ここで、統合部２７から入力された統合推定方向情報をＨ（θ，ｎ）とする。波源方向算出部２８は、以下の式３９を用いて、統合推定方向情報Ｈ（θ，ｎ）が最大値を取るような統合推定方向情報Ｈ（θ，ｎ）の引数を要素とする集合を波源方向Θとして算出できる。なお、式３９において、θは全ての波源方向または波源方向候補を表す。
Figure JPOXMLDOC01-appb-I000039
　また、統合推定方向情報のピークが閾値を超える場合には、波源方向算出部２８は、閾値を超えるピークを有する方向を音源とみなし、閾値を超える方向を推定方向として出力することもできる。
　また、本実施形態の波源方向推定装置は、一定時間Ｔごとに統合推定方向情報が最大になる時刻に対応する方向を音源方向として推定することもできる。ただし、音源の方向は、一定時間Ｔの間で変化しない、あるいは変化の大きさが無視できるほど小さいと仮定する。このように仮定することによって、波源方向の推定精度を向上させることができる。
　以上のように、本実施形態の波源方向推定装置は、統合手段によって算出される統合推定方向情報に基づいて波動の波源方向を算出する波源方向算出手段を備える。例えば、波源方向算出手段は、一定時間ごとに統合推定方向情報が最大になる時刻に対応する方向を波源方向として算出する。本実施形態の波源方向推定装置によれば、虚像音源を誤推定することなく、高精度に音源の方向を推定することができる。
　（ハードウェア）
　ここで、各実施形態に係る波源方向推定装置の処理を実行するハードウェア構成について、図１１の情報処理装置９０を一例として挙げて説明する。なお、図１１の情報処理装置９０は、各実施形態の波源方向推定装置の処理を実行するための構成例であって、本発明の範囲を限定するものではない。
　図１１のように、情報処理装置９０は、プロセッサ９１、主記憶装置９２、補助記憶装置９３、入出力インターフェース９５および通信インターフェース９６を備える。図１１においては、インターフェースをＩ／Ｆ（Interface）と略して表記する。プロセッサ９１、主記憶装置９２、補助記憶装置９３、入出力インターフェース９５および通信インターフェース９６は、バス９９を介して互いにデータ通信可能に接続される。また、プロセッサ９１、主記憶装置９２、補助記憶装置９３および入出力インターフェース９５は、通信インターフェース９６を介して、インターネットやイントラネットなどのネットワークに接続される。
　プロセッサ９１は、補助記憶装置９３等に格納されたプログラムを主記憶装置９２に展開し、展開されたプログラムを実行する。本実施形態においては、情報処理装置９０にインストールされたソフトウェアプログラムを用いる構成とすればよい。プロセッサ９１は、本実施形態に係る波源方向推定装置による処理を実行する。
　主記憶装置９２は、プログラムが展開される領域を有する。主記憶装置９２は、例えばＤＲＡＭ（Dynamic Random Access Memory）などの揮発性メモリとすればよい。また、ＭＲＡＭ（Magnetoresistive Random Access Memory）などの不揮発性メモリを主記憶装置９２として構成・追加してもよい。
　補助記憶装置９３は、種々のデータを記憶する。補助記憶装置９３は、ハードディスクやフラッシュメモリなどのローカルディスクによって構成される。なお、種々のデータを主記憶装置９２に記憶させる構成とし、補助記憶装置９３を省略することも可能である。
　入出力インターフェース９５は、情報処理装置９０と周辺機器とを接続するためのインターフェースである。通信インターフェース９６は、規格や仕様に基づいて、インターネットやイントラネットなどのネットワークを通じて、外部のシステムや装置に接続するためのインターフェースである。入出力インターフェース９５および通信インターフェース９６は、外部機器と接続するインターフェースとして共通化してもよい。
　情報処理装置９０には、必要に応じて、キーボードやマウス、タッチパネルなどの入力機器を接続するように構成してもよい。それらの入力機器は、情報や設定の入力に使用される。なお、タッチパネルを入力機器として用いる場合は、表示機器の表示画面が入力機器のインターフェースを兼ねる構成とすればよい。プロセッサ９１と入力機器との間のデータ通信は、入出力インターフェース９５に仲介させればよい。
　また、情報処理装置９０には、情報を表示するための表示機器を備え付けてもよい。表示機器を備え付ける場合、情報処理装置９０には、表示機器の表示を制御するための表示制御装置（図示しない）が備えられていることが好ましい。表示機器は、入出力インターフェース９５を介して情報処理装置９０に接続すればよい。
　また、情報処理装置９０には、必要に応じて、ディスクドライブを備え付けてもよい。ディスクドライブは、バス９９に接続される。ディスクドライブは、プロセッサ９１と図示しない記録媒体（プログラム記録媒体）との間で、記録媒体からのデータ・プログラムの読み出し、情報処理装置９０の処理結果の記録媒体への書き込みなどを仲介する。記録媒体は、例えば、ＣＤ（Compact Disc）やＤＶＤ（Digital Versatile Disc）などの光学記録媒体で実現できる。また、記録媒体は、ＵＳＢ（Universal Serial Bus）メモリやＳＤ（Secure Digital）カードなどの半導体記録媒体や、フレキシブルディスクなどの磁気記録媒体、その他の記録媒体によって実現してもよい。
　以上が、各実施形態に係る波源方向推定装置を可能とするためのハードウェア構成の一例である。なお、図１１のハードウェア構成は、各実施形態に係る波源方向推定装置の演算処理を実行するためのハードウェア構成の一例であって、本発明の範囲を限定するものではない。また、各実施形態に係る波源方向推定装置に関する処理をコンピュータに実行させるプログラムも本発明の範囲に含まれる。さらに、各実施形態に係るプログラムを記録したプログラム記録媒体も本発明の範囲に含まれる。
　各実施形態の波源方向推定装置の構成要素は、任意に組み合わせることができる。また、各実施形態の波源方向推定装置の構成要素は、ソフトウェアによって実現してもよいし、回路によって実現してもよい。
　以上、実施形態を参照して本発明を説明してきたが、本発明は上記実施形態に限定されるものではない。本発明の構成や詳細には、本発明のスコープ内で当業者が理解し得る様々な変更をすることができる。
　上記の実施形態の一部又は全部は、以下の付記のようにも記載されうるが、以下には限られない。
（付記１）
　複数のセンサによって取得された波動から変換される電気信号を入力信号として取得する複数の入力手段と、
　複数の前記入力信号から少なくとも二つの前記入力信号を組み合わせたペアを少なくとも二つ選択する信号選択手段と、
　前記入力信号のペアを構成する少なくとも二つの前記入力信号間において波源探索方向ごとの波動の到達時間差を相対遅延時間として算出する相対遅延時間計算手段と、
　前記入力信号のペアと前記相対遅延時間とを用いて前記波動の波源の推定方向情報を周波数ごとに生成する少なくとも一つの周波数別推定方向情報生成手段と、
　前記周波数別推定方向情報生成手段によって生成される周波数ごとの前記推定方向情報を統合する統合手段とを備える波源方向推定装置。
（付記２）
　前記信号選択手段は、
　前記センサの間隔に基づいて、複数の前記入力信号から少なくとも二つの前記入力信号を組み合わせたペアを選択する付記１に記載の波源方向推定装置。
（付記３）
　前記相対遅延時間計算手段は、
　一つの前記入力信号のペアの供給元である前記センサのペアに対する前記波源探索方向を基準として、前記信号選択手段によって選択された全ての入力信号のペアの相対遅延時間を、基準である前記波源探索方向の関数として計算する付記１または２に記載の波源方向推定装置。
（付記４）
　前記周波数別推定方向情報生成手段は、
　ペアをなす少なくとも２つの前記入力信号を周波数領域の変換信号に変換する変換手段と、
　前記変換手段によって変換された前記変換信号を用いてクロススペクトルを計算するクロススペクトル計算手段と、
　前記クロススペクトル計算手段によって算出される前記クロススペクトルを用いて平均クロススペクトルを計算する平均計算手段と、
　前記平均計算手段によって算出される平均クロススペクトルを用いて分散を計算する分散計算手段と、
　前記平均計算手段によって算出される前記平均クロススペクトルと、前記分散計算手段によって算出される分散とを用いて周波数別クロススペクトルを計算する周波数別クロススペクトル生成手段と、
　前記周波数別クロススペクトル生成手段によって算出される前記周波数別クロススペクトルを逆変換して周波数別相互相関関数を計算する逆変換手段と、
　前記逆変換手段によって算出される前記周波数別相互相関関数と前記相対遅延時間とを用いて周波数ごとに前記推定方向情報を計算する周波数別推定方向情報計算手段とを有する付記１乃至３のいずれか一項に記載の波源方向推定装置。
（付記５）
　前記周波数別クロススペクトル生成手段は、
　前記平均計算手段から前記平均クロススペクトルを取得し、取得した前記平均クロススペクトルを用いて周波数別基本クロススペクトルを計算する周波数別基本クロススペクトル計算手段と、
　前記分散計算手段から前記分散を取得し、取得した前記分散を用いてカーネル関数スペクトルを計算するカーネル関数スペクトル生成手段と、
　前記周波数別基本クロススペクトル計算手段によって算出される前記周波数別基本クロススペクトルと、前記カーネル関数スペクトル生成手段によって算出される前記カーネル関数スペクトルとの積を計算して前記周波数別クロススペクトルを計算する乗算手段とを含む付記４に記載の波源方向推定装置。
（付記６）
　前記統合手段は、
　複数の周波数ごとの前記推定方向情報を複数の前記入力信号のペアについて統合した周波数別統合推定方向情報を算出し、算出した前記周波数別統合推定方向情報を全ての周波数について統合して統合推定方向情報を算出する付記１乃至５のいずれか一項に記載の波源方向推定装置。
（付記７）
　前記統合手段は、
　複数の周波数ごとの前記推定方向情報を全ての周波数について統合した入力信号組み合わせ別統合推定方向情報を算出し、算出した前記入力信号組み合わせ別統合推定方向情報を全ての入力信号の組み合わせについて統合して統合推定方向情報を算出する付記１乃至５のいずれか一項に記載の波源方向推定装置。
（付記８）
　前記統合手段によって算出される前記統合推定方向情報に基づいて前記波動の波源方向を算出する波源方向算出手段を備える付記１乃至７のいずれか一項に記載の波源方向推定装置。
（付記９）
　前記波源方向算出手段は、
　一定時間ごとに前記統合推定方向情報が最大になる時刻に対応する方向を波源方向として算出する付記８に記載の波源方向推定装置。
（付記１０）
　複数の前記入力手段のそれぞれに対応して配置される前記センサを備える付記１乃至９のいずれか一項に記載の波源方向推定装置。
（付記１１）
　情報処理装置が、
　複数のセンサによって取得された波動から変換される電気信号を入力信号として取得し、
　複数の前記入力信号から少なくとも二つの前記入力信号を組み合わせたペアを少なくとも二つ選択し、
　前記入力信号のペアを構成する少なくとも二つの前記入力信号間において波源探索方向ごとの波動の到達時間差を相対遅延時間として算出し、
　前記入力信号のペアと前記相対遅延時間とを用いて前記波動の波源の推定方向情報を周波数ごとに少なくとも一つ生成し、
　周波数ごとの前記推定方向情報を統合する波源方向推定方法。
（付記１２）
　複数のセンサによって取得された波動から変換される電気信号を入力信号として取得する処理と、
　複数の前記入力信号から少なくとも二つの前記入力信号を組み合わせたペアを少なくとも二つ選択する処理と、
　前記入力信号のペアを構成する少なくとも二つの前記入力信号間において波源探索方向ごとの波動の到達時間差を相対遅延時間として算出する処理と、
　前記入力信号のペアと前記相対遅延時間とを用いて前記波動の波源の推定方向情報を周波数ごとに少なくとも一つ生成する処理と、
　周波数ごとの前記推定方向情報を統合する処理とコンピュータに実行させるプログラムを記録させたプログラム記録媒体。
　１０、２０　　波源方向推定装置
　１１、２１　　入力端子
　１２、２２　　信号選択部
　１３、２３　　相対遅延時間計算部
　１５、２５　　周波数別推定方向情報生成部
　１７、２７　　統合部
　２８　　波源方向算出部
　１５１　　変換部
　１５２　　クロススペクトル計算部
　１５３　　平均計算部
　１５４　　分散計算部
　１５５　　周波数別クロススペクトル生成部
　１５６　　逆変換部
　１５７　　周波数別推定方向情報計算部
　５５１　　周波数別基本クロススペクトル計算部
　５５２　　カーネル関数スペクトル生成部
　５５３　　乗算部

Claims (12)

  1. 　複数のセンサによって取得された波動から変換される電気信号を入力信号として取得する複数の入力手段と、
    　複数の前記入力信号から少なくとも二つの前記入力信号を組み合わせたペアを少なくとも二つ選択する信号選択手段と、
    　前記入力信号のペアを構成する少なくとも二つの前記入力信号間において波源探索方向ごとの波動の到達時間差を相対遅延時間として算出する相対遅延時間計算手段と、
    　前記入力信号のペアと前記相対遅延時間とを用いて前記波動の波源の推定方向情報を周波数ごとに生成する少なくとも一つの周波数別推定方向情報生成手段と、
    　前記周波数別推定方向情報生成手段によって生成される周波数ごとの前記推定方向情報を統合する統合手段とを備える波源方向推定装置。
  2. 　前記信号選択手段は、
    　前記センサの間隔に基づいて、複数の前記入力信号から少なくとも二つの前記入力信号を組み合わせたペアを選択する請求項１に記載の波源方向推定装置。
  3. 　前記相対遅延時間計算手段は、
    　一つの前記入力信号のペアの供給元である前記センサのペアに対する前記波源探索方向を基準として、前記信号選択手段によって選択された全ての入力信号のペアの相対遅延時間を、基準である前記波源探索方向の関数として計算する請求項１または２に記載の波源方向推定装置。
  4. 　前記周波数別推定方向情報生成手段は、
    　ペアをなす少なくとも２つの前記入力信号を周波数領域の変換信号に変換する変換手段と、
    　前記変換手段によって変換された前記変換信号を用いてクロススペクトルを計算するクロススペクトル計算手段と、
    　前記クロススペクトル計算手段によって算出される前記クロススペクトルを用いて平均クロススペクトルを計算する平均計算手段と、
    　前記平均計算手段によって算出される平均クロススペクトルを用いて分散を計算する分散計算手段と、
    　前記平均計算手段によって算出される前記平均クロススペクトルと、前記分散計算手段によって算出される分散とを用いて周波数別クロススペクトルを計算する周波数別クロススペクトル生成手段と、
    　前記周波数別クロススペクトル生成手段によって算出される前記周波数別クロススペクトルを逆変換して周波数別相互相関関数を計算する逆変換手段と、
    　前記逆変換手段によって算出される前記周波数別相互相関関数と前記相対遅延時間とを用いて周波数ごとに前記推定方向情報を計算する周波数別推定方向情報計算手段とを有する請求項１乃至３のいずれか一項に記載の波源方向推定装置。
  5. 　前記周波数別クロススペクトル生成手段は、
    　前記平均計算手段から前記平均クロススペクトルを取得し、取得した前記平均クロススペクトルを用いて周波数別基本クロススペクトルを計算する周波数別基本クロススペクトル計算手段と、
    　前記分散計算手段から前記分散を取得し、取得した前記分散を用いてカーネル関数スペクトルを計算するカーネル関数スペクトル生成手段と、
    　前記周波数別基本クロススペクトル計算手段によって算出される前記周波数別基本クロススペクトルと、前記カーネル関数スペクトル生成手段によって算出される前記カーネル関数スペクトルとの積を計算して前記周波数別クロススペクトルを計算する乗算手段とを含む請求項４に記載の波源方向推定装置。
  6. 　前記統合手段は、
    　複数の周波数ごとの前記推定方向情報を複数の前記入力信号のペアについて統合した周波数別統合推定方向情報を算出し、算出した前記周波数別統合推定方向情報を全ての周波数について統合して統合推定方向情報を算出する請求項１乃至５のいずれか一項に記載の波源方向推定装置。
  7. 　前記統合手段は、
    　複数の周波数ごとの前記推定方向情報を全ての周波数について統合した入力信号組み合わせ別統合推定方向情報を算出し、算出した前記入力信号組み合わせ別統合推定方向情報を全ての入力信号の組み合わせについて統合して統合推定方向情報を算出する請求項１乃至５のいずれか一項に記載の波源方向推定装置。
  8. 　前記統合手段によって算出される前記統合推定方向情報に基づいて前記波動の波源方向を算出する波源方向算出手段を備える請求項７に記載の波源方向推定装置。
  9. 　前記波源方向算出手段は、
    　一定時間ごとに前記統合推定方向情報が最大になる時刻に対応する方向を波源方向として算出する請求項８に記載の波源方向推定装置。
  10. 　複数の前記入力手段のそれぞれに対応して配置される前記センサを備える請求項１乃至９のいずれか一項に記載の波源方向推定装置。
  11. 　情報処理装置が、
    　複数のセンサによって取得された波動から変換される電気信号を入力信号として取得し、
    　複数の前記入力信号から少なくとも二つの前記入力信号を組み合わせたペアを少なくとも二つ選択し、
    　前記入力信号のペアを構成する少なくとも二つの前記入力信号間において波源探索方向ごとの波動の到達時間差を相対遅延時間として算出し、
    　前記入力信号のペアと前記相対遅延時間とを用いて前記波動の波源の推定方向情報を周波数ごとに少なくとも一つ生成し、
    　周波数ごとの前記推定方向情報を統合する波源方向推定方法。
  12. 　複数のセンサによって取得された波動から変換される電気信号を入力信号として取得する処理と、
    　複数の前記入力信号から少なくとも二つの前記入力信号を組み合わせたペアを少なくとも二つ選択する処理と、
    　前記入力信号のペアを構成する少なくとも二つの前記入力信号間において波源探索方向ごとの波動の到達時間差を相対遅延時間として算出する処理と、
    　前記入力信号のペアと前記相対遅延時間とを用いて前記波動の波源の推定方向情報を周波数ごとに少なくとも一つ生成する処理と、
    　周波数ごとの前記推定方向情報を統合する処理とコンピュータに実行させるプログラムを記録させたプログラム記録媒体。
PCT/JP2018/023970 2018-06-25 2018-06-25 波源方向推定装置、波源方向推定方法、およびプログラム記録媒体 Ceased WO2020003342A1 (ja)

Priority Applications (3)

Application Number Priority Date Filing Date Title
PCT/JP2018/023970 WO2020003342A1 (ja) 2018-06-25 2018-06-25 波源方向推定装置、波源方向推定方法、およびプログラム記録媒体
JP2020526724A JP7056739B2 (ja) 2018-06-25 2018-06-25 波源方向推定装置、波源方向推定方法、およびプログラム
US17/252,391 US20210263125A1 (en) 2018-06-25 2018-06-25 Wave-source-direction estimation device, wave-source-direction estimation method, and program storage medium

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
PCT/JP2018/023970 WO2020003342A1 (ja) 2018-06-25 2018-06-25 波源方向推定装置、波源方向推定方法、およびプログラム記録媒体

Publications (1)

Publication Number Publication Date
WO2020003342A1 true WO2020003342A1 (ja) 2020-01-02

Family

ID=68986225

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/JP2018/023970 Ceased WO2020003342A1 (ja) 2018-06-25 2018-06-25 波源方向推定装置、波源方向推定方法、およびプログラム記録媒体

Country Status (3)

Country Link
US (1) US20210263125A1 (ja)
JP (1) JP7056739B2 (ja)
WO (1) WO2020003342A1 (ja)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2023054780A (ja) * 2021-10-04 2023-04-14 ノキア テクノロジーズ オサケユイチア 空間オーディオキャプチャ

Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2008089312A (ja) * 2006-09-29 2008-04-17 Kddi Corp 信号到来方向推定装置及び方法、並びに信号分離装置及び方法、コンピュータプログラム
JP2010175431A (ja) * 2009-01-30 2010-08-12 Nippon Telegr & Teleph Corp <Ntt> 音源方向推定装置とその方法と、プログラム
JP2010193323A (ja) * 2009-02-19 2010-09-02 Casio Hitachi Mobile Communications Co Ltd 録音装置、再生装置、録音方法、再生方法、及びコンピュータプログラム
JP2011071686A (ja) * 2009-09-25 2011-04-07 Nec Corp 映像音響処理装置、映像音響処理方法及びプログラム
JP2015161551A (ja) * 2014-02-26 2015-09-07 株式会社東芝 音源方向推定装置、音源方向推定方法およびプログラム
WO2018003158A1 (ja) * 2016-06-29 2018-01-04 日本電気株式会社 相関関数生成装置、相関関数生成方法、相関関数生成プログラムおよび波源方向推定装置
WO2018131099A1 (ja) * 2017-01-11 2018-07-19 日本電気株式会社 相関関数生成装置、相関関数生成方法、相関関数生成プログラムおよび波源方向推定装置

Family Cites Families (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US4791672A (en) * 1984-10-05 1988-12-13 Audiotone, Inc. Wearable digital hearing aid and method for improving hearing ability
JP3536996B2 (ja) * 1994-09-13 2004-06-14 ソニー株式会社 パラメータ変換方法及び音声合成方法
US6266003B1 (en) * 1998-08-28 2001-07-24 Sigma Audio Research Limited Method and apparatus for signal processing for time-scale and/or pitch modification of audio signals
JP2005172760A (ja) * 2003-12-15 2005-06-30 Mitsubishi Electric Corp 方位探知装置
US7733224B2 (en) * 2006-06-30 2010-06-08 Bao Tran Mesh network personal emergency response appliance
WO2017039632A1 (en) * 2015-08-31 2017-03-09 Nunntawi Dynamics Llc Passive self-localization of microphone arrays
US10957338B2 (en) * 2018-05-16 2021-03-23 Synaptics Incorporated 360-degree multi-source location detection, tracking and enhancement

Patent Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2008089312A (ja) * 2006-09-29 2008-04-17 Kddi Corp 信号到来方向推定装置及び方法、並びに信号分離装置及び方法、コンピュータプログラム
JP2010175431A (ja) * 2009-01-30 2010-08-12 Nippon Telegr & Teleph Corp <Ntt> 音源方向推定装置とその方法と、プログラム
JP2010193323A (ja) * 2009-02-19 2010-09-02 Casio Hitachi Mobile Communications Co Ltd 録音装置、再生装置、録音方法、再生方法、及びコンピュータプログラム
JP2011071686A (ja) * 2009-09-25 2011-04-07 Nec Corp 映像音響処理装置、映像音響処理方法及びプログラム
JP2015161551A (ja) * 2014-02-26 2015-09-07 株式会社東芝 音源方向推定装置、音源方向推定方法およびプログラム
WO2018003158A1 (ja) * 2016-06-29 2018-01-04 日本電気株式会社 相関関数生成装置、相関関数生成方法、相関関数生成プログラムおよび波源方向推定装置
WO2018131099A1 (ja) * 2017-01-11 2018-07-19 日本電気株式会社 相関関数生成装置、相関関数生成方法、相関関数生成プログラムおよび波源方向推定装置

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
KATO, MASANORI ET AL.: "TDOA Estimation Based on Phase-Voting Cross Correlation and Circular Standard Deviation", 2017 25TH EUROPEAN SIGNAL PROCESSING CONFERENCE(EUSIPCO, 2017, pages 1230 - 1234, XP033236133, ISBN: 978-0-9928626-7-1, DOI: 10.23919/EUSIPCO.2017.8081404 *
KATO, MASANORI ET AL.: "TDOA Estimation Based on Phase-Voting Cross Correlation and Circular Standard Deviation", IEICE TECHNICAL REPORT, vol. 6, no. 477, 11 March 2017 (2017-03-11), pages 1230 - 1234, XP033236133, ISSN: 0913-5685 *

Cited By (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2023054780A (ja) * 2021-10-04 2023-04-14 ノキア テクノロジーズ オサケユイチア 空間オーディオキャプチャ
US12323762B2 (en) 2021-10-04 2025-06-03 Nokia Technologies Oy Spatial audio capture
JP7708730B2 (ja) 2021-10-04 2025-07-15 ノキア テクノロジーズ オサケユイチア 空間オーディオキャプチャ

Also Published As

Publication number Publication date
JPWO2020003342A1 (ja) 2021-06-24
JP7056739B2 (ja) 2022-04-19
US20210263125A1 (en) 2021-08-26

Similar Documents

Publication Publication Date Title
JP6933303B2 (ja) 波源方向推定装置、波源方向推定方法、およびプログラム
KR101025163B1 (ko) 진동 및 소음 전달경로 해석 시스템과 진동 및 소음 전달경로 해석 방법
JP6289936B2 (ja) 音源方向推定装置、音源方向推定方法およびプログラム
JP6392289B2 (ja) レーダーおよびソナー用途における固定小数点高速フーリエ変換のスケーリング
CN109923430B (zh) 用于进行相位差展开的装置及方法
US9858946B2 (en) Signal processing apparatus, signal processing method, and signal processing program
JPWO2018116943A1 (ja) ノイズ抑圧装置、ノイズ抑圧方法、及びこれらを用いた受信装置、受信方法
JP6862799B2 (ja) 信号処理装置、方位算出方法及び方位算出プログラム
Padois Acoustic source localization based on the generalized cross-correlation and the generalized mean with few microphones
JP2017151076A (ja) 音源探査装置、音源探査方法およびそのプログラム
Zhang et al. Transient nearfield acoustic holography based on an interpolated time-domain equivalent source method
JP2010127771A (ja) 合成開口ソーナー、合成開口ソーナーの位相誤差補正方法及びプログラム
Padois et al. On the use of geometric and harmonic means with the generalized cross-correlation in the time domain to improve noise source maps
JP2010185682A (ja) 一般調和解析装置および周波数分析装置
JP7056739B2 (ja) 波源方向推定装置、波源方向推定方法、およびプログラム
JP2007288463A (ja) 特性取得装置、方法およびプログラム
JP2648110B2 (ja) 信号検出方法及び装置
JP2010169644A (ja) 目標信号検出装置、方法及びプログラム
JP2012063152A (ja) 速度推定装置及びコンピュータプログラム及び速度推定方法
JP5034920B2 (ja) エコーキャンセル装置、カラオケ装置、エコーキャンセル方法およびプログラム
WO2015053068A1 (ja) 音場測定装置、音場測定方法および音場測定プログラム
JP7196504B2 (ja) 音響特性計測システム、音響特性計測方法、およびプログラム
JP6119604B2 (ja) 信号処理装置、信号処理方法、および信号処理プログラム
KR101294405B1 (ko) 위상 변환된 잡음 신호를 이용한 음성 영역 검출 방법 및 그 장치
JP7023324B2 (ja) マルチメディア信号認識のための電子装置およびその動作方法

Legal Events

Date Code Title Description
121 Ep: the epo has been informed by wipo that ep was designated in this application

Ref document number: 18923984

Country of ref document: EP

Kind code of ref document: A1

ENP Entry into the national phase

Ref document number: 2020526724

Country of ref document: JP

Kind code of ref document: A

NENP Non-entry into the national phase

Ref country code: DE

122 Ep: pct application non-entry in european phase

Ref document number: 18923984

Country of ref document: EP

Kind code of ref document: A1