WO2018003158A1 - 相関関数生成装置、相関関数生成方法、相関関数生成プログラムおよび波源方向推定装置 - Google Patents

相関関数生成装置、相関関数生成方法、相関関数生成プログラムおよび波源方向推定装置 Download PDF

Info

Publication number
WO2018003158A1
WO2018003158A1 PCT/JP2017/004028 JP2017004028W WO2018003158A1 WO 2018003158 A1 WO2018003158 A1 WO 2018003158A1 JP 2017004028 W JP2017004028 W JP 2017004028W WO 2018003158 A1 WO2018003158 A1 WO 2018003158A1
Authority
WO
WIPO (PCT)
Prior art keywords
frequency
correlation function
cross spectrum
spectrum
calculating
Prior art date
Application number
PCT/JP2017/004028
Other languages
English (en)
French (fr)
Inventor
正徳 加藤
裕三 仙田
Original Assignee
日本電気株式会社
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 日本電気株式会社 filed Critical 日本電気株式会社
Priority to JP2018524877A priority Critical patent/JPWO2018003158A1/ja
Priority to US16/309,542 priority patent/US20190250240A1/en
Publication of WO2018003158A1 publication Critical patent/WO2018003158A1/ja

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/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
    • 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
    • G—PHYSICS
    • G10—MUSICAL INSTRUMENTS; ACOUSTICS
    • G10L—SPEECH ANALYSIS OR SYNTHESIS; SPEECH RECOGNITION; SPEECH OR VOICE PROCESSING; 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/48—Speech or voice analysis techniques not restricted to a single one of groups G10L15/00 - G10L21/00 specially adapted for particular use
    • G10L25/51—Speech or voice analysis techniques not restricted to a single one of groups G10L15/00 - G10L21/00 specially adapted for particular use for comparison or discrimination
    • H—ELECTRICITY
    • H04—ELECTRIC COMMUNICATION TECHNIQUE
    • H04R—LOUDSPEAKERS, MICROPHONES, GRAMOPHONE PICK-UPS OR LIKE ACOUSTIC ELECTROMECHANICAL TRANSDUCERS; DEAF-AID SETS; PUBLIC ADDRESS SYSTEMS
    • H04R3/00—Circuits for transducers, loudspeakers or microphones
    • H04R3/005—Circuits for transducers, loudspeakers or microphones for combining the signals of two or more microphones
    • 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 OR SYNTHESIS; SPEECH RECOGNITION; SPEECH OR VOICE PROCESSING; 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
    • H—ELECTRICITY
    • H04—ELECTRIC COMMUNICATION TECHNIQUE
    • H04R—LOUDSPEAKERS, MICROPHONES, GRAMOPHONE PICK-UPS OR LIKE ACOUSTIC ELECTROMECHANICAL TRANSDUCERS; DEAF-AID SETS; 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]
    • H—ELECTRICITY
    • H04—ELECTRIC COMMUNICATION TECHNIQUE
    • H04R—LOUDSPEAKERS, MICROPHONES, GRAMOPHONE PICK-UPS OR LIKE ACOUSTIC ELECTROMECHANICAL TRANSDUCERS; DEAF-AID SETS; PUBLIC ADDRESS SYSTEMS
    • H04R2499/00—Aspects covered by H04R or H04S not otherwise provided for in their subgroups
    • H04R2499/10—General applications
    • H04R2499/11—Transducers incorporated or for use in hand-held devices, e.g. mobile phones, PDA's, camera's
    • H—ELECTRICITY
    • H04—ELECTRIC COMMUNICATION TECHNIQUE
    • H04R—LOUDSPEAKERS, MICROPHONES, GRAMOPHONE PICK-UPS OR LIKE ACOUSTIC ELECTROMECHANICAL TRANSDUCERS; DEAF-AID SETS; PUBLIC ADDRESS SYSTEMS
    • H04R2499/00—Aspects covered by H04R or H04S not otherwise provided for in their subgroups
    • H04R2499/10—General applications
    • H04R2499/15—Transducers incorporated in visual displaying devices, e.g. televisions, computer displays, laptops

