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

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

Info

Publication number
WO2018131099A1
WO2018131099A1 PCT/JP2017/000680 JP2017000680W WO2018131099A1 WO 2018131099 A1 WO2018131099 A1 WO 2018131099A1 JP 2017000680 W JP2017000680 W JP 2017000680W WO 2018131099 A1 WO2018131099 A1 WO 2018131099A1
Authority
WO
WIPO (PCT)
Prior art keywords
correlation function
frequency
spectrum
cross spectrum
variance
Prior art date
Application number
PCT/JP2017/000680
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 JP2018561142A priority Critical patent/JP6769495B2/ja
Priority to PCT/JP2017/000680 priority patent/WO2018131099A1/ja
Priority to US16/476,897 priority patent/US11336997B2/en
Publication of WO2018131099A1 publication Critical patent/WO2018131099A1/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/02—Direction-finders for determining the direction from which infrasonic, sonic, ultrasonic, or electromagnetic waves, or particle emission, not having a directional significance, are being received using radio waves
    • G01S3/04—Details
    • G01S3/043—Receivers
    • G—PHYSICS
    • G10—MUSICAL INSTRUMENTS; ACOUSTICS
    • G10L—SPEECH ANALYSIS TECHNIQUES OR SPEECH SYNTHESIS; SPEECH RECOGNITION; SPEECH OR VOICE PROCESSING TECHNIQUES; SPEECH OR AUDIO CODING OR DECODING
    • G10L25/00—Speech or voice analysis techniques not restricted to a single one of groups G10L15/00 - G10L21/00
    • G10L25/03—Speech or voice analysis techniques not restricted to a single one of groups G10L15/00 - G10L21/00 characterised by the type of extracted parameters
    • G10L25/18—Speech or voice analysis techniques not restricted to a single one of groups G10L15/00 - G10L21/00 characterised by the type of extracted parameters the extracted parameters being spectral information of each sub-band
    • G—PHYSICS
    • G10—MUSICAL INSTRUMENTS; ACOUSTICS
    • G10L—SPEECH ANALYSIS TECHNIQUES OR SPEECH SYNTHESIS; SPEECH RECOGNITION; SPEECH OR VOICE PROCESSING TECHNIQUES; SPEECH OR AUDIO CODING OR DECODING
    • G10L25/00—Speech or voice analysis techniques not restricted to a single one of groups G10L15/00 - G10L21/00
    • G10L25/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
    • H04R1/00—Details of transducers, loudspeakers or microphones
    • H04R1/20—Arrangements for obtaining desired frequency or directional characteristics
    • H04R1/32—Arrangements for obtaining desired frequency or directional characteristics for obtaining desired directional characteristic only
    • H04R1/40—Arrangements for obtaining desired frequency or directional characteristics for obtaining desired directional characteristic only by combining a number of identical transducers
    • H04R1/406—Arrangements for obtaining desired frequency or directional characteristics for obtaining desired directional characteristic only by combining a number of identical transducers microphones
    • H—ELECTRICITY
    • H04—ELECTRIC COMMUNICATION TECHNIQUE
    • H04R—LOUDSPEAKERS, MICROPHONES, GRAMOPHONE PICK-UPS OR LIKE ACOUSTIC ELECTROMECHANICAL TRANSDUCERS; 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
    • G10—MUSICAL INSTRUMENTS; ACOUSTICS
    • G10L—SPEECH ANALYSIS TECHNIQUES OR SPEECH SYNTHESIS; SPEECH RECOGNITION; SPEECH OR VOICE PROCESSING TECHNIQUES; SPEECH OR AUDIO CODING OR DECODING
    • G10L25/00—Speech or voice analysis techniques not restricted to a single one of groups G10L15/00 - G10L21/00
    • G10L25/03—Speech or voice analysis techniques not restricted to a single one of groups G10L15/00 - G10L21/00 characterised by the type of extracted parameters
    • G10L25/06—Speech or voice analysis techniques not restricted to a single one of groups G10L15/00 - G10L21/00 characterised by the type of extracted parameters the extracted parameters being correlation coefficients
    • 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
    • H—ELECTRICITY
    • H04—ELECTRIC COMMUNICATION TECHNIQUE
    • H04S—STEREOPHONIC SYSTEMS 
    • H04S7/00—Indicating arrangements; Control arrangements, e.g. balance control
    • H04S7/40—Visual indication of stereophonic sound image

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.
  • Patent Document 1 discloses that a cross spectrum is calculated from two acoustic signals, a frequency component having a small time variation of a phase component of the cross spectrum is a speech component, and a portion having a large time variation is a noise component. A technique for suppressing a component is disclosed. Further, in the above technical field, Patent Document 2 detects a direction of arrival of radio waves by calculating a time difference between incoming radio waves from a frequency gradient of a phase in a frequency band having a high coherence level indicated by a cross spectrum of two received signals. Is disclosed.
  • Non-Patent Document 1 and Non-Patent Document 2 disclose a method of estimating the direction of a sound source (sound wave generation source, generation location) using sound reception signals of two microphones. Yes. Specifically, a cross correlation function between the received sound signals is obtained from the two received sound signals. And the arrival direction of a sound wave is estimated by calculating the time difference which a cross correlation function gives the maximum value as the arrival time difference of a sound wave.
  • 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 plurality of input signal acquisition means to obtain a plurality of frequency domain signals; Cross spectrum calculation means for calculating a cross spectrum based on a plurality of the frequency domain signals; A dispersion calculating means for calculating the dispersion of the cross spectrum; A correlation function calculating means for calculating a correlation function based on the cross spectrum and the variance; Is provided.
  • 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 a plurality of the input signals acquired by a plurality of the input signal acquisition means to obtain a plurality of frequency domain signals; A cross spectrum calculating step of calculating a cross spectrum based on the plurality of frequency domain signals; A variance calculating step for calculating the variance of the cross spectrum; A correlation function calculation step of calculating and generating a correlation function based on the cross spectrum and the variance; 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 a plurality of the input signals acquired by a plurality of the input signal acquisition means to obtain a plurality of frequency domain signals; A cross spectrum calculating step of calculating a cross spectrum based on the plurality of frequency domain signals; A variance calculating step for calculating the variance of the cross spectrum; A correlation function calculation step of calculating and generating a correlation function based on the cross spectrum and the variance; Is executed on the computer.
  • a wave source direction estimation device Having the correlation function generating device, Based on the correlation function generated by the correlation function generator, the direction of the wave source is estimated.
  • a correlation function having a clear peak can be generated even in an environment where the ambient noise level is high.
  • 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 having a clear peak.
  • the correlation function generation device 100 includes a plurality of input signal acquisition units 101, a conversion unit 102, a cross spectrum calculation unit 103, a variance calculation unit 104, and a correlation function calculation unit 105. .
  • the input signal acquisition unit 101 acquires a wave generated at the wave source as an input signal.
  • the conversion unit 102 converts the plurality of input signals acquired by the plurality of input signal acquisition units 101 to obtain a plurality of frequency domain signals.
  • the cross spectrum calculation unit 103 calculates a cross spectrum based on a plurality of frequency domain signals.
  • the variance calculation unit 104 calculates the variance of the cross spectrum.
  • the correlation function calculation unit 105 calculates and generates a correlation function based on the cross spectrum and the variance.
  • a correlation function having a clear peak can be generated even in an environment where the ambient noise level is high.
  • FIG. 2 is a block diagram for explaining the configuration of the wave source direction estimation apparatus 200 according to the present embodiment.
  • the wave source direction estimation apparatus 200 of this embodiment functions as a part of an apparatus such as a digital video camera, a smartphone, a mobile phone, a notebook 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.
  • the wave source direction estimation apparatus 200 includes an input terminal 20 1 , an input terminal 20 2 , a conversion unit 201, a cross spectrum calculation unit 202, a variance calculation unit 203, a correlation function calculation unit 204, an estimated direction information generation unit 205, and a relative delay time calculation unit. 206.
  • the conversion unit 201, the cross spectrum calculation unit 202, the variance calculation unit 203, and the correlation function calculation unit 204 constitute a correlation function generation unit 210.
  • the input terminal 20 1 and the input terminal 20 2 a microphone is sound and the sound collecting device of the target sound source (hereinafter, microphone) as sound and digital signal various noises are mixed that occurs around the (sample value sequence) Entered.
  • a sound signal inputted 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 and the input terminal 20 2 at time t x 1 (t), x 2 (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”.
  • the cross spectrum calculation unit 202 calculates a cross spectrum based on the converted signal supplied from the conversion unit 201 and transmits the cross spectrum to the variance calculation unit 203 and the correlation function calculation unit 204.
  • 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 is calculated as follows.
  • conj (X 2 (k, n)) represents a complex conjugate of X 2 (k, n).
  • a cross spectrum normalized by the amplitude component may be used instead of the cross spectrum calculated by Equation (2). Such a cross spectrum is calculated as follows.
  • the variance calculation unit 203 calculates the variance of the cross spectrum using the cross spectrum supplied from the cross spectrum calculation unit 202 and transmits the calculated variance to the correlation function calculation unit 204.
  • variance is calculated for each frequency bin from a cross spectrum input in the past. You may calculate not using a frequency bin unit but the subband unit which bundled several frequency bins, or all the frequency bins. Further, the variance may be calculated by combining them. If the calculation is performed in units of subbands, the number of bands for calculating the variance can be reduced, and the amount of calculation can be reduced.
  • the variance calculated by the variance calculation unit 203 is not the cross spectrum but the phase variance of the cross spectrum.
  • the phase of the cross spectrum is a declination of the cross spectrum when the cross spectrum is regarded as a complex number.
  • the dispersion is an index indicating how far the time series of the phase is scattered with respect to a certain reference value, and is the spread of the phase distribution. Therefore, the variance in the present embodiment is not limited to the definition of variance in probability theory, and the standard deviation is also an example of variance. A typical example of variance follows the definition of probability theory.
  • phase of the cross spectrum S 12 (k, n) in the frequency bin k of the nth frame is arg (S 12 (k, n)
  • the variance V 12 (k, n) of the frequency bin k is calculated as follows.
  • the variance is calculated as follows.
  • V 12 (k, n) is calculated as follows.
  • arg (S 0 ) is a reference phase.
  • an absolute value error may be used instead of the square error.
  • the dispersion V12 (k, n) is calculated as follows.
  • R (k, n) is an average of the cross spectrum and is called a circumferential average in the direction statistics.
  • R (k, n) is calculated as follows.
  • S 12 (k, n) in Expression (9) is limited to “cross spectrum normalized by amplitude component” in Expression (3).
  • leak integration may be used as follows instead of the moving average of the equation (8).
  • ⁇ is a real number and satisfies 0 ⁇ ⁇ 1.
  • V 12 (k, n) is calculated as follows.
  • ln (x) represents the natural logarithm of x.
  • angular deviation known as a statistic similar to the circumferential standard deviation may be used.
  • the method for calculating the variance is not limited to the above-described mathematical formula, and a circumferential average multidimensional function, polynomial, nonlinear function, and a combination thereof may be used.
  • the corrected variance may be transmitted to the correlation function calculation unit 204.
  • the correction is performed using dispersion existing in the periphery on the time-frequency plane (dispersion existing in frames and frequencies around the correction target).
  • a typical example is a method of calculating an average value of dispersion.
  • the corrected variance VV 12 (k, n) is calculated as follows.
  • P 1 and P 2 are predetermined constants.
  • the possibility of the target signal being included in the frame / frequency is low, so the dispersion is increased. Also, if a low value is found due to a prominent dispersion value in the adjacent frequency / frame, the dispersion is lowered.
  • the correlation function calculation unit 204 calculates a correlation function based on the cross spectrum supplied from the cross spectrum calculation unit 202 and the variance supplied from the variance calculation unit 203 and transmits the correlation function to the estimated direction information generation unit 205. Next, the details of the correlation function calculation unit 204 will be described with reference to FIG.
  • FIG. 3 is a block diagram illustrating a configuration of the correlation function calculation unit 204 included in the wave source direction estimation apparatus 200 according to the present embodiment.
  • the correlation function calculation unit 204 includes a weighting unit 241 and an inverse conversion unit 242.
  • Weighting unit 241 weights the cross spectrum supplied from the cross spectrum calculation unit 202 based on the variance supplied from the variance calculation unit 203, and transmits the weighted cross spectrum to the inverse conversion unit 242.
  • the weighting unit 241 reduces the value of the cross spectrum because it is unlikely that the target sound source information is included in the cross spectrum.
  • a typical method is a method of calculating a weight from dispersion using a mapping function prepared in advance and multiplying the weight.
  • the weighted cross spectrum WS 12 (k, n) is calculated as follows.
  • G (k, n) is a weight calculated based on the variance.
  • G (k, n) is calculated as follows.
  • V 12 (k, n) is the variance supplied from the variance calculation unit 203.
  • VV 12 (k, n) which is a corrected dispersion may be used.
  • a and b are real numbers, respectively, and satisfy a> 0. Further, as in the following equation, a method of simply switching values based on threshold determination may be used.
  • V th is a threshold value and is a positive real number.
  • a dispersion function expressed in other forms such as a linear mapping function, a higher-order polynomial function, and a nonlinear function for the calculation of the gain.
  • the inverse transform unit 242 obtains the inverse transform of the cross spectrum supplied from the weighting unit 241.
  • the transform unit 201 uses Fourier transform, a method using inverse Fourier transform for inverse transform will be described. Assuming that the cross spectrum supplied from the weighting unit 241 is WS 12 (k, n), the correlation function s 12 ( ⁇ , n) obtained by the inverse transformation of WS 12 (k, n) is calculated as follows. .
  • 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 the relative delay time 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, if the direction search range is 0 degree to 90 degrees in 10 degree increments, that is, 0 degree, 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.
  • Estimated direction information generation unit 205 Based on the correlation function supplied from the correlation function calculation unit 204 and the relative delay time supplied from the relative delay time calculation unit 206, 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 s 12 ( ⁇ , 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.
  • FIG. 4A is a diagram illustrating an example of a correlation function table 401 included in the wave source direction estimation device 200 according to the present embodiment.
  • Correlation function table 401 stores frequency domain signal 412, cross spectrum 413, variance 414, relative delay time 415, and correlation function 416 in association with input signal 411.
  • the wave source direction estimation apparatus 200 calculates a correlation function 416 for the input signal 411 with reference to the correlation function table 401, for example, and generates estimated direction information.
  • FIG. 4B is a diagram illustrating an example of the conversion function table 402 provided in the wave source direction estimation apparatus 200 according to the present embodiment.
  • the conversion function table 402 stores a conversion function 421 and an inverse conversion function 422 in association with the input signal 411.
  • the wave source direction estimation apparatus 200 converts the input signal 411 with reference to the conversion function table 402, for example, and calculates the frequency domain signal 412, the cross spectrum 413, and the like.
  • FIG. 4C is a diagram illustrating an example of the dispersion table 403 included in the wave source direction estimation device 200 according to the present embodiment.
  • the dispersion table 403 stores the dispersion technique 431 in association with the cross spectrum 424.
  • the wave source direction estimation apparatus 200 calculates the variance 414 of the cross spectrum 424 with reference to the variance table 403, for example.
  • FIG. 5 is a diagram illustrating 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 implements a functional configuration unit of the wave source direction estimation apparatus 200 of FIGS. 2 and 3 by executing a program.
  • a ROM (Read Only Memory) 520 stores fixed data such as initial data and programs and other programs.
  • the network interface 530 communicates with other devices via the network. Note that the number of CPUs 510 is not limited to one, and a plurality of CPUs may be included, or a GPU (Graphics Processing Unit) for image processing may be included.
  • GPU Graphics Processing Unit
  • the network interface 530 preferably includes a CPU independent of the CPU 510 and writes or reads transmission / reception data in a RAM (Random Access Memory) 540 area. Also, it is desirable to provide a DMAC (Direct Memory Access Controller) that transfers data between the RAM 540 and the storage 550 (not shown). Furthermore, 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 network interface 530, the DMAC, or the input / output interface 560.
  • a CPU independent of the CPU 510 and writes or reads transmission / reception data in a RAM (Random Access Memory) 540 area.
  • DMAC Direct Memory Access Controller
  • the RAM 540 is a random access memory used by the CPU 510 as a work area for temporary storage. In the RAM 540, an area for storing data necessary for realizing the present embodiment is secured.
  • Input signal 541, sound input to the input terminal 20 1 and the input terminal 20 2 is a signal such as.
  • the frequency domain signal 542 is a signal obtained by converting the input signal by a predetermined conversion method.
  • Cross spectrum 543 is calculated based on frequency domain signal 542.
  • the variance 544 is the variance of the cross spectrum 543 calculated from the cross spectrum 543.
  • the relative delay time 545 is a difference in arrival time of sound waves and the like determined based on the microphone interval and the sound source direction.
  • Correlation function 546 is calculated from cross spectrum 543 and variance 544. These data are, for example, data developed from the correlation function table 401.
  • the RAM 540 may include data expanded from the conversion function table 402 and the distribution table 403.
  • the input / output data 547 is data input / output via the input / output interface 560.
  • the transmission / reception data 548 is data transmitted / received via the network interface 530.
  • the RAM 540 has an application execution area 549 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 a correlation function table 401, a conversion function table 402, and a distribution table 403.
  • the correlation function table 401 is a table for managing the relationship between the input signal 411 and the correlation function 416 shown in FIG. 4A.
  • the conversion function table 402 is a table for managing the relationship between the input signal 411 and the conversion function 421 shown in FIG. 4B.
  • the dispersion table 403 is a table for managing the relationship between the cross spectrum 424 and the dispersion method 431 shown in FIG. 4C.
  • the storage 550 further stores a conversion module 551, a cross spectrum calculation module 552, a variance calculation module 553, a correlation function calculation module 554, a relative delay time calculation module 555, and an estimated direction information generation module 556.
  • the conversion module 551 is a module that converts an input signal to obtain 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 variance calculation module 553 is a module for calculating the variance of the cross spectrum.
  • the correlation function calculation module 554 is a module that calculates a correlation function from the cross spectrum and the variance.
  • the relative delay time calculation module 555 is a module that calculates the relative delay time between the microphone pairs.
  • the estimated direction information generation module 556 is a module that generates a correspondence relationship between a direction and a correlation value as estimated direction information from the correlation function and the relative delay time.
  • 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 for explaining the processing procedure of the wave source direction estimation apparatus 200 according to the present embodiment. This flowchart is executed by the CPU 510 using the RAM 540, and implements the functional components of the wave source direction estimation apparatus 200 shown in FIGS.
  • step S601 the wave source direction estimation apparatus 200 acquires an input signal.
  • step S603 the wave source direction estimation device 200 converts the input signal to obtain a frequency domain signal.
  • step S605 the wave source direction estimation apparatus 200 calculates a cross spectrum from the frequency domain signal.
  • step S ⁇ b> 607 the wave source direction estimation apparatus 200 calculates the dispersion of the phase of the cross spectrum.
  • step S609 the wave source direction estimation apparatus 200 weights the cross spectrum based on the calculated variance.
  • step S611 the wave source direction estimation apparatus 200 inversely transforms the weighted cross spectrum to calculate a correlation function.
  • step S613, the wave source direction estimation apparatus 200 calculates a relative delay time.
  • step S615 estimated direction information is generated from the 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 and the direction of the wave source can be estimated with high accuracy. 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 where the environmental noise level is high, it is possible to accurately specify the occurrence position of a scream / anger, the appearance position of a suspicious ship / drone, and the like.
  • the difference in arrival time of sound waves to the two microphones, that is, the sound source direction can be obtained from the phase of the cross spectrum. Therefore, the target sound and the environmental sound can be distinguished based on the magnitude of dispersion obtained from the time series of phases (in this embodiment, the phase of the past L-1 frame). Since the time change in the direction of the target sound source is small, the phase dispersion is small. Even if the target sound generation time is about several seconds, for example, when the phase is obtained every 10 milliseconds, the dispersion of the phase obtained from the several seconds becomes small. On the other hand, since it is difficult to specify the direction of the environmental noise, the time change of the direction becomes large and the dispersion of the phase also becomes large. Since the variance is calculated for each frequency, it is possible to specify the frequency band in which the target sound source exists.
  • FIG. 7 is a block diagram for explaining the configuration of the wave source direction estimation apparatus 700 according to this embodiment.
  • the wave source direction estimation apparatus 700 according to the present embodiment is different from the second embodiment in that an average calculation unit 701 is provided. 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 conversion unit 201, the cross spectrum calculation unit 202, the variance calculation unit 203, the correlation function calculation unit 204, and the average calculation unit 701 constitute a correlation function generation unit 710.
  • the average calculation unit 701 calculates the average of the cross spectrum supplied from the cross spectrum calculation unit 202 and transmits the average to the correlation function calculation unit 204.
  • the average calculation unit 701 calculates the average of the cross spectrum supplied from the cross spectrum calculation unit 202 and transmits the average to the correlation function calculation unit 204.
  • an average is calculated for each frequency bin from a cross spectrum input in the past. You may calculate not the frequency bin unit but the subband unit which bundled the several frequency bin. If the cross spectrum in the frequency bin k of the nth frame is S 12 (k, n), the average cross spectrum SS 12 (k, n) obtained from the past L ⁇ 1 frames is calculated as follows.
  • leak integration may be used as follows.
  • ⁇ is a real number and satisfies 0 ⁇ ⁇ 1.
  • the average cross spectrum is transmitted to the correlation function calculation unit 204 instead of the cross spectrum.
  • the target sound source whose position temporal variation is small
  • the peak position of the correlation function and the temporal variation of the estimated direction information are small, so the accuracy of direction estimation is improved.
  • the cross spectrum averaging effect will be high, so the correlation function Peak position and estimated direction information become clear.
  • FIG. 8 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 is different from the third embodiment in that a dispersion calculation unit 801 is provided instead of the dispersion calculation unit 203. 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.
  • the conversion unit 201, the cross spectrum calculation unit 202, the variance calculation unit 801, the average calculation unit 701, and the correlation function calculation unit 204 constitute a correlation function generation unit 810.
  • the variance calculation unit 801 calculates the variance using the average cross spectrum supplied from the average calculation unit 301 and transmits the variance to the correlation function calculation unit 204.
  • the difference is that an average cross spectrum is used instead of a cross spectrum for calculation of dispersion.
  • the calculation of the circumferential average can be omitted in the calculation process. Assuming that the average cross spectrum is SS 12 (k, n), the dispersion V 12 (k, n) is calculated as follows when circumferential dispersion is used in the calculation of “dispersion of phase of cross spectrum”.
  • V 12 (k, n) is calculated as follows.
  • the variance calculation unit calculates the variance using the average cross spectrum calculated by the average cross spectrum calculation unit. For this reason, it is not necessary to calculate the average cross spectrum in the variance calculation unit. Therefore, the variance can be calculated with a smaller calculation amount than in the third embodiment.
  • FIG. 9 is a block diagram for explaining the configuration of the correlation function calculation unit 904 according to this embodiment.
  • the correlation function calculation unit 904 according to the present embodiment replaces the weighting unit 241 and the inverse conversion unit 242 with the frequency-specific cross spectrum generation units 941 1 and 941 2. ,... 941 K and the integrated correlation function calculation unit 942. 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.
  • FIG. 9 is a block diagram of the correlation function calculation unit 904.
  • the correlation function calculation unit 904 includes frequency-specific cross spectrum generation units 941 1 , 941 2 ,..., 941 K and an integrated correlation function calculation unit 942.
  • the frequency-specific cross spectrum generation units 941 1 , 941 2 ,..., 941 K include the cross spectrum S 12 (k, n) supplied from the cross spectrum calculation unit 202 and the variance p ( k, n) is used to calculate a cross spectrum corresponding to each frequency k of S 12 (k, n), and is transmitted to the integrated correlation function calculation unit 942 as a frequency-specific cross spectrum.
  • the cross spectrum for each frequency is generated in order to calculate a correlation function for each frequency component. That is, the 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) in the subsequent stage.
  • Correlation function calculation section 904 integrates the correlation function calculating unit 942 and the frequency-cross spectrum generation unit 941 1, 941 2, ..., and calculates an integrated correlation function based on frequency-cross spectrum supplied from 941 K, the estimated This is transmitted to the direction information generation unit 205.
  • the frequency-specific cross spectrum generation unit 941 k that calculates the frequency-specific cross spectrum of a certain frequency k will be described in detail with reference to FIG. Then, with reference to FIG. 11, the detail of the integrated correlation function calculation part 942 is demonstrated.
  • FIG. 10 is a block diagram of the frequency-specific cross spectrum generation unit 941 k .
  • the frequency-specific cross spectrum generation unit 941 k includes a frequency-specific cross spectrum generation unit 1011, a kernel function spectrum generation unit 1012, and a multiplication unit 1013.
  • the frequency-specific cross spectrum generation unit 941 k uses the cross spectrum S 12 (k, n) supplied from the cross spectrum calculation unit 202 and the variance supplied from the variance calculation unit 203 to calculate S 12 (k, n).
  • the cross spectrum corresponding to the frequency k is calculated, and the cross spectrum for each frequency is transmitted to the integrated correlation function calculation unit 942.
  • the frequency-specific cross spectrum generation units 941 1 , 941 2 ,..., 941 K use the cross spectrum S 12 (k, n) supplied from the cross spectrum calculation unit 202 to generate S 12 (k, n).
  • a cross spectrum corresponding to each frequency k is calculated and transmitted to the multiplier 1013 as a frequency-specific cross spectrum.
  • the cross spectrum by frequency is performed in order to calculate a correlation function for each frequency component. That is, the 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) in the subsequent stage.
  • the frequency-specific cross spectrum generation unit 941 k that calculates the frequency-specific cross spectrum of a certain frequency k will be described in detail.
  • the frequency-specific cross spectrum generation unit 941 k obtains the phase component and the amplitude component separately and then integrates them. If the cross spectrum Uk (w, n) for each frequency k is assumed, its amplitude component is
  • 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 cross spectrum by frequency from the cross spectrum S 12 (k, n) of the frequency k explain.
  • 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
  • phase component arg U k (w, n)
  • a 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)), 2arg (S 12 (k, n)), 3arg (S 12 (k, n)), 4arg (S 12 (k, n)) is used.
  • phase component of the frequency that is a non-constant multiple of the frequency k is set to zero. Therefore, the phase component arg (U k (w, n)) of the cross spectrum by frequency 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 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.
  • the kernel function spectrum generation unit 1012 obtains a kernel function spectrum based on the variance supplied from the variance calculation unit 203 and outputs the kernel function spectrum to the multiplication unit 1013.
  • 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. Alternatively, the absolute value may be squared.
  • the kernel function spectrum is described as G (w) and the kernel function is expressed as g ( ⁇ ).
  • a Gaussian function is used as the kernel function. At this time, the Gaussian function is given by the following equation.
  • g 1 , g 2 , and g 3 are positive real numbers.
  • 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 (31), spread of the Gaussian function is greater the greater the g 3.
  • g 1 and g 2 are positive 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.
  • g 3 and g 5 that affect the spread of the kernel function are determined based on the variance supplied from the variance calculation unit 203.
  • the basic calculation method of the spread control parameter is a method of converting a dispersion value into a spread control parameter using a mapping function set in advance. For example, if the variance exceeds a certain threshold, the spread control parameter is set to a large value (for example, 10), and if the variance is below the threshold, the small value (for example, 0.01) is set. In this case, when the variance is V 12 (k, n) and the threshold is p th , the spread control parameter q (k, n) in the frequency bin k of the nth frame is calculated as follows.
  • q 1 and q 2 are positive real numbers and satisfy q 1 > q 2 . Moreover, you may calculate using a linear function as follows.
  • q 3 and q 4 are real numbers, respectively, and satisfy q 3 > 0.
  • dispersion functions represented in other forms such as a linear mapping function, a higher-order polynomial function, and a nonlinear function can also be used for calculation of dispersion.
  • the function for obtaining the spread control parameter may be a function of frequency k as well as dispersion. For example, a function that decreases as the frequency increases. An example in which the inverse of k is used as such a representative example will be described. In this case, the spread control parameter q (k, n) is calculated using the following function instead of the equation (33).
  • the spread control parameter q (k, n) is calculated using the following function.
  • the multiplication unit 1013 calculates the product of the frequency-specific cross spectrum supplied from the frequency-specific cross spectrum generation unit 1011 and the kernel function spectrum supplied from the kernel function spectrum generation unit 1012, and obtains a new value obtained by the calculation.
  • the frequency cross spectrum is transmitted to the integrated correlation function calculation unit 942.
  • the frequency-specific spectrum supplied from the frequency-specific cross spectrum generation unit 1011 is U k (w, n) and the kernel function spectrum supplied from the kernel function spectrum generation unit 1012 is G (w)
  • a new one obtained by calculation is obtained.
  • the frequency-specific cross spectrum UM k (w, n) is calculated as follows.
  • FIG. 11 is a block diagram of the integrated correlation function calculation unit 942.
  • the integrated correlation function calculation unit 942 includes inverse transformation units 1121 1 , 1121 2 ,..., 1121 K, and an integration unit 1122.
  • Inverse transform unit 1121 1, 1121 2, ..., 1121 K performs frequency-cross spectrum generation unit 941 1, 941 2, ..., an inverse transform of the frequency-cross spectrum supplied from 941 K, the frequency Each is transmitted to the integration unit 1122 as another correlation function.
  • 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 generation unit 941 k is UM k (w, n)
  • the frequency-specific correlation function u k ( ⁇ , n) obtained by the inverse transformation of UM k (w, n).
  • the integrating unit 1122 integrates the frequency-specific correlation functions supplied from the inverse converting units 1121 1 , 1121 2 ,..., 1121 K , and transmits them to the estimated direction information generating 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 1122 calculates the total 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 the 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 frequency-specific cross spectrum obtained by the frequency-specific cross spectrum generation unit 1011 is based on the cross spectrum of a certain frequency k, “the phase component arg (S 12 (K, n)) multiplied by p is defined as "assigned.”
  • p is an integer of 1 or more. That is, the frequency-specific cross spectrum is defined such that the phase component arg (Uk (w, n)) satisfies at least the following expression.
  • p 1 and 2 and 3
  • p 2 and 3.
  • the direction estimation accuracy is equivalent to that of the prior art, and high accuracy of direction estimation cannot be achieved.
  • p 1, 2, 3, 4,...
  • the number of p increases, that is, if the allocation to the phase component of frequency pk increases, the peak of the correlation function by frequency becomes Since it becomes sharper, the accuracy of direction estimation is improved.
  • FIG. 12B shows the relationship between the frequency-specific cross spectrum multiplied by the kernel function spectrum and the frequency-specific correlation function. For comparison, the cross spectrum by frequency before multiplication by the kernel function spectrum is also shown. As shown in the left side view of FIG. 12B, if the kernel function spectrum is not multiplied, components exist up to a high frequency, and thus the peak of the correlation function for each frequency becomes sharp. On the other hand, when 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 sharpness of the peak of the correlation function for each frequency becomes small. That is, as the peak of the kernel function spectrum becomes sharper (the bottom of the kernel function spectrum becomes narrower), the peak of the correlation function for each frequency becomes smaller. Further, as shown in the right side of FIG. 12B, when the skirt of the correlation function for each frequency is greatly widened, the skirts of adjacent peaks overlap and a correlation function for each frequency having a shallow valley is obtained.
  • FIG. 12C is a diagram illustrating a relationship between the presence / absence of the kernel function and the integrated correlation function.
  • the peak positions of the correlation functions u 1 ( ⁇ , n) to u 3 ( ⁇ , n) by frequency are close, but u 1 ( ⁇ , n) to u Since 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) have a large peak due to integration. Can be formed. For this reason, the position of the peak becomes clear as compared with the case without the kernel function (1210).
  • FIG. 12D is a diagram showing the relationship between the difference in the width of the kernel function spectrum and the integrated correlation function.
  • FIG. 12B when a wide kernel function spectrum is used, a correlation function for each frequency having a shallow valley is formed due to the periodicity of the correlation function. Therefore, as shown by 1230 in FIG. 12D, 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. 12D is a diagram showing the relationship between the difference in the width of the kernel function spectrum and the integrated correlation function.
  • the product of the kernel function spectrum obtained by the Fourier transform of the kernel function and the cross spectrum by frequency is calculated.
  • the product 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 inverse conversion unit 1121 k in the integrated correlation function calculation unit 942 is provided for each frequency supplied from the inverse conversion unit 1121 k.
  • a kernel function may be convoluted with the correlation function.
  • 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 peak of the correlation function by frequency obtained by inversely transforming the cross spectrum by frequency becomes sharp, and the peak position of the correlation function becomes clear.
  • the accuracy of the sound source direction estimation is improved.
  • the shape of the frequency-specific correlation function can be changed by multiplying the frequency-specific cross spectrum by the kernel function spectrum.
  • the width of the kernel function is changed according to the magnitude of the variance, the peak of the correlation function for each frequency of the frequency where the target sound source exists can be emphasized. Therefore, the peak of the correlation function can be enhanced compared with the second embodiment in which the correlation function for each frequency is not obtained and the peak enhancement of the correlation function for each frequency is not performed, so that the direction estimation accuracy is improved.
  • FIG. 13 is a block diagram for explaining the configuration of the frequency-specific cross spectrum generation unit 1341 k according to the present embodiment.
  • the frequency-specific cross spectrum generation unit 1341 k according to the present embodiment replaces the kernel function spectrum generation unit 1012 with a kernel function spectrum storage unit 1314.
  • the difference is that it has a kernel function spectrum selector 1312. Since other configurations and operations are the same as those of the fifth embodiment, the same configurations and operations are denoted by the same reference numerals, and detailed description thereof is omitted.
  • FIG. 13 is a block diagram of the frequency-specific cross spectrum generation unit 1341 k .
  • the frequency-specific cross spectrum generation unit 1341 k includes a frequency-specific cross spectrum generation unit 1011, a kernel function spectrum selection unit 1312, a kernel function spectrum storage unit 1314, and a multiplication unit 1013.
  • the frequency-specific cross spectrum generation unit 1011 k uses the cross spectrum S 12 (k, n) supplied from the cross spectrum calculation unit 202 and the variance supplied from the dispersion calculation unit 203 to calculate S 12 (k, n).
  • the cross spectrum corresponding to the frequency k is calculated, and the cross spectrum for each frequency is transmitted to the integrated correlation function calculation unit 942.
  • the kernel function spectrum storage unit 1314 stores a kernel function spectrum and transmits the kernel function spectrum to the kernel function spectrum selection unit 1312.
  • the kernel function spectrum storage unit 1314 stores a plurality of kernel function spectra corresponding to different variances. For example, kernel function spectra corresponding to 10 types of dispersion (0.1, 0.2,..., 1.0) are stored in increments of 0.1 to 0.1. That is, kernel function spectra are calculated in advance for all variance values from 0.1 to 1.0, and stored in the kernel function spectrum storage unit 1314. In this case, 10 types of kernel function spectra are stored. In the present embodiment, G 0.1 (w), G 0.2 (w),... G 1.0 (kernel function spectra having variances of 0.1 , 0.2 ,. w).
  • a kernel function spectrum corresponding to many dispersion values is prepared, for example, 10,000 kinds of dispersions (0.001, 0.002, 0.003,...) In increments of 0.001 to 0.001. It is desirable to store a kernel function spectrum corresponding to. However, there arises a problem that the storage capacity increases. Therefore, in order to generate an integrated correlation function with a clear peak while suppressing an increase in storage capacity, a kernel corresponding to about 5 to 20 types of variance in increments of 0.05 to 0.05 or 0.1 It is necessary to prepare a function spectrum. Note that the step size does not necessarily have to be uniform.
  • a kernel function spectrum corresponding to one type of variance value having a sufficiently large value may be prepared.
  • the kernel function spectrum selection unit 1312 selects one kernel function spectrum from a plurality of kernel function spectra stored in the kernel function spectrum storage unit 1314 based on the variance supplied from the variance calculation unit 203, This is transmitted to the multiplication unit 1013.
  • the kernel function spectrum selection unit 1312 selects a kernel function spectrum generated with a variance value closest to the variance value supplied from the variance calculation unit 203. For example, if the supplied dispersion value is 0 to less than 0.15, the kernel function spectrum G 0.1 (w) is selected. Similarly, if the supplied dispersion value is not less than 0.15 and less than 0.25, the kernel function spectrum G 0.2 (w) is selected.
  • FIG. 14 is a block diagram for explaining the configuration of the integrated correlation function calculation unit 1442 according to this embodiment.
  • the integrated correlation function calculation unit 1442 according to the present embodiment includes an inverse transformation unit 1121 1 , 1121 2 ,..., 1121 K and an integration unit 1122 compared to the integrated correlation function calculation unit 942 of the fifth embodiment. Instead, it is different in that it has an integration unit 1421 and an inverse conversion unit 1422. Since other configurations and operations are the same as those of the fifth embodiment, the same configurations and operations are denoted by the same reference numerals, and detailed description thereof is omitted.
  • FIG. 14 is a block diagram of the integrated correlation function calculation unit 1442.
  • the integrated correlation function calculation unit 1442 includes an integration unit 1421 and an inverse conversion unit 1422.
  • the integration unit 1421 integrates the frequency-specific cross spectrums supplied from the frequency-specific cross spectrum generation units 941 1 , 941 2 ,..., 941 K , and transmits them to the inverse conversion unit 1422 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 1122 of the fifth embodiment.
  • the frequency-specific cross spectrum supplied from the frequency-specific cross spectrum generation unit 941 k is expressed as UM k (0, n), UM k (1, n),..., UM k (W ⁇ 1, n), the integrated cross spectrum U (k, n) is calculated as follows.
  • the integrated cross spectrum U (k, n) is calculated as follows.
  • the integrated cross spectrum U (k, n) is generated. You may correct to.
  • the degree of influence is controlled in the form of frequency selection and weighting as in the fifth embodiment. For example, when the set of frequencies in which the target sound exists is ⁇ , the calculation is performed as follows when obtaining the integrated cross spectrum U (k, n) 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 inverse transform unit 1422 performs inverse transform of the integrated cross spectrum supplied from the integration unit 1421 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. If the integrated cross spectrum supplied from the integration unit 1421 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 5th 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 fifth embodiment.
  • FIG. 15A is a diagram for describing a configuration of a wave source direction estimation system 1500 according to the present embodiment.
  • the wave source direction estimation system 1500 according to the present embodiment uses the wave source direction estimation device 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 1500 includes microphones 150 1 and 150 2 , an AD conversion unit 1501, and a display unit 1502.
  • wave source direction estimation devices 700 and 800 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.
  • the microphones 150 1 and 150 2 convert sound around the device including sound generated from the target object that is an estimation target into an electric signal, and transmits the electric signal to the AD conversion unit 1501.
  • 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 converter 1501 converts the sound electrical signal supplied from the microphones 150 1 and 150 2 into a digital signal and transmits the digital signal to the input terminals 20 1 and 20 2 .
  • the display unit 1502 converts the estimated direction information supplied from the wave source direction estimation device 200 into visualization data such as an image and displays the visualization data on a display device such as a display.
  • a 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. 15B is a diagram illustrating an example of an image 1510 displayed on the display unit 1502 of the wave source direction estimation system 1500 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.
  • a sound in a situation where a drone flying sound is generated at an orientation of 40 degrees in an outdoor environment was used. Sound was collected using two microphones installed at intervals of several centimeters.
  • FIG. 15B shows that the correlation value is higher as the color is blacker.
  • the range of azimuth is 0 to 180 degrees.
  • the horizontal axis represents time. Referring to FIG. 15B, it can be seen that the correlation value at the orientation of 40 degrees is high. From this, it can be seen that the arrival direction of the drone is about 40 degrees.
  • the direction of the wave source can be estimated with high accuracy. Moreover, since the estimated direction information is displayed as visualization data such as an image, 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 part or all of the above-described embodiment can be described as in the following supplementary notes, but is not limited thereto.
  • (Appendix 1) 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 plurality of input signal acquisition means to obtain a plurality of frequency domain signals; Cross spectrum calculation means for calculating a cross spectrum based on a plurality of the frequency domain signals; A dispersion calculating means for calculating the dispersion of the cross spectrum; A correlation function calculating means for calculating a correlation function based on the cross spectrum and the variance; A correlation function generation device comprising: (Appendix 2) The correlation function generation device according to attachment 1, wherein the variance calculation means calculates variance based on the phase of the cross spectrum.
  • the correlation function generation device (Appendix 3) The correlation function generation device according to attachment 1, wherein the variance calculation means calculates a circumferential variance of the cross spectrum. (Appendix 4) The correlation function generating device according to claim 1, wherein the variance calculating unit calculates a variance based on a cross spectrum calculated based on a frequency domain signal converted from an input signal acquired in the past by the input signal acquiring unit. . (Appendix 5) The correlation function calculating means includes Weighting means for weighting the cross spectrum based on the variance; 5. The correlation function generation device according to any one of supplementary notes 1 to 4, wherein a correlation function is calculated based on the weighted cross spectrum. (Appendix 6) Average calculating means for calculating an average cross spectrum by averaging the cross spectra; 6.
  • the correlation function generation device includes A frequency-specific cross spectrum generation means for generating a plurality of first frequency-specific cross spectra based on the dispersion and the cross spectrum; An integrated correlation function calculation means for calculating a correlation function by integrating a plurality of cross spectra by the first frequency;
  • the correlation function generation device according to any one of appendices 1 to 6, further comprising: (Appendix 8)
  • the integrated correlation function calculating means includes A first inverse transform means for calculating a plurality of first frequency-specific cross spectra and obtaining a plurality of frequency-specific correlation functions; First integration means for calculating a correlation function by integrating the plurality of frequency-specific correlation functions calculated by the first inverse transform means;
  • the correlation function generator according to appendix 7, comprising: (Appendix 9)
  • the frequency-specific cross spectrum generation means includes: Kernel function spectrum generating means for generating a
  • the integrated correlation function calculating means includes A second integration means for integrating a plurality of first frequency-specific cross spectra to obtain one integrated cross spectrum; A second inverse transform means for calculating an inverse transform of the one integrated cross spectrum to obtain one correlation function;
  • the correlation function generation device according to any one of appendices 7 to 10, comprising: (Appendix 12) A plurality of input signal acquisition steps for acquiring waves generated by the wave source as input signals; A conversion step of converting a plurality of the input signals acquired by a plurality of the input signal acquisition means to obtain a plurality of frequency domain signals; A cross spectrum calculating step of calculating a cross spectrum based on the plurality of frequency domain signals; A variance calculating step for calculating the variance of the cross spectrum; A correlation function calculation step of calculating and generating a correlation function based on the cross spectrum and the variance;
  • a correlation function generation method including: (Appendix 13) A plurality of input signal acquisition steps for acquiring waves generated by the wave source as input signals; A conversion step of converting converting