Definitions

  • the present invention relates to a correlation function generation device, a correlation function generation method, a correlation function generation program, and a wave source direction estimation device.
  • Non-Patent Document 1 and Non-Patent Document 2 describe a method of estimating the direction of a sound source (sound wave generation source, generation location) using sound reception signals of two microphones. Specifically, a cross correlation function between the received sound signals is obtained from the two received sound signals. And the technique which estimates the arrival direction of a sound wave by calculating the time difference which a cross correlation function gives the maximum value as the arrival time difference of a sound wave is disclosed.
  • An object of the present invention is to provide a technique for solving the above-described problems.
  • a correlation function generation device provides: A plurality of input signal acquisition means for acquiring a wave generated by the wave source as an input signal; Conversion means for converting the plurality of input signals acquired by the input signal acquisition means into a plurality of frequency domain signals; Cross spectrum calculation means for calculating a cross spectrum based on the frequency domain signal; A frequency-specific cross-spectrum calculating means for calculating a frequency-specific cross-spectrum based on the cross-spectrum; An integrated correlation function calculating means for calculating an integrated correlation function based on the frequency-specific cross spectrum; Equipped with.
  • a correlation function generation method includes: A plurality of input signal acquisition steps for acquiring waves generated by the wave source as input signals; A conversion step of converting the plurality of input signals acquired in the input signal acquisition step into a plurality of frequency domain signals; A cross spectrum calculating step of calculating a cross spectrum based on the frequency domain signal; A frequency-specific cross spectrum calculation step of calculating a frequency-specific cross spectrum based on the cross spectrum; An integrated correlation function calculating step of calculating an integrated correlation function based on the cross spectrum by frequency; including.
  • a correlation function generation program provides: A plurality of input signal acquisition steps for acquiring waves generated by the wave source as input signals; A conversion step of converting the plurality of input signals acquired in the input signal acquisition step into a plurality of frequency domain signals; A cross spectrum calculating step of calculating a cross spectrum based on the frequency domain signal; A frequency-specific cross spectrum calculation step of calculating a frequency-specific cross spectrum based on the cross spectrum; An integrated correlation function calculating step of calculating an integrated correlation function based on the cross spectrum by frequency; Is executed on the computer.
  • a wave source direction estimation device In order to achieve the above object, a wave source direction estimation device according to the present invention is provided.
  • a correlation function having a clear peak can be generated even in an environment where the ambient noise level is high.
  • the direction of the wave source can be estimated with high accuracy.
  • the estimation target of the wave source direction estimation device is not limited to the generation source of sound waves that are vibration waves of air or water. It can also be applied to sources of vibration waves using earth and solids such as earthquakes and landslides. In that case, a vibration sensor is used instead of a microphone as a device that converts vibration waves into electrical signals.
  • the wave source direction estimation apparatus can be applied to the case where the direction is estimated using radio waves as well as gas / liquid / solid vibration waves. In that case, an antenna is used as a device that converts radio waves into electrical signals.
  • the wave source is assumed to be a sound source.
  • a correlation function generation device 100 as a first embodiment of the present invention will be described with reference to FIG.
  • the correlation function generation device 100 is a device that generates a correlation function based on an input signal.
  • the correlation function generation device 100 includes an input signal acquisition unit 101, a conversion unit 102, a cross spectrum calculation unit 103, a frequency-specific cross spectrum calculation unit 104, and an integrated correlation function calculation unit 105. Including.
  • the plurality of input signal acquisition units 101 acquire waves generated by the wave source as input signals.
  • the conversion unit 102 converts the plurality of input signals acquired by the input signal acquisition unit into a plurality of frequency domain signals.
  • the cross spectrum calculation unit 103 calculates a cross spectrum based on the frequency domain signal.
  • the frequency-specific cross spectrum calculation unit 104 calculates the frequency-specific cross spectrum based on the cross spectrum.
  • the integrated correlation function calculation unit 105 calculates an integrated correlation function based on the frequency-specific cross spectrum.
  • a correlation function having a clear peak can be generated even in an environment where the ambient noise level is high.
  • the direction of the wave source can be estimated with high accuracy.
  • Non-Patent Document 1 and Non-Patent Document 2 it is difficult to accurately estimate the direction of a sound source that exists far away in an environment with a high ambient noise level such as outdoors.
  • a high ambient noise level such as outdoors.
  • the sound source to be estimated target sound source
  • the volume of sound radiated from the target sound source is significantly reduced when reaching the microphone.
  • the sound of the target sound source is buried in the ambient noise, and a correlation function having a clear peak cannot be generated. For this reason, the direction estimation accuracy of the target sound source may be reduced.
  • FIG. 2A is a block diagram showing the configuration of the wave source direction estimating apparatus according to the present embodiment.
  • FIG. 2B is a block diagram illustrating a configuration of an integrated correlation function calculation unit included in the wave source direction estimation apparatus according to the present embodiment.
  • the wave source direction estimation apparatus 200 functions as part of a device such as a digital video camera, a smartphone, a mobile phone, a laptop computer, or a passive sonar. It is also installed in an abnormal sound detection device that detects abnormalities based on voices and sounds such as suspicious drone detection, scream detection, and vehicle accident detection.
  • a device such as a digital video camera, a smartphone, a mobile phone, a laptop computer, or a passive sonar. It is also installed in an abnormal sound detection device that detects abnormalities based on voices and sounds such as suspicious drone detection, scream detection, and vehicle accident detection.
  • application examples of the wave source direction estimation apparatus 200 according to the present embodiment are not limited to these, and can be applied to any wave source direction estimation apparatus that is required to estimate the direction of a target sound source from received sound.
  • Wave source direction estimation apparatus 200 includes an input terminal 20 1, the input terminal 20 2, a converting unit 201, a cross spectrum calculation section 202, 1 and ⁇ 203 k frequency-cross spectrum calculation unit 203, a.
  • the wave source direction estimation apparatus 200 further includes an integrated correlation function calculation unit 204, an estimation direction information generation unit 205, and a relative delay time calculation unit 206.
  • a microphone (hereinafter, microphone) a sound collector sound and digital signal varying noise from ambient are mixed in (the sample value series) Is entered as A sound signal input to the input terminal 20 1 and the input terminal 20 2 in the present embodiment referred to as an input signal. Then, representing the input signal at the input terminal 20 1 at time t x 1 (t), the input signal of the input terminal 20 2 x 2 and (t).
  • the sound input to the input terminal is collected by a microphone that is a sound collector. Since there are a plurality of input terminals, when collecting the sound of the target sound source, two microphones as many as the number of terminals are used simultaneously. In the present embodiment, it is assumed that the input terminal and the microphone correspond one-to-one, and the sound collected by the mth microphone is supplied to the mth input terminal. Therefore, the input signal supplied to the mth input terminal is also referred to as “mth microphone input signal”.
  • the wave source direction estimation apparatus 200 estimates the direction of the sound source using the time difference at which the sound of the target sound source reaches the two microphones. For this reason, since the microphone interval is also important information, not only the input signal but also the microphone position information is supplied to the wave source direction estimation apparatus 200.
  • Conversion unit 201 converts an input signal supplied from the input terminal 20 1 and the input terminal 20 2 which supplies to the cross spectrum calculation section 202. The conversion is performed for the purpose of decomposing the input signal into a plurality of frequency components.
  • a case where a representative Fourier transform is used will be described.
  • Two types of input signals x m (t) are input to the conversion unit 201.
  • m is an input terminal number.
  • the conversion unit 201 cuts out a waveform having an appropriate length from the input signal supplied from the input terminal while shifting the waveform with a certain period.
  • the signal section cut out in this way is called a frame
  • the length of the cut-out waveform is called a frame length
  • the period of shifting the frame is called a frame period.
  • the signal cut out using Fourier transformation is converted into a frequency domain signal.
  • n the frame number
  • j represents an imaginary unit (-1 square root), and exp represents an exponential function.
  • K represents a frequency bin number and is an integer from 0 to K-1.
  • k is not simply a frequency bin number but simply called “frequency”.
  • Cross spectrum calculation unit 202 calculates a cross spectral, frequency-cross spectrum calculation unit 203 1, 203 2, ..., and transferred to 203 K.
  • the cross spectrum calculation unit 202 calculates the product of the complex conjugate of the converted signal X 2 (k, n) and the converted signal X 1 (k, n).
  • the cross spectrum of the converted signal is S 12 (k, n)
  • the cross spectrum is calculated as follows.
  • conj (X 2 (k, n)) represents the complex conjugate of X 2 (k, n).
  • the frequency-specific cross spectrum calculation sections 203 1 , 203 2 ,..., 203 K use the cross spectrum S 12 (k, n) supplied from the cross spectrum calculation section 202 to calculate S 12 (k, n).
  • a cross spectrum corresponding to each frequency k is calculated and transmitted to the integrated correlation function calculation unit 204 as a frequency-specific cross spectrum.
  • the cross spectrum for each frequency is calculated to calculate a correlation function for each frequency component. That is, a cross spectrum for each frequency is calculated in order to obtain a correlation function corresponding to a certain frequency k (referred to as a correlation function for each frequency) at a later stage.
  • FIG. 2B is a block diagram of the frequency-specific cross spectrum calculation unit 203 k .
  • the frequency-specific cross spectrum calculation unit 203 k includes a frequency-specific basic cross spectrum calculation unit 2031 k .
  • the frequency-specific cross spectrum calculation unit 203 k uses the cross spectrum S 12 (k, n) supplied from the cross spectrum calculation unit 202 to calculate a frequency-specific basic cross spectrum, and calculates an integrated correlation function as a frequency-specific cross spectrum.
  • To the unit 204 uses the cross spectrum S 12 (k, n) supplied from the cross spectrum calculation unit 202 to calculate a frequency-specific basic cross spectrum, and calculates an integrated correlation function as a frequency-specific cross spectrum.
  • the frequency-based basic cross spectrum calculation unit 2031 k obtains the phase component and the amplitude component separately in advance when calculating the frequency-based basic cross spectrum based on the cross spectrum S 12 (k, n) of the frequency k. To integrate. Assuming that the fundamental cross spectrum U k (w, n) by frequency of frequency k, its amplitude component
  • w represents a frequency and is an integer of 0 or more and W ⁇ 1 or less.
  • and the phase component arg (U k (w, n)) of the frequency-specific basic cross spectrum from the cross spectrum S 12 (k, n) of the frequency k. Will be described.
  • 1.0 is used as a frequency obtained by multiplying k by an integer.
  • the phase component of the frequency that is a non-constant multiple of the frequency k is set to zero.
  • p is an integer of 1 or more and P or less. Since the important information when performing the wave source direction estimation is the phase component, an appropriate constant is used for the amplitude component in this way.
  • may be used instead of 1.0. That is, the amplitude component
  • the frequency obtained by multiplying k by an integer is obtained by multiplying the cross spectrum S 12 (k, n) of the frequency k by a constant.
  • the phase components of the frequencies k, 2k, 3k, and 4k are obtained by multiplying the phase component arg (S 12 (k, n)) of the frequency k by an integer multiple, that is, arg (S 12 (k, n, n)), 2 arg (S 12 (k, n)), 3 arg (S 12 (k, n)), 4 arg (S 12 (k, n)) are used.
  • phase component of the frequency that is a non-constant multiple of the frequency k is set to zero. Accordingly, the phase component arg (U k (w, n)) of the frequency-specific basic cross spectrum corresponding to the frequency k is calculated as follows.
  • p is an integer of 1 or more and P or less.
  • P is an integer greater than 1.
  • the amplitude component and the phase component obtained by the above method are integrated using the above equation (3) to obtain a frequency-based basic cross spectrum U k (w, n).
  • the spectrum for each frequency is obtained after the amplitude component and the phase component are obtained separately.
  • the power of the cross spectrum is used as shown in the following formula, the frequency-specific spectrum U k (w, n) can be obtained without obtaining the amplitude component and the phase component.
  • Integrated correlation function calculation section 204, frequency-cross spectrum calculation unit 203 1, 203 2, ..., the integrated correlation function calculated based on the frequency-cross spectrum supplied from 203 K, the estimated direction-information generating unit 205 introduce.
  • FIG. 2C is a block diagram illustrating a configuration of the integrated correlation function calculation unit 204 included in the wave source direction estimation apparatus 200 according to the present embodiment.
  • the integrated correlation function calculation unit 204 includes frequency-specific correlation function generation units 241 1 , 241 2 ,..., 241 K, and an integration unit 242.
  • the transform unit 201 uses Fourier transform, a method using inverse Fourier transform for inverse transform will be described.
  • the frequency-specific cross spectrum supplied from the frequency-specific cross spectrum calculation unit 203 k is U k (w, n)
  • the frequency-specific correlation function u k ( ⁇ , n) obtained by inverse transformation of U k (w, n).
  • the integration unit 242 integrates the frequency-specific correlation functions supplied from the frequency-specific correlation function generation units 241 1 , 241 2 ,..., 241 K , and transmits them to the estimated direction information generation unit 205 as an integrated correlation function.
  • a single correlation function is obtained by mixing or superposing a plurality of individually obtained correlation functions by frequency.
  • the integration unit 242 calculates the sum of the correlation functions for each frequency. If the integrated correlation function is u ( ⁇ , n), u ( ⁇ , n) is calculated as follows.
  • u ( ⁇ , n) is calculated as follows.
  • the integrated correlation function may be obtained using only the correlation function for each frequency corresponding to the frequency. Moreover, you may control the influence degree of the correlation function classified by frequency in integration in the form of weighting. For example, when a set of frequencies where the target sound exists is ⁇ , u ( ⁇ , n) is calculated as follows when the frequency is selected.
  • u ( ⁇ , n) is calculated as follows.
  • a and b are real numbers and satisfy a>b> 0.
  • a correlation function that is less influenced by non-target sounds such as noise can be generated, so that the direction estimation accuracy is improved.
  • the relative delay time calculation unit 206 obtains the relative delay time between the microphone pairs from the input microphone position information and the sound source search target direction, and transmits it to the estimated direction information generation unit 205 as a set with the sound source search target direction.
  • the relative delay time is a difference in arrival time of sound waves that is uniquely determined based on the microphone interval and the sound source direction. Assuming that the speed of sound is c and the interval between two microphones is d and the direction of the sound source, that is, the direction of sound arrival, is ⁇ , the relative delay time ⁇ ( ⁇ ) with respect to the sound source direction ⁇ is calculated by the following equation.
  • the relative delay time is calculated for all sound source search target directions. For example, when the direction search range is 0 degree to 90 degrees in increments of 10 degrees, that is, 0 degrees, 10 degrees, 20 degrees,..., 90 degrees, 10 types of relative delay times are calculated. Then, the direction to be searched and the relative delay time are supplied to the estimated direction information generation unit 205 as a pair.
  • the estimated direction information generation unit 205 estimates the correspondence between the direction and the correlation value. Output as information. If the correlation function is u ( ⁇ , n) and the relative delay time ⁇ ( ⁇ ), the estimated direction information H ( ⁇ , n) is given by the following equation.
  • the correlation value is determined for each direction, if the correlation value is basically high, it can be determined that there is a high possibility that a sound source exists in that direction.
  • Such estimated direction information is used in various forms. For example, when the function has a plurality of peaks, it is considered that there are a plurality of sound sources having each peak as an arrival direction. Therefore, not only can the direction of each sound source be estimated simultaneously, but it can also be used to estimate the number of sound sources.
  • the possibility of the sound source based on the difference between the peak and non-peak of the correlation function. If the difference between the peak and the non-peak is large, it can be determined that the possibility of the sound source is high. At the same time, it can be determined that the reliability of the estimated direction is high. If the number of sound sources can be assumed to be one in advance, the direction with the maximum correlation value may be output as estimated direction information. In this case, the estimated direction information is not the correspondence between the direction and the correlation value, but the direction itself.
  • the frequency-specific cross spectrum is based on the cross spectrum of a certain frequency k: “the phase component arg (S 12 (k, n)) of the frequency k is added to the phase component of the frequency pk obtained by multiplying the frequency k by an integer. It is defined as "assigned p times". Here, p is an integer of 1 or more. That is, the frequency-specific cross spectrum is defined as the phase component arg (U k (w, n)) satisfying at least the following expression.
  • p 1 and 2 and 3
  • p 2 and 3.
  • the frequency-specific cross spectrum is generated by extracting only the component of frequency k, but the direction estimation accuracy is equivalent to that of the prior art, and high accuracy of direction estimation cannot be achieved.
  • the peak of the correlation function by frequency appears periodically, and the peak interval is inversely proportional to the frequency k. If the frequency k increases, the two adjacent frequency function correlation peaks approach each other, and the peaks become inconspicuous due to the overlap of the correlation functions.
  • FIG. 4 is a diagram illustrating an example of the configuration of the integrated correlation function table 401 included in the wave source direction estimation apparatus 200 according to the present embodiment.
  • the integrated correlation function table 401 stores a frequency domain signal 412, a cross spectrum 413, a frequency-specific cross spectrum 414, and an integrated correlation function 415 in association with the input signal 411.
  • the wave source direction estimation apparatus 200 may calculate an integrated correlation function every time an input signal is acquired, or obtain an integrated correlation function corresponding to the input signal in advance and refer to the integrated correlation function table 401 to determine the integrated correlation function. A function may be calculated.
  • FIG. 5 is a block diagram showing a hardware configuration of the wave source direction estimation apparatus 200 according to the present embodiment.
  • a CPU (Central Processing Unit) 510 is a processor for arithmetic control, and realizes a functional component of the wave source direction estimation apparatus 200 in FIG. 2A by executing a program.
  • a ROM (Read Only Memory) 520 stores fixed data and programs such as initial data and programs.
  • the communication control unit 530 communicates with other devices via a network. Note that the number of CPUs 510 is not limited to one, and may be a plurality of CPUs or may include a GPU (Graphics Processing Unit) for image processing.
  • the communication control unit 530 preferably includes a CPU independent of the CPU 510 and writes or reads transmission / reception data in an area of a RAM (Random Access Memory) 540.
  • the input / output interface 560 preferably has a CPU independent of the CPU 510 and writes or reads input / output data in the RAM 540 area. Therefore, the CPU 510 recognizes that the data has been received or transferred to the RAM 540 and processes the data. Further, the CPU 510 prepares the processing result in the RAM 540 and leaves the subsequent transmission or transfer to the communication control unit 530, the DMAC, or the input / output interface 560.
  • DMAC Direct Memory Access Controller
  • the RAM 540 is a random access memory used by the CPU 510 as a temporary storage work area. In the RAM 540, an area for storing data necessary for realizing the present embodiment is secured.
  • the input signal 541 is sound signal data collected by a sound collection device such as a microphone or signal data input to and obtained by an input signal acquisition device.
  • the frequency domain signal 542 is a signal obtained by converting the input signal 541 by the conversion unit 201.
  • a cross spectrum 543 is a spectrum calculated by the cross spectrum calculation unit 202.
  • the frequency-specific cross spectrum 544 is a spectrum calculated by the frequency-specific cross spectrum calculation unit 203 k .
  • the integrated correlation function 545 is a function calculated by the integrated correlation function calculation unit 204.
  • the input / output data 546 is data input / output via the input / output interface 560.
  • Transmission / reception data 547 is data transmitted / received via the network interface 530.
  • the RAM 540 includes an application execution area 548 for executing various application modules.
  • the storage 550 stores a database, various parameters, or the following data or programs necessary for realizing the present embodiment.
  • the storage 550 stores an integrated correlation function table 401.
  • the integrated correlation function table 401 is a table for managing the relationship between the input signal and the integrated correlation function shown in FIG.
  • the storage 550 further stores a conversion module 551, a cross spectrum calculation module 552, a frequency-specific cross spectrum calculation module 553, and an integrated correlation function calculation module 554.
  • the storage 550 stores an estimated direction information generation module 555 and a relative delay time calculation module 556.
  • the conversion module 551 is a module that converts an input signal into a frequency domain signal.
  • the cross spectrum calculation module 552 is a module that calculates a cross spectrum based on the frequency domain signal.
  • the frequency-specific cross spectrum calculation module 553 is a module that calculates a frequency-specific cross spectrum using the cross spectrum.
  • the integrated correlation function calculation module 554 is a module that calculates an integrated correlation function based on the frequency-specific cross spectrum.
  • the estimated direction information generation module 555 is a module that generates the estimated direction information of the wave source based on the integrated envelope function.
  • the relative delay time calculation module 556 is a module that calculates a relative delay time. These modules 551 to 556 are read by the CPU 510 into the application execution area 548 of the RAM 540 and executed.
  • the control program 557 is a program for controlling the entire wave source direction estimation apparatus 200.
  • the input / output interface 560 interfaces input / output data with input / output devices.
  • a display unit 561 and an operation unit 562 are connected to the input / output interface 560.
  • a storage medium 564 may be further connected to the input / output interface 560.
  • a speaker 563 that is an audio output unit, a microphone that is an audio input unit, or a GPS position determination unit may be connected.
  • the RAM 540 and the storage 550 shown in FIG. 5 do not show programs and data related to general-purpose functions and other realizable functions that the wave source direction estimation apparatus 200 has.
  • FIG. 6 is a flowchart showing a processing procedure of the wave source direction estimation apparatus 200 according to the present embodiment. This flowchart is executed by the CPU 510 in FIG. 5 using the RAM 540, and realizes a functional component of the wave source direction estimation apparatus 200 in FIG.
  • step S601 the wave source direction estimation device 200 acquires an input signal.
  • step S603 the conversion unit 201 of the wave source direction estimation apparatus 200 converts an input signal supplied from the input terminal 20 1 and the input terminal 20 2.
  • the conversion unit 201 supplies the frequency domain signal obtained by the conversion to the cross spectrum calculation unit 202.
  • step S604 the cross spectrum calculation unit 202 calculates a cross spectrum based on the supplied converted signal.
  • the cross spectrum calculation unit 202 transmits the calculated cross spectrum to the frequency-specific cross spectrum calculation units 203 1 , 203 k ,... 203 K.
  • step S607 the frequency-specific cross spectrum calculation units 203 1 , 203 k ,... 203 K calculate a cross spectrum corresponding to each frequency k of the cross spectrum. That is, the frequency-specific cross spectrum calculation units 203 1 , 203 k ,... 203 K calculate the frequency-specific cross spectrum. Then, the frequency-specific cross spectrum calculation units 203 1 , 203 k ,... 203 K transmit the frequency-specific cross spectrum to the integrated correlation function calculation unit 204.
  • step S609 the frequency-specific correlation function generators 241 1 , 241 2 ,..., 241 K inversely transform the frequency-specific cross spectrum to calculate the frequency-specific correlation function.
  • step S611 the integration unit 242 integrates the frequency-specific correlation functions and calculates an integrated correlation function.
  • step S613 the relative delay time calculation unit 206 calculates the relative delay time between the microphone pairs from the microphone position information and the sound source search target direction.
  • step S615 the estimated direction information generation unit 205 generates estimated direction information from the integrated correlation function and the relative delay time.
  • the arrival direction of the target sound included in the input signal that is, the direction in which the target object exists is estimated. This is effective in estimating the direction in which the target exists in the environment where the environmental noise level is high, using the sound generated by the target as a clue.
  • environmental noise include busy streets, streets, along highways, and places where many people and cars gather.
  • the target include humans, animals, automobiles, aircraft, ships, water bikes, and drones (small drones).
  • the position of the target sound source can be specified by performing sound source direction estimation at a plurality of locations. As a result, even in an environment with a high environmental noise level, it is possible to accurately specify a scream, a gunshot, a location where an automobile collision sound is generated, or the like.
  • FIG. 7 is a block diagram for explaining the configuration of the integrated correlation function generation unit 704 provided in the wave source direction estimation apparatus according to the present embodiment.
  • the integrated correlation function generation unit 704 included in the wave source direction estimation apparatus according to the present embodiment has frequency-specific correlation function generation units 241 1 , 241 2 ,. , 241 K and the integration unit 242 are replaced with an integration unit 741 and an integrated correlation function generation unit 742. Since other configurations and operations are the same as those of the second embodiment, the same configurations and operations are denoted by the same reference numerals, and detailed description thereof is omitted.
  • the integration unit 741 integrates the frequency-specific cross spectra supplied from the frequency-specific cross spectrum calculation units 203 1 , 203 2 ,..., 203 K , and transmits them to the integrated correlation function generation unit 742 as an integrated cross spectrum.
  • a single integrated cross spectrum is obtained by mixing or overlapping a plurality of frequency-specific cross spectra obtained individually.
  • a summation or a sum of powers is used as in the integration unit 242 of the second embodiment.
  • the integrated cross spectrum U (k, n) is calculated as follows.
  • the integrated cross spectrum U (k, n) is generated. You may correct to. Similar to the second embodiment, the degree of influence is controlled in the form of frequency selection and weighting. For example, when the set of frequencies in which the target sound exists is ⁇ , the calculation is performed as follows when the integrated cross spectrum U (k, n) is obtained by selecting a band.
  • U (k, n) is calculated as follows.
  • a and b are real numbers and satisfy a>b> 0.
  • a correlation function that is less influenced by non-target sounds such as noise can be generated, so that the direction estimation accuracy is improved.
  • the integrated correlation function generation unit 742 performs inverse transformation of the integrated cross spectrum supplied from the integration unit 741, and transmits the result to the estimated direction information generation unit 205 as an integrated correlation function. Also in this embodiment, a method using inverse Fourier transform for inverse transform will be described.
  • the integrated cross spectrum supplied from the integration unit 741 is U (k, n)
  • the integrated correlation function u ( ⁇ , n) obtained by the inverse transformation of U (k, n) is calculated as follows.
  • the cross correlation by frequency is integrated and then inverse transformation is performed to obtain an integrated correlation function.
  • count of reverse conversion decreases compared with 2nd Embodiment which performed reverse conversion for every cross spectrum according to frequency. Therefore, the integrated correlation function can be obtained with a smaller calculation amount than in the second embodiment.
  • FIG. 8A is a block diagram for explaining the configuration of the wave source direction estimation apparatus 800 according to this embodiment.
  • the wave source direction estimation apparatus 800 according to the present embodiment replaces the frequency-specific cross spectrum calculation units 203 1 , 203 2 ,..., 203 K with a frequency-specific cross spectrum calculation unit 803. 1 , 803 2 ,..., 803 K. Since other configurations and operations are the same as those in the first embodiment, the same configurations and operations are denoted by the same reference numerals, and detailed description thereof is omitted.
  • FIG. 8B is a block diagram of the frequency-specific cross spectrum calculation unit 803 k .
  • the frequency-specific cross spectrum calculation unit 803 k includes a frequency-specific basic cross spectrum calculation unit 2031 k , a kernel function spectrum storage unit 831, and a multiplication unit 832.
  • the frequency-based basic cross spectrum calculation unit 2031 k uses the cross spectrum S 12 (k, n) supplied from the cross spectrum calculation unit 202 to calculate a cross spectrum corresponding to the frequency k of S 12 (k, n). Then, it is transmitted to the multiplier 832 as a frequency-specific basic cross spectrum. Operation of the frequency-specific basic cross-spectrum calculation unit 2031 k, since with the exception of its output destination is the same as the frequency-specific basic cross-spectrum calculation unit 2031 k of the second embodiment, detailed description is omitted.
  • the kernel function spectrum storage unit 831 stores a kernel function spectrum, and outputs the kernel function spectrum to the multiplication unit 832.
  • the kernel function spectrum is obtained by Fourier transforming a kernel function and taking its absolute value. You may square instead of taking an absolute value.
  • a Gaussian function is used as the kernel function.
  • the Gaussian function is given by the following mathematical expression using three real numbers g 1 , g 2 , and g 3 given in advance.
  • g 1 is the magnitude of the Gaussian function
  • g 2 is the peak position of the Gaussian function
  • g 3 is the spread of the Gaussian function.
  • g 3 for adjusting the spread of the Gaussian function is important because it greatly affects the sharpness of the peak of the correlation function for each frequency. As seen from equation (21), spread of the Gaussian function is greater the greater the g 3.
  • g 1 and g 2 are real numbers.
  • the logistic function 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 logistic function is an important parameter that greatly affects the sharpness of the peak of the correlation function for each frequency, as in the case of g 3 in the Gaussian function.
  • a cosine function or a uniform function may be used.
  • the parameters g 1 to g 5 used in the kernel function are not constants and may be different values depending on the frequency k. That is, it may be a function of the frequency k as g 1 (k) to g 5 (k).
  • g 3 is a function g 3 (k) of the frequency k, and a function having a smaller value as the frequency becomes higher.
  • g 3 (k) is given as follows.
  • G 3 is a real number.
  • the kernel function G (k) becomes a function with a sharp peak and a narrow tail as the frequency k increases.
  • the multiplication unit 832 calculates the product of the frequency-specific basic cross spectrum supplied from the frequency-specific basic cross spectrum calculation unit 2031 k and the kernel function spectrum supplied from the kernel function spectrum storage unit 831 and integrates the product as a frequency-specific cross spectrum. This is transmitted to the correlation function calculation unit 204.
  • the frequency-basis fundamental spectrum supplied from the frequency-basis fundamental spectrum calculation unit 2031 k is U k (w, n) and the kernel function spectrum supplied from the kernel function spectrum storage unit 831 is G (w)
  • the frequency-specific cross spectrum The spectrum UM k (w, n) is calculated as follows.
  • FIG. 9 shows the relationship between the frequency-specific cross spectrum multiplied by the kernel function spectrum and the frequency-specific correlation function.
  • the cross spectrum by frequency before multiplication by the kernel function spectrum is also shown.
  • the kernel function spectrum is not multiplied, components exist up to a high frequency, so that the peak of the frequency-specific correlation function becomes sharp.
  • the kernel function spectrum is multiplied as shown in the center diagram and the right diagram in FIG.
  • the high frequency components are attenuated, so that the peak sharpness of the correlation function for each frequency is reduced. That is, as the peak of the kernel function spectrum becomes sharper (the bottom of the kernel function spectrum becomes narrower), the sharpness of the peak of the correlation function by frequency becomes smaller. Further, as shown in the right side of FIG. 9, when the skirt of the correlation function for each frequency is widened, the skirts of adjacent peaks overlap and a correlation function for each frequency having a shallow valley is obtained.
  • FIG. 10 is a diagram showing the relationship between the presence / absence of the kernel function and the integrated correlation function.
  • the peak positions of the frequency-specific correlation functions u 1 ( ⁇ , n) to u 3 ( ⁇ , n) are close, but u 1 ( ⁇ , n) to Since the width of u 3 ( ⁇ , n) is narrow, a large peak cannot be formed during integration. For this reason, the position of the peak is not clear.
  • the kernel function as shown in FIG.
  • the width of the correlation function for each frequency is wide, so that u 1 ( ⁇ , n) to u 3 ( ⁇ , n) are large due to integration.
  • a peak can be formed. For this reason, the position of the peak is clearer than in the case of (a) without the kernel function.
  • FIG. 11 is a diagram showing the relationship between the difference in the kernel function spectrum width and the integrated correlation function.
  • a correlation function for each frequency having a shallow valley is formed due to the periodicity of the correlation function. Therefore, as shown in FIG. 11C, when the correlation functions by frequency with shallow valleys are integrated, an integrated correlation function with shallow valleys, that is, inconspicuous peaks, is generated.
  • FIG. 11D when a narrow kernel function spectrum is used, a correlation function by frequency having a deeper valley than the right diagram of FIG. 9 is formed. Therefore, as shown in FIG. 11D, an integrated correlation function with a clear peak is generated.
  • the product of the kernel function spectrum obtained by the Fourier transform of the kernel function and the basic cross spectrum for each frequency is calculated, but it can also be realized in the time domain due to the nature of the Fourier transform.
  • a “convolution operation unit” that convolves the kernel function in the subsequent stage of the frequency-specific correlation function generation unit 241 k in the integrated correlation function calculation unit 204 is provided.
  • the kernel function may be convolved with the correlation function for each frequency supplied from 241 k .
  • the convolution operation requires a large amount of calculation, it is more efficient to calculate the product in the frequency domain as in this embodiment.
  • the frequency-specific cross spectrum is generated by multiplying the frequency-specific basic cross spectrum by the kernel function spectrum. For this reason, the width of the correlation function for each frequency obtained by the inverse transformation becomes wide, and the peak of the integrated correlation function becomes clear. In particular, when the peak positions of the individual frequency-related correlation functions are close to each other and each has a sharp peak, the effect of clarifying the peak of the integrated correlation function is enhanced by performing correction.
  • FIG. 12 is a diagram for explaining a configuration of a frequency-specific cross spectrum calculation unit included in the wave source direction estimation apparatus according to the present embodiment.
  • the frequency-specific cross spectrum calculation unit 1203 k included in the wave source direction estimation apparatus according to the present embodiment has a kernel function spectrum generation unit 1231 instead of the kernel function spectrum storage unit 831.
  • the same configurations and operations are denoted by the same reference numerals, and detailed description thereof is omitted.
  • the kernel function spectrum generation unit 1231 generates a kernel function spectrum using the cross spectrum supplied from the cross spectrum calculation unit 202, and transmits the generated kernel function spectrum to the multiplication unit 832.
  • the kernel function spectrum generation unit 1231 analyzes the supplied cross spectrum to obtain the possibility that the target sound exists in the input signal, and generates a kernel function spectrum having a shape reflecting the possibility of existence. Basically, if the possibility of existence is low, a kernel function spectrum having a narrow width and a small spread is generated. Thereby, since the peak of the correlation function classified by frequency becomes low, the possibility that an erroneous peak appears in the integrated correlation function can be reduced.
  • a method for estimating the SNR (Signal-to-Noise Ratio) of the input signal will be described as a method for obtaining the possibility of the target sound.
  • the absolute value of the supplied cross spectrum is calculated.
  • the square of the Fourier transform obtained by the conversion unit 201 is referred to as an input signal power spectrum, but in this embodiment, the absolute value of the cross spectrum is handled as the input signal power spectrum.
  • the power spectrum of the noise component (non-target sound component) included in the input signal is estimated based on the input signal power spectrum. Assuming that the input signal power spectrum is P X (k, n), P X (k, n) is calculated as follows.
  • the power spectrum of the noise component is estimated based on the input signal power spectrum.
  • the method described in Non-Patent Document 3 is used.
  • the estimated noise power spectrum is an average of the power spectrum at the initial stage of estimation when the input signal power spectrum is supplied. In this case, it is necessary to satisfy the condition that the target sound is not included immediately after the estimation is started.
  • the estimated noise power spectrum is P N (k, n)
  • P N (k, n) is calculated as follows.
  • N 0 is a predetermined integer.
  • Non-Patent Document 4 discloses a method for obtaining an estimated noise power spectrum from a minimum value (minimum statistic) of an input signal power spectrum.
  • the minimum value of the input signal power spectrum within a fixed time is held for each frequency, and the noise component is estimated from the minimum value. Since the minimum value of the input signal power spectrum is similar to the noise power spectrum and the spectrum shape, it can be used as an estimate of the noise power spectrum.
  • the estimated noise power spectrum is obtained, the ratio with the input signal power spectrum is taken as the estimated SNR value.
  • the estimated SNR ⁇ (k, n) is calculated as follows.
  • this estimated SNR ⁇ (k, n) is used as it is for the existence possibility q (k, n) of the target sound.
  • the estimated SNR obtained in this way is called an estimated a-posteriori SNR in Non-Patent Document 3.
  • an estimated a-priori SNR obtained by the method described in Non-Patent Document 3 may be used instead of the estimated a-posteriori SNR.
  • the SNR is estimated after suppressing the noise component. Therefore, although the calculation amount increases, higher estimation accuracy than the a-posteriori SNR can be achieved.
  • the method for calculating the presence possibility of the target sound using the input signal power spectrum and the estimated noise power spectrum is not limited to the ratio of both, such as the estimated SNR.
  • the difference between the two may be used instead of the ratio.
  • a simple magnitude relationship may be used.
  • the method for obtaining the possibility of the target sound by analyzing the cross spectrum is not limited to the method using the power spectrum.
  • Another typical example is a method of analyzing the phase component of the cross spectrum.
  • a method for analyzing the phase component a method using a cross spectrum group delay (a phase component differentiated in the frequency direction) will be described.
  • the group delay of the cross spectrum is obtained.
  • the group delay is gd (k, n)
  • the group delay of the cross spectrum S 12 (k, n) can be calculated as follows.
  • the average value of gd (k, n) is calculated, and the degree of divergence from the average value is assumed to exist.
  • the possibility of existence q (k, n) can be calculated as follows.
  • q 0 is a positive real number.
  • the gd (k, n) bar is obtained by averaging gd (k, n) in the frequency direction. There are various methods for averaging. For example, the following arithmetic average can be used.
  • q (k, n) approaches 1 if gd (k, n) is close to gd (k, n) bar. On the other hand, q (k, n) approaches 0 as gd (k, n) moves away from the gd (k, n) bar.
  • a kernel function spectrum is generated using the obtained existence possibility.
  • an example of controlling the parameters of the kernel function that is the basis of the kernel function spectrum will be described.
  • An example of using a Gaussian function as a kernel function will be described.
  • At high possible presence of the target sound is to reduce the g 3.
  • the higher the possibility of existence the narrower the width of g ( ⁇ ), and the closer to the shape in which the g ( ⁇ ) peak is emphasized.
  • a linear function with the reciprocal of the possibility of existence as a variable is used. At this time, assuming that the existence possibility is q (k, n), g 3 is calculated as follows.
  • a 1 and b 1 are real numbers and satisfy a 1 > 0.0 and b 1 > 0.0.
  • the function for obtaining g 3 from the existence possibility q (k, n) of the target sound is not limited to a linear function. It is also possible to use functions expressed in other forms such as sigmoid functions, higher-order polynomial functions, and non-linear functions instead of linear functions.
  • g 5 may be calculated by the same method as g 3 .
  • the width of the kernel function g (tau) is narrow, approaches the shape of the peak is emphasized.
  • the possibility of the target sound is obtained, and the kernel function is calculated based on the possibility.
  • the probability is high, the width of the kernel function spectrum widens and the shape approaches a flat one.
  • the possibility is low, the width of the kernel function spectrum becomes narrow.
  • the peak of the correlation function classified by frequency of the frequency where the target sound exists becomes high, and the peak of the correlation function classified by frequency of the frequency where the target sound does not exist becomes low.
  • the peak of the integrated correlation function is emphasized more than in the fourth embodiment, and the direction estimation accuracy of the target sound is improved.
  • the possibility that an erroneous peak appears in the integrated correlation function can be reduced.
  • FIG. 13A is a block diagram for explaining a configuration of a wave source direction estimation apparatus 1300 according to the present embodiment.
  • the wave source direction estimation apparatus 1300 according to the present embodiment is different from the third embodiment in that an integrated correlation function calculation unit 1304 is provided instead of the integrated correlation function calculation unit 204. Since other configurations and operations are the same as those of the third embodiment, the same configurations and operations are denoted by the same reference numerals, and detailed description thereof is omitted.
  • FIG. 13B is a block diagram for explaining a configuration of a frequency-specific cross spectrum included in the wave source direction estimation device according to the sixth embodiment of the present invention.
  • the frequency-specific cross spectrum calculation unit 203 k according to the present embodiment is different from the third embodiment in that an integrated cross spectrum generation unit 1341 is provided instead of the integration unit 741. Since other configurations and operations are the same as those of the third embodiment, the same configurations and operations are denoted by the same reference numerals, and detailed description thereof is omitted.
  • Integrated cross spectrum generation unit 1341 based on the cross spectrum supplied from the cross spectrum calculation unit 202, integrated frequency-cross spectrum calculation unit 203 1, 203 2, ..., a frequency-cross spectrum supplied from 203 K And transmitted to the integrated correlation function generation unit 742 as an integrated cross spectrum.
  • the third embodiment the case where the frequency where the target sound exists or the frequency where the power of the target sound is large is known in advance has been described.
  • the possibility that the target sound exists in the input signal is obtained, and integration is performed based on the existence possibility.
  • the possibility of the target sound is obtained based on the supplied cross spectrum.
  • the method described in the fifth embodiment can be similarly used for calculating the possibility of existence.
  • frequency cross spectrum integration is performed using the obtained existence possibility.
  • the possibility of the target sound is q (k, n)
  • a set ⁇ of frequencies with high possibility of the target sound being found is obtained based on q (k, n). If q (k, n) for a certain frequency k exceeds a predetermined threshold value ⁇ q , that frequency is made an element of the set ⁇ . This can be expressed as a mathematical expression as follows.
  • the weight may be calculated using the existence possibility q (k, n), and the integration by the weighted sum may be performed using the weight. If the weight function is ⁇ (q (k, n)), the integrated cross spectrum U (k, n) is calculated as follows.
  • the weighting function ⁇ (q (w, n)) is a monotonically increasing function that takes a large value for a large q (w, n).
  • the integrated cross spectrum is calculated using the existence possibility. For this reason, even when the possibility of the target sound existing in advance is unknown, band selection and weighting at the time of generating the integrated cross spectrum can be performed appropriately, and high estimation accuracy can be achieved.
  • FIG. 14 is a diagram for explaining a configuration of a wave source direction estimation system 1400 according to the present embodiment.
  • the wave source direction estimation system 1400 according to the present embodiment uses the wave source direction estimation apparatus 200 according to the second embodiment. Therefore, the same configurations and operations as those of the second embodiment are denoted by the same reference numerals, and detailed description thereof is omitted.
  • the wave source direction estimation system 1400 includes a microphone 140 1 , a microphone 140 2 , an AD conversion unit 1401, and a display unit 1402.
  • a microphone 140 1 it is possible to use the wave source direction estimation device 800 or the wave source direction estimation device 1300 instead of the wave source direction estimation device 200.
  • an example using a microphone will be described.
  • various sensors that can receive a wave radiated from the wave source and convert it into an electrical signal are available. Used in place of a microphone.
  • Microphone 140 1 and microphone 140 2 converts the sound of peripheral devices, including a sound generated from the target object is estimated subject to an electrical signal and transmitted to the AD conversion unit 1401.
  • the medium through which the sound is transmitted is an air medium
  • the sound reaches the microphone as air vibration.
  • the microphone converts the vibration of the air that has arrived into an electrical signal.
  • the AD conversion unit 1401 converts the electrical signal of the sound supplied from the microphone 140 1 and the microphone 140 2 into a digital signal and transmits the digital signal to the input terminal 20 1 and the input terminal 20 2 .
  • the display unit 1402 converts the estimated direction information supplied from the wave source direction estimating apparatus 200 into visualization data such as an image, and displays it on a display apparatus such as a display.
  • the most basic visualization method is a method of displaying a correlation function at a certain time in a two-dimensional graph. At that time, the direction is displayed on the horizontal axis and the correlation value is displayed on the vertical axis. It is also effective to display not only a certain time but also a temporal change of the correlation function in three dimensions. By displaying the time change, it becomes possible to clarify the appearance of the target sound source, to predict the movement pattern of the target sound source, and to predict the movement direction of the target sound source.
  • a method of projecting onto a two-dimensional plane instead of three-dimensional is also effective. There is a problem that it is difficult to see the back side when displayed in 3D. If displayed on a plane projected from above, the blind spot disappears and the listability is improved.
  • the correlation value may be expressed by contour lines instead of color shading.
  • FIG. 15 is a diagram illustrating an example of an image displayed on the display unit 1402 of the wave source direction estimation system 1400 according to the present embodiment, which is obtained from the estimated direction information supplied from the wave source direction estimation apparatus 200. This was acquired for the purpose of confirming the effect of this embodiment.
  • FIG. 15 shows that the correlation value is higher as the color is blacker.
  • the range of azimuth is 0 to 180 degrees.
  • the vertical axis represents time. Referring to FIG. 15, it can be seen that the correlation value of the azimuth of 30 degrees is high around the time from 20 seconds to 25 seconds. From this, it can be seen that a scream occurred at 20 to 25 seconds and the direction of scream occurrence was about 30 degrees.
  • the user can visually grasp the direction estimation information of the wave source.
  • the present invention may be applied to a system composed of a plurality of devices, or may be applied to a single device. Furthermore, the present invention can also be applied to a case where an information processing program that implements the functions of the embodiments is supplied directly or remotely to a system or apparatus. Therefore, in order to realize the functions of the present invention on a computer, a program installed on the computer, a medium storing the program, and a WWW (World Wide Web) server that downloads the program are also included in the scope of the present invention. . In particular, at least a non-transitory computer readable medium storing a program for causing a computer to execute the processing steps included in the above-described embodiments is included in the scope of the present invention.
  • a plurality of input signal acquisition means for acquiring a wave generated by the wave source as an input signal; Conversion means for converting the plurality of input signals acquired by the input signal acquisition means into a plurality of frequency domain signals; Cross spectrum calculation means for calculating a cross spectrum based on the frequency domain signal; A frequency-specific cross-spectrum calculating means for calculating a frequency-specific cross-spectrum based on the cross-spectrum; An integrated correlation function calculating means for calculating an integrated correlation function based on the frequency-specific cross spectrum; A correlation function generation device.
  • the integrated correlation function calculating means includes A frequency-specific correlation function generating means for generating a frequency-specific correlation function by inverse transformation of the frequency-specific cross spectrum; An integrated correlation function generating means for generating a single integrated correlation function by integrating the frequency-specific correlation functions; The correlation function generator according to appendix 1, wherein (Appendix 3)
  • the integrated correlation function calculating means includes An integrated cross spectrum generating means for generating an integrated cross spectrum by integrating the cross spectra by frequency; Integrated correlation function generating means for generating an integrated correlation function by inverse transformation of the integrated cross spectrum;
  • the frequency-specific cross spectrum calculation means includes: A frequency-basis basic cross spectrum calculating means for calculating a frequency-basic basic cross spectrum based on the cross spectrum; 4.
  • the frequency-specific cross spectrum calculation means includes: A frequency-basis basic cross spectrum calculating means for calculating a frequency-basic basic cross spectrum based on the cross spectrum; Kernel function storage means for storing a kernel function spectrum; Multiplying means for multiplying the fundamental cross spectrum by frequency and the kernel function spectrum to obtain the cross spectrum by frequency;
  • the correlation function generation device includes: A frequency-basis basic cross spectrum calculating means for calculating a frequency-basic basic cross spectrum based on the cross spectrum; Kernel function spectrum calculating means for calculating a kernel function spectrum based on the cross spectrum; Multiplying means for multiplying the fundamental cross spectrum by frequency and the kernel function spectrum to obtain the cross spectrum by frequency;
  • the correlation function generation device according to any one of appendices 1 to 3, wherein, Appendix 6)
  • the frequency-specific cross spectrum calculation means includes: A frequency-basis basic cross spectrum calculating means for calculating a frequency-basic basic cross spectrum based on the cross spectrum; Kernel function spectrum calculating means for calculating a kernel

Abstract

周囲のノイズレベルが高い環境の場合であっても、明瞭なピークを持つ相関関数を生成する。相関関数生成装置であって、波源で発生した波を入力信号として取得する複数の入力信号取得部と、入力信号取得部で取得した複数の入力信号を複数の周波数領域信号に変換する変換部と、周波数領域信号に基づいて、クロススペクトルを計算するクロススペクトル計算部と、クロススペクトルに基づいて、周波数別クロススペクトルを計算する周波数別クロススペクトル計算部と、周波数別クロススペクトルに基づいて、統合相関関数を計算する統合相関関数計算部と、を備えた。

Description

相関関数生成装置、相関関数生成方法、相関関数生成プログラムおよび波源方向推定装置
　本発明は、相関関数生成装置、相関関数生成方法、相関関数生成プログラムおよび波源方向推定装置に関する。
　上記技術分野において、非特許文献１および非特許文献２には、２つのマイクロフォンの受音信号を用いて、音源(音波の発生源、発生場所)の方向を推定する方法が記載されている。具体的には、２つの受音信号から、受音信号間の相互相関関数を求める。そして、相互相関関数が最大値を与える時間差を音波の到来時間差として算出することで、音波の到来方向を推定する技術が開示されている。
C.H.Knapp and G.C.Carter, "The generalized correlation method for estimation of time delay," IEEE Trans. Acoustics, Speech and Signal Processing, vol.24, no. 4, pp. 320-327, Aug. 1976. J.P. Ianniello, "Time delay estimation via cross-correlation in the presence of large estimation errors," IEEE Trans. Acoustics, Speech and Signal Processing, vol.30, no. 6, pp. 998-1003, Dec. 1982. Y. Ephraim and D. Malah, "Speech enhancement using a minimum mean-square error short-time spectral amplitude estimator," IEEE Trans. Acoust., Speech, Signal Processing, vol. ASSP-32, no. 6, pp. 1109-1121, Dec. 1984. R. Martin, "Spectral subtraction based on minimum statistics," Proc. of EUSPICO-94, pp.1182-1185, Sept. 1994
　しかしながら、上記文献に記載の技術では、周囲のノイズレベルが高い環境の場合には、明瞭なピークを持つ相関関数が生成できなかった。また、波源の方向を高精度に推定することができなかった。
　本発明の目的は、上述の課題を解決する技術を提供することにある。
　上記目的を達成するため、本発明に係る相関関数生成装置は、
　波源で発生した波を入力信号として取得する複数の入力信号取得手段と、
　前記入力信号取得手段で取得した複数の前記入力信号を複数の周波数領域信号に変換する変換手段と、
　前記周波数領域信号に基づいて、クロススペクトルを計算するクロススペクトル計算手段と、
　前記クロススペクトルに基づいて、周波数別クロススペクトルを計算する周波数別クロススペクトル計算手段と、
　前記周波数別クロススペクトルに基づいて、統合相関関数を計算する統合相関関数計算手段と、
　を備えた。
　上記目的を達成するため、本発明に係る相関関数生成方法は、
　波源で発生した波を入力信号として取得する複数の入力信号取得ステップと、
　前記入力信号取得ステップにおいて取得した複数の前記入力信号を複数の周波数領域信号に変換する変換ステップと、
　前記周波数領域信号に基づいて、クロススペクトルを計算するクロススペクトル計算ステップと、
　前記クロススペクトルに基づいて、周波数別クロススペクトルを計算する周波数別クロススペクトル計算ステップと、
　前記周波数別クロススペクトルに基づいて、統合相関関数を計算する統合相関関数計算ステップと、
　を含む。
　上記目的を達成するため、本発明に係る相関関数生成プログラムは、
　波源で発生した波を入力信号として取得する複数の入力信号取得ステップと、
　前記入力信号取得ステップにおいて取得した複数の前記入力信号を複数の周波数領域信号に変換する変換ステップと、
　前記周波数領域信号に基づいて、クロススペクトルを計算するクロススペクトル計算ステップと、
　前記クロススペクトルに基づいて、周波数別クロススペクトルを計算する周波数別クロススペクトル計算ステップと、
　前記周波数別クロススペクトルに基づいて、統合相関関数を計算する統合相関関数計算ステップと、
　をコンピュータに実行させる。
　上記目的を達成するため、本発明に係る波源方向推定装置は、
　上記相関関数生成装置と、
　統合相関関数に基づいて、波源の推定方向情報を生成する推定方向情報生成手段と、
　を備えた。
　本発明によれば、周囲のノイズレベルが高い環境の場合であっても、明瞭なピークを持つ相関関数を生成できる。また、波源の方向を高精度に推定することができる。
本発明の第１実施形態に係る情報処理装置の構成を示すブロック図である。 本発明の第２実施形態に係る波源方向推定装置の構成を示すブロック図である。 本発明の第２実施形態に係る波源方向推定装置の備える周波数別クロススペクトル計算部の構成を示すブロック図である。 本発明の第２実施形態に係る波源方向推定装置の備える統合相関関数計算部の構成を示すブロック図である。 本発明の第２実施形態に係る波源方向推定装置により得られた周波数別相関関数の一例について説明する図である。 本発明の第２実施形態に係る波源方向推定装置により得られた周波数別相関関数を統合した統合相関関数の一例について説明する図である。 本発明の第２実施形態に係る波源方向推定装置の備える統合相関関数テーブルの構成の一例を示す図である。 本発明の第２実施形態に係る波源方向推定装置のハードウェア構成を示すブロック図である。 本発明の第２実施形態に係る波源方向推定装置の処理手順を説明するフローチャートである。 本発明の第３実施形態に係る波源方向推定装置の備える統合相関関数生成部の構成を示すブロック図である。 本発明の第４実施形態に係る波源方向推定装置の構成を示すブロック図である。 本発明の第４実施形態に係る波源方向推定装置の備える周波数別クロススペクトル計算部の構成を示すブロック図である。 本発明の第４実施形態に係る波源方向推定装置の周波数別クロススペクトル計算部における、カーネル関数スペクトルを乗じた周波数別クロススペクトルと、周波数別相関関数との関係を示す図である。 本発明の第４実施形態に係る波源方向推定装置の周波数別クロススペクトル計算部における、カーネル関数による周波数別相関関数の大きさ制御の効果について説明する図である。 本発明の第４実施形態に係る波源方向推定装置の周波数別クロススペクトル計算部における、カーネル関数スペクトル幅の違いと、統合相関関数との関係を示す図である。 本発明の第５実施形態に係る波源方向推定装置の備える周波数別クロススペクトル計算部の構成を説明するための図である。 本発明の第６実施形態に係る波源方向推定装置の構成を説明するためのブロック図である。 本発明の第６実施形態に係る波源方向推定装置の備える周波数別クロススペクトルの構成を説明するためのブロック図である。 本発明の第７実施形態に係る波源方向推定システムの構成を説明するための図である。 本発明の第７実施形態に係る波源方向推定システムの表示部に表示される画像の一例を示す図である。
　以下に、本発明を実施するための形態について、図面を参照して、例示的に詳しく説明記載する。ただし、以下の実施の形態に記載されている、構成、数値、処理の流れ、機能要素などは一例に過ぎず、その変形や変更は自由であって、本発明の技術範囲を以下の記載に限定する趣旨のものではない。
　また、以下の実施形態に係る波源方向推定装置の推定対象は、空気や水の振動波である音波の発生源に限定されない。地震や地滑りなどの土や固体を媒質とする振動波の発生源にも適用できる。その場合、振動波を電気信号に変換する装置には、マイクロフォンではなく振動センサが用いられる。さらに、気体・液体・固体の振動波だけでなく、電波を用いて方向を推定する場合にも以下の実施形態に係る波源方向推定装置を適用できる。その場合、電波を電気信号に変換する装置にはアンテナが用いられる。以下の実施形態においては、波源は音源と仮定して説明する。
　［第１実施形態］
　本発明の第１実施形態としての相関関数生成装置１００について、図１を用いて説明する。相関関数生成装置１００は、入力信号に基づいて相関関数を生成する装置である。
　図１に示すように、相関関数生成装置１００は、入力信号取得部１０１と、変換部１０２と、クロススペクトル計算部１０３と、周波数別クロススペクトル計算部１０４と、統合相関関数計算部１０５とを含む。
　複数の入力信号取得部１０１は、波源で発生した波を入力信号として取得する。変換部１０２は、入力信号取得手段で取得した複数の前記入力信号を複数の周波数領域信号に変換する。クロススペクトル計算部１０３は、周波数領域信号に基づいて、クロススペクトルを計算する。周波数別クロススペクトル計算部１０４は、クロススペクトルに基づいて、周波数別クロススペクトルを計算する。統合相関関数計算部１０５は、周波数別クロススペクトルに基づいて、統合相関関数を計算する。
　本実施形態によれば、周囲のノイズレベルが高い環境の場合であっても、明瞭なピークを持つ相関関数を生成できる。また、波源の方向を高精度に推定することができる。
　［第２実施形態］
　次に本発明の第２実施形態に係る波源方向推定装置について、図２Ａ乃至図６を用いて説明する。
　＜前提技術＞
　上記非特許文献１および非特許文献２に記載の技術では、屋外などの周囲騒音レベルが高い環境において、遠方に存在する音源の方向を高精度に推定することが困難であった。例えば、推定対象の音源（目標音源）が、マイクロフォンから遠く離れた場所に存在する場合、目標音源から放射される音の音量は、マイクロフォンに到達する時点で大幅に小さくなる。このため、目標音源の音が周囲環境雑音に埋もれて、明瞭なピークを持つ相関関数を生成できなかった。このため、目標音源の方向推定精度が低下することがあった。
　＜本実施形態の技術＞図２Ａは、本実施形態に係る波源方向推定装置の構成を示すブロック図である。図２Ｂは、本実施形態に係る波源方向推定装置の備える統合相関関数計算部の構成を示すブロック図である。
　本実施形態に係る波源方向推定装置２００は、例えば、デジタルビデオカメラ、スマートフォン、携帯電話、ノートパソコン、パッシブソーナーなどといった装置の一部として機能する。また、不審ドローン検知、悲鳴検知、車両事故検知などの声や音に基づいて異常を検知する異常音検知装置にも搭載される。しかし、本実施形態に係る波源方向推定装置２００の適用例はこれらに限定されるものではなく、受信音から目標音源の方向推定を要求されるあらゆる波源方向推定装置に適用可能である。
　波源方向推定装置２００は、入力端子２０１と、入力端子２０２と、変換部２０１と、クロススペクトル計算部２０２と、周波数別クロススペクトル計算部２０３１～２０３ｋと、を備える。波源方向推定装置２００は、さらに、統合相関関数計算部２０４と、推定方向情報生成部２０５と、相対遅延時間計算部２０６とを備える。
　入力端子２０１と入力端子２０２とには、目標音源の音と、集音装置であるマイクロフォン(以下、マイク)の周囲で生じる様々な雑音が混在した音とがデジタル信号(サンプル値系列)として入力される。入力端子２０１と入力端子２０２とに入力された音信号を本実施形態では入力信号と呼ぶ。そして、時刻ｔにおける入力端子２０１の入力信号をｘ１（ｔ）と、入力端子２０２の入力信号をｘ２（ｔ）と表す。
　入力端子に入力される音は、集音装置であるマイクで集音される。入力端子は複数存在するので、目標音源の音を集音する場合には、端子数と同じ２個のマイクが同時に使用される。本実施形態では、入力端子とマイクとは一対一に対応しており、ｍ番目のマイクが集音した音は、ｍ番目の入力端子に供給されると仮定する。したがって、ｍ番目の入力端子に供給された入力信号のことを「ｍ番目のマイクの入力信号」とも呼ぶ。
　波源方向推定装置２００は、目標音源の音が２つのマイクに到達する時間差を利用して音源の方向を推定する。このため、マイク間隔も重要な情報となるので、波源方向推定装置２００には、入力信号だけでなく、マイク位置情報も供給される。
　変換部２０１は、入力端子２０１と入力端子２０２とから供給された入力信号を変換し、クロススペクトル計算部２０２へ供給する。変換は、入力信号を複数の周波数成分に分解する目的で実施される。ここでは、代表的なフーリエ変換を用いた場合について説明する。
　変換部２０１には２種類の入力信号ｘｍ（ｔ）が入力される。ここで、ｍは入力端子番号である。変換部２０１は、入力端子から供給された入力信号から、適当な長さの波形を一定の周期でずらしながら切り出す。こうして切り出した信号区間をフレーム、切り出した波形の長さをフレーム長、フレームをずらす周期をフレーム周期と呼ぶ。そして、フーリエ変換を用いて切り出された信号を周波数領域信号に変換する。ｎをフレーム番号とし、切り出した入力信号ｘｍ（ｔ，ｎ）（ｔ＝０，１，・・・，Ｋ－１）とすると、ｘｍ（ｔ，ｎ）のフーリエ変換Ｘｍ（ｋ，ｎ）は次のように計算される。
Figure JPOXMLDOC01-appb-M000001
 　ここで、ｊは虚数単位（－１の平方根）、ｅｘｐは指数関数を表す。また、ｋは周波数ビン番号を表し、０以上Ｋ－１以下の整数である。以下、簡略化のため、ｋのことを周波数ビン番号ではなく、単に「周波数」と呼ぶ。