Landscapes

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

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.
　しかしながら、上記文献に記載の技術では、周囲のノイズレベルが高い環境の場合には、明瞭なピークを持つ相関関数が生成できなかった。
　本発明の目的は、上述の課題を解決する技術を提供することにある。
　上記目的を達成するため、本発明に係る相関関数生成装置は、
　波源で発生した波を入力信号として取得する複数の入力信号取得手段と、
　複数の前記入力信号取得手段で取得した複数の前記入力信号を変換して、複数の周波数領域信号を得る変換手段と、
　複数の前記周波数領域信号に基づいて、クロススペクトルを計算するクロススペクトル計算手段と、
　前記クロススペクトルの分散を計算する分散計算手段と、
　前記クロススペクトルと前記分散とに基づいて、相関関数を計算する相関関数計算手段と、
　を備える。
　上記目的を達成するため、本発明に係る相関関数生成方法は、
　波源で発生した波を入力信号として取得する複数の入力信号取得ステップと、
　複数の前記入力信号取得手段で取得した複数の前記入力信号を変換して、複数の周波数領域信号を得る変換ステップと、
　複数の前記周波数領域信号に基づいて、クロススペクトルを計算するクロススペクトル計算ステップと、
　前記クロススペクトルの分散を計算する分散計算ステップと、
　前記クロススペクトルと前記分散とに基づいて、相関関数を計算して生成する相関関数計算ステップと、
　を含む。
　上記目的を達成するため、本発明に係る相関関数生成プログラムは、
　波源で発生した波を入力信号として取得する複数の入力信号取得ステップと、
　複数の前記入力信号取得手段で取得した複数の前記入力信号を変換して、複数の周波数領域信号を得る変換ステップと、
　複数の前記周波数領域信号に基づいて、クロススペクトルを計算するクロススペクトル計算ステップと、
　前記クロススペクトルの分散を計算する分散計算ステップと、
　前記クロススペクトルと前記分散とに基づいて、相関関数を計算して生成する相関関数計算ステップと、
　をコンピュータに実行させる。
　上記目的を達成するため、本発明に係る波源方向推定装置は、
　上記相関関数生成装置を有し、
　前記相関関数生成装置により生成された相関関数に基づいて、前記波源の方向を推定する。
　本発明によれば、周囲のノイズレベルが高い環境の場合であっても、明瞭なピークを持つ相関関数を生成できる。
本発明の第１実施形態に係る相関関数生成装置の構成を示すブロック図である。 本発明の第２実施形態に係る相関関数生成部を含む波源方向推定装置の構成を説明するためのブロック図である。 本発明の第２実施形態に係る波源方向推定装置の備える相関関数計算部の構成を示すブロック図である。 本発明の第２実施形態に係る波源方向推定装置の備える相関関数テーブルの一例を説明する図である。 本発明の第２実施形態に係る波源方向推定装置の備える変換関数テーブルの一例を説明する図である。 本発明の第２実施形態に係る波源方向推定装置の備える分散テーブルの一例を説明する図である。 本発明の第２実施形態に係る波源方向推定装置のハードウェア構成を説明するブロック図である。 本発明の第２実施形態に係る波源方向推定装置の処理手順を説明するフローチャートである。 本発明の第３実施形態に係る相関関数生成部を含む波源方向推定装置の構成を説明するためのブロック図である。 本発明の第４実施形態に係る相関関数生成部を含む波源方向推定装置の構成を説明するためのブロック図である。 本発明の第５実施形態に係る波源方向推定装置の備える相関関数計算部の構成を説明するためのブロック図である。 本発明の第５実施形態に係る波源方向推定装置の備える周波数別クロススペクトル生成の構成を説明するためのブロック図である。 本発明の第５実施形態に係る波源方向推定装置の備える統合相関関数計算部の構成を説明するためのブロック図である。 本発明の第５実施形態に係る波源方向推定装置により計算された周波数別クロススペクトルおよび周波数別相関関数を説明する図である。 本発明の第５実施形態に係る波源方向推定装置によりカーネル関数スペクトルが乗じられた周波数別クロススペクトルと周波数別相関関数との関係を説明する図である。 本発明の第５実施形態に係る波源方向推定装置においてカーネル関数の有無と統合相関関数との関係を説明する図である。 本発明の第５実施形態に係る波源方向推定装置においてカーネル関数スペクトルの幅の違いと統合相関関数との関係を説明する図である。 本発明の第６実施形態に係る波源方向推定装置の備える周波数別クロススペクトル生成の構成を説明するブロック図である。 本発明の第７実施形態に係る統合相関関数計算部の構成を説明するためのブロック図である。 本発明の第８実施形態に係る波源方向推定システムの構成を説明するための図である。 本発明の第８実施形態に係る波源方向推定システムの表示部で表示した画像の一例を示す図である。
　以下に、本発明を実施するための形態について、図面を参照して、例示的に詳しく説明記載する。ただし、以下の実施の形態に記載されている、構成、数値、処理の流れ、機能要素などは一例に過ぎず、その変形や変更は自由であって、本発明の技術範囲を以下の記載に限定する趣旨のものではない。
　また、以下の実施形態に係る波源方向推定装置の推定対象は、空気や水の振動波である音波の発生源に限定されない。地震や地滑りなどの土や固体を媒質とする振動波の発生源にも適用できる。その場合、振動波を電気信号に変換する装置には、マイクロフォンではなく振動センサが用いられる。さらに、気体・液体・固体の振動波だけでなく、電波を用いて方向を推定する場合にも以下の実施形態に係る波源方向推定装置を適用できる。その場合、電波を電気信号に変換する装置にはアンテナが用いられる。以下の実施形態においては、波源は音源と仮定して説明する。
　［第１実施形態］
　本発明の第１実施形態としての相関関数生成装置１００について、図１を用いて説明する。相関関数生成装置１００は、明瞭なピークを持つ相関関数を生成する装置である。
　図１に示すように、相関関数生成装置１００は、複数の入力信号取得部１０１と、変換部１０２と、クロススペクトル計算部１０３と、分散計算部１０４と、相関関数計算部１０５と、を含む。
　入力信号取得部１０１は、波源で発生した波を入力信号として取得する。変換部１０２は、複数の入力信号取得部１０１で取得した複数の入力信号を変換して、複数の周波数領域信号を得る。クロススペクトル計算部１０３は、複数の周波数領域信号に基づいて、クロススペクトルを計算する。分散計算部１０４は、クロススペクトルの分散を計算する。相関関数計算部１０５は、クロススペクトルと分散とに基づいて、相関関数を計算して生成する。
　本実施形態によれば、周囲のノイズレベルが高い環境の場合であっても、明瞭なピークを持つ相関関数を生成できる。
　［第２実施形態］
　＜前提技術＞
　上記非特許文献１および非特許文献２に記載の技術では、屋外などの周囲騒音レベルが高い環境において、遠方に存在する音源の方向を高精度に推定することが困難であった。例えば、推定対象の音源（目標音源）が、マイクロフォンから遠く離れた場所に存在する場合、目標音源から放射される音の音量は、マイクロフォンに到達する時点で大幅に小さくなる。このため、目標音源の音が周囲環境雑音に埋もれて、明瞭なピークを持つ相関関数を生成できなかった。よって、目標音源の方向推定精度が低下することがあった。
　＜本実施形態の技術＞
　次に本発明の第２実施形態としての波源方向推定装置２００について、図２乃至図６を用いて説明する。
　図２は、本実施形態に係る波源方向推定装置２００の構成を説明するためのブロック図である。本実施形態の波源方向推定装置２００は、例えば、デジタルビデオカメラ、スマートフォン、携帯電話、ノートパソコン、パッシブソーナーなどといった装置の一部として機能する。また、不審ドローン検知、悲鳴検知、車両事故検知などの声や音に基づいて異常を検知する異常音検知装置にも搭載される。しかし、本実施形態に係る波源方向推定装置２００の適用例はこれらに限定されるものではなく、受信音から目標音源の方向推定を要求されるあらゆる波源方向推定装置に適用可能である。
　波源方向推定装置２００は、入力端子２０１、入力端子２０２、変換部２０１、クロススペクトル計算部２０２、分散計算部２０３、相関関数計算部２０４、推定方向情報生成部２０５および相対遅延時間計算部２０６を備える。なお、変換部２０１、クロススペクトル計算部２０２、分散計算部２０３および相関関数計算部２０４は、相関関数生成部２１０を構成する。
　（入力端子２０）
　入力端子２０１と入力端子２０２とには、目標音源の音と集音装置であるマイクロフォン(以下、マイク)の周囲で生じる様々な雑音が混在した音とがデジタル信号(サンプル値系列)として入力される。入力端２０１と入力端子２０２とに入力された音信号を本実施形態では入力信号と呼ぶ。そして、時刻tにおける入力端子２０１と入力端子２０２との入力信号をｘ１（ｔ）、ｘ２（ｔ）と表す。
　入力端子に入力される音は、集音装置であるマイクで集音される。入力端子は複数存在するので、目標音源の音を集音する場合には、端子数と同じ２個のマイクが同時に使用される。本実施形態では、入力端子とマイクとは一対一に対応しており、ｍ番目のマイクが集音した音は、ｍ番目の入力端子に供給されると仮定する。したがって、ｍ番目の入力端子に供給された入力信号のことを「ｍ番目のマイクの入力信号」とも呼ぶ。
　波源方向推定装置２００は、目標音源の音が２つのマイクに到達する時間差を利用して音源の方向を推定する。このため、マイク間隔も重要な情報となるので、波源方向推定装置２００には、入力信号だけでなく、マイク位置情報も供給される。
　（変換部２０１）
　変換部２０１は、入力端子２０１と入力端子２０２とから供給された入力信号を変換し、クロススペクトル計算部２０２へ供給する。変換は、入力信号を複数の周波数成分に分解する目的で実施される。ここでは、代表的なフーリエ変換を用いた場合について説明する。
　変換部２０１には２種類の入力信号ｘｍ（ｔ）が入力される。ここで、ｍは入力端子番号である。変換部２０１は、入力端子から供給された入力信号から、適当な長さの波形を一定の周期でずらしながら切り出す。こうして切り出した信号区間をフレーム、切り出した波形の長さをフレーム長、フレームをずらす周期をフレーム周期と呼ぶ。そして、フーリエ変換を用いて切り出された信号を周波数領域信号に変換する。ｎをフレーム番号とし、切り出した入力信号ｘｍ（ｔ,ｎ）（ｔ＝０，１,・・・,Ｋ－１）とすると、ｘｍ（ｔ,ｎ）のフーリエ変換Ｘｍ（ｋ,ｎ）は次のように計算される。
Figure JPOXMLDOC01-appb-M000001
 　ここで、ｊは虚数単位（－１の平方根）、ｅｘｐは指数関数を表す。また、ｋは周波数ビン番号を表し、０以上Ｋ－１以下の整数である。以下、簡略化のため、ｋのことを周波数ビン番号ではなく、単に「周波数」と呼ぶ。
　（クロススペクトル計算部２０２）
　クロススペクトル計算部２０２は、変換部２０１から供給される変換信号に基づいて、クロススペクトルを計算し、分散計算部２０３と相関関数計算部２０４へ伝達する。クロススペクトル計算部２０２は、変換信号Ｘ２（ｋ,ｎ）の複素共役と変換信号Ｘ１（ｋ,ｎ）の積を計算する。変換信号のクロススペクトルをＳ１２（ｋ,ｎ）とすると、クロススペクトルは次のように計算される。
Figure JPOXMLDOC01-appb-M000002
 　ここで、ｃｏｎｊ（Ｘ２（ｋ,ｎ））は、Ｘ２（ｋ,ｎ）の複素共役を表す。また、式（２）で計算したクロススペクトルの代わりに、振幅成分で正規化したクロススペクトルを用いてもよい。このようなクロススペクトルは以下のように計算される。
Figure JPOXMLDOC01-appb-M000003
 　（分散計算部２０３）
　分散計算部２０３は、クロススペクトル計算部２０２から供給されたクロススペクトルを用いて、クロススペクトルの分散を計算し、相関関数計算部２０４へ伝達する。本実施形態では、過去に入力されたクロススペクトルから周波数ビンごとに分散を計算する例を説明する。周波数ビン単位ではなく、複数の周波数ビンを束ねたサブバンド単位や、全周波数ビンを用いて計算してもよい。また、それらを組み合わせて分散を計算してもよい。サブバンド単位で計算すると、分散を計算する帯域数が削減できるので、演算量が低減できる。
　分散計算部２０３で計算される分散とは、クロススペクトルではなく、クロススペクトルの位相の分散である。ここで、クロススペクトルの位相とは、クロススペクトルを複素数と捉えたとき、クロススペクトルの偏角のことである。また、分散とは、位相の時系列がある基準値に対してどの程度散らばっているかを示す指標、位相の分布の広がりのことである。したがって、本実施形態における分散とは、確率論における分散の定義に限定されるものではなく、標準偏差も分散の一例であるとする。分散の代表例は、確率論の定義に従ったものである。この場合、第ｎフレームの周波数ビンｋにおけるクロススペクトルＳ１２（ｋ,ｎ）の位相をａｒｇ（Ｓ１２（ｋ,ｎ））とすると、過去Ｌ－１フレーム分のクロススペクトルの位相ａｒｇ（Ｓ１２（ｋ,ｎ））,ａｒｇ（Ｓ１２（ｋ,ｎ－１））,ａｒｇ（Ｓ１２（ｋ,ｎ－２））,・・・,ａｒｇ（Ｓ１２（ｋ,ｎーＬ＋１））から、周波数ビンｋの分散Ｖ１２（ｋ,ｎ）は次のように計算される。
Figure JPOXMLDOC01-appb-M000004
 　また、標準偏差の定義に基づいて計算した場合、分散は次のように計算される。