　クロススペクトル計算部２０２は、変換部２０１から供給される変換信号に基づいて、クロススペクトルを計算し、周波数別クロススペクトル計算部２０３１，２０３２，・・・，２０３Ｋへ伝達する。クロススペクトル計算部２０２は、変換信号Ｘ２（ｋ，ｎ）の複素共役と変換信号Ｘ１（ｋ，ｎ）との積を計算する。変換信号のクロススペクトルをＳ１２（ｋ，ｎ）とすると、クロススペクトルは次のように計算される。
Figure JPOXMLDOC01-appb-M000002
 　ここで、ｃｏｎｊ（Ｘ２（ｋ，ｎ））は、Ｘ２（ｋ，ｎ）の複素共役を表す。
　＜周波数別クロススペクトル計算部＞
　周波数別クロススペクトル計算部２０３１，２０３２，・・・，２０３Ｋは、クロススペクトル計算部２０２から供給されるクロススペクトルＳ１２（ｋ，ｎ）を用いて、Ｓ１２（ｋ，ｎ）の各周波数ｋに対応するクロススペクトルを計算し、周波数別クロススペクトルとして統合相関関数計算部２０４に伝達する。周波数別クロススペクトルの計算は、周波数成分ごとに相関関数を計算するために行われる。つまり、ある周波数ｋに対応する相関関数(周波数別相関関数と呼ぶ)を後段で求めるために、周波数別クロススペクトルを計算する。
　次に、ある周波数kの周波数別クロススペクトルを計算する周波数別クロススペクトル計算部２０３ｋについて詳細に説明する。図２Ｂは、周波数別クロススペクトル計算部２０３ｋのブロック図である。周波数別クロススペクトル計算部２０３ｋは、周波数別基本クロススペクトル計算部２０３１ｋを有する。周波数別クロススペクトル計算部２０３ｋは、クロススペクトル計算部２０２から供給されるクロススペクトルＳ１２（ｋ，ｎ）を用いて、周波数別基本クロススペクトルを計算し、周波数別クロススペクトルとして統合相関関数計算部２０４へ伝達する。周波数別基本クロススペクトル計算部２０３１ｋでは、周波数ｋのクロススペクトルＳ１２（ｋ，ｎ）に基づいて、周波数別基本クロススペクトルを計算する際、位相成分と振幅成分とをあらかじめ別々に求めたのちに統合する。周波数ｋの周波数別基本クロススペクトルＵｋ（ｗ，ｎ）、その振幅成分を｜Ｕｋ（ｗ，ｎ）｜、位相成分ａｒｇ（Ｕｋ（ｗ，ｎ））とすると、次のような関係が成立する。
Figure JPOXMLDOC01-appb-M000003
 　ここで、ｗは周波数を表し、０以上Ｗ－１以下の整数である。以下、周波数ｋのクロススペクトルＳ１２（ｋ，ｎ）から、周波数別基本クロススペクトルの振幅成分｜Ｕｋ（ｗ，ｎ）｜と位相成分ａｒｇ（Ｕｋ（ｗ，ｎ））とを求める方法について説明する。