Figure JPOXMLDOC01-appb-M000005
 　また、あらかじめ定めた基準値との２乗誤差の総和も分散の例であるとする。この場合、分散Ｖ１２（ｋ,ｎ）は以下のように計算される。
Figure JPOXMLDOC01-appb-M000006
 　ここで、ａｒｇ（Ｓ０）は基準となる位相である。他には、２乗誤差の代わりに絶対値誤差を用いてもよい。
Figure JPOXMLDOC01-appb-M000007
 　ここで、｜ｘ｜はｘの絶対値を表す。
　位相が角度情報であるという点に着目すると、方向統計学（Directional Statistics）で定義されている分散(円周分散、circular variance)や標準偏差(円周標準偏差、circular standard deviation)を用いてもよい。円周分散を用いた場合、分散Ｖ１２（ｋ,ｎ）は以下のように計算される。
Figure JPOXMLDOC01-appb-M000008
 　ここでＲ（ｋ,ｎ）は、クロススペクトルの平均であり、方向統計学では円周平均と呼ばれる。Ｒ（ｋ,ｎ）は、次のように計算される。
Figure JPOXMLDOC01-appb-M000009
 　ただし、式（９）のＳ１２（ｋ,ｎ）は、式（３）の「振幅成分で正規化したクロススペクトル」に限定される。クロススペクトルの平均は、式（８）の移動平均の代わりに、次のようにリーク積分を用いてもよい。
Figure JPOXMLDOC01-appb-M000010
 　ここで、αは実数で、０＜α＜１を満足する。
　式（８）の変形型である次の分散も有効である。
Figure JPOXMLDOC01-appb-M000011
 　また、円周標準偏差を用いた場合、分散Ｖ１２（ｋ,ｎ）は以下のように計算される。
Figure JPOXMLDOC01-appb-M000012
 　ここで、ｌｎ（ｘ）はｘの自然対数を表す。この他には、円周標準偏差に類する統計量として知られる角状変位（angular deviation）を用いてもよい。
Figure JPOXMLDOC01-appb-M000013
 　以上のとおり、方向統計学を用いた場合には、円周平均に基づいて各種分散を求めることが可能である。したがって、分散の計算方法としては、上記に示した数式に限定されることなく、円周平均の多次元関数、多項式、非線形関数、およびそれらの組み合わせを用いてもよい。
　以上のとおりに求めた分散に対して補正処理を追加することも、分散の精度を向上させる上では有効である。つまり、分散の代わりに、補正された分散を相関関数計算部２０４へ伝達してもよい。補正は、時間周波数平面上で周辺に存在する分散（補正対象の周辺のフレームや周波数に存在する分散）を用いて行う。代表的な例は、分散の平均値を計算する方法である。この場合、補正された分散ＶＶ１２（ｋ,ｎ）は次のように計算される。
Figure JPOXMLDOC01-appb-M000014
 　ここで、Ｐ１とＰ２とは、それぞれあらかじめ定めた定数である。この他の例としては、隣接する複数の周波数やフレームの分散にばらつきがある場合、そのフレーム・周波数に目標信号が含まれている可能性が低いので、分散を大きくする。また、隣接周波数・フレームに分散値が突出して低い値が見つかった場合、分散を低くする。
　（相関関数計算部２０４）
　相関関数計算部２０４は、クロススペクトル計算部２０２から供給されたクロススペクトルと、分散計算部２０３から供給された分散に基づき、相関関数を計算し、推定方向情報生成部２０５に伝達する。次に図３を参照して、相関関数計算部２０４の詳細を説明する。
　図３は、本実施形態に係る波源方向推定装置２００の備える相関関数計算部２０４の構成を示すブロック図である。相関関数計算部２０４は、重み付け部２４１および逆変換部２４２を備える。
　（重み付け部２４１）
　重み付け部２４１は、分散計算部２０３から供給された分散に基づいて、クロススペクトル計算部２０２から供給されたクロススペクトルに重み付けを行い、重み付けされたクロススペクトルを逆変換部２４２に伝達する。
　重み付け部２４１では、分散が大きい場合には、クロススペクトルに目標音源の情報が含まれている可能性が低いので、クロススペクトルの値を小さくする。代表的な方法は、あらかじめ用意した写像関数を用いて、分散から重みを計算し、その重みを乗じる方法である。このとき、重み付けされたクロススペクトルＷＳ１２（ｋ,ｎ）は次のように計算される。
Figure JPOXMLDOC01-appb-M000015
 　ここで、Ｇ（ｋ,ｎ）は分散に基づいて算出された重みである。Ｇ（ｋ,ｎ）は、次のように計算される。
Figure JPOXMLDOC01-appb-M000016
 　ここで、Ｖ１２（ｋ,ｎ）は分散計算部２０３から供給された分散である。Ｖ１２（ｋ,ｎ）の代わりに補正された分散であるＶＶ１２（ｋ,ｎ）を用いてもよい。また、ａとｂはそれぞれ実数であり、ａ＞０を満足する。また、次式のように、単純に閾値判定に基づいて値を切り替える方法でもよい。
Figure JPOXMLDOC01-appb-M000017
 　ここで、Ｖｔｈは閾値であり、正の実数である。この他にも、線形写像関数や高次の多項式関数、非線形関数など、他の形で表される分散の関数をゲインの計算に用いることも可能である。
　（逆変換部２４２）
　逆変換部２４２は、重み付け部２４１から供給されたクロススペクトルの逆変換を求める。本実施形態では、変換部２０１でフーリエ変換を用いたので、逆変換には逆フーリエ変換を用いる方法について説明する。重み付け部２４１から供給されたクロススペクトルをＷＳ１２（ｋ,ｎ）とすると、ＷＳ１２（ｋ,ｎ）の逆変換により得られる相関関数ｓ１２（τ,ｎ）は次のように計算される。
Figure JPOXMLDOC01-appb-M000018
 　（相対遅延時間計算部２０６）
　相対遅延時間計算部２０６は、入力されたマイク位置情報と音源探索対象方向とから、マイクペア間の相対遅延時間を求め、音源探索対象方向とセットで推定方向情報生成部２０５に伝達する。相対遅延時間とは、マイク間隔と音源方向とに基づいて一意に定まる音波の到達時間差のことである。音速をｃと仮定し、２つのマイクの間隔をｄ、音源方向、つまり音の到来方向をθとした場合、音源方向θに対する相対遅延時間τ(θ)は次の式で計算される。
Figure JPOXMLDOC01-appb-M000019
 　相対遅延時間は、全ての音源探索対象方向に対して計算される。例えば、方向の探索範囲が１０度刻みで０度から９０度まで、つまり０度、１０度、２０度、・・・９０度である場合、１０種類の相対遅延時間を計算する。そして、探索対象の方向と相対遅延時間とをペアで推定方向情報生成部２０５へ供給する。
　（推定方向情報生成部２０５）
　推定方向情報生成部２０５は、相関関数計算部２０４から供給される相関関数と、相対遅延時間計算部２０６から供給される相対遅延時間と、に基づき、方向と相関値との対応関係を推定方向情報として出力する。相関関数をｓ１２（τ,ｎ）、相対遅延時間τ(θ)とすると、推定方向情報Ｈ（θ,ｎ）は、次式で与えられる。
Figure JPOXMLDOC01-appb-M000020
 　各方向に対して相関値が定まるので、基本的に相関値が高ければ、その方向に音源が存在する可能性が高いと判断できる。