　振幅成分｜Ｕｋ（ｗ，ｎ）｜は、ｋを整数倍した周波数には、１．０を用いる。一方、周波数kの非定数倍した周波数の位相成分はゼロにする。これを数式で表現すると、振幅成分｜Ｕｋ（ｗ，ｎ）｜は以下のように与えられる。
Figure JPOXMLDOC01-appb-M000004
 　ここで、ｐは１以上Ｐ以下の整数である。波源方向推定を行うときに重要な情報は位相成分であるため、振幅成分はこのように適当な定数を用いる。その他には、１．０の代わりに｜Ｓ１２（ｋ，ｎ）｜を用いてもよい。つまり、振幅成分｜Ｕｋ（ｗ，ｎ）｜を次式のように求めてもよい。
Figure JPOXMLDOC01-appb-M000005
 　位相成分ａｒｇ（Ｕｋ（ｗ，ｎ））は、ｋを整数倍した周波数には、周波数ｋのクロススペクトルＳ１２（ｋ，ｎ）を定数倍したものを用いる。例えば、周波数ｋ、２ｋ、３ｋ、４ｋの位相成分には、周波数ｋの位相成分ａｒｇ（Ｓ１２（ｋ，ｎ））をそれぞれ同一の倍率で整数倍したもの、つまりａｒｇ（Ｓ１２（ｋ，ｎ））、２ａｒｇ（Ｓ１２（ｋ，ｎ））、３ａｒｇ（Ｓ１２（ｋ，ｎ））、４ａｒｇ（Ｓ１２（ｋ，ｎ））を用いる。一方、周波数ｋの非定数倍した周波数の位相成分はゼロにする。したがって、周波数ｋに対応する周波数別基本クロススペクトルの位相成分ａｒｇ（Ｕｋ（ｗ，ｎ））は、以下のように計算される。
Figure JPOXMLDOC01-appb-M000006
 　ここで、ｐは１以上Ｐ以下の整数である。また、Ｐは１よりも大きい整数である。
　以上の方法で求めた振幅成分と位相成分とを、上述の式（３）を用いて統合し、周波数ｋの周波数別基本クロススペクトルＵｋ（ｗ，ｎ）を得る。
　これまで説明してきた方法では、振幅成分と位相成分とを別々に求めてから周波数別スペクトルを得ていた。しかし、次に示す数式に示すようにクロススペクトルのべき乗を用いれば、振幅成分と位相成分とを求めることなく周波数別スペクトルＵｋ（ｗ，ｎ）を求めることが可能である。
Figure JPOXMLDOC01-appb-M000007
 　統合相関関数計算部２０４は、周波数別クロススペクトル計算部２０３１，２０３２，・・・，２０３Ｋから供給される周波数別クロススペクトルに基づき統合相関関数を計算し、推定方向情報生成部２０５へ伝達する。
　＜統合相関関数計算部＞
　図２Ｃは、本実施形態に係る波源方向推定装置２００の備える統合相関関数計算２０４部の構成を示すブロック図である。統合相関関数計算部２０４は、周波数別相関関数生成部２４１１，２４１２，・・・，２４１Ｋと、統合部２４２を備える。
　周波数別相関関数生成部２４１１，２４１２，・・・，２４１Ｋは、周波数別クロススペクトル計算部２０３１，２０３２，・・・，２０３Ｋから供給される周波数別クロススペクトルの逆変換を行い、周波数別相関関数として統合部２４２にそれぞれ伝達する。本実施形態では、変換部２０１でフーリエ変換を用いたので、逆変換には逆フーリエ変換を用いる方法について説明する。周波数別クロススペクトル計算部２０３ｋから供給された周波数別クロススペクトルをＵｋ（ｗ，ｎ）とすると、Ｕｋ（ｗ，ｎ）の逆変換により得られる周波数別相関関数ｕｋ（τ，ｎ）は次のように計算される。
Figure JPOXMLDOC01-appb-M000008
 　統合部２４２は、周波数別相関関数生成部２４１１，２４１２，・・・，２４１Ｋから供給される周波数別相関関数を統合し、統合相関関数として推定方向情報生成部２０５へ伝達する。個別に求めた複数の周波数別相関関数を混合したり、重ね合わせたりすることにより、一つの相関関数を求める。統合方法に単純な和を用いた場合、統合部２４２は、周波数別相関関数の総和を計算する。統合相関関数をｕ（τ，ｎ）とすると、ｕ（τ，ｎ）は次のように計算される。
Figure JPOXMLDOC01-appb-M000009
 　また、総和では無く総乗を用いてもよい。この場合、ｕ（τ，ｎ）は次のように計算される。
Figure JPOXMLDOC01-appb-M000010
 　目標音の存在する周波数や、目的音のパワーが大きい周波数があらかじめ判明している場合には、その周波数に対応する周波数別相関関数だけ用いて、統合された相関関数を求めてもよい。また、重み付けという形で、統合における周波数別相関関数の影響度を制御してもよい。例えば目標音の存在する周波数の集合をΩとすると、周波数の選択によりｕ（τ，ｎ）を求める場合は、次のように計算される。
Figure JPOXMLDOC01-appb-M000011
 　また、重み付けを用いる場合には、ｕ（τ，ｎ）は次のように計算される。
Figure JPOXMLDOC01-appb-M000012
 　ここで、ａとｂとは実数で、ａ＞ｂ＞０を満足する。このように、目的音が存在する周波数の周波数別相関関数を重点的に用いて統合すると、雑音などの非目的音の影響が小さい相関関数を生成できるため、方向推定精度が向上する。
　相対遅延時間計算部２０６は、入力されたマイク位置情報と音源探索対象方向とから、マイクペア間の相対遅延時間を求め、音源探索対象方向とセットで推定方向情報生成部２０５に伝達する。相対遅延時間とは、マイク間隔と音源方向とに基づいて一意に定まる音波の到達時間差のことである。音速をｃと仮定し、２つのマイクの間隔をｄ、音源方向、つまり音の到来方向をθとした場合、音源方向θに対する相対遅延時間τ(θ)は次の式で計算される。
Figure JPOXMLDOC01-appb-M000013
 　相対遅延時間は、全ての音源探索対象方向に対して計算される。例えば、方向の探索範囲が１０度刻みで０度から９０度まで、つまり０度、１０度、２０度、・・・、９０度である場合、１０種類の相対遅延時間を計算する。そして、探索対象の方向と相対遅延時間をペアで推定方向情報生成部２０５へ供給する。
　推定方向情報生成部２０５は、統合相関関数計算部２０４から供給される統合相関関数と、相対遅延時間計算部２０６から供給される相対遅延時間に基づき、方向と相関値との対応関係を推定方向情報として出力する。相関関数をｕ（τ，ｎ）、相対遅延時間τ（θ）をとすると、推定方向情報Ｈ（θ，ｎ）は、次式で与えられる。
Figure JPOXMLDOC01-appb-M000014
 　各方向に対して相関値が定まるので、基本的に相関値が高ければ、その方向に音源が存在する可能性が高いと判断できる。
　このような推定方向情報は、様々な形で利用される。例えば、関数が複数のピークを有する場合には、各ピークを到来方向とする音源が複数存在すると考えられる。したがって、各音源の方向を同時に推定できるだけでなく、音源数の推定にも用いることが可能である。
　また、相関関数のピークと非ピークとの差分に基づき、音源の存在可能性を判定することも可能である。もし、ピークと非ピークとの差分が大きければ音源の存在可能性が高いと判定できる。同時に、推定方向の信頼性も高いと判断できる。もし、音源数が１つと事前に仮定できる場合には、相関値が最大となる方向を推定方向情報として出力してもよい。この場合、推定方向情報は、方向と相関値との対応関係ではなく、方向そのものとなる。
　＜周波数別クロススペクトルの説明＞
　以上の方法により周波数別クロススペクトルを計算すると、周波数別クロススペクトルを逆変換して得られる周波数別相関関数のピークが鋭くなり、相関関数のピーク位置が明確になる。相関関数のピーク位置に基づいて波源方向推定を行う本実施形態の場合には、ピークが鋭くなれば音源方向推定の精度が向上する。さらに、Ｐの値が大きくなるにつれて、すなわちｋを整数倍した周波数の成分が増えるにつれて、相関関数のピークは鋭くなる。図３Ａに、この様子を示す。ここで、図中のＱは３よりも大きい整数である。Ｐ＝１の場合、つまり位相成分が一つしかない場合、これを逆変換して得られる相関関数は、このようにピーク位置が不明確な相関関数になる。Ｐが大きくなると、図３Ａに示す通り、相関関数のピークは鋭くなる。
　従って本発明では、周波数別クロススペクトルは、ある周波数kのクロススペクトルを基に「周波数kを整数倍した周波数ｐｋの位相成分に、周波数kの位相成分ａｒｇ（Ｓ１２（ｋ，ｎ））をｐ倍したものを割り当てたもの」として定義される。ここで、ｐは１以上の整数である。つまり、周波数別クロススペクトルは、その位相成分ａｒｇ（Ｕｋ（ｗ，ｎ））が少なくとも以下の式を満足するものとして定義される。
Figure JPOXMLDOC01-appb-M000015
 　加えてｐは、ｐ＝１と２、ｐ＝１と３、ｐ＝２と３など、２つ以上に限定される。ｐが１のみの場合には、周波数kの成分だけを抜き出して周波数別クロススペクトルが生成されることになるが、方向推定精度は従来技術と同等となり、方向推定の高精度を達成できない。なお、図３Ａで説明した通り、ｐ＝１，２，３，４・・・とｐの数が増えれば、つまり、周波数ｐｋの位相成分への割り当てが増えれば、周波数別相関関数のピークは鋭くなるため、方向推定の精度は向上する。