　このような推定方向情報は、様々な形で利用される。例えば、関数が複数のピークを有する場合には、各ピークを到来方向とする音源が複数存在すると考えられる。したがって、各音源の方向を同時に推定できるだけでなく、音源数の推定にも用いることが可能である。
　また、相関関数のピークと非ピークとの差分に基づき、音源の存在可能性を判定することも可能である。もし、ピークと非ピークとの差分が大きければ音源の存在可能性が高いと判定できる。同時に、推定方向の信頼性も高いと判断できる。もし、音源数が１つと事前に仮定できる場合には、相関値が最大となる方向を推定方向情報として出力してもよい。この場合、推定方向情報は、方向と相関値の対応関係ではなく、方向そのものとなる。
　図４Ａは、本実施形態に係る波源方向推定装置２００の備える相関関数テーブル４０１の一例を説明する図である。相関関数テーブル４０１は、入力信号４１１に関連付けて、周波数領域信号４１２、クロススペクトル４１３、分散４１４、相対遅延時間４１５および相関関数４１６を記憶する。波源方向推定装置２００は、例えば、相関関数テーブル４０１を参照して、入力信号４１１に対する相関関数４１６を算出し、推定方向情報を生成する。
　図４Ｂは、本実施形態に係る波源方向推定装置２００の備える変換関数テーブル４０２の一例を説明する図である。変換関数テーブル４０２は、入力信号４１１に関連付けて変換関数４２１および逆変換関数４２２を記憶する。波源方向推定装置２００は、例えば、変換関数テーブル４０２を参照して、入力信号４１１を変換して、周波数領域信号４１２やクロススペクトル４１３などを算出する。
　図４Ｃは、本実施形態に係る波源方向推定装置２００の備える分散テーブル４０３の一例を説明する図である。分散テーブル４０３は、クロススペクトル４２４に関連付けて分散手法４３１を記憶する。波源方向推定装置２００は、例えば、分散テーブル４０３を参照してクロススペクトル４２４の分散４１４を算出する。
　図５は、本実施形態に係る波源方向推定装置２００のハードウェア構成を示す図である。ＣＰＵ(Central Processing Unit)５１０は演算制御用のプロセッサであり、プログラムを実行することで図２および図３の波源方向推定装置２００の機能構成部を実現する。ＲＯＭ(Read Only Memory)５２０は、初期データおよびプログラムなどの固定データおよびその他のプログラムを記憶する。また、ネットワークインタフェース５３０は、ネットワークを介して他の装置などと通信する。なお、ＣＰＵ５１０は１つに限定されず、複数のＣＰＵであっても、あるいは画像処理用のＧＰＵ(Graphics Processing Unit)を含んでもよい。また、ネットワークインタフェース５３０は、ＣＰＵ５１０とは独立したＣＰＵを有して、ＲＡＭ(Random Access Memory)５４０の領域に送受信データを書き込みあるいは読み出しするのが望ましい。また、ＲＡＭ５４０とストレージ５５０との間でデータを転送するＤＭＡＣ(Direct Memory Access Controller)を設けるのが望ましい（図示なし）。さらに、入出力インタフェース５６０は、ＣＰＵ５１０とは独立したＣＰＵを有して、ＲＡＭ５４０の領域に入出力データを書き込みあるいは読み出しするのが望ましい。したがって、ＣＰＵ５１０は、ＲＡＭ５４０にデータが受信あるいは転送されたことを認識してデータを処理する。また、ＣＰＵ５１０は、処理結果をＲＡＭ５４０に準備し、後の送信あるいは転送はネットワークインタフェース５３０やＤＭＡＣ、あるいは入出力インタフェース５６０に任せる。
　ＲＡＭ５４０は、ＣＰＵ５１０が一時記憶のワークエリアとして使用するランダムアクセスメモリである。ＲＡＭ５４０には、本実施形態の実現に必要なデータを記憶する領域が確保されている。入力信号５４１は、入力端子２０１および入力端子２０２に入力された音などの信号である。周波数領域信号５４２は、入力信号を所定の変換方法で変換して得られた信号である。クロススペクトル５４３は、周波数領域信号５４２に基づいて計算される。分散５４４は、クロススペクトル５４３から計算されたクロススペクトル５４３の分散である。相対遅延時間５４５は、マイク間隔と音源方向とに基づいて決まる音波などの到達時間差である。相関関数５４６は、クロススペクトル５４３と分散５４４とから計算される。これらのデータは、例えば、相関関数テーブル４０１から展開されたデータである。なお、図示していないが、ＲＡＭ５４０は、変換関数テーブル４０２および分散テーブル４０３から展開されたデータを含んでいてもよい。
　入出力データ５４７は、入出力インタフェース５６０を介して入出力されるデータである。送受信データ５４８は、ネットワークインタフェース５３０を介して送受信されるデータである。また、ＲＡＭ５４０は、各種アプリケーションモジュールを実行するためのアプリケーション実行領域５４９を有する。
　ストレージ５５０には、データベースや各種のパラメータ、あるいは本実施形態の実現に必要な以下のデータまたはプログラムが記憶されている。ストレージ５５０は、相関関数テーブル４０１、変換関数テーブル４０２および分散テーブル４０３を格納する。相関関数テーブル４０１は、図４Ａに示した、入力信号４１１と相関関数４１６などとの関係を管理するテーブルである。変換関数テーブル４０２は、図４Ｂに示した、入力信号４１１と変換関数４２１などとの関係を管理するテーブルである。分散テーブル４０３は、図４Ｃに示した、クロススペクトル４２４と分散手法４３１などとの関係を管理するテーブルである。ストレージ５５０は、さらに、変換モジュール５５１、クロススペクトル計算モジュール５５２、分散計算モジュール５５３、相関関数計算モジュール５５４、相対遅延時間計算モジュール５５５および推定方向情報生成モジュール５５６を格納する。
　変換モジュール５５１は、入力信号を変換して周波数領域信号を取得するモジュールである。クロススペクトル計算モジュール５５２は、周波数領域信号に基づいてクロススペクトルを計算するモジュールである。分散計算モジュール５５３は、クロススペクトルの分散を計算するモジュールである。相関関数計算モジュール５５４は、クロススペクトルと分散とから相関関数を計算するモジュールである。相対遅延時間計算モジュール５５５は、マイクペア間の相対遅延時間を計算するモジュールである。推定方向情報生成モジュール５５６は、相関関数と相対遅延時間とから方向と相関値との対応関係を推定方向情報として生成するモジュールである。これらのモジュール５５１～５５６は、ＣＰＵ５１０によりＲＡＭ５４０のアプリケーション実行領域５４９に読み出され、実行される。制御プログラム５５７は、波源方向推定装置２００の全体を制御するためのプログラムである。
　入出力インタフェース５６０は、入出力機器との入出力データをインタフェースする。入出力インタフェース５６０には、表示部５６１、操作部５６２、が接続される。また、入出力インタフェース５６０には、さらに、記憶媒体５６４が接続されてもよい。さらに、音声出力部であるスピーカ５６３や、音声入力部であるマイク、あるいは、ＧＰＳ位置判定部が接続されてもよい。なお、図５に示したＲＡＭ５４０やストレージ５５０には、波源方向推定装置２００が有する汎用の機能や他の実現可能な機能に関するプログラムやデータは図示されていない。
　図６は、本実施形態に係る波源方向推定装置２００の処理手順を説明するフローチャートである。このフローチャートは、ＣＰＵ５１０がＲＡＭ５４０を使用して実行し、図２および図３の波源方向推定装置２００の機能構成部を実現する。
　ステップＳ６０１において、波源方向推定装置２００は、入力信号を取得する。ステップＳ６０３において、波源方向推定装置２００は、入力信号を変換して周波数領域信号を取得する。ステップＳ６０５において、波源方向推定装置２００は、周波数領域信号からクロススペクトルを計算する。ステップＳ６０７において、波源方向推定装置２００は、クロススペクトルの位相の分散を計算する。ステップＳ６０９において、波源方向推定装置２００は、計算した分散に基づいて、クロススペクトルを重み付けする。ステップＳ６１１において、波源方向推定装置２００は、重み付けしたクロススペクトルを逆変換して、相関関数を計算する。ステップＳ６１３において、波源方向推定装置２００は、相対遅延時間を計算する。ステップＳ６１５において、相関関数と相対遅延時間とから推定方向情報を生成する。
　＜本実施形態の効果＞
　本実施形態によれば、入力信号に含まれる目標音の到来方向、すなわち目標物体が存在する方向や波源の方向を高精度に推定できる。環境雑音レベルが高い環境において、目標物が発する音を手掛かりに、目標物が存在する方向を推定する場合に有効である。環境雑音の例としては、繁華街や街頭、街道沿い、人や自動車が多く集まる場所が挙げられる。また、目標物の例としては、人間や動物、自動車、航空機、船舶、水上バイク、ドローン（小型無人機）が挙げられる。
　例えば、屋外のテーマパークや展示会場などに接近する不審な自動車・船舶・ドローンなどを検知し、それらの方向を推定することで、不審者や不審物を効率的に取り締ることが可能である。また、音源方向推定を複数箇所で実施することで、目標音源の位置を特定できる。これにより、環境雑音レベルが高い環境でも、悲鳴・怒号の発生位置や、不審な船舶・ドローンの出現位置などを正確に特定することが可能となる。
　クロススペクトルの位相から、２つのマイクへの音波の到達時間差、すなわち音源方向が得られる。したがって、位相の時系列（本実施形態では過去Ｌ－１フレームの位相）から得られた分散の大きさに基づいて、目標音と環境音との区別が可能になる。目標音源の方向の時間変化は小さいため、位相の分散は小さくなる。もし、目標音の発生時間が数秒程度であっても、例えば１０ミリ秒ごとに位相を求める場合には、その数秒から得られる位相の分散は小さくなる。一方で、環境雑音は方向の特定が困難なため、方向の時間変化は大きくなり、位相の分散も大きくなる。分散は周波数ごとに計算するので、目標音源が存在する周波数帯域を特定することが可能になる。
　［第３実施形態］
　次に本発明の第３実施形態に係る波源方向推定装置について、図７を用いて説明する。図７は、本実施形態に係る波源方向推定装置７００の構成を説明するためのブロック図である。本実施形態に係る波源方向推定装置７００は、上記第２実施形態と比べると、平均計算部７０１を有する点で異なる。その他の構成および動作は、第２実施形態と同様であるため、同じ構成および動作については同じ符号を付してその詳しい説明を省略する。なお、変換部２０１、クロススペクトル計算部２０２、分散計算部２０３、相関関数計算部２０４および平均計算部７０１は、相関関数生成部７１０を構成する。
　（平均計算部７０１）
　平均計算部７０１は、クロススペクトル計算部２０２から供給されたクロススペクトルの平均を計算し、相関関数計算部２０４へ伝達する。本実施形態では、過去に入力されたクロススペクトルから周波数ビンごとに平均を計算する例を説明する。周波数ビン単位ではなく、複数の周波数ビンを束ねたサブバンド単位で計算してもよい。第ｎフレームの周波数ビンｋにおけるクロススペクトルをＳ１２（ｋ,ｎ）とすると、過去Ｌ－１フレームから求めた平均クロススペクトルＳＳ１２（ｋ,ｎ）は次のように計算される。
Figure JPOXMLDOC01-appb-M000021
 　また、次のようにリーク積分を用いてもよい。
Figure JPOXMLDOC01-appb-M000022
 　ここで、αは実数で、０＜α＜１を満足する。
　＜本実施形態の効果＞
　本実施形態によれば、クロススペクトルの代わりに平均クロススペクトルを相関関数計算部２０４に伝達する。このため、第２実施形態と比べると、位置の時間変動が小さな目標音源に対しては、相関関数のピーク位置、および推定方向情報の時間変動が小さくなるので、方向推定の精度が向上する。特に、自動車の走行音、航空機やドローンの飛行音、船舶の航行音(スクリュー音)など、目標物が発する音が継続的に生じる場合、クロススペクトルの平均化効果は高くなるので、相関関数のピーク位置、および推定方向情報が明確になる。
　［第４実施形態］
　次に本発明の第４実施形態に係る波源方向推定装置について、図８を用いて説明する。図８は、本実施形態に係る波源方向推定装置８００の構成を説明するためのブロック図である。本実施形態に係る波源方向推定装置８００は、上記第３実施形態と比べると、分散計算部２０３に代えて、分散計算部８０１を有する点で異なる。その他の構成および動作は、第３実施形態と同様であるため、同じ構成および動作については同じ符号を付してその詳しい説明を省略する。なお、変換部２０１、クロススペクトル計算部２０２、分散計算部８０１、平均計算部７０１および相関関数計算部２０４は、相関関数生成部８１０を構成する。
　（分散計算部８０１）
　分散計算部８０１は、平均計算部３０１から供給された平均クロススペクトルを用いて分散を計算し、相関関数計算部２０４へ伝達する。第３実施形態と比較すると、分散の計算にクロススペクトルではなく、平均クロススペクトルを用いている点が異なる。特に、方向統計学で知られている分散や標準偏差を計算する場合には、その計算過程において円周平均の計算を省略できる。平均クロススペクトルをＳＳ１２（ｋ,ｎ）とすると、「クロススペクトルの位相の分散」の計算において円周分散を用いた場合、分散Ｖ１２（ｋ,ｎ）は以下のように計算される。
Figure JPOXMLDOC01-appb-M000023
 式（２３）の変形型である次の分散も有効である。
Figure JPOXMLDOC01-appb-M000024
 　また、円周標準偏差を用いた場合、分散Ｖ１２（ｋ,ｎ）は以下のように計算される。
Figure JPOXMLDOC01-appb-M000025
 　＜本実施形態の効果＞
　本実施形態によれば、分散計算部は、平均クロススペクトル計算部で計算された平均クロススペクトルを利用して、分散を計算する。このため、分散計算部において平均クロススペクトルの計算が不要になる。したがって、第３実施形態よりも少ない計算量で分散を計算することができる。
　［第５実施形態］
　次に本発明の第５実施形態に係る波源方向推定装置について、図９乃至図１２Ｄを用いて説明する。本実施形態に係る波源方向推定装置は、上記第２実施形態の波源方向推定装置と比べると、相関関数計算部２０４に代えて、相関関数計算部９０４を有する点で異なる。図９は、本実施形態に係る相関関数計算部９０４の構成を説明するためのブロック図である。本実施形態に係る相関関数計算部９０４は、上記第２実施形態の相関関数計算部２０４と比べると、重み付け部２４１と逆変換部２４２に代えて、周波数別クロススペクトル生成部９４１１，９４１２,・・・，９４１Ｋと統合相関関数計算部９４２とを有する点で異なる。その他の構成および動作は、第２実施形態と同様であるため、同じ構成および動作については同じ符号を付してその詳しい説明を省略する。
　図９は、相関関数計算部９０４のブロック図である。相関関数計算部９０４は、周波数別クロススペクトル生成部９４１１，９４１２,・・・，９４１Ｋと統合相関関数計算部９４２とを備える。周波数別クロススペクトル生成部９４１１，９４１２,・・・，９４１Ｋは、クロススペクトル計算部２０２から供給されるクロススペクトルＳ１２（ｋ,ｎ）と分散計算部２０３から供給される分散ｐ（ｋ,ｎ）とを用いて、Ｓ１２（ｋ,ｎ）の各周波数ｋに対応するクロススペクトルを計算し、周波数別クロススペクトルとして統合相関関数計算部９４２に伝達する。周波数別クロススペクトルの生成は、周波数成分ごとに相関関数を計算するために行われる。つまり、ある周波数ｋに対応する相関関数（周波数別相関関数と呼ぶ）を後段で求めるために、周波数別クロススペクトルを計算する。
　相関関数計算部９０４は、統合相関関数計算部９４２および周波数別クロススペクトル生成部９４１１，９４１２,・・・，９４１Ｋから供給される周波数別クロススペクトルに基づき統合相関関数を計算し、推定方向情報生成部２０５へ伝達する。
　次に、図１０を参照して、ある周波数ｋの周波数別クロススペクトルを計算する周波数別クロススペクトル生成部９４１ｋについて詳細に説明する。その後、図１１を参照して、統合相関関数計算部９４２の詳細について説明する。
　＜周波数別クロススペクトル生成部９４１ｋ＞
　図１０は、周波数別クロススペクトル生成部９４１ｋのブロック図である。周波数別クロススペクトル生成部９４１ｋは、周波数別クロススペクトル生成部１０１１、カーネル関数スペクトル生成部１０１２、乗算部１０１３を有する。周波数別クロススペクトル生成部９４１ｋは、クロススペクトル計算部２０２から供給されるクロススペクトルＳ１２（ｋ,ｎ）と分散計算部２０３から供給される分散を用いて、Ｓ１２（ｋ,ｎ）の周波数ｋに対応するクロススペクトルを計算し、周波数別クロススペクトルを統合相関関数計算部９４２へ伝達する。
　（周波数別クロススペクトル生成部９４１）
　周波数別クロススペクトル生成部９４１１，９４１２,・・・，９４１Ｋは、クロススペクトル計算部２０２から供給されるクロススペクトルＳ１２（ｋ,ｎ）を用いて、Ｓ１２（ｋ,ｎ）の各周波数ｋに対応するクロススペクトルを計算し、周波数別クロススペクトルとして乗算部１０１３に伝達する。周波数別クロススペクトルは、周波数成分ごとに相関関数を計算するために行われる。つまり、ある周波数ｋに対応する相関関数（周波数別相関関数と呼ぶ）を後段で求めるために、周波数別クロススペクトルを計算する。
　次に、ある周波数ｋの周波数別クロススペクトルを計算する周波数別クロススペクトル生成部９４１ｋについて詳細に説明する。周波数別クロススペクトル生成部９４１ｋでは、周波数ｋのクロススペクトルＳ１２（ｋ,ｎ）に基づいて周波数別クロススペクトルを計算する際、位相成分と振幅成分をあらかじめ別々に求めたのちに統合する。周波数ｋの周波数別クロススペクトルＵｋ（ｗ,ｎ）をとし、その振幅成分を｜Ｕｋ（ｗ,ｎ）｜、位相成分をａｒｇ（Ｕｋ（ｗ,ｎ））とすると、次のような関係が成立する。
Figure JPOXMLDOC01-appb-M000026
 　ここで、ｗは周波数を表し、０以上Ｗ－１以下の整数である。以下、周波数ｋのクロススペクトルＳ１２（ｋ,ｎ）から、周波数別クロススペクトルの振幅成分|Ｕｋ（ｗ,ｎ）|と位相成分ａｒｇ（Ｕｋ（ｗ,ｎ））とを求める方法について説明する。