　周波数別スペクトルを計算する際にPを大きくする効果について、図３Ｂを用いて説明する。図３Ｂは、異なるＰ＝１とＰ＝５との場合に対する、２つの周波数別相関関数と、それらを統合した統合相関関数の例を示した図である。図３Ｂの下段に示す通り、周波数別相関関数のピークは周期的に出現し、そのピーク間隔は周波数ｋに反比例する。周波数ｋが高くなれば、隣接する２つの周波数別相関関数のピークが接近し、相関関数の重なりからピークは目立たなくなる。Ｐ＝１の場合、周波数別相関関数のピークは鋭くないため、統合して得られた統合関数のピークも滑らかである。このため、推定方向は不明確となる。一方、Ｐ＝５の場合、周波数別相関関数、およびその統合相関関数のピークは鋭くなる。このため、推定方向は明確になり、推定精度の向上につながる。
　図４は、本実施形態に係る波源方向推定装置２００の備える統合相関関数テーブル４０１の構成の一例を示す図である。統合相関関数テーブル４０１は、入力信号４１１に関連付けて周波数領域信号４１２と、クロススペクトル４１３と、周波数別クロススペクトル４１４と、統合相関関数４１５とを記憶する。波源方向推定装置２００は、入力信号を取得するたびに統合相関関数を計算してもよいし、入力信号に対応する統合相関関数をあらかじめ求めておき、統合相関関数テーブル４０１を参照して統合相関関数を算出してもよい。
　図５は、本実施形態に係る波源方向推定装置２００のハードウェア構成を示すブロック図である。
　ＣＰＵ（Central Processing Unit）５１０は演算制御用のプロセッサであり、プログラムを実行することで図２Ａの波源方向推定装置２００の機能構成部を実現する。ＲＯＭ（Read Only Memory）５２０は、初期データおよびプログラムなどの固定データおよびプログラムを記憶する。また、通信制御部５３０は、ネットワークを介してその他の装置などと通信する。なお、ＣＰＵ５１０は１つに限定されず、複数のＣＰＵであっても、あるいは画像処理用のＧＰＵ（Graphics Processing Unit）を含んでもよい。また、通信制御部５３０は、ＣＰＵ５１０とは独立したＣＰＵを有して、ＲＡＭ（Random Access Memory）５４０の領域に送受信データを書き込みあるいは読み出しするのが望ましい。また、ＲＡＭ５４０とストレージ５５０との間でデータを転送するＤＭＡＣ（Direct Memory Access Controller）を設けるのが望ましい（図示なし）。さらに、入出力インタフェース５６０は、ＣＰＵ５１０とは独立したＣＰＵを有して、ＲＡＭ５４０の領域に入出力データを書き込みあるいは読み出しするのが望ましい。したがって、ＣＰＵ５１０は、ＲＡＭ５４０にデータが受信あるいは転送されたことを認識してデータを処理する。また、ＣＰＵ５１０は、処理結果をＲＡＭ５４０に準備し、後の送信あるいは転送は通信制御部５３０やＤＭＡＣ、あるいは入出力インタフェース５６０に任せる。
　ＲＡＭ５４０は、ＣＰＵ５１０が一時記憶のワークエリアとして使用するランダムアクセスメモリである。ＲＡＭ５４０には、本実施形態の実現に必要なデータを記憶する領域が確保されている。入力信号５４１は、マイクなどの集音装置で集音した音信号データや入力信号取得装置などに入力され、取得された信号データなどである。
　周波数領域信号５４２は、入力信号５４１を変換部２０１で変換した信号である。クロススペクトル５４３は、クロススペクトル計算部２０２で計算したスペクトルである。周波数別クロススペクトル５４４は、周波数別クロススペクトル計算部２０３ｋで計算したスペクトルである。統合相関関数５４５は、統合相関関数計算部２０４で計算した関数である。
　入出力データ５４６は、入出力インタフェース５６０を介して入出力されるデータである。送受信データ５４７は、ネットワークインタフェース５３０を介して送受信されるデータである。また、ＲＡＭ５４０は、各種アプリケーションモジュールを実行するためのアプリケーション実行領域５４８を有する。
　ストレージ５５０には、データベースや各種のパラメータ、あるいは本実施形態の実現に必要な以下のデータまたはプログラムが記憶されている。ストレージ５５０は、統合相関関数テーブル４０１を格納する。統合相関関数テーブル４０１は、図４に示した、入力信号と統合相関関数との関係を管理するテーブルである。
　ストレージ５５０は、さらに、変換モジュール５５１と、クロススペクトル計算モジュール５５２と、周波数別クロススペクトル計算モジュール５５３と、統合相関関数計算モジュール５５４とを格納する。また、ストレージ５５０は、推定方向情報生成モジュール５５５と、相対遅延時間計算モジュール５５６とを格納する。
　変換モジュール５５１は、入力信号を周波数領域信号に変換するモジュールである。クロススペクトル計算モジュール５５２は、周波数領域信号に基づいてクロススペクトルを計算するモジュールである。周波数別クロススペクトル計算モジュール５５３は、クロススペクトルを用いて、周波数別クロススペクトルを計算するモジュールである。統合相関関数計算モジュール５５４は、周波数別クロススペクトルに基づいて、統合相関関数を計算するモジュールである。
　推定方向情報生成モジュール５５５は、統合包絡関数に基づいて波源の推定方向情報を生成するモジュールである。相対遅延時間計算モジュール５５６は、相対遅延時間を計算するモジュールである。これらのモジュール５５１～５５６は、ＣＰＵ５１０によりＲＡＭ５４０のアプリケーション実行領域５４８に読み出され、実行される。制御プログラム５５７は、波源方向推定装置２００の全体を制御するためのプログラムである。
　入出力インタフェース５６０は、入出力機器との入出力データをインタフェースする。入出力インタフェース５６０には、表示部５６１、操作部５６２、が接続される。また、入出力インタフェース５６０には、さらに、記憶媒体５６４が接続されてもよい。さらに、音声出力部であるスピーカ５６３や、音声入力部であるマイク、あるいは、ＧＰＳ位置判定部が接続されてもよい。なお、図５に示したＲＡＭ５４０やストレージ５５０には、波源方向推定装置２００が有する汎用の機能や他の実現可能な機能に関するプログラムやデータは図示されていない。
　図６本実施形態に係る波源方向推定装置２００の処理手順を示すフローチャートである。このフローチャートは、図５のＣＰＵ５１０がＲＡＭ５４０を使用して実行し、図２の波源方向推定装置２００の機能構成部を実現する。
　ステップＳ６０１において、波源方向推定装置２００は、入力信号を取得する。ステップＳ６０３において、波源方向推定装置２００の変換部２０１は、入力端子２０１および入力端子２０２から供給された入力信号を変換する。変換部２０１は、変換して得られた周波数領域信号をクロススペクトル計算部２０２へ供給する。ステップＳ６０４において、クロススペクトル計算部２０２は、供給された変換信号に基づいて、クロススペクトルを計算する。クロススペクトル計算部２０２は、計算したクロススペクトルを周波数別クロススペクトル計算部２０３１，２０３ｋ，・・・２０３Ｋへ伝達する。
　ステップＳ６０７において、周波数別クロススペクトル計算部２０３１，２０３ｋ，・・・２０３Ｋは、クロススペクトルの各周波数ｋに対応するクロススペクトルを計算する。すなわち、周波数別クロススペクトル計算部２０３１，２０３ｋ，・・・２０３Ｋは、周波数別クロススペクトルを計算する。そして、周波数別クロススペクトル計算部２０３１，２０３ｋ，・・・２０３Ｋは、周波数別クロススペクトルを統合相関関数計算部２０４へ伝達する。
　ステップＳ６０９において、周波数別相関関数生成部２４１１，２４１２，・・・，２４１Ｋは、周波数別クロススペクトルを逆変換して、周波数別相関関数を計算する。ステップＳ６１１において、統合部２４２は、周波数別相関関数を統合して、統合相関関数を計算する。
　ステップＳ６１３において、相対遅延時間計算部２０６は、マイク位置情報と音源探索対象方向とから、マイクペア間の相対遅延時間を計算する。ステップＳ６１５において、推定方向情報生成部２０５は、統合相関関数と相対遅延時間とから推定方向情報を生成する。
　本実施形態によれば、入力信号に含まれる目標音の到来方向、すなわち目標物体が存在する方向を推定する。環境雑音レベルが高い環境において、目標物が発する音を手掛かりに、目標物が存在する方向を推定する場合に有効である。環境雑音の例としては、繁華街や街頭、街道沿い、人や自動車が多く集まる場所が挙げられる。また、目標物の例としては、人間や動物、自動車、航空機、船舶、水上バイク、ドローン（小型無人機）が挙げられる。
　例えば、屋外のテーマパークや展示会場などに接近する不審な自動車・船舶・ドローンなどを検知し、それらの方向を推定することで、不審者や不審物を効率的に取り締ることが可能である。また、音源方向推定を複数箇所で実施することで、目標音源の位置を特定できる。これにより、環境雑音レベルが高い環境でも、悲鳴や銃声、自動車の衝突音の発生箇所などを正確に特定することが可能となる。
　［第３実施形態］
　次に本発明の第３実施形態に係る波源方向推定装置について、図７を用いて説明する。図７は、本実施形態に係る波源方向推定装置の備える統合相関関数生成部７０４の構成を説明するためのブロック図である。本実施形態に係る波源方向推定装置の備える統合相関関数生成部７０４は、上記第２実施形態の統合相関関数生成部２０４と比べると、周波数別相関関数生成部２４１１，２４１２，・・・，２４１Ｋと統合部２４２とに代えて、統合部７４１と統合相関関数生成部７４２とを有する点で異なる。その他の構成および動作は、第２実施形態と同様であるため、同じ構成および動作については同じ符号を付してその詳しい説明を省略する。
　統合部７４１は、周波数別クロススペクトル計算部２０３１，２０３２，・・・，２０３Ｋから供給される周波数別クロススペクトルを統合し、統合クロススペクトルとして統合相関関数生成部７４２へ伝達する。個別に求めた複数の周波数別クロススペクトルを混合したり、重ね合わせたりすることにより、一つの統合クロススペクトルを求める。統合には、第２実施形態の統合部２４２と同様に総和や総乗が用いられる。統合に総和を用いた場合、統合クロススペクトルＵ（ｋ，ｎ）は次のように計算される。
Figure JPOXMLDOC01-appb-M000016
 　また、総乗を用いた場合、統合クロススペクトルＵ（ｋ，ｎ）は次のように計算される。
Figure JPOXMLDOC01-appb-M000017
 　第２実施形態の統合部２４２と同様に、目標音源の存在する周波数や、目的音源のパワーが大きい周波数があらかじめ判明している場合には、統合クロススペクトルＵ（ｋ，ｎ）を生成する際に補正してもよい。第２実施形態と同様に、周波数の選択や重み付けという形で影響度を制御する。例えば目標音の存在する周波数の集合をΩとすると、帯域の選択により統合クロススペクトルＵ（ｋ，ｎ）を求める場合は、次のように計算する。
Figure JPOXMLDOC01-appb-M000018
 　また、重み付けを用いる場合には、Ｕ（ｋ，ｎ）は次のように計算される。
Figure JPOXMLDOC01-appb-M000019
 　ここで、ａとｂとは実数で、ａ＞ｂ＞０を満足する。このように、目的音が存在する周波数の周波数別相関関数を重点的に用いて統合すると、雑音などの非目的音の影響が小さい相関関数を生成できるため、方向推定精度が向上する。
　統合相関関数生成部７４２は、統合部７４１から供給される統合クロススペクトルの逆変換を行い、統合相関関数として推定方向情報生成部２０５に伝達する。本実施形態でも、逆変換には逆フーリエ変換を用いる方法について説明する。統合部７４１から供給された統合クロススペクトルをＵ（ｋ，ｎ）とすると、Ｕ（ｋ，ｎ）の逆変換により得られる統合相関関数ｕ（τ，ｎ）は次のように計算される。
Figure JPOXMLDOC01-appb-M000020
 　本実施形態によれば、周波数別クロススペクトルを統合してから逆変換を行い、統合相関関数を得る。このため、周波数別クロススペクトルごとに逆変換を行っていた第２実施形態と比べると、逆変換の回数が少なくなる。したがって、第２実施形態よりも少ない計算量で統合相関関数を求めることができる。
　［第４実施形態］
　次に本発明の第４実施形態に係る波源方向推定装置について、図８Ａ乃至図１１を用いて説明する。図８Ａは、本実施形態に係る波源方向推定装置８００の構成を説明するためのブロック図である。本実施形態に係る波源方向推定装置８００は、上記第２実施形態と比べると、周波数別クロススペクトル計算部２０３１，２０３２，・・・，２０３Ｋに代えて、周波数別クロススペクトル計算部８０３１，８０３２，・・・，８０３Ｋを有する点で異なる。その他の構成および動作は、第1実施形態と同様であるため、同じ構成および動作については同じ符号を付してその詳しい説明を省略する。
　図８Ｂは、周波数別クロススペクトル計算部８０３ｋのブロック図である。周波数別クロススペクトル計算部８０３ｋは、周波数別基本クロススペクトル計算部２０３１ｋ、カーネル関数スペクトル記憶部８３１、乗算部８３２を有する。周波数別基本クロススペクトル計算部２０３１ｋは、クロススペクトル計算部２０２から供給されるクロススペクトルＳ１２（ｋ，ｎ）を用いて、Ｓ１２（ｋ，ｎ）の周波数ｋに対応するクロススペクトルを計算し、周波数別基本クロススペクトルとして乗算部８３２へ伝達する。周波数別基本クロススペクトル計算部２０３１ｋの動作は、その出力先を除いて第２実施形態の周波数別基本クロススペクトル計算部２０３１ｋと同様であるため、詳しい説明は省略する。
　カーネル関数スペクトル記憶部８３１は、カーネル関数スペクトルを記憶しており、乗算部８３２へカーネル関数スペクトルを出力する。カーネル関数スペクトルとは、カーネル関数をフーリエ変換し、その絶対値を取ったものである。絶対値を取る代わりに二乗してもよい。カーネル関数としては、ガウス関数を用いる。ガウス関数は、あらかじめ与えられた３つの実数ｇ１，ｇ２，ｇ３を用いて、次のような数式で与えられる。
Figure JPOXMLDOC01-appb-M000021
 　ｇ１はガウス関数の大きさ、ｇ２はガウス関数のピークの位置、ｇ３はガウス関数の広がりを制御する。特に、ガウス関数の広がりを調整するｇ３は、周波数別相関関数のピークの鋭さに大きな影響を与えるので重要である。式（２１）から分かるように、ｇ３が大きくなればガウス関数の広がりは大きくなる。
　他には、以下に示すロジスティック関数を用いてもよい。
Figure JPOXMLDOC01-appb-M000022
 　ここで、ｇ１とｇ２とは実数である。ロジスティック関数は、ガウス関数と同様の形状をしているが、ガウス関数よりも裾が長いという性質を持つ。特に、ロジスティック関数の広がりを調整するｇ５は、ガウス関数におけるｇ３の場合と同様、周波数別相関関数のピークの鋭さに大きな影響を与える重要なパラメータである。この他にも、コサイン関数や一様関数を用いてもよい。
　カーネル関数で用いられるパラメータｇ１～ｇ５は、定数ではなく、周波数ｋに応じて異なる値を用いてもよい。つまり、ｇ１（ｋ）～ｇ５（ｋ）のように周波数kの関数としてもよい。例えば、ｇ３を周波数ｋの関数ｇ３（ｋ）とし、周波数が高くなるにつれて小さな値を持つ関数とする。このような代表例としてｋの逆数を関数ｇ３（ｋ）関数とする場合、ｇ３（ｋ）は次のように与えられる。
Figure JPOXMLDOC01-appb-M000023
 　ここで、Ｇ３は実数である。このとき、カーネル関数Ｇ（ｋ）は、周波数ｋが高くなるにつれて、ピークが鋭く裾が狭い関数になる。
　乗算部８３２は、周波数別基本クロススペクトル計算部２０３１ｋから供給される周波数別基本クロススペクトルと、カーネル関数スペクトル記憶部８３１から供給されるカーネル関数スペクトルの積を計算し、周波数別クロススペクトルとして統合相関関数計算部２０４に伝達する。周波数別基本クロススペクトル計算部２０３１ｋから供給される周波数別基本スペクトルをＵｋ（ｗ，ｎ）、カーネル関数スペクトル記憶部８３１から供給されるカーネル関数スペクトルをＧ（ｗ）とすると、周波数別クロススペクトルＵＭｋ（ｗ，ｎ）は次のように計算される。
Figure JPOXMLDOC01-appb-M000024
 　このように、周波数別基本クロススペクトルにカーネル関数スペクトルを乗じることで、統合相関関数計算部２０４にある周波数別相関関数生成部２４１ｋで得られる周波数別相関関数の大きさを変更できる。図９は、カーネル関数スペクトルが乗じられた周波数別クロススペクトルと、周波数別相関関数の関係を示している。比較のため、カーネル関数スペクトルを乗じる前の周波数別クロススペクトルも掲載する。図９の左側図に示すように、カーネル関数スペクトルを乗じていなければ、高い周波数まで成分が存在するため、周波数別相関関数のピークは鋭くなる。一方で、図９の中央図、右側図のようにカーネル関数スペクトルを乗じると、高い周波数の成分が減衰するので、周波数別相関関数のピークの鋭さは小さくなる。つまり、カーネル関数スペクトルのピークが鋭くなる(カーネル関数スペクトルの裾が狭くなる)につれて、周波数別相関関数のピークの鋭さは小さくなる。また、図９の右側図のように、周波数別相関関数の裾が大きく広がると、隣接する山の裾が重なり、谷が浅い周波数別相関関数が得られる。