　振幅成分|Ｕｋ（ｗ,ｎ）|は、ｋを整数倍した周波数には、１．０を用いる。一方、周波数ｋの非定数倍した周波数の位相成分はゼロにする。これを数式で表現すると、振幅成分|Ｕｋ（ｗ,ｎ）|は以下のように与えられる。
Figure JPOXMLDOC01-appb-M000027
 　ここで、ｐは１以上Ｐ以下の整数である。波源方向推定を行うときに重要な情報は位相成分であるため、振幅成分はこのように適当な定数を用いる。その他には、１．０の代わりに|Ｓ１２（ｋ,ｎ）|を用いてもよい。つまり、振幅成分|Ｕｋ（ｗ,ｎ）|を次式のように求めてもよい。
Figure JPOXMLDOC01-appb-M000028
 　位相成分ａｒｇ（Ｕｋ（ｗ,ｎ））は、ｋを整数倍した周波数には、周波数ｋのクロススペクトルＳ１２（ｋ,ｎ）を定数倍したものを用いる。例えば、周波数ｋ、２ｋ、３ｋ、４ｋの位相成分には、周波数ｋの位相成分ａｒｇ（Ｓ１２（ｋ,ｎ））をそれぞれ同一の倍率で整数倍したもの、つまりａｒｇ（Ｓ１２（ｋ,ｎ））、２ａｒｇ（Ｓ１２（ｋ,ｎ））、３ａｒｇ（Ｓ１２（ｋ,ｎ））、４ａｒｇ（Ｓ１２（ｋ,ｎ））を用いる。一方、周波数ｋの非定数倍した周波数の位相成分はゼロにする。したがって、周波数ｋに対応する周波数別クロススペクトルの位相成分ａｒｇ（Ｕｋ（ｗ,ｎ））は、以下のように計算される。
Figure JPOXMLDOC01-appb-M000029
 　ここで、ｐは１以上Ｐ以下の整数である。また、Ｐは１よりも大きい整数である。
　以上の方法で求めた振幅成分と位相成分とを、先に記載した式（２６）を用いて統合し、周波数ｋの周波数別クロススペクトルＵｋ（ｗ,ｎ）を得る。
　これまで説明してきた方法では、振幅成分と位相成分とを別々に求めてから周波数別スペクトルを得ていた。しかし、次に示す数式に示すようにクロススペクトルのべき乗を用いれば、振幅成分と位相成分とを求めることなく周波数別スペクトルＵｋ（ｗ,ｎ）を求めることが可能である。
Figure JPOXMLDOC01-appb-M000030
 　(カーネル関数スペクトル生成部１０１２)
　カーネル関数スペクトル生成部１０１２は、分散計算部２０３から供給された分散に基づいてカーネル関数スペクトルを求め、乗算部１０１３へカーネル関数スペクトルを出力する。カーネル関数スペクトルとは、カーネル関数をフーリエ変換し、その絶対値を取ったものである。絶対値を取る代わりに二乗してもよい。また、絶対値の二乗としてもよい。本実施形態では、カーネル関数スペクトルをＧ（ｗ）、カーネル関数をｇ（τ）と表記して説明する。カーネル関数としては、ガウス関数を用いる。この時、ガウス関数は、次のような数式で与えられる。
Figure JPOXMLDOC01-appb-M000031
 　ここで、ｇ１,ｇ２,ｇ３は正の実数である。ｇ１はガウス関数の大きさ、ｇ２はガウス関数のピークの位置、ｇ３はガウス関数の広がりを制御する。特に、ガウス関数の広がりを調整するｇ３は、周波数別相関関数のピークの鋭さに大きな影響を与えるので重要である。式（３１）から分かるように、ｇ３が大きくなればガウス関数の広がりは大きくなる。
　他には、以下に示すロジスティック関数を用いてもよい。
Figure JPOXMLDOC01-appb-M000032
 　ここで、ｇ１とｇ２とは正の実数である。ロジスティック関数は、ガウス関数と同様の形状をしているが、ガウス関数よりも裾が長いという性質を持つ。特に、ロジスティック関数の広がりを調整するｇ５は、ガウス関数におけるｇ３の場合と同様、周波数別相関関数のピークの鋭さに大きな影響を与える重要なパラメータである。この他にも、コサイン関数や一様関数を用いてもよい。
　カーネル関数のパラメータのうち、カーネル関数の広がりに影響を与えるｇ３やｇ５は、分散計算部２０３から供給された分散に基づいて決定する。本実施形態では、これらのパラメータのことを広がり制御パラメータと呼び、ｑ（ｋ,ｎ）で表現する。したがって、カーネル関数がガウス関数の場合には、ｇ３＝ｑ（ｋ,ｎ）である。分散が小さければ、周波数別相関関数のピークが鋭く、裾が狭くなるようにパラメータを変化させる。したがって、広がり制御パラメータを小さくする。
　広がり制御パラメータの基本的な算出方法は、あらかじめ設定した写像関数を用いて、分散の値を広がり制御パラメータに変換する方法である。例えば、分散がある閾値を上回れば広がり制御パラメータを大きな値（例えば１０）に、閾値を下回れば小さな値（例えば０．０１）に設定する。この場合、分散をＶ１２（ｋ,ｎ）、閾値をｐｔｈとすると、第ｎフレームの周波数ビンｋにおける広がり制御パラメータｑ（ｋ,ｎ）は次のように計算される。
Figure JPOXMLDOC01-appb-M000033
 　ここで、ｑ１とｑ２とは正の実数であり、ｑ１＞ｑ２を満足する。また、次のように一次関数を用いて計算してもよい。
Figure JPOXMLDOC01-appb-M000034
 　ここで、ｑ３とｑ４とはそれぞれ実数であり、ｑ３＞０を満足する。この他にも、線形写像関数や高次の多項式関数、非線形関数など、他の形で表される分散の関数を分散の計算に用いることも可能である。また、分散をそのまま広がり制御パラメータにしてもよい（ｑ（ｋ,ｎ）＝Ｖ１２（ｋ,ｎ）とすると同じ)。
　広がり制御パラメータを求める関数は、分散だけでなく、周波数ｋの関数としてもよい。例えば、周波数が高くなるにつれて小さくなる関数とする。このような代表例としてｋの逆数を用いる例を説明する。この場合、式（３３）の代わりに、次の関数を用いて広がり制御パラメータｑ（ｋ,ｎ）を計算する。
Figure JPOXMLDOC01-appb-M000035
 　また、式（３４）の代わりとしては、次の関数を用いて広がり制御パラメータｑ（ｋ,ｎ）を計算する。
Figure JPOXMLDOC01-appb-M000036
 　（乗算部１０１３）
　乗算部１０１３は、周波数別クロススペクトル生成部１０１１から供給される周波数別クロススペクトルと、カーネル関数スペクトル生成部１０１２から供給されるカーネル関数スペクトルと、の積を計算し、計算によって得られた新たな周波数クロススペクトルを統合相関関数計算部９４２に伝達する。周波数別クロススペクトル生成部１０１１から供給される周波数別スペクトルをＵｋ（ｗ,ｎ）、カーネル関数スペクトル生成部１０１２から供給されるカーネル関数スペクトルをＧ（ｗ）とすると、計算により得られる新たな周波数別クロススペクトルＵＭｋ（ｗ,ｎ）は次のように計算される。
Figure JPOXMLDOC01-appb-M000037
 　＜統合相関関数計算部＞
　図１１は、統合相関関数計算部９４２のブロック図である。統合相関関数計算部９４２は、逆変換部１１２１１,１１２１２,・・・,１１２１Ｋと、統合部１１２２を備える。
　逆変換部１１２１１,１１２１２,・・・,１１２１Ｋは、周波数別クロススペクトル生成部９４１１,９４１２,・・・,９４１Ｋから供給される周波数別クロススペクトルの逆変換を行い、周波数別相関関数として統合部１１２２にそれぞれ伝達する。本実施形態では、変換部２０１でフーリエ変換を用いたので、逆変換には逆フーリエ変換を用いる方法について説明する。周波数別クロススペクトル生成部９４１ｋから供給された周波数別クロススペクトルをＵＭｋ（ｗ,ｎ）とすると、ＵＭｋ（ｗ,ｎ）の逆変換により得られる周波数別相関関数ｕｋ（τ,ｎ）は次のように計算される。
Figure JPOXMLDOC01-appb-M000038
 　統合部１１２２は、逆変換部１１２１１,１１２１２,・・・,１１２１Ｋから供給される周波数別相関関数を統合し、統合相関関数として推定方向情報生成部２０５へ伝達する。個別に求めた複数の周波数別相関関数を混合したり、重ね合わせたりすることにより、一つの相関関数を求める。統合方法に単純な和を用いた場合、統合部１１２２は、周波数別相関関数の総和を計算する。統合相関関数をｕ（τ,ｎ）とすると、ｕ（τ,ｎ）は次のように計算される。
Figure JPOXMLDOC01-appb-M000039
 　また、総和ではなく総乗を用いてもよい。この場合、ｕ（τ,ｎ）は次のように計算される。
Figure JPOXMLDOC01-appb-M000040
 　目標音の存在する周波数や、目的音のパワーが大きい周波数があらかじめ判明している場合には、その周波数に対応する周波数別相関関数だけ用いて、統合された相関関数を求めてもよい。また、重み付けという形で、統合における周波数別相関関数の影響度を制御してもよい。例えば目標音の存在する周波数の集合をΩとすると、周波数の選択によりｕ（τ,ｎ）を求める場合は、次のように計算される。
Figure JPOXMLDOC01-appb-M000041
 　また、重み付けを用いる場合には、ｕ（τ,ｎ）は次のように計算される。
Figure JPOXMLDOC01-appb-M000042
 　ここで、ａとｂとは実数で、ａ＞ｂ＞０を満足する。このように、目的音が存在する周波数の周波数別相関関数を重点的に用いて統合すると、雑音などの非目的音の影響が小さい相関関数を生成できるため、方向推定精度が向上する。
　＜周波数別クロススペクトルの説明＞
　以上の方法により周波数別クロススペクトルを計算すると、周波数別クロススペクトルを逆変換して得られる周波数別相関関数のピークが鋭くなり、相関関数のピーク位置が明確になる。相関関数のピーク位置に基づいて波源方向推定を行う本発明の場合には、ピークが鋭くなれば音源方向推定の精度が向上する。さらに、Ｐの値が大きくなるにつれて、すなわちｋを整数倍した周波数の成分が増えるにつれて、相関関数のピークは鋭くなる。図１２Ａに、この様子を示す。ここで、図中のＱは３よりも大きい整数である。Ｐ＝１の場合、つまり位相成分が一つしかない場合、これを逆変換して得られる相関関数は、このようにピーク位置が不明確な相関関数になる。Ｐが大きくなると、図１２Ａに示す通り、相関関数のピークは鋭くなる。
　したがって、周波数別クロススペクトル生成部１０１１で得られる周波数別クロススペクトルは、ある周波数ｋのクロススペクトルを基に「周波数ｋを整数倍した周波数ｐｋの位相成分に、周波数ｋの位相成分ａｒｇ（Ｓ１２（ｋ,ｎ））をｐ倍したものを割り当てたもの」として定義される。ここでｐは１以上の整数である。つまり、周波数別クロススペクトルは、その位相成分ａｒｇ（Ｕｋ（ｗ,ｎ））が少なくとも以下の式を満足するものとして定義される。
Figure JPOXMLDOC01-appb-M000043
 　加えてｐは、ｐ＝１と２、ｐ＝１と３、ｐ＝２と３など、２つ以上に限定される。ｐが１のみの場合には、周波数ｋの成分だけを抜き出して周波数別クロススペクトルが生成されることになるが、方向推定精度は従来技術と同等となり、方向推定の高精度を達成できない。なお、図１２Ａで説明した通り、ｐ＝１，２，３，４,・・・とｐの数が増えれば、つまり周波数ｐｋの位相成分への割り当てが増えれば、周波数別相関関数のピークは鋭くなるため、方向推定の精度は向上する。