　ここで、カーネル関数とカーネル関数スペクトルの形状の関係について補足する。フーリエ変換の性質により、形状の関係は逆になる。カーネル関数のピークが鋭く、裾が狭くなるにつれて、カーネル関数スペクトルのピークは平坦に近づき、裾が広くなる。ガウス関数の広がりを調整するｇ３との関係も含めて説明すると、ｇ３が大きくなれば、ガウス関数の広がりは大きくなるが、そのスペクトルの広がりは小さくなる。
　カーネル関数による周波数別相関関数の大きさ制御の効果について、図１０において説明する。図１０は、カーネル関数の有無と統合相関関数の関係を示した図である。図１０の（ａ）に示すようなカーネル関数無しの場合、周波数別相関関数ｕ１（τ，ｎ）～ｕ３（τ，ｎ）のピーク位置は近いが、ｕ１（τ，ｎ）～ｕ３（τ，ｎ）の幅が細いため、統合時に大きなピークを形成できない。このため、ピークの位置が明確にならない。一方、図１０の（ｂ）に示すようなカーネル関数有りの場合、周波数別相関関数の幅は太くなっているため、統合によりｕ１（τ，ｎ）～ｕ３（τ，ｎ）は大きなピークを形成できる。このため、（ａ）のカーネル関数無しの場合に比べて、ピークの位置が明確になる。
　さらに、カーネル関数による周波数別相関関数の大きさ制御の別の効果について、図１１で説明する。図１１は、カーネル関数スペクトルの幅の違いと、統合相関関数との関係を示した図である。図９の右側図に示すように、幅の広いカーネル関数スペクトルを用いると、相関関数の周期性により、谷が浅い周波数別相関関数が形成される。したがって、図１１の（ｃ）に示す通り、谷が浅い周波数別相関関数を統合すると、谷が浅い、つまりピークの目立たない統合相関関数が生成される。一方、図９の中央図に示すように、幅の狭いカーネル関数スペクトルを用いた場合には、図９の右図よりも谷が深い周波数別相関関数が形成される。したがって、図１１の（ｄ）に示す通り、ピークが明確な統合相関関数が生成される。
　本実施形態では、カーネル関数のフーリエ変換により得られたカーネル関数スペクトルと、周波数別基本クロススペクトルとの積を計算しているが、フーリエ変換の性質により時間領域でも実現できる。周波数別クロススペクトル計算部８０３Ｋでは無く、統合相関関数計算部２０４の中にある周波数別相関関数生成部２４１ｋの後段にカーネル関数を畳み込む「畳み込み演算部」を設け、周波数別相関関数生成部２４１ｋから供給される周波数別相関関数にカーネル関数を畳み込んでもよい。ただし、畳み込み演算は計算量が多いため、本実施形態のように周波数領域で積を計算するほうが効率的である。
　本実施形態によれば、周波数別基本クロススペクトルにカーネル関数スペクトルを乗じて周波数別クロススペクトルを生成する。このため、逆変換により得られる周波数別相関関数の幅が太くなり、統合相関関数のピークが明確になる。特に、個別の周波数別相関関数のピーク位置が近いにも関わらず、それぞれ鋭いピークを持つ場合、補正を行うことにより統合相関関数のピークの明確化効果が高くなる。
　［第５実施形態］
　次に本発明の第５実施形態に係る波源方向推定装置について、図１２を用いて説明する。図１２は、本実施形態に係る波源方向推定装置の備える周波数別クロススペクトル計算部の構成を説明するための図である。本実施形態に係る波源方向推定装置の備える周波数別クロススペクトル計算部１２０３ｋは、上記第４実施形態と比べると、カーネル関数スペクトル記憶部８３１に代えて、カーネル関数スペクトル生成部１２３１を有する点で異なる。その他の構成および動作は、第４実施形態と同様であるため、同じ構成および動作については同じ符号を付してその詳しい説明を省略する。
　カーネル関数スペクトル生成部１２３１は、クロススペクトル計算部２０２から供給されるクロススペクトルに用いてカーネル関数スペクトルを生成し、生成したカーネル関数スペクトルを乗算部８３２に伝達する。カーネル関数スペクトル生成部１２３１は、供給されたクロススペクトルを分析することで、入力信号に目標音が存在する可能性を求め、存在可能性を反映した形状を有するカーネル関数スペクトルを生成する。基本的に、存在可能性が低ければ、幅が狭い、広がりが小さいカーネル関数スペクトルを生成する。これにより、周波数別相関関数のピークは低くなるため、統合相関関数に誤ったピークが出現する可能性を低減できる。
　目標音の存在可能性を求める方法として、入力信号のＳＮＲ（目的信号対雑音比：Signal-to-Noise Ratio）を推定する方法について説明する。まず、供給されたクロススペクトルの絶対値を算出する。一般的には、変換部２０１で得られたフーリエ変換を二乗したものを入力信号パワースペクトルと呼ぶが、本実施形態ではクロススペクトルの絶対値を入力信号パワースペクトルとして扱う。次に、入力信号パワースペクトルに基づき、入力信号に含まれる雑音成分(非目的音成分)のパワースペクトルを推定する。入力信号パワースペクトルをＰＸ（ｋ，ｎ）とすると、ＰＸ（ｋ，ｎ）は以下のように計算される。
Figure JPOXMLDOC01-appb-M000025
 　次に、入力信号パワースペクトルに基づいて、雑音成分のパワースペクトルを推定する。ここでは、非特許文献３に記載の方法を用いる。推定雑音パワースペクトルを、入力信号パワースペクトルが供給され始めた推定初期のパワースペクトルを平均化したものとする。この場合、推定が開始された直後には目的音が含まれないという条件を満たす必要がある。推定雑音パワースペクトルをＰＮ（ｋ，ｎ）とすると、ＰＮ（ｋ，ｎ）は以下のように計算される。
Figure JPOXMLDOC01-appb-M000026
 　ここで、Ｎ０はあらかじめ定めた整数である。
　他の方法としては、推定雑音パワースペクトルを、入力信号パワースペクトルの最小値（最小統計量）から求める方法が非特許文献４に開示されている。この方法では、一定時間内における入力信号パワースペクトルの最小値を周波数ごとに保持し、その最小値から雑音成分を推定する。入力信号パワースペクトルの最小値は、雑音のパワースペクトルとスペクトル形状と似ているため、雑音パワースペクトルの推定値として用いることができる。
　推定雑音パワースペクトルが得られたあと、入力信号パワースペクトルとの比をとってＳＮＲの推定値とする。入力信号パワースペクトルをＰＸ（ｋ，ｎ）、推定雑音パワースペクトルをＰＮ（ｋ，ｎ）とすると、推定ＳＮＲγ（ｋ，ｎ）は、次のように計算される。
Figure JPOXMLDOC01-appb-M000027
 　そして、この推定ＳＮＲであるγ（ｋ，ｎ）をそのまま目的音の存在可能性ｑ（ｋ，ｎ）に用いる。
　このようにして得られた推定ＳＮＲは、非特許文献３において、推定a-posteriori　ＳＮＲと呼ばれている。推定ＳＮＲには、推定a-posteriori　ＳＮＲの代わりに、非特許文献３に記載の方法で得られる推定a-priori　ＳＮＲを用いてもよい。a-priori　ＳＮＲの推定には、雑音成分の抑圧を行ってからＳＮＲを推定するので、計算量は多くなるもののa-posteriori　ＳＮＲよりも高い推定精度を達成できる。
　入力信号パワースペクトルと推定雑音パワースペクトルを用いて目標音の存在可能性を計算する方法は、推定ＳＮＲのような両者の比に限定されない。例えば、比の代わりに、両者の差分を用いてもよい。また、単純な大小関係を用いてもよい。
　クロススペクトルの分析により目標音が存在する可能性を求める方法としては、パワースペクトルを用いる方法に限定されない。他の代表的な例は、クロススペクトルの位相成分を分析する方法が挙げられる。位相成分の分析方法として、クロススペクトルの群遅延（位相成分を周波数方向に微分したもの）を用いた方法について説明する。はじめに、クロススペクトルの群遅延を求める。群遅延をｇｄ（ｋ，ｎ）とすると、クロススペクトルＳ１２（ｋ，ｎ）の群遅延は次のように算出できる。
Figure JPOXMLDOC01-appb-M000028
 　そして、ｇｄ（ｋ，ｎ）の平均値を計算し、その平均値からの乖離度を存在可能性とする。例えば、ガウス関数を用いて目標音の存在可能性を計算する場合、存在可能性ｑ（ｋ，ｎ）は次のように計算できる。
Figure JPOXMLDOC01-appb-M000029
 　ここで、ｑ０は正の実数である。また、ｇｄ（ｋ，ｎ）バーは、ｇｄ（ｋ，ｎ）を周波数方向に平均化したものである。平均化には様々な方法があるが、例えば、以下のような相加平均を用いることができる。
Figure JPOXMLDOC01-appb-M000030
 　式（２９）を参照すると、ｇｄ（ｋ，ｎ）がｇｄ（ｋ，ｎ）バーに近ければ、ｑ（ｋ，ｎ）は１に近づく。一方、ｇｄ（ｋ，ｎ）がｇｄ（ｋ，ｎ）バーから遠ざかるにつれて、ｑ（ｋ，ｎ）は０に近づく。
　次に得られた存在可能性を用いて、カーネル関数スペクトルを生成する。ここでは、カーネル関数スペクトルの元となるカーネル関数のパラメータを制御する例について説明する。また、カーネル関数としては、ガウス関数を用いる例を説明する。目的音の存在可能性が高いときは、ｇ３を小さくする。これにより、存在可能性が高くなるほど、ｇ（τ）の幅は狭くなり、ｇ（τ）ピークが強調された形状に近づく。目的音の存在可能性からｇ３を求めるには、存在可能性の逆数を変数とする一次関数を用いる。このとき、存在可能性をｑ（ｋ，ｎ）とすると、ｇ３は次のように計算される。
Figure JPOXMLDOC01-appb-M000031
 　ここで、ａ１とｂ１とは実数であり、ａ１＞０．０、ｂ１＞０．０を満足する。目的音の存在可能性ｑ（ｋ，ｎ）からｇ３を求める関数は、１次関数に限定されない。シグモイド関数や高次の多項式関数、非線形関数など、他の形で表される関数を1次関数の代わりに用いることも可能である。
　カーネル関数として、ロジスティック関数を用いる場合には、ｇ３と同様の方法でｇ５を計算すればよい。その結果、目的音の存在可能性が高ければ、ｇ５は小さくなるので、カーネル関数ｇ（τ）の幅は狭く、ピークが強調された形状に近づく。
　このように存在可能性からパラメータを生成した後、カーネル関数、およびカーネル関数スペクトルを生成する。
　本実施形態によれば、目的音の存在可能性を求め、その可能性に基づいてカーネル関数を算出する。可能性が高い場合には、カーネル関数スペクトルの幅は広がり、形状は平坦なものに近づく。逆に可能性が低い場合には、カーネル関数スペクトルの幅は狭くなる。これにより、目的音が存在する周波数の周波数別相関関数のピークが高くなり、目的音が存在しない周波数の周波数別相関関数のピークが低くなる。以上から、第４実施形態よりも統合相関関数のピークが強調され、目的音の方向推定精度が向上する。特に、非目的音の周波数別相関関数のピークが低くなるので、統合相関関数に誤ったピークが出現する可能性を低減できる。
　［第６実施形態］
　次に本発明の第６実施形態に係る波源方向推定装置について、図１３Ａおよび図１３Ｂを用いて説明する。図１３Ａは、本実施形態に係る波源方向推定装置１３００の構成を説明するためのブロック図である。本実施形態に係る波源方向推定装置１３００は、上記第３実施形態と比べると、統合相関関数計算部２０４に代えて、統合相関関数計算部１３０４を有する点で異なる。その他の構成および動作は、第３実施形態と同様であるため、同じ構成および動作については同じ符号を付してその詳しい説明を省略する。
　図１３Ｂは、本発明の第６実施形態に係る波源方向推定装置の備える周波数別クロススペクトルの構成を説明するためのブロック図である。本実施形態に係る周波数別クロススペクトル計算部２０３ｋは、上記第３実施形態と比べると、統合部７４１に代えて、統合クロススペクトル生成部１３４１を有する点で異なる。その他の構成および動作は、第３実施形態と同様であるため、同じ構成および動作については同じ符号を付してその詳しい説明を省略する。
　統合クロススペクトル生成部１３４１は、クロススペクトル計算部２０２から供給されるクロススペクトルに基づき、周波数別クロススペクトル計算部２０３１，２０３２，・・・，２０３Ｋから供給される周波数別クロススペクトルを統合し、統合クロススペクトルとして統合相関関数生成部７４２へ伝達する。第３実施形態では、目標音の存在する周波数や、目的音のパワーが大きい周波数があらかじめ判明している場合について説明した。本実施形態では、供給されたクロススペクトルを分析することで、入力信号に目標音が存在する可能性を求め、存在可能性に基づいて統合する。
　まず、供給されたクロススペクトルに基づき、目標音の存在可能性を求める。存在可能性の計算には、第５実施形態で説明した方法を同様に用いることができる。次に、求めた存在可能性を用いて、周波数クロススペクトルの統合を行う。まず、目標音の存在可能性をｑ（ｋ，ｎ）とすると、目標音の存在する可能性が高い周波数の集合Ωをｑ（ｋ，ｎ）に基づいて求める。ある周波数kに対するｑ（ｋ，ｎ）があらかじめ定めた閾値θｑを超えれば、その周波数を集合Ωの要素とする。これを数式で表現すると、次のようになる。
Figure JPOXMLDOC01-appb-M000032
 　集合Ωが定まれば、第３実施形態で説明した方法を用いればよい。具体的には、式（１７）や式（１８）に示した計算式で求めればよい。
　また、存在可能性ｑ（ｋ，ｎ）を用いて重みを計算し、その重みを用いて重み付け和による統合を行ってもよい。重み関数をη（ｑ（ｋ，ｎ））とすると、統合クロススペクトルＵ（ｋ，ｎ）を次のように計算する。
Figure JPOXMLDOC01-appb-M000033
 　ただし、重み関数η（ｑ（ｗ，ｎ））は、大きなｑ（ｗ，ｎ）に対して大きな値を取る単調増加関数とする。
　本実施形態によれば、クロススペクトルに基づき目標音の存在可能性を求めてから、その存在可能性を用いて統合クロススペクトルを計算する。このため、あらかじめ目標音の存在可能性が不明の状態でも、統合クロススペクトル生成時の帯域選択や重み付けが適切に行えるようになり、高い推定精度を達成できる。
　［第７実施形態］
　次に本発明の第７実施形態に係る波源方向推定システムについて、図１４および図１５を用いて説明する。図１４は、本実施形態に係る波源方向推定システム１４００の構成を説明するための図である。本実施形態に係る波源方向推定システム１４００は、上記第２実施形態に係る波源方向推定装置２００を用いている。したがって、第２実施形態と同じ構成および動作については同じ符号を付してその詳しい説明を省略する。
　本実施形態に係る波源方向推定システム１４００は、マイク１４０１、マイク１４０２、ＡＤ変換部１４０１および表示部１４０２を有する。なお、本実施形態では、波源方向推定装置２００の代わりに波源方向推定装置８００、または波源方向推定装置１３００を用いることが可能である。また、波源が音源であるという仮定で説明するため、マイクを用いた例を説明するが、音源以外の場合には、その波源から放射される波動を受信して電気信号に変換できる各種センサが、マイクの代わりに用いられる。
　マイク１４０１およびマイク１４０２は、推定対象である目標物体から生じる音を含めた装置周辺の音を電気信号に変換し、ＡＤ変換部１４０１へ伝達する。音が伝わる媒質が空気媒質である場合、音は空気の振動としてマイクに到達する。マイクは、到達した空気の振動を電気信号に変換する。
　ＡＤ変換部１４０１は、マイク１４０１、およびマイク１４０２から供給された音の電気信号をデジタル信号に変換し、入力端子２０１、および入力端子２０２に伝達する。
　表示部１４０２は、波源方向推定装置２００から供給された推定方向情報を画像などの可視化データに変換し、ディスプレイなどの表示装置に表示する。最も基本的な可視化方法は、ある時刻における相関関数を２次元グラフで表示する方法である。その際、横軸に方向、縦軸に相関値を表示する。ある時刻だけでなく、相関関数の時間変化を３次元で表示する方法も有効である。時間変化を表示することにより、目標音源の出現の明確化、目標音源の移動パタン、目標音源の移動方向の予測などが可能になる。３次元ではなく、２次元平面に投影する方法も有効である。３次元だと表示されたときに裏側が見づらいという問題がある。上から投影した平面上に表示すれば、死角が無くなり一覧性が向上する。相関値は、色の濃淡ではなく、等高線で表現してもよい。
　図１５は、本実施形態に係る波源方向推定システム１４００の表示部１４０２で表示した画像の一例を示す図であり、波源方向推定装置２００から供給された推定方向情報から得られたものである。これは、本実施形態の効果を確認する目的で取得した。例の作成には、街頭環境において、時刻２０秒～２５秒に方位３０度で悲鳴が発生する状況の音を用いた。数センチ間隔で設置した二つのマイクを用いて集音した。
　図１５は、色が黒いほど相関値が高いことを表している。方位角の範囲は０から１８０度である。縦軸は時刻を表している。図１５を参照すると、時刻２０秒～２５秒あたりで方位３０度の相関値が高くなっていることが分かる。このことから、悲鳴が時刻２０秒～２５秒に発生していること、悲鳴の発生方向が約３０度であることが分かる。
　本実施形態によれば、推定方向情報を画像などの可視化データとして表示するので、ユーザが波源の方向推定情報を視覚的に把握することができる。
　［他の実施形態］
　以上、実施形態を参照して本願発明を説明したが、本願発明は上記実施形態に限定されるものではない。本願発明の構成や詳細には、本願発明のスコープ内で当業者が理解し得る様々な変更をすることができる。また、それぞれの実施形態に含まれる別々の特徴を如何様に組み合わせたシステムまたは装置も、本発明の範疇に含まれる。
　また、本発明は、複数の機器から構成されるシステムに適用されてもよいし、単体の装置に適用されてもよい。さらに、本発明は、実施形態の機能を実現する情報処理プログラムが、システムあるいは装置に直接あるいは遠隔から供給される場合にも適用可能である。したがって、本発明の機能をコンピュータで実現するために、コンピュータにインストールされるプログラム、あるいはそのプログラムを格納した媒体、そのプログラムをダウンロードさせるＷＷＷ(World Wide Web)サーバも、本発明の範疇に含まれる。特に、少なくとも、上述した実施形態に含まれる処理ステップをコンピュータに実行させるプログラムを格納した非一時的コンピュータ可読媒体（non-transitory computer readable medium）は本発明の範疇に含まれる。
　［実施形態の他の表現］
　上記の実施形態の一部または全部は、以下の付記のようにも記載されうるが、以下には限られない。（付記１）
　波源で発生した波を入力信号として取得する複数の入力信号取得手段と、
　前記入力信号取得手段で取得した複数の前記入力信号を複数の周波数領域信号に変換する変換手段と、
　前記周波数領域信号に基づいて、クロススペクトルを計算するクロススペクトル計算手段と、
　前記クロススペクトルに基づいて、周波数別クロススペクトルを計算する周波数別クロススペクトル計算手段と、
　前記周波数別クロススペクトルに基づいて、統合相関関数を計算する統合相関関数計算手段と、
　を備えた相関関数生成装置。（付記２）
　前記統合相関関数計算手段は、
　前記周波数別クロススペクトルの逆変換により周波数別相関関数を生成する周波数別相関関数生成手段と、
　前記周波数別相関関数を統合して一つの統合相関関数を生成する統合相関関数生成手段と、
　を有する付記１に記載の相関関数生成装置。（付記３）
　前記統合相関関数計算手段は、
　前記周波数別クロススペクトルを統合して統合クロススペクトルを生成する統合クロススペクトル生成手段と、
　前記統合クロススペクトルの逆変換により統合相関関数を生成する統合相関関数生成手段と、
　を有する付記１に記載の相関関数生成装置。（付記４）
　前記周波数別クロススペクトル計算手段は、
　前記クロススペクトルに基づいて、周波数別基本クロススペクトルを計算する周波数別基本クロススペクトル計算手段を有し、
　前記周波数別基本クロススペクトルを前記周波数別クロススペクトルとして求める付記１乃至３のいずれか１項に記載の相関関数生成装置。（付記５）
　前記周波数別クロススペクトル計算手段は、
　前記クロススペクトルに基づいて、周波数別基本クロススペクトルを計算する周波数別基本クロススペクトル計算手段と、
　カーネル関数スペクトルを記憶するカーネル関数記憶手段と、
　前記周波数別基本クロススペクトルと前記カーネル関数スペクトルとを乗じて、前記周波数別クロススペクトルを求める乗算手段と、
　を有する付記１乃至３のいずれか１項に記載の相関関数生成装置。（付記６）
　前記周波数別クロススペクトル計算手段は、
　前記クロススペクトルに基づいて、周波数別基本クロススペクトルを計算する周波数別基本クロススペクトル計算手段と、
　前記クロススペクトルに基づいて、カーネル関数スペクトルを計算するカーネル関数スペクトル計算手段と、
　前記周波数別基本クロススペクトルと前記カーネル関数スペクトルとを乗じて、前記周波数別クロススペクトルを求める乗算手段と、
　を有する付記１乃至３のいずれか１項に記載の相関関数生成装置。（付記７）
　波源で発生した波を入力信号として取得する複数の入力信号取得ステップと、
　前記入力信号取得ステップにおいて取得した複数の前記入力信号を複数の周波数領域信号に変換する変換ステップと、
　前記周波数領域信号に基づいて、クロススペクトルを計算するクロススペクトル計算ステップと、
　前記クロススペクトルに基づいて、周波数別クロススペクトルを計算する周波数別クロススペクトル計算ステップと、
　前記周波数別クロススペクトルに基づいて、統合相関関数を計算する統合相関関数計算ステップと、
　を含む相関関数生成方法。（付記８）
　波源で発生した波を入力信号として取得する複数の入力信号取得ステップと、
　前記入力信号取得ステップにおいて取得した複数の前記入力信号を複数の周波数領域信号に変換する変換ステップと、
　前記周波数領域信号に基づいて、クロススペクトルを計算するクロススペクトル計算ステップと、
　前記クロススペクトルに基づいて、周波数別クロススペクトルを計算する周波数別クロススペクトル計算ステップと、
　前記周波数別クロススペクトルに基づいて、統合相関関数を計算する統合相関関数計算ステップと、
　をコンピュータに実行させる相関関数生成プログラム。（付記９）
　付記１乃至６のいずれか１項に記載の相関関数生成装置と、
　統合相関関数に基づいて、波源の推定方向情報を生成する推定方向情報生成手段と、
　を備えた波源方向推定装置。
　この出願は、２０１６年６月２９日に出願された日本出願特願２０１６－１２８４８６を基礎とする優先権を主張し、その開示の全てをここに取り込む。