　＜周波数別クロススペクトルにカーネル関数スペクトルを乗じる効果の説明＞
　周波数別クロススペクトルにカーネル関数スペクトルを乗じることで、周波数別相関関数の形状を変更できる。図１２Ｂは、カーネル関数スペクトルが乗じられた周波数別クロススペクトルと、周波数別相関関数の関係を示している。比較のため、カーネル関数スペクトルを乗じる前の周波数別クロススペクトルも掲載する。図１２Ｂの左側図に示すように、カーネル関数スペクトルを乗じていなければ、高い周波数まで成分が存在するため、周波数別相関関数のピークは鋭くなる。一方で、図１２Ｂの中央図、右側図のようにカーネル関数スペクトルを乗じると、高い周波数の成分が減衰するので、周波数別相関関数のピークの鋭さは小さくなる。つまり、カーネル関数スペクトルのピークが鋭くなる（カーネル関数スペクトルの裾が狭くなる）につれて、周波数別相関関数のピークの鋭さは小さくなる。また、図１２Ｂの右側図のように、周波数別相関関数の裾が大きく広がると、隣接する山の裾が重なり、谷が浅い周波数別相関関数が得られる。
　ここで、カーネル関数とカーネル関数スペクトルとの形状の関係について補足する。フーリエ変換の性質により、形状の関係は逆になる。カーネル関数のピークが鋭く、裾が狭くなるにつれて、カーネル関数スペクトルのピークは平坦に近づき、裾が広くなる。ガウス関数の広がりを調整するｇ３との関係も含めて説明すると、ｇ３が大きくなれば、ガウス関数の広がりは大きくなるが、そのスペクトルの広がりは小さくなる。
　カーネル関数による周波数別相関関数の大きさ制御の効果について、図１２Ｃにおいて説明する。図１２Ｃは、カーネル関数の有無と統合相関関数との関係を示した図である。図１２Ｃに示すようなカーネル関数無しの場合（１２１０）、周波数別相関関数ｕ１（τ,ｎ）～ｕ３（τ,ｎ）のピーク位置は近いが、ｕ１（τ,ｎ）～ｕ３（τ,ｎ）の幅が細いため、統合時に大きなピークを形成できない。このため、ピークの位置が明確にならない。一方、図１２Ｃに示すようなカーネル関数有りの場合（１２２０）、周波数別相関関数の幅は太くなっているため、統合によりｕ１（τ,ｎ）～ｕ３（τ,ｎ）は大きなピークを形成できる。このため、カーネル関数無しの場合（１２１０）に比べて、ピークの位置が明確になる。
　さらに、カーネル関数による周波数別相関関数の大きさ制御の別の効果について、図１２Ｄを参照して説明する。図１２Ｄは、カーネル関数スペクトルの幅の違いと、統合相関関数と、の関係を示した図である。図１２Ｂの右側図に示すように、幅の広いカーネル関数スペクトルを用いると、相関関数の周期性により、谷が浅い周波数別相関関数が形成される。したがって、図１２Ｄの１２３０に示す通り、谷が浅い周波数別相関関数を統合すると、谷が浅い、つまりピークの目立たない統合相関関数が生成される。一方、図１２Ｂの中央図に示すように、幅の狭いカーネル関数スペクトルを用いた場合には、図１２Ｂの右側図よりも谷が深い周波数別相関関数が形成される。したがって、図１２Ｄの１２４０に示す通り、ピークが明確な統合相関関数が生成される。
　本実施形態では、カーネル関数のフーリエ変換により得られたカーネル関数スペクトルと、周波数別クロススペクトルと、の積を計算しているが、フーリエ変換の性質により時間領域でも実現できる。周波数別クロススペクトル計算部ではなく、統合相関関数計算部９４２の中にある逆変換部１１２１ｋの後段にカーネル関数を畳み込む「畳み込み演算部」を設け、逆変換部１１２１ｋから供給される周波数別相関関数にカーネル関数を畳み込んでもよい。ただし、畳み込み演算は計算量が多いため、本実施形態のように周波数領域で積を計算するほうが効率的である。
　＜本実施形態の効果＞
　本実施形態によれば、周波数別クロススペクトルを逆変換して得られる周波数別相関関数のピークが鋭くなり、相関関数のピーク位置が明確になる。相関関数のピーク位置に基づいて波源方向推定を行う本実施形態の場合には、音源方向推定の精度が向上する。また、周波数別クロススペクトルにカーネル関数スペクトルを乗じることで、周波数別相関関数の形状を変更できる。その際、カーネル関数の広さは、分散の大きさに応じて変更するので、目標音源が存在する周波数の周波数別相関関数のピークを強調できる。したがって、周波数ごとの相関関数を求めず、また周波数別相関関数のピーク強調を行わない第２実施形態と比較すると、相関関数のピークが強調できるので、方向推定精度が向上する。
　［第６実施形態］
　次に本発明の第６実施形態に係る波源方向推定装置について、図１３を用いて説明する。本実施形態に係る波源方向推定装置は、上記第５実施形態の波源方向推定装置と比べると、周波数別クロススペクトル生成部９４１ｋに代えて、周波数別クロススペクトル生成部１３４１ｋを有する点で異なる。図１３は、本実施形態に係る周波数別クロススペクトル生成部１３４１ｋの構成を説明するためのブロック図である。本実施形態に係る周波数別クロススペクトル生成部１３４１ｋは、上記第５実施形態の周波数別クロススペクトル生成部９４１ｋと比べると、カーネル関数スペクトル生成部１０１２に代えて、カーネル関数スペクトル記憶部１３１４とカーネル関数スペクトル選択部１３１２とを有する点で異なる。その他の構成および動作は、第５実施形態と同様であるため、同じ構成および動作については同じ符号を付してその詳しい説明を省略する。
　図１３は、周波数別クロススペクトル生成部１３４１ｋのブロック図である。周波数別クロススペクトル生成部１３４１ｋは、周波数別クロススペクトル生成部１０１１、カーネル関数スペクトル選択部１３１２、カーネル関数スペクトル記憶部１３１４、乗算部１０１３を有する。周波数別クロススペクトル生成部１０１１ｋは、クロススペクトル計算部２０２から供給されるクロススペクトルＳ１２（ｋ,ｎ）と分散計算部２０３から供給される分散を用いて、Ｓ１２（ｋ,ｎ）の周波数ｋに対応するクロススペクトルを計算し、周波数別クロススペクトルを統合相関関数計算部９４２へ伝達する。
　カーネル関数スペクトル記憶部１３１４は、カーネル関数スペクトルを記憶しており、カーネル関数スペクトル選択部１３１２へカーネル関数スペクトルを伝達する。カーネル関数スペクトル記憶部１３１４は、異なる分散に対応する複数のカーネル関数スペクトルを記憶する。例えば、０．１から０．１刻みで１０種類の分散（０．１,０．２,．．．,１．０）に対応するカーネル関数スペクトルを記憶しておく。すなわち、０．１から１．０までの全て分散値に対してカーネル関数スペクトルをあらかじめ計算しておき、それらをカーネル関数スペクトル記憶部１３１４に記憶しておく。この場合、１０種類のカーネル関数スペクトルが記憶される。本実施形態では、分散が０．１、０．２、・・・１．０のカーネル関数スペクトルをＧ０．１（ｗ）、Ｇ０．２（ｗ）、・・・Ｇ１．０（ｗ）と表記する。
　第５実施形態で説明した通り、ピークが明確な統合相関関数を生成する上では、様々な幅を有するカーネル関数を用意することが望ましい。したがって、多くの分散値に対応するカーネル関数スペクトルを用意する、例えば、０．００１から０．００１刻みで１０，０００種類の分散（０．００１,０．００２,０．００３,・・・）に対応するカーネル関数スペクトルを記憶することが望ましい。しかし、記憶容量が増大する問題が生じる。したがって、記憶容量の増大を抑えつつ、ピークが明確な統合相関関数を生成するためには、０．０５から０．０５刻み、もしくは０．１刻みで約５～２０種類の分散に対応するカーネル関数スペクトルを用意する必要がある。なお、刻み幅は必ずしも均等である必要はない。統合相関関数のピークを強調する上では、０に近い分散値では細かく、０から遠ざかるにつれて刻み幅を広くしたほうがよい。また、分散値がある閾値を超えた場合には、十分に大きな値を有する１種類の分散値に対応するカーネル関数スペクトルを用意すればよい。
　カーネル関数スペクトル選択部１３１２は、分散計算部２０３から供給された分散に基づいて、カーネル関数スペクトル記憶部１３１４に記憶されている複数のカーネル関数スペクトルの中から、一つのカーネル関数スペクトルを選択し、乗算部１０１３へ伝達する。カーネル関数スペクトル選択部１３１２は、分散計算部２０３から供給された分散値に最も近い分散値で生成されたカーネル関数スペクトルを選択する。例えば、供給された分散値が０から０．１５未満の値であれば、カーネル関数スペクトルＧ０．１（ｗ）を選択する。同様に、供給された分散値が０．１５以上から０．２５未満の値であれば、カーネル関数スペクトルＧ０．２（ｗ）を選択する。
　＜本実施形態の効果＞
　本実施形態によれば、あらかじめカーネル関数を計算しておくので、方向推定処理中にはカーネル関数の計算が不要になる。したがって、第５実施形態よりも少ない計算量で方向推定を行うことができる。
　［第７実施形態］
　次に本発明の第７実施形態に係る波源方向推定装置について、図１４を用いて説明する。本実施形態に係る波源方向推定装置は、上記第５実施形態の波源方向推定装置と比べると、統合相関関数計算部９４２に代えて、統合相関関数計算部１４４２を有する点で異なる。図１４は、本実施形態に係る統合相関関数計算部１４４２の構成を説明するためのブロック図である。本実施形態に係る統合相関関数計算部１４４２は、上記第５実施形態の統合相関関数計算部９４２と比べると、逆変換部１１２１１，１１２１２,・・・，１１２１Ｋと統合部１１２２とに代えて、統合部１４２１と逆変換部１４２２とを有する点で異なる。その他の構成および動作は、第５実施形態と同様であるため、同じ構成および動作については同じ符号を付してその詳しい説明を省略する。
　図１４は、統合相関関数計算部１４４２のブロック図である。統合相関関数計算部１４４２は、統合部１４２１と逆変換部１４２２とを備える。
　統合部１４２１は、周波数別クロススペクトル生成部９４１１，９４１２,・・・，９４１Ｋから供給される周波数別クロススペクトルを統合し、統合クロススペクトルとして逆変換部１４２２へ伝達する。個別に求めた複数の周波数別クロススペクトルを混合したり、重ね合わせたりすることにより、一つの統合クロススペクトルを求める。統合には、第５実施形態の統合部１１２２と同様に総和や総乗が用いられる。統合に総和を用いた場合、周波数別クロススペクトル生成部９４１ｋから供給された周波数別クロススペクトルをＵＭｋ（０,ｎ）,ＵＭｋ（１,ｎ）,・・・,ＵＭｋ（Ｗ・１,ｎ）とすると、統合クロススペクトルＵ（ｋ,ｎ）は次のように計算される。
Figure JPOXMLDOC01-appb-M000044
 　また、総乗を用いた場合、統合クロススペクトルＵ（ｋ,ｎ）は次のように計算される。
Figure JPOXMLDOC01-appb-M000045
 　第５実施形態の統合部１１２２と同様に、目標音源の存在する周波数や、目的音源のパワーが大きい周波数があらかじめ判明している場合には、統合クロススペクトルＵ（ｋ,ｎ）を生成する際に補正してもよい。この場合、第５実施形態と同様に、周波数の選択や重み付けという形で影響度を制御する。例えば目標音の存在する周波数の集合をΩとすると、帯域の選択により統合クロススペクトルＵ（ｋ,ｎ）を求める場合は、次のように計算する。
Figure JPOXMLDOC01-appb-M000046
 　また、重み付けを用いる場合には、Ｕ（ｋ,ｎ）は次のように計算される。
Figure JPOXMLDOC01-appb-M000047
 　ここで、ａとｂとは実数で、ａ＞ｂ＞０を満足する。このように、目的音が存在する周波数の周波数別相関関数を重点的に用いて統合すると、雑音などの非目的音の影響が小さい相関関数を生成できるため、方向推定精度が向上する。
　逆変換部１４２２は、統合部１４２１から供給される統合クロススペクトルの逆変換を行い、統合相関関数として推定方向情報生成部２０５に伝達する。本実施形態でも、逆変換には逆フーリエ変換を用いる方法について説明する。統合部１４２１から供給された統合クロススペクトルをＵ（ｋ,ｎ）とすると、Ｕ（ｋ,ｎ）の逆変換により得られる統合相関関数ｕ（τ,ｎ）は次のように計算される。
Figure JPOXMLDOC01-appb-M000048
 　＜本実施形態の効果＞
　本実施形態によれば、周波数別クロススペクトルを統合してから逆変換を行い、統合相関関数を得る。このため、周波数別クロススペクトルごとに逆変換を行っていた第５実施形態と比べると、逆変換の回数が少なくなる。したがって、第５実施形態よりも少ない計算量で統合相関関数を求めることができる。
　［第８実施形態］
　次に本発明の第８実施形態に係る波源方向推定システムについて、図１５Ａおよび図１５Ｂを用いて説明する。図１５Ａは、本実施形態に係る波源方向推定システム１５００の構成を説明するための図である。本実施形態に係る波源方向推定システム１５００は、上記第２実施形態に係る波源方向推定装置を用いている。したがって、第２実施形態と同じ構成および動作については同じ符号を付してその詳しい説明を省略する。
　本実施形態に係る波源方向推定システム１５００は、マイク１５０１、１５０２、ＡＤ変換部１５０１および表示部１５０２を有する。なお、本実施形態では、波源方向推定装置２００の代わりに波源方向推定装置７００、８００を用いることが可能である。また、波源が音源であるという仮定で説明するため、マイクを用いた例を説明するが、音源以外の場合には、その波源から放射される波動を受信して電気信号に変換できる各種センサが、マイクの代わりに用いられる。
　マイク１５０１、１５０２は、推定対象である目標物体から生じる音を含めた装置周辺の音を電気信号に変換し、ＡＤ変換部１５０１へ伝達する。音が伝わる媒質が空気媒質である場合、音は空気の振動としてマイクに到達する。マイクは、到達した空気の振動を電気信号に変換する。
　ＡＤ変換部１５０１は、マイク１５０１、１５０２から供給された音の電気信号をデジタル信号に変換し、入力端子２０１、２０２に伝達する。
　表示部１５０２は、波源方向推定装置２００から供給された推定方向情報を画像などの可視化データに変換し、ディスプレイなどの表示装置に表示する。基本的な可視化方法は、ある時刻における相関関数を２次元グラフで表示する方法である。その際、横軸に方向、縦軸に相関値を表示する。ある時刻だけでなく、相関関数の時間変化を３次元で表示する方法も有効である。時間変化を表示することにより、目標音源の出現の明確化、目標音源の移動パタン、目標音源の移動方向の予測などが可能になる。３次元ではなく、２次元平面に投影する方法も有効である。３次元だと表示されたときに裏側が見づらいという問題がある。上から投影した平面上に表示すれば、死角が無くなり一覧性が向上する。相関値は、色の濃淡ではなく、等高線で表現してもよい。
　図１５Ｂは、本実施形態に係る波源方向推定システム１５００の表示部１５０２で表示した画像１５１０の一例を示す図であり、波源方向推定装置２００から供給された推定方向情報から得られたものである。これは、本実施形態の効果を確認する目的で取得した。例の作成には、屋外環境において、方位４０度にドローンの飛翔音が発生する状況の音を用いた。数センチ間隔で設置した二つのマイクを用いて集音した。
　図１５Ｂは、色が黒いほど相関値が高いことを表している。方位角の範囲は０から１８０度である。横軸は時刻を表している。図１５Ｂを参照すると、方位４０度の相関値が高くなっていることが分かる。このことから、ドローンの到来方向が約４０度であることが分かる。
　＜本実施形態の効果＞
　本実施形態によれば、波源の方向を高精度に推定することができる。また、推定方向情報を画像などの可視化データとして表示するので、ユーザが波源の方向推定情報を視覚的に把握することができる。
　［他の実施形態］
　以上、実施形態を参照して本願発明を説明したが、本願発明は上記実施形態に限定されるものではない。本願発明の構成や詳細には、本願発明のスコープ内で当業者が理解し得る様々な変更をすることができる。また、それぞれの実施形態に含まれる別々の特徴を如何様に組み合わせたシステムまたは装置も、本発明の範疇に含まれる。
　また、本発明は、複数の機器から構成されるシステムに適用されてもよいし、単体の装置に適用されてもよい。さらに、本発明は、実施形態の機能を実現する情報処理プログラムが、システムあるいは装置に直接あるいは遠隔から供給される場合にも適用可能である。したがって、本発明の機能をコンピュータで実現するために、コンピュータにインストールされるプログラム、あるいはそのプログラムを格納した媒体、そのプログラムをダウンロードさせるＷＷＷ(World Wide Web)サーバも、本発明の範疇に含まれる。特に、少なくとも、上述した実施形態に含まれる処理ステップをコンピュータに実行させるプログラムを格納した非一時的コンピュータ可読媒体（non-transitory computer readable medium）は本発明の範疇に含まれる。
　［実施形態の他の表現］
　上記の実施形態の一部又は全部は、以下の付記のようにも記載されうるが、以下には限られない。
（付記１）
　波源で発生した波を入力信号として取得する複数の入力信号取得手段と、
　複数の前記入力信号取得手段で取得した複数の前記入力信号を変換して、複数の周波数領域信号を得る変換手段と、
　複数の前記周波数領域信号に基づいて、クロススペクトルを計算するクロススペクトル計算手段と、
　前記クロススペクトルの分散を計算する分散計算手段と、
　前記クロススペクトルと前記分散とに基づいて、相関関数を計算する相関関数計算手段と、
　を備える相関関数生成装置。
（付記２）
　前記分散計算手段は、前記クロススペクトルの位相に基づいて分散を計算する付記１に記載の相関関数生成装置。
（付記３）
　前記分散計算手段は、前記クロススペクトルの円周分散を計算する付記１に記載の相関関数生成装置。
（付記４）
　前記分散計算手段は、前記入力信号取得手段で過去に取得された入力信号から変換された周波数領域信号に基づいて計算されたクロススペクトルに基づいて分散を計算する付記１に記載の相関関数生成装置。
（付記５）
　前記相関関数計算手段は、
　　前記分散に基づいて前記クロススペクトルを重み付けする重み付け手段を備え、
　　重み付されたクロススペクトルに基づいて、相関関数を計算する、付記１乃至４のいずれか１項に記載の相関関数生成装置。
（付記６）
　前記クロススペクトルを平均して平均クロススペクトルを計算する平均計算手段をさらに備え、
　前記分散計算手段は、前記平均クロススペクトルから分散を計算する付記１乃至５のいずれか１項に記載の相関関数生成装置。
（付記７）
　前記相関関数計算手段は、
　　前記分散と前記クロススペクトルとに基づいて、複数の第１周波数別クロススペクトルを生成する周波数別クロススペクトル生成手段と、
　　複数の前記第１周波数別クロススペクトルを統合して一つの相関関数を計算する統合相関関数計算手段と、
　を備える付記１乃至６のいずれか１項に記載の相関関数生成装置。
（付記８）
　前記統合相関関数計算手段は、
　　複数の前記第１周波数別クロススペクトルの逆変換を計算し、複数の周波数別相関関数を得る第１逆変換手段と、
　　前記第１逆変換手段により計算された複数の前記周波数別相関関数を統合して一つの相関関数を計算する第１統合手段と、
　を備える付記７に記載の相関関数生成装置。
（付記９）
　前記周波数別クロススペクトル生成手段は、
　　前記分散に基づいて、カーネル関数スペクトルを生成するカーネル関数スペクトル生成手段と、
　　複数の前記第１周波数別クロススペクトルに、前記カーネル関数スペクトルを乗じた複数の第２周波数別クロススペクトルを得る乗算手段と、
　を備える付記７または８に記載の相関関数生成装置。
（付記１０）
　前記周波数別クロススペクトル生成手段は、
　　複数のカーネル関数スペクトルを記憶するカーネル関数スペクトル記憶手段と、
　　前記分散に基づいて、前記カーネル関数スペクトル記憶手段に記憶された複数の前記カーネル関数スペクトルから一つのカーネル関数スペクトルを選択するカーネル関数スペクトル選択手段と、
　をさらに備え、
　前記乗算手段は、複数の前記第１周波数別クロススペクトルに、選択された前記一つのカーネル関数スペクトルを乗じた複数の第２周波数別クロススペクトルを得る、付記９に記載の相関関数生成装置。
（付記１１）
　前記統合相関関数計算手段は、
　　複数の前記第１周波数別クロススペクトルを統合して、一つの統合クロススペクトルを得る第２統合手段と、
　　前記一つの統合クロススペクトルの逆変換を計算して一つの相関関数を得る第２逆変換手段と、
　を備える付記７乃至１０のいずれか１項に記載の相関関数生成装置。
（付記１２）
　波源で発生した波を入力信号として取得する複数の入力信号取得ステップと、
　複数の前記入力信号取得手段で取得した複数の前記入力信号を変換して、複数の周波数領域信号を得る変換ステップと、
　複数の前記周波数領域信号に基づいて、クロススペクトルを計算するクロススペクトル計算ステップと、
　前記クロススペクトルの分散を計算する分散計算ステップと、
　前記クロススペクトルと前記分散とに基づいて、相関関数を計算して生成する相関関数計算ステップと、
　を含む相関関数生成方法。
（付記１３）
　波源で発生した波を入力信号として取得する複数の入力信号取得ステップと、
　複数の前記入力信号取得手段で取得した複数の前記入力信号を変換して、複数の周波数領域信号を得る変換ステップと、
　複数の前記周波数領域信号に基づいて、クロススペクトルを計算するクロススペクトル計算ステップと、
　前記クロススペクトルの分散を計算する分散計算ステップと、
　前記クロススペクトルと前記分散とに基づいて、相関関数を計算して生成する相関関数計算ステップと、
　をコンピュータに実行させる相関関数生成プログラム。
（付記１４）
　付記１乃至１１のいずれか１項に記載の相関関数生成装置を有し、
　前記相関関数生成装置により生成された相関関数に基づいて、前記波源の方向を推定する波源方向推定装置。

Claims (14)

  1. 　波源で発生した波を入力信号として取得する複数の入力信号取得手段と、
    　複数の前記入力信号取得手段で取得した複数の前記入力信号を変換して、複数の周波数領域信号を得る変換手段と、
    　複数の前記周波数領域信号に基づいて、クロススペクトルを計算するクロススペクトル計算手段と、
    　前記クロススペクトルの分散を計算する分散計算手段と、
    　前記クロススペクトルと前記分散とに基づいて、相関関数を計算する相関関数計算手段と、
    　を備える相関関数生成装置。
  2. 　前記分散計算手段は、前記クロススペクトルの位相に基づいて分散を計算する請求項１に記載の相関関数生成装置。
  3. 　前記分散計算手段は、前記クロススペクトルの円周分散を計算する請求項１に記載の相関関数生成装置。
  4. 　前記分散計算手段は、前記入力信号取得手段で過去に取得された入力信号から変換された周波数領域信号に基づいて計算されたクロススペクトルに基づいて分散を計算する請求項１に記載の相関関数生成装置。
  5. 　前記相関関数計算手段は、
    　　前記分散に基づいて前記クロススペクトルを重み付けする重み付け手段を備え、
    　　重み付されたクロススペクトルに基づいて、相関関数を計算する、請求項１乃至４のいずれか１項に記載の相関関数生成装置。
  6. 　前記クロススペクトルを平均して平均クロススペクトルを計算する平均計算手段をさらに備え、
    　前記分散計算手段は、前記平均クロススペクトルから分散を計算する請求項１乃至５のいずれか１項に記載の相関関数生成装置。
  7. 　前記相関関数計算手段は、
    　　前記分散と前記クロススペクトルとに基づいて、複数の第１周波数別クロススペクトルを生成する周波数別クロススペクトル生成手段と、
    　　複数の前記第１周波数別クロススペクトルを統合して一つの相関関数を計算する統合相関関数計算手段と、
    　を備える請求項１乃至６のいずれか１項に記載の相関関数生成装置。
  8. 　前記統合相関関数計算手段は、
    　　複数の前記第１周波数別クロススペクトルの逆変換を計算し、複数の周波数別相関関数を得る第１逆変換手段と、
    　　前記第１逆変換手段により計算された複数の前記周波数別相関関数を統合して一つの相関関数を計算する第１統合手段と、
    　を備える請求項７に記載の相関関数生成装置。
  9. 　前記周波数別クロススペクトル生成手段は、
    　　前記分散に基づいて、カーネル関数スペクトルを生成するカーネル関数スペクトル生成手段と、
    　　複数の前記第１周波数別クロススペクトルに、前記カーネル関数スペクトルを乗じた複数の第２周波数別クロススペクトルを得る乗算手段と、
    　を備える請求項７または８に記載の相関関数生成装置。
  10. 　前記周波数別クロススペクトル生成手段は、
    　　複数のカーネル関数スペクトルを記憶するカーネル関数スペクトル記憶手段と、
    　　前記分散に基づいて、前記カーネル関数スペクトル記憶手段に記憶された複数の前記カーネル関数スペクトルから一つのカーネル関数スペクトルを選択するカーネル関数スペクトル選択手段と、
    　をさらに備え、
    　前記乗算手段は、複数の前記第１周波数別クロススペクトルに、選択された前記一つのカーネル関数スペクトルを乗じた複数の第２周波数別クロススペクトルを得る、請求項９に記載の相関関数生成装置。
  11. 　前記統合相関関数計算手段は、
    　　複数の前記第１周波数別クロススペクトルを統合して、一つの統合クロススペクトルを得る第２統合手段と、
    　　前記一つの統合クロススペクトルの逆変換を計算して一つの相関関数を得る第２逆変換手段と、
    　を備える請求項７乃至１０のいずれか１項に記載の相関関数生成装置。
  12. 　波源で発生した波を入力信号として取得する複数の入力信号取得ステップと、
    　複数の前記入力信号取得手段で取得した複数の前記入力信号を変換して、複数の周波数領域信号を得る変換ステップと、
    　複数の前記周波数領域信号に基づいて、クロススペクトルを計算するクロススペクトル計算ステップと、
    　前記クロススペクトルの分散を計算する分散計算ステップと、
    　前記クロススペクトルと前記分散とに基づいて、相関関数を計算して生成する相関関数計算ステップと、
    　を含む相関関数生成方法。
  13. 　波源で発生した波を入力信号として取得する複数の入力信号取得ステップと、
    　複数の前記入力信号取得手段で取得した複数の前記入力信号を変換して、複数の周波数領域信号を得る変換ステップと、
    　複数の前記周波数領域信号に基づいて、クロススペクトルを計算するクロススペクトル計算ステップと、
    　前記クロススペクトルの分散を計算する分散計算ステップと、
    　前記クロススペクトルと前記分散とに基づいて、相関関数を計算して生成する相関関数計算ステップと、
    　をコンピュータに実行させる相関関数生成プログラム。
  14. 　請求項１乃至１１のいずれか１項に記載の相関関数生成装置を有し、
    　前記相関関数生成装置により生成された相関関数に基づいて、前記波源の方向を推定する波源方向推定装置。
PCT/JP2017/000680 2017-01-11 2017-01-11 相関関数生成装置、相関関数生成方法、相関関数生成プログラムおよび波源方向推定装置 WO2018131099A1 (ja)

Priority Applications (3)

Application Number Priority Date Filing Date Title
JP2018561142A JP6769495B2 (ja) 2017-01-11 2017-01-11 相関関数生成装置、相関関数生成方法、相関関数生成プログラムおよび波源方向推定装置
PCT/JP2017/000680 WO2018131099A1 (ja) 2017-01-11 2017-01-11 相関関数生成装置、相関関数生成方法、相関関数生成プログラムおよび波源方向推定装置
US16/476,897 US11336997B2 (en) 2017-01-11 2017-01-11 Correlation function generation apparatus, correlation function generation method, correlation function generation program, and wave source direction estimation apparatus

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
PCT/JP2017/000680 WO2018131099A1 (ja) 2017-01-11 2017-01-11 相関関数生成装置、相関関数生成方法、相関関数生成プログラムおよび波源方向推定装置

Publications (1)

Publication Number Publication Date
WO2018131099A1 true WO2018131099A1 (ja) 2018-07-19

Family

ID=62840122

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/JP2017/000680 WO2018131099A1 (ja) 2017-01-11 2017-01-11 相関関数生成装置、相関関数生成方法、相関関数生成プログラムおよび波源方向推定装置

Country Status (3)

Country Link
US (1) US11336997B2 (ja)
JP (1) JP6769495B2 (ja)
WO (1) WO2018131099A1 (ja)

Cited By (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2020003342A1 (ja) * 2018-06-25 2020-01-02 日本電気株式会社 波源方向推定装置、波源方向推定方法、およびプログラム記録媒体
WO2021044470A1 (ja) * 2019-09-02 2021-03-11 日本電気株式会社 波源方向推定装置、波源方向推定方法、およびプログラム記録媒体
WO2023157159A1 (ja) * 2022-02-17 2023-08-24 日本電信電話株式会社 位相差スペクトル推定方法、チャネル間関係情報推定方法、信号符号化方法、信号処理方法、これらの装置、プログラム
WO2023181144A1 (ja) * 2022-03-23 2023-09-28 三菱電機株式会社 ノイズ除去装置及び方法
WO2024053353A1 (ja) * 2022-09-08 2024-03-14 パナソニック インテレクチュアル プロパティ コーポレーション オブ アメリカ 信号処理装置、及び、信号処理方法

Families Citing this family (3)

* 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
US10629226B1 (en) * 2018-10-29 2020-04-21 Bestechnic (Shanghai) Co., Ltd. Acoustic signal processing with voice activity detector having processor in an idle state
US11392462B2 (en) 2019-09-26 2022-07-19 Gulfstream Aerospace Corporation Voting of triple redundant circular data

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2011033717A (ja) * 2009-07-30 2011-02-17 Secom Co Ltd 雑音抑圧装置
JP2012185009A (ja) * 2011-03-04 2012-09-27 Mitsubishi Electric Corp 方位探知装置
WO2013042201A1 (ja) * 2011-09-20 2013-03-28 トヨタ自動車株式会社 音源検出装置

Family Cites Families (3)

* 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
JP5018643B2 (ja) 2008-05-27 2012-09-05 三菱電機株式会社 方位探知装置
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

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2011033717A (ja) * 2009-07-30 2011-02-17 Secom Co Ltd 雑音抑圧装置
JP2012185009A (ja) * 2011-03-04 2012-09-27 Mitsubishi Electric Corp 方位探知装置
WO2013042201A1 (ja) * 2011-09-20 2013-03-28 トヨタ自動車株式会社 音源検出装置

Cited By (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2020003342A1 (ja) * 2018-06-25 2020-01-02 日本電気株式会社 波源方向推定装置、波源方向推定方法、およびプログラム記録媒体
WO2021044470A1 (ja) * 2019-09-02 2021-03-11 日本電気株式会社 波源方向推定装置、波源方向推定方法、およびプログラム記録媒体
JPWO2021044470A1 (ja) * 2019-09-02 2021-03-11
JP7276469B2 (ja) 2019-09-02 2023-05-18 日本電気株式会社 波源方向推定装置、波源方向推定方法、およびプログラム
WO2023157159A1 (ja) * 2022-02-17 2023-08-24 日本電信電話株式会社 位相差スペクトル推定方法、チャネル間関係情報推定方法、信号符号化方法、信号処理方法、これらの装置、プログラム
WO2023181144A1 (ja) * 2022-03-23 2023-09-28 三菱電機株式会社 ノイズ除去装置及び方法
WO2024053353A1 (ja) * 2022-09-08 2024-03-14 パナソニック インテレクチュアル プロパティ コーポレーション オブ アメリカ 信号処理装置、及び、信号処理方法

Also Published As

Publication number Publication date
US20190335273A1 (en) 2019-10-31
JPWO2018131099A1 (ja) 2019-11-07
JP6769495B2 (ja) 2020-10-14
US11336997B2 (en) 2022-05-17

Similar Documents

Publication Publication Date Title
WO2018131099A1 (ja) 相関関数生成装置、相関関数生成方法、相関関数生成プログラムおよび波源方向推定装置
US10515650B2 (en) Signal processing apparatus, signal processing method, and signal processing program
WO2018003158A1 (ja) 相関関数生成装置、相関関数生成方法、相関関数生成プログラムおよび波源方向推定装置
US20090141908A1 (en) Distance based sound source signal filtering method and apparatus
WO2017216999A1 (ja) 波源方向推定装置、波源方向推定システム、波源方向推定方法および波源方向推定プログラム
CN102147458B (zh) 一种针对宽带声源的波达方向估计方法及其装置
EP3227704B1 (en) Method for tracking a target acoustic source
EP3260875A1 (en) An automotive testing method, system and computer program product
EP2878515B1 (en) Generating an audio signal with a configurable distance cue
CN110047507B (zh) 一种声源识别方法及装置
KR101581885B1 (ko) 복소 스펙트럼 잡음 제거 장치 및 방법
WO2016119388A1 (zh) 一种基于语音信号构造聚焦协方差矩阵的方法及装置
JP2007336232A (ja) 特定方向収音装置、特定方向収音プログラム、記録媒体
Al-Aboosi et al. Improved underwater signal detection using efficient time–frequency de-noising technique and Pre-whitening filter
CN105144290A (zh) 信号处理装置、信号处理方法和信号处理程序
Hosseini et al. Time difference of arrival estimation of sound source using cross correlation and modified maximum likelihood weighting function
JP6933303B2 (ja) 波源方向推定装置、波源方向推定方法、およびプログラム
JP4984668B2 (ja) ソーナーシステムおよびその位相誤差補正方法
JP5713933B2 (ja) 音源距離測定装置、音響直間比推定装置、雑音除去装置、それらの方法、及びプログラム
Pici et al. Design of spectrally adaptive noise radar waveforms
CN114814728A (zh) 一种声源定位方法、系统、电子设备及介质
CN117198276A (zh) 基于毫米波信号非接触式耳机语音感知方法及系统
KR20230108435A (ko) 인공지능 신경망 학습방법 및 시스템
Chen et al. A Fast Estimation Method for 3-D Acoustic Source Localization
Mizumachi Statistical confidence measure for direction-of-arrival estimate

Legal Events

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

Ref document number: 17891785

Country of ref document: EP

Kind code of ref document: A1

ENP Entry into the national phase

Ref document number: 2018561142

Country of ref document: JP

Kind code of ref document: A

NENP Non-entry into the national phase

Ref country code: DE

122 Ep: pct application non-entry in european phase

Ref document number: 17891785

Country of ref document: EP

Kind code of ref document: A1