Claims (9)

  1. 　波源で発生した波を入力信号として取得する複数の入力信号取得手段と、
    　前記入力信号取得手段で取得した複数の前記入力信号を複数の周波数領域信号に変換する変換手段と、
    　前記周波数領域信号に基づいて、クロススペクトルを計算するクロススペクトル計算手段と、
    　前記クロススペクトルに基づいて、周波数別クロススペクトルを計算する周波数別クロススペクトル計算手段と、
    　前記周波数別クロススペクトルに基づいて、統合相関関数を計算する統合相関関数計算手段と、
    　を備えた相関関数生成装置。
  2. 　前記統合相関関数計算手段は、
    　前記周波数別クロススペクトルの逆変換により周波数別相関関数を生成する周波数別相関関数生成手段と、
    　前記周波数別相関関数を統合して一つの統合相関関数を生成する統合相関関数生成手段と、
    　を有する請求項１に記載の相関関数生成装置。
  3. 　前記統合相関関数計算手段は、
    　前記周波数別クロススペクトルを統合して統合クロススペクトルを生成する統合クロススペクトル生成手段と、
    　前記統合クロススペクトルの逆変換により統合相関関数を生成する統合相関関数生成手段と、
    　を有する請求項１に記載の相関関数生成装置。
  4. 　前記周波数別クロススペクトル計算手段は、
    　前記クロススペクトルに基づいて、周波数別基本クロススペクトルを計算する周波数別基本クロススペクトル計算手段を有し、
    　前記周波数別基本クロススペクトルを前記周波数別クロススペクトルとして求める請求項１乃至３のいずれか１項に記載の相関関数生成装置。
  5. 　前記周波数別クロススペクトル計算手段は、
    　前記クロススペクトルに基づいて、周波数別基本クロススペクトルを計算する周波数別基本クロススペクトル計算手段と、
    　カーネル関数スペクトルを記憶するカーネル関数記憶手段と、
    　前記周波数別基本クロススペクトルと前記カーネル関数スペクトルとを乗じて、前記周波数別クロススペクトルを求める乗算手段と、
    　を有する請求項１乃至３のいずれか１項に記載の相関関数生成装置。
  6. 　前記周波数別クロススペクトル計算手段は、
    　前記クロススペクトルに基づいて、周波数別基本クロススペクトルを計算する周波数別基本クロススペクトル計算手段と、
    　前記クロススペクトルに基づいて、カーネル関数スペクトルを計算するカーネル関数スペクトル計算手段と、
    　前記周波数別基本クロススペクトルと前記カーネル関数スペクトルとを乗じて、前記周波数別クロススペクトルを求める乗算手段と、
    　を有する請求項１乃至３のいずれか１項に記載の相関関数生成装置。
  7. 　波源で発生した波を入力信号として取得する複数の入力信号取得ステップと、
    　前記入力信号取得ステップにおいて取得した複数の前記入力信号を複数の周波数領域信号に変換する変換ステップと、
    　前記周波数領域信号に基づいて、クロススペクトルを計算するクロススペクトル計算ステップと、
    　前記クロススペクトルに基づいて、周波数別クロススペクトルを計算する周波数別クロススペクトル計算ステップと、
    　前記周波数別クロススペクトルに基づいて、統合相関関数を計算する統合相関関数計算ステップと、
    　を含む相関関数生成方法。
  8. 　波源で発生した波を入力信号として取得する複数の入力信号取得ステップと、
    　前記入力信号取得ステップにおいて取得した複数の前記入力信号を複数の周波数領域信号に変換する変換ステップと、
    　前記周波数領域信号に基づいて、クロススペクトルを計算するクロススペクトル計算ステップと、
    　前記クロススペクトルに基づいて、周波数別クロススペクトルを計算する周波数別クロススペクトル計算ステップと、
    　前記周波数別クロススペクトルに基づいて、統合相関関数を計算する統合相関関数計算ステップと、
    　をコンピュータに実行させる相関関数生成プログラム。
  9. 　請求項１乃至６のいずれか１項に記載の相関関数生成装置と、
    　統合相関関数に基づいて、波源の推定方向情報を生成する推定方向情報生成手段と、
    　を備えた波源方向推定装置。
PCT/JP2017/004028 2016-06-29 2017-02-03 相関関数生成装置、相関関数生成方法、相関関数生成プログラムおよび波源方向推定装置 WO2018003158A1 (ja)

Priority Applications (2)

Application Number Priority Date Filing Date Title
JP2018524877A JPWO2018003158A1 (ja) 2016-06-29 2017-02-03 相関関数生成装置、相関関数生成方法、相関関数生成プログラムおよび波源方向推定装置
US16/309,542 US20190250240A1 (en) 2016-06-29 2017-02-03 Correlation function generation device, correlation function generation method, correlation function generation program, and wave source direction estimation device

Applications Claiming Priority (2)

Application Number Priority Date Filing Date Title
JP2016128486 2016-06-29
JP2016-128486 2016-06-29

Publications (1)

Publication Number Publication Date
WO2018003158A1 true WO2018003158A1 (ja) 2018-01-04

Family

ID=60786280

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/JP2017/004028 WO2018003158A1 (ja) 2016-06-29 2017-02-03 相関関数生成装置、相関関数生成方法、相関関数生成プログラムおよび波源方向推定装置

Country Status (3)

Country Link
US (1) US20190250240A1 (ja)
JP (1) JPWO2018003158A1 (ja)
WO (1) WO2018003158A1 (ja)

Cited By (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2020003343A1 (ja) * 2018-06-25 2020-01-02 日本電気株式会社 波源方向推定装置、波源方向推定方法、およびプログラム記録媒体
WO2020003342A1 (ja) * 2018-06-25 2020-01-02 日本電気株式会社 波源方向推定装置、波源方向推定方法、およびプログラム記録媒体
US20220342026A1 (en) * 2019-09-02 2022-10-27 Nec Corporation Wave source direction estimation device, wave source direction estimation method, and program recording medium

Families Citing this family (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US10334360B2 (en) * 2017-06-12 2019-06-25 Revolabs, Inc Method for accurately calculating the direction of arrival of sound at a microphone array

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JPH11344408A (ja) * 1998-06-02 1999-12-14 Hitachi Ltd 音源探査装置
JP2011033717A (ja) * 2009-07-30 2011-02-17 Secom Co Ltd 雑音抑圧装置
JP2012149906A (ja) * 2011-01-17 2012-08-09 Mitsubishi Electric Corp 音源位置推定装置、音源位置推定方法および音源位置推定プログラム
JP2014525037A (ja) * 2011-07-14 2014-09-25 マイクロソフト コーポレーション 位相スペクトルを使った音源定位

Family Cites Families (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US7363172B2 (en) * 2006-01-05 2008-04-22 United States Of America As Represented By The Secretary Of The Navy Method and apparatus for detecting damage in structures
JP4828308B2 (ja) * 2006-05-31 2011-11-30 三菱電機株式会社 位相変調系列再生装置
US8306132B2 (en) * 2009-04-16 2012-11-06 Advantest Corporation Detecting apparatus, calculating apparatus, measurement apparatus, detecting method, calculating method, transmission system, program, and recording medium
JP2012244846A (ja) * 2011-05-23 2012-12-10 Mitsubishi Electric Engineering Co Ltd 環境発電装置、環境発電装置システムおよびセンサ装置
JP2013213739A (ja) * 2012-04-02 2013-10-17 Nippon Telegr & Teleph Corp <Ntt> 音源位置推定装置、音源位置推定方法及びそのプログラム

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JPH11344408A (ja) * 1998-06-02 1999-12-14 Hitachi Ltd 音源探査装置
JP2011033717A (ja) * 2009-07-30 2011-02-17 Secom Co Ltd 雑音抑圧装置
JP2012149906A (ja) * 2011-01-17 2012-08-09 Mitsubishi Electric Corp 音源位置推定装置、音源位置推定方法および音源位置推定プログラム
JP2014525037A (ja) * 2011-07-14 2014-09-25 マイクロソフト コーポレーション 位相スペクトルを使った音源定位

Cited By (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2020003343A1 (ja) * 2018-06-25 2020-01-02 日本電気株式会社 波源方向推定装置、波源方向推定方法、およびプログラム記録媒体
WO2020003342A1 (ja) * 2018-06-25 2020-01-02 日本電気株式会社 波源方向推定装置、波源方向推定方法、およびプログラム記録媒体
JPWO2020003343A1 (ja) * 2018-06-25 2021-06-24 日本電気株式会社 波源方向推定装置、波源方向推定方法、およびプログラム
JPWO2020003342A1 (ja) * 2018-06-25 2021-06-24 日本電気株式会社 波源方向推定装置、波源方向推定方法、およびプログラム
JP7056739B2 (ja) 2018-06-25 2022-04-19 日本電気株式会社 波源方向推定装置、波源方向推定方法、およびプログラム
US11408963B2 (en) 2018-06-25 2022-08-09 Nec Corporation Wave-source-direction estimation device, wave-source-direction estimation method, and program storage medium
US20220342026A1 (en) * 2019-09-02 2022-10-27 Nec Corporation Wave source direction estimation device, wave source direction estimation method, and program recording medium

Also Published As

Publication number Publication date
JPWO2018003158A1 (ja) 2019-05-09
US20190250240A1 (en) 2019-08-15

Similar Documents

Publication Publication Date Title
JP6769495B2 (ja) 相関関数生成装置、相関関数生成方法、相関関数生成プログラムおよび波源方向推定装置
CN107122738B (zh) 基于深度学习模型的无线电信号识别方法及其实现系统
EP2380033B1 (en) Estimating a sound source location using particle filtering
CN106251877B (zh) 语音声源方向估计方法及装置
WO2018003158A1 (ja) 相関関数生成装置、相関関数生成方法、相関関数生成プログラムおよび波源方向推定装置
JP6780644B2 (ja) 信号処理装置、信号処理方法、および信号処理プログラム
Dorfan et al. Tree-based recursive expectation-maximization algorithm for localization of acoustic sources
CN102147458B (zh) 一种针对宽带声源的波达方向估计方法及其装置
EP3227704B1 (en) Method for tracking a target acoustic source
JP7214798B2 (ja) 音声信号処理方法、音声信号処理装置、電子機器及び記憶媒体
WO2017216999A1 (ja) 波源方向推定装置、波源方向推定システム、波源方向推定方法および波源方向推定プログラム
Al-Aboosi et al. Improved underwater signal detection using efficient time–frequency de-noising technique and Pre-whitening filter
Song et al. High-resolution time delay estimation algorithms through cross-correlation post-processing
Hosseini et al. Time difference of arrival estimation of sound source using cross correlation and modified maximum likelihood weighting function
CN103784164A (zh) 超声信号的预处理方法及系统
Zhong et al. Direction of arrival tracking of an underwater acoustic source using particle filtering: Real data experiments
CN114495974B (zh) 音频信号处理方法
CN114814728A (zh) 一种声源定位方法、系统、电子设备及介质
JP5713933B2 (ja) 音源距離測定装置、音響直間比推定装置、雑音除去装置、それらの方法、及びプログラム
CN117198276A (zh) 基于毫米波信号非接触式耳机语音感知方法及系统
Chen et al. A Fast Estimation Method for 3-D Acoustic Source Localization
CN117956376A (en) Audio judgment method and device, electronic equipment and storage medium
Llerena et al. Synchronizing Speech Mixtures in Speech Separation Problems under Reverberant Conditions
CN117809683A (zh) 一种基于注意力机制的单麦克风声学成像方法
Mizumachi Statistical confidence measure for direction-of-arrival estimate

Legal Events

Date Code Title Description
ENP Entry into the national phase

Ref document number: 2018524877

Country of ref document: JP

Kind code of ref document: A

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

Ref document number: 17819528

Country of ref document: EP

Kind code of ref document: A1

NENP Non-entry into the national phase

Ref country code: DE

122 Ep: pct application non-entry in european phase

Ref document number: 17819528

Country of ref document: EP

Kind code of ref document: A1