JP6183067B2 - Data analysis apparatus and method, program, and recording medium - Google Patents

Description

  The present invention relates to a data analysis apparatus and method, a program, and a recording medium, and more particularly to detection of non-Gaussian or non-linearity of time series data.

The non-Gaussian property of time series data means that the amplitude distribution of the time series data is not symmetrical with respect to the average value. On the other hand, the non-linearity of time series data means that a non-linear term is included in the generation process of time series data. Nonlinear data has non-Gaussian properties.
As a data analysis method for confirming non-Gaussianity and nonlinearity of acquired data (observation data), it has been proposed to use a normalized bispectrum related to the skewness of the data (Non-Patent Documents 1 to 3). ).

  The problem with the conventional method is that even if the observation data is Gaussian, if finite data is extracted, the data may appear to have non-Gaussian or non-linearity depending on the extraction method, and erroneous determination may occur. That is.

Hereinafter, a conventional method will be described with reference to FIG. In FIG. 1, time series data obtained as a result of observation is represented as x (t), and a discretized version of x (t) is represented as x [n]. Here, n is an index representing time.
As shown in FIG. 2, a plurality of data blocks DB _{1 to} DB _P are generated from the time series data x [n]. Each of the data blocks DB _i (i is any one of 1 to P) includes N sample values.
N-sample discrete Fourier transform is performed on each data block DB _i by the Fourier transform unit 14 to obtain a component X _i [k] of each frequency. The Fourier transform is represented by the following formula (1).

ここで、ｋは周波数を表すインデックスである。

Here, k is an index representing a frequency.

The bispectrum calculation unit 22 and the power spectrum calculation unit 24 are supplied with frequency components X _i [k] for a plurality of frequencies constituting Q frequency pairs. For example, Q is 49, and the arrangement of the 49 frequency pairs on the two-dimensional frequency plane is as indicated by “◯” in FIG.

但し、Ｘ_ｉ ^＊［ｋ_１＋ｋ_２］はＸ_ｉ［ｋ_１＋ｋ_２］の複素共役を表す。 The bispectrum calculation unit 22 calculates the bispectrum Ba [k ₁ , k ₂ ] using the frequency component X _i [k] for each of the data blocks DB _{1 to} DB _P.
For this purpose, first, a bispectrum B _i [k ₁ , k ₂ ] for each data block (i-th data block) DB _i is calculated by the following equation (2).

^{_{However, X i * [k 1 +}} k 2] represents the complex conjugate of _{X i [k 1 + k 2} ].

The bispectrum B _i [k ₁ , k ₂ ] obtained by Expression (2) is an index indicating the degree of correlation between the frequencies of the observation data x [n].

The bispectrum calculation unit 22 performs the calculation of Expression (2) for each of the _P data blocks DB _{1 to} DB _P , and calculates the bispectrum B _i [k ₁ , k ₂ ] obtained as a result of the calculation. By averaging, the bispectrum Ba [k ₁ , k ₂ ] for the time-series data DS composed of P data blocks is obtained. The calculation for obtaining the bispectrum Ba [k ₁ , k ₂ ] is expressed by the following equation (3).

  If the observation data contains only white noise, the bispectrum has a small value because there is no correlation between the frequencies of the observation data. On the other hand, if the observation data includes a non-stationary component, a correlation appears between the frequencies of the observation data, and the bispectrum expressed by Equation (3) has a large value.

  The bispectrum represented by Expression (3) has a value depending on the magnitude (power) of each frequency component included in the observation data. For this reason, it is difficult to determine the magnitude of correlation between frequencies in the bispectrum itself. Therefore, normally, a normalized value is used in which the magnitude of the bispectrum does not depend on the magnitude of each frequency component of the observation data.

ここで、
Ｘ_ｉ ^＊［ｋ］は、Ｘ_ｉ［ｋ］の複素共役である。
式（４）の計算は、周波数ｋ_１、ｋ_２、ｋ_３（＝ｋ_１＋ｋ_２）の各々について行われ、従って、パワースペクトル計算部２４からは、Ｓａ［ｋ_１］、Ｓａ［ｋ_２］、Ｓａ［ｋ_１＋ｋ_２］が出力される。 For normalization, the power spectrum calculation unit 24 calculates the power spectrum Sa [k] of the observation data as follows.

here,
X _i ^* [k] is a complex conjugate of X _i [k].
The calculation of Expression (4) is performed for each of the frequencies k ₁ , k ₂ , and k ₃ (= k ₁ + k ₂ ). Therefore, the power spectrum calculation unit 24 receives Sa [k ₁ ] and Sa [k ₂ ], Sa [k ₁ + k ₂ ] are output.

The normalization unit 27 includes the bispectrum Ba [k ₁ , k ₂ ] calculated by the bispectrum calculation unit 22, the power spectra Sa [k ₁ ] and Sa [k ₂ ] calculated by the power spectrum calculation unit 24, and Using Sa [k ₁ + k ₂ ], a normalized bispectrum Bn [k ₁ , k ₂ ] is calculated by the following equation (5).

By having a power spectrum term in the denominator, the size of each frequency component of the observation data is normalized. Therefore, the normalized bispectrum Bn [k ₁ , k ₂ ] in Expression (5) does not depend on the size (power) of the observation data, and shows a value corresponding to the correlation between the frequencies.
However, in the actual calculation, the corrected standard obtained by correcting as in the following equation (6) in order to set the variance of the bispectrum Bn [k ₁ , k ₂ ] in the equation (5) to 1. Bispectrum Bnc [k ₁ , k ₂ ] is used.

The normalized bispectrum Bnc [k ₁ , k ₂ ] represented by the equation (6) is a value calculated for a specific frequency pair (k ₁ , k ₂ ), and the processing described above for Q frequency pairs. Is repeated to obtain the normalized bispectrum Bnc [k ₁ , k ₂ ] expressed by the equation (6) for each of the Q frequency pairs.

The verification amount calculation unit 29 detects non-Gaussianity and nonlinearity of the observation data by using the normalized bispectrum Bnc [k ₁ , k ₂ ] obtained for the Q frequency pairs as described above.

If the observation data is white noise, the normalized bispectrum Bnc [k ₁ , k ₂ ] obtained for each of Q frequency pairs (k ₁ , k ₂ ) follows a normal distribution with an average of 0 and a variance of 1. In order to detect non-stationarity of the observation data using this property, the following test amount z is calculated.

In Equation (7), Bnc _q (where q is any one of 1 to Q) represents the same as Bnc [k ₁ , k ₂ ] for each of the Q frequency pairs (k ₁ , k ₂ ).
The numerator of Equation (7) represents adding Q normalized bispectral values. The denominator of Equation (7) is a term for normalizing the variance of the test amount z.
If the observation data contains only white noise and Q is a sufficiently large value, the test amount z follows a normal distribution with an average of 0 and a variance of 1. If the observed data contains non-Gaussian or non-linear components, the value of the normalized bispectrum increases, and as a result, the test quantity z does not follow a normal distribution with mean 0 and variance 1. Large value. By utilizing this fact, non-Gaussian or non-linearity of the observation data can be detected.
The above is a conventional method for detecting non-Gaussianity or non-linearity from observation data using a normalized bispectrum.

Melvin J.M. Hinich, “Bisspectrum of ship-radiated noise”, J. Am. Acoustic. Soc. Am. 85 (4), April 1989 Melvin J.M. Hinich, “Detecting a Transient Signal by Bidirectional Analysis”, IEEE Trans. On Acoustics, Speech and Signal Processing, Vol. 36, no. 7, July 1990 Melvin J.M. Hinich et. al, “On the Principal Domain of the Discrete Biscopy of a Stationary Signal”, IEEE Trans. On Signal Processing, Vol. 43, no. 9, September 1995

  The above-described conventional method is processing focusing on the fact that the observation data has a certain degree of distortion if the observation data has non-Gaussian or non-linearity. However, even if the observation data has Gaussian or linearity, if finite data is extracted, it may appear to have non-Gaussian or non-linearity depending on how it is extracted, and therefore an error may occur in the determination.

  The present invention is for solving the above-described problems, and an object of the present invention is to make it possible to obtain a processing result in which the influence of the size and distortion of observation data is suppressed.

（但し、
Ｑは前記周波数ペアの数、
γ _ｑ は、前記周波数ペアの各々について、前記バイコヒーレンス計算部で計算されたバイコヒーレンス、
ｂｉａｓは、前記検定量の平均を０とするための調整量、
σ ^２ は、前記検定量の分散を補正するための調整量）
により前記検定量ｚ _γ を計算するものである。 The data analysis device according to the first aspect of the present invention includes:
A Fourier transform unit that Fourier-transforms a plurality of data blocks each composed of successive sample values in time-series data representing observation values, and calculates a plurality of frequency components for each data block;
A bicoherence calculator that calculates the bicoherence of each frequency pair consisting of two frequencies of the frequency components for each data block calculated by the Fourier transform unit;
A test amount calculator for calculating a test amount for the time series data based on the bicoherence for each of the frequency pairs;
A comparison unit that determines whether or not the verification amount calculated by the verification amount calculation unit is greater than a predetermined threshold;
The determination result by the comparison unit is output as the analysis result ,
The test amount calculation unit, the average value and variance of the test amount, the number of sample values constituting each data block, the number of data blocks constituting the time series data, and
The verification amount is calculated using an adjustment amount for keeping the data block constant regardless of the overlap ratio of the data block .
The data analysis apparatus according to the second aspect of the present invention includes:
A Fourier transform unit that Fourier-transforms a plurality of data blocks each composed of successive sample values in time-series data representing observation values, and calculates a plurality of frequency components for each data block;
A bicoherence calculator that calculates the bicoherence of each frequency pair consisting of two frequencies of the frequency components for each data block calculated by the Fourier transform unit;
A test amount calculator for calculating a test amount for the time series data based on the bicoherence for each of the frequency pairs;
A comparison unit that determines whether or not the verification amount calculated by the verification amount calculation unit is greater than a predetermined threshold;
The determination result by the comparison unit is output as the analysis result,
The test statistic calculation unit, when representing the test amount z _gamma,

(However,
Q is the number of frequency pairs,
γ _q is the bicoherence calculated by the bicoherence calculator for each of the frequency pairs,
bias is an adjustment amount for setting the average of the test amount to 0,
σ ² is an adjustment amount for correcting the variance of the test amount)
_{Is used} to calculate the test amount zγ .

  According to the present invention, non-Gaussianity and nonlinearity of observation data are detected using bicoherence, which is an index representing the phase-wise correlation between frequencies of observation data, and using a test quantity based on bicoherence. Therefore, it is possible to accurately determine non-Gaussian or non-linearity while suppressing the influence of the size and skewness of the observation data.

It is a block diagram which shows the outline of the apparatus for enforcing the conventional data analysis method. It is a figure which shows the production | generation method of the data block from time series data. It is a figure which shows an example of arrangement | positioning on the two-dimensional frequency plane of the frequency pair used for a data analysis. It is a block diagram which shows the outline of the data analyzer of this invention. FIG. 5 is a block diagram illustrating a configuration example of a detection processing unit of the data analysis apparatus in FIG. 4. It is a block diagram which shows the structural example of the coherence calculation part of the detection process part of FIG. It is a block diagram which shows the structural example of the apparatus which produces | generates the parameter used by the detection process part of FIG. (A) And (b) is a figure which shows the waveform which added the white noise to the simulated transient signal and the simulation signal. (A) And (b) is a figure which shows the result of having applied the conventional method and the method of this invention with respect to the signal shown by FIG.8 (b).

  FIG. 4 shows a time-series data analysis apparatus according to the present invention. Hereinafter, as an example, a case where an acoustic signal in an audible range is received and analyzed will be described. The illustrated time-series data analysis apparatus includes an input signal processing unit 2, a time-series data memory 4, and a detection processing unit 6.

  The input signal processing unit 2 includes an A / D converter 3 and is connected to a microphone 1 that converts an acoustic signal into an electric signal. The analog signal x (t) given from the microphone 1 is sampled and digitized. Thus, time-series data (time-series data representing observation values) x [n] is generated. The sampling frequency of the A / D converter 3 is, for example, 50 kHz, and the time series data output from the A / D converter 3 is represented by x [n]. n is an index indicating time.

  The time series data memory 4 stores time series data x [n] output from the input signal processing unit 2.

As shown in FIG. 5, the detection processing unit 6 includes a data block generation unit 12, a Fourier transform unit 14, a frequency component memory 16, a frequency component selection unit 18, a bicoherence calculation unit 20, and a test amount calculation unit. 32, a comparison unit 34, a parameter memory 36, and an adjustment amount calculation unit 38. The same reference numerals as those in FIG. 1 denote similar blocks. Although the data block generation unit 12, the frequency component memory 16, and the frequency component selection unit 18 are not shown in FIG. 1, the same methods are actually used in the conventional method.

The data block generation unit 12 includes a plurality of data blocks DB _i (i = 1 to 1) each composed of a plurality of sample values from the time series data x [n] stored in the time series data memory 4. P) is generated and supplied to the Fourier transform unit 14.
For example, as shown in FIG. 2 , the data block DB _i is composed of N sample values that are continuous with each other, and the overlap rate L is 50%. That is, the N / 2 sample values in the second half of each data block (i-th data block) and the N / 2 sample values in the first half of the next data block (i + 1-th data block) are the same. Is.

The Fourier transform unit 14 performs a discrete Fourier transform on the N sample values (observation data) x [0] to x [N−1] of each of the _P data blocks DB _{1 to} DB _P , and performs transformation. The resulting frequency components X ₁ [k] to X _P [k] are output. k is an index representing a frequency. The result (Frequency k component) of the Fourier transform for each data block (i-th data block) DB _i (i is any one of 1 to P) is represented by X _i [k].

The frequency component X _i [k] obtained as a result of the Fourier transform is represented by the following equation (1). Equation (1) is the same as equation (1) described above with respect to the conventional method, but will be described again. The same applies to other mathematical expressions below.

Here, k is an index representing one of M frequency values f ₁ , f ₂ ,... F _M. f _{1 to} f _M are generalized and represented by f _m (m is _{1 to M} ).
Frequency _{f m} has the following values.
When m = 1, f _m = 0 (8A)
When 2 ≦ m <M, f _m = f _(m−1) + ΔF (8B)
However, ΔF = 1562.5 Hz, which is equal to 50 kHz (sampling frequency) divided by N = 32 (number of sample values in each data block).
For example, M is 21.

The frequency component X _i [k] (k = f _{1 to} f _M ) for each of the plurality of data blocks DB _i (i = 1 to P) calculated by the Fourier transform unit 14 is stored in the frequency component memory 16. Is done.

The frequency component selection unit 18 uses a frequency component X _i [k] (i = _{1 to} P, k = f ₁ ) for each of the plurality of data blocks DB _i (i = 1 to P) stored in the frequency component memory 16. ˜F _M ) that are used for bicoherence calculation are sequentially read out and supplied to the bicoherence calculation unit 20.

  The bicoherence calculation in the bicoherence calculation unit 20 is sequentially performed for a plurality of, that is, Q frequency pairs. Q is 49, for example, and the arrangement of the 49 frequency pairs on the two-dimensional frequency plane is, for example, as shown by “◯” in FIG.

When calculating the bicoherence for the frequency pair (k ₁ , k ₂ ), the frequency components X _i [k ₁ ], X _i of the frequencies k ₁ , k ₂ for each DB _i of a plurality, that is, P data blocks. Not only [k ₂ ] but also a frequency component X _i [k ₃ ] having a frequency k ₃ equal to the sum of k ₁ and k ₂ is used. Therefore, the frequency component selection unit 18 bicoherences the frequency components X _i [k ₁ ], X _i [k ₂ ], and X _i [k ₃ ] of the frequencies k ₁ , k ₂ , and k ₃ for each data block DB _i. This is supplied to the calculation unit 20.
As described above,
k ₃ = k ₁ + k ₂
In the following, “k ₃ ” is expressed as “k ₁ + k ₂ ”.
On the other hand, the symbol “k” is used in the description of an arbitrary frequency or the description applicable to all of k ₁ , k ₂ , and k ₃ .

The bicoherence calculation unit 20 uses the frequency components X _i [k ₁ ], X _i [k ₂ ], and X _i [k ₁ + k ₂ ] supplied from the frequency component selection unit 18 for each frequency pair (k ₁ , k ₂ ). ] Is used to calculate the bicoherence γ [k ₁ , k ₂ ] of the frequency pair (k ₁ , k ₂ ). As will be described in detail below, Q bicoherence values are obtained for time-series data DS composed of P data blocks. The verification amount calculation unit 32 obtains a verification amount z _γ for the time series data DS based on Q bicoherence values. The comparison unit 34 compares the obtained test amount z _γ with a predetermined threshold value zth.

  As shown in FIG. 6, for example, the bicoherence calculation unit 20 includes a bispectrum calculation unit 22, a power spectrum calculation unit 24, a product root mean square calculation unit 26, and a normalization unit 28.

  The bispectrum calculation unit 22, the power spectrum calculation unit 24, and the product root mean square calculation unit 26 receive the frequency component supplied from the frequency component selection unit 18 and calculate the bispectrum, power spectrum, and product square mean, respectively. .

When calculating the bicoherence of the frequency pair (k ₁ , k ₂ ), the bispectrum calculation unit 22 includes frequency components X _i [k ₁ ], X _i [k ₂ ], and X _i [k ₁ + k ₂ ]. (I = 1 to P) is supplied, the frequency component X _i [k ₁ + k ₂ ] (i = ₁ to P) is supplied to the power spectrum calculation unit 24, and the frequency square mean calculation unit 26 has a frequency Components X _i [k ₁ ], X _i [k ₂ ] (i = 1 to P) are supplied.

ここで、
Ｘ_ｉ［ｋ_１］はｉ番目のデータブロックの周波数ｋ_１の周波数成分であり、
Ｘ_ｉ［ｋ_２］はｉ番目のデータブロックの周波数ｋ_２の周波数成分であり、
Ｘ_ｉ ^＊［ｋ_１＋ｋ_２］はＸ_ｉ［ｋ_１＋ｋ_２］の複素共役を表し、
Ｘ_ｉ［ｋ_１＋ｋ_２］はｉ番目のデータブロックの周波数ｋ_１＋ｋ_２の周波数成分である。 The bispectrum calculation unit 22 calculates the bispectrum of the frequency pair (k ₁ , k ₂ ) for each of the P data blocks, and further obtains an average of the bispectrum for the P data blocks, The average is output as a bispectrum for time-series data DS composed of P data blocks. That is, first, a bispectrum B _i [k ₁ , k ₂ ] for each data block (i-th data block) DB _i is calculated by the following equation (2).

here,
X _i [k ₁ ] is a frequency component of the frequency k ₁ of the i-th data block,
X _i [k ₂ ] is a frequency component of the frequency k ₂ of the i-th data block,
X _i ^* [k ₁ + k ₂ ] represents the complex conjugate of X _i [k ₁ + k ₂ ],
X _i [k ₁ + k ₂ ] is a frequency component of the frequency k ₁ + k ₂ of the i-th data block.

The bispectrum B _i [k ₁ , k ₂ ] obtained by Expression (2) is an index indicating the degree of correlation between the frequencies of the observation data x [n].

The bispectrum calculation unit 22 performs the calculation of Expression (2) for each of the P data blocks DB _i (i = 1 to P), and the P bispectrum B _i [k] obtained as a result of the calculation. ₁ , k ₂ ] is averaged to obtain a bispectrum Ba [k ₁ , k ₂ ] for time-series data DS composed of P data blocks. A calculation formula for obtaining the bispectrum Ba [k ₁ , k ₂ ] is expressed by the following formula (3).

When the formula (2) and the formula (3) are combined and rewritten, the following formula (9) is obtained.

ここで、
Ｘ_ｉ［ｋ_１＋ｋ_２］はｉ番目のデータブロックの、周波数ｋ_１＋ｋ_２の周波数成分であり、
Ｘ_ｉ ^＊［ｋ_１＋ｋ_２］は、Ｘ_ｉ［ｋ_１＋ｋ_２］の複素共役である。 The power spectrum calculation unit 24 has a frequency (k ₁ + k ₂ ) equal to the sum of _two frequencies k ₁ and k ₂ constituting each frequency pair (k ₁ , k ₂ ) for each of the P data blocks. By calculating the power spectrum of the frequency component X _i [k ₁ + k ₂ ] and obtaining the average of the power spectrum for P data blocks, the power for time-series data DS composed of P data blocks is calculated. A spectrum Sa [k ₁ + k ₂ ] is obtained. That is, the calculation for obtaining the power spectrum Sa [k ₁ + k ₂ ] for the time series data DS is expressed by the following equation (10).

here,
X _i [k ₁ + k ₂ ] is a frequency component of the frequency k ₁ + k ₂ of the i-th data block,
X _i ^* [k ₁ + k ₂ ] is a complex conjugate of X _i [k ₁ + k ₂ ].

ここで、
Ｘ_ｉ［ｋ_１］はｉ番目のデータブロックの周波数ｋ_１の周波数成分であり、
Ｘ_ｉ［ｋ_２］はｉ番目のデータブロックの周波数ｋ_２の周波数成分である。 The product root mean square calculation unit 26, for each of the P data blocks, frequency components X _i [k ₁ ], X of the _two frequencies k ₁ , k ₂ constituting each frequency pair (k ₁ , k ₂ ). The square of the absolute value of the product of _i [k ₂ ] is obtained, the mean of the square is obtained for P data blocks, and the mean is the product square mean Ca for time-series data DS composed of P data blocks. Output as [k ₁ , k ₂ ]. The calculation for obtaining the product root mean square Ca [k ₁ , k ₂ ] is expressed by the following equation (11).

here,
X _i [k ₁ ] is a frequency component of the frequency k ₁ of the i-th data block,
X _i [k ₂ ] is a frequency component of the frequency k ₂ of the i-th data block.

この計算は、バイスペクトルＢａ［ｋ_１，ｋ_２］を、積自乗平均Ｃａ［ｋ_１，ｋ_２］とパワースペクトルＳａ［ｋ_１＋ｋ_２］の積の平方根で割ることで規格化することを意味する。
規格化部２８で算出されたバイコヒーレンスγ［ｋ_１，ｋ_２］は、バイコヒーレンス計算部２０の出力となる。
式（１２）で求められるバイコヒーレンスγ［ｋ_１，ｋ_２］は、周波数ペアごとに計算され、周波数ペアを構成する２つの周波数ｋ_１及びｋ_２の成分間の位相的な相関関係を表す複素数であり、その振幅が０〜１の間になるように規格化されており、周波数ｋ_１とｋ_２の関数である。式（５）または（６）で求められる規格化バイスペクトルは、振幅の大きさが観測データの歪度に依存している。歪度は、観測データの抽出の仕方によって大きさが変化する場合があるので、結果的に検定量の大きさに影響を及ぼす。
一方、式（１２）で求められるバイコヒーレンスは、周波数ｋ_１及びｋ_２の成分間の位相的な相関関係のみを抽出し、観測データの抽出の仕方によって大きさが変化する可能性のある歪度の影響を抑制している。 The normalization unit 28 includes the bispectrum Ba [k ₁ , k ₂ ] calculated by the bispectrum calculation unit 22, the power spectrum Sa [k ₁ + k ₂ ] calculated by the power spectrum calculation unit 24, and the product root mean square. Bicoherence γ [k ₁ , k ₂ ] is calculated based on the product root mean square Ca [k ₁ , k ₂ ] calculated by the calculation unit 26. The calculation is represented by the following equation (12).

This calculation normalizes the bispectrum Ba [k ₁ , k ₂ ] by dividing it by the square root of the product of the product square mean Ca [k ₁ , k ₂ ] and the power spectrum Sa [k ₁ + k ₂ ]. means.
The bicoherence γ [k ₁ , k ₂ ] calculated by the normalization unit 28 becomes the output of the bicoherence calculation unit 20.
The bicoherence γ [k ₁ , k ₂ ] obtained by the equation (12) is calculated for each frequency pair, and represents a phase correlation between the components of the _two frequencies k ₁ and k ₂ constituting the frequency pair. It is a complex number, normalized so that its amplitude is between 0 and _1, and is a function of the frequencies k ₁ and k ₂ . In the normalized bispectrum obtained by the equation (5) or (6), the magnitude of the amplitude depends on the skewness of the observation data. Since the magnitude of the skewness may vary depending on how the observation data is extracted, the result affects the magnitude of the test amount.
On the other hand, the bicoherence obtained by Equation (12) is a distortion that extracts only the topological correlation between the components of the frequencies k ₁ and k ₂ , and the magnitude of which may change depending on how the observation data is extracted. The influence of the degree is suppressed.

The verification amount calculation unit 32 calculates the bicoherence γ [k] for the time series data DS composed of P data blocks calculated by the bicoherence calculation unit 20 for each of the Q frequency pairs (k ₁ , k ₂ ). ₁ , k ₂ ] is integrated for all of the Q frequency pairs to calculate the test amount z _γ . In this calculation, the average value and variance of the calculated test amount z _γ are the number N of sample values constituting each data block, the number P of data blocks constituting the time series data DS, and the overlap ratio of the data blocks. This is performed using an adjustment amount (bias, σ ² below) for keeping constant regardless of L, and is expressed by the following equation (13).

In Expression (13), γ _q (where q is any one of 1 to Q) represents the same as γ [k ₁ , k ₂ ] for each of the Q frequency pairs (k ₁ , k ₂ ). .
In the bias term of the numerator of formula (13), the average value of the numerator is zero. Therefore, an adjustment term for setting the average value of z _γ to 0, and σ ² of the denominator are adjustment terms for correcting the variance of the test amount z _γ , and the adjustment amount calculation unit 38 uses the following formulas (14 ) And Equation (15).

In the equations (14) and (15), α and β are constants and are stored in the parameter memory 36.
For example, when calculating the N = 32 sample discrete Fourier transform of P = 32 data blocks with an overlap rate L of 50%, setting α = 2 and β = 6 will make the observed data Gaussian. And having a linearity, the test amount z _γ follows a normal distribution with a mean of 0 and a variance of 1.
This sets the threshold value with respect to test statistic z _gamma, when test statistic z _gamma is greater than the threshold, the observed data has a non-Gaussian, or observed data is determined to have a nonlinearity that Is possible.

The comparison unit 34 determines whether or not the verification amount z _γ calculated by the verification amount calculation unit 32 is larger than a predetermined threshold value zth, and outputs a determination result DJ. That is, if it is larger than the threshold value zth, a determination result indicating that there is non-Gaussianity or non-linearity is output, and if it is equal to or less than the threshold value zth, a determination result indicating that there is no non-Gaussian property or non-linearity is output. The determination result DJ is output as a result of analysis by the data analysis apparatus, and the result of analysis is displayed on the display unit 8 shown in FIG. 4, for example.

As described above, the bias term of the numerator of Equation (13) is an adjustment term for setting the average value of the numerator to 0, and σ ² of the denominator is an adjustment term for correcting the variance of the test amount z _γ. The two adjustment terms are used to keep the magnitude of the test quantity z _γ constant for Gaussian and linear signals even when N (number of samples in the time direction) and P (number of data blocks) are changed. Yes, it is obtained by the above equations (14) and (15).
That is, the parameters α and β need to be determined according to P, N, and L.

  Α and β in the equations (14) and (15) are determined as follows. That is, the parameters α and β are obtained experimentally by using pseudo white noise and stored in the parameter memory 36.

Specifically, bicoherence γw [k ₁ , k ₂ ] is calculated with respect to pseudo white noise using equation (12), and a test amount (test) calculated using the same equation as equation (13) using this is calculated. Test amount) Parameters α and β are determined from the average value zwa and variance zwv of zw.

  The parameter generator shown in FIG. 7 is used for determining the parameters α and β. Among the parameter generation devices, the same reference numerals as those in FIGS. 4 and 5 have the same functions, and can be shared with the devices in FIGS. 4 and 5.

(A) First, the white noise generation unit 42 generates time-series data representing the pseudo white noise (represented by the same code x [n] as the observation data). The generation of white noise can be performed by software, for example. The generated time-series data x [n] representing white noise is stored in the time-series data memory 4 instead of the observation data (time-series data representing observation values) described with reference to FIG.

(B) The data block generator 12 reads out the P data blocks each having N sample values from the time series data x [n] stored in the time series data memory 4 with the overlap ratio L, and Fourier This is supplied to the conversion unit 14.
The Fourier transform unit 14 and the bicoherence calculation unit 20 perform the same processing as the observation data described with reference to FIG. 4 on the data x [n] representing the white noise read from the time series data memory 4. The bicoherence calculation unit 20 outputs the calculated bicoherence. The output bicoherence is represented by a symbol γw [k ₁ , k ₂ ] different from that in the case of processing the observation data.

式（１６）でγｗ_ｑの添え字ｑは、それぞれ異なる周波数ペア（ｋ_１，ｋ_２）に対応するものであり、Ｑ個の周波数ペアを用いる場合、ｑは１からＱまでの値を取る。例えば、図３に示されるのと同じ配置の周波数ペアを用いる場合、Ｑは例えば４９である。 (C) The verification amount calculation unit 32 calculates a verification amount (experimental verification amount) zw based on the bicoherence γw [k ₁ , k ₂ ]. This calculation is expressed by the following equation (16).

In Expression (16), the suffixes _q of γw _q correspond to different frequency pairs (k ₁ , k ₂ ). When Q frequency pairs are used, q takes a value from 1 to Q. . For example, when using frequency pairs with the same arrangement as shown in FIG.

(D) Repeat the processes (a) to (c) a plurality of times, and average the mean and variance of the test amount zw _r (r is any one of 1 to R) obtained a plurality of times, that is, R times. This is supplied to the variance calculation unit 44. The average variance calculator 44, the average of multiple tests amount zw _r obtained in trials Zwa, calculating the variance Zwv. This calculation is expressed by the following equations (17) and (18).

The average variance calculation unit 44 supplies the average value zwa and the variance zwv thus obtained to the parameter calculation unit 46. The parameter calculation unit 46 calculates the parameters α and β based on the average value zwa and the variance zwv. This calculation is expressed by the following equations (19) and (20).

  Here, N and P are the number of data blocks and the number of samples when time series data representing white noise is read and supplied to the detection processing unit 6, respectively.

(E) The calculations from steps (a) to (d) are performed for different values of the number of data blocks P, the number of samples N, and the overlap rate L, and values of α and β are obtained and stored in the parameter memory 36. To do. For example, it is stored in association with P, N, and L, for example, in the form of a lookup table.

  When analysis is performed on the observation data, α and β corresponding to P, N, and L at that time are read out and used for calculation of Expressions (14) and (15).

As described above, for each setting condition (P, N, L), instead of storing α and β in a table format, bias, σ ² is set for each setting condition P, N, L. It may be calculated and stored in a memory. In that case, the adjustment amount calculation unit 38 of FIG. 5 is not necessary.

  As described above, in the present invention, as in the conventional method, data of N samples × P data blocks are extracted from observation data, and P N sample discrete Fourier transforms are calculated. Thereafter, the processing up to the calculation of the power spectrum and the bispectrum using the output of the Fourier transform of each of the P data blocks is the same as the conventional method.

であるのに対し、バイコヒーレンスγ［ｋ_１，ｋ_２］の分母は

であり、両者の違いは、Ｓａ［ｋ_１］Ｓａ［ｋ_２］と、Ｃａ［ｋ_１，ｋ_２］の部分にある。この違いによって、バイコヒーレンスγ［ｋ_１，ｋ_２］は、式の上では信号Ｘ［ｋ_１］Ｘ［ｋ_２］と信号Ｘ［ｋ_１＋ｋ_２］間のコヒーレンスを求めているが、結果的には観測データの周波数ｋ_１及びｋ_２の成分間の位相的な相関関係を評価していることになる。 The difference between the bicoherence γ [k ₁ , k ₂ ] (formula (12)) and the conventional normalized bispectrum Bn [k ₁ , k ₂ ] (formula (6)) is in the denominator in the formula. The denominator of the conventional normalized bispectrum Bn [k ₁ , k ₂ ] is

Whereas the denominator of bicoherence γ [k ₁ , k ₂ ] is

The difference between the two is in the portions of Sa [k ₁ ] Sa [k ₂ ] and Ca [k ₁ , k ₂ ]. Due to this difference, the bicoherence γ [k ₁ , k ₂ ] obtains the coherence between the signal X [k ₁ ] X [k ₂ ] and the signal X [k ₁ + k ₂ ] in the equation. In other words, the phase correlation between the components of the frequencies k ₁ and k ₂ of the observation data is evaluated.

  The bispectrum itself also represents the phase relationship between the frequencies, but the bispectrum also depends on the amplitude of the observation data. On the other hand, the bicoherence used in the present invention is normalized as shown in Equation (12) so that the amplitude is between 0 and 1, and the two frequency components are Only the information about the phase is extracted.

If the observation data is non-Gaussian and non-linear, there is a phase connection between the signals of the frequencies k ₁ and k ₂ , so that the bicoherence has a high value. From the above, it is possible to confirm whether or not the observation data has non-Gaussian and non-linearity by bicoherence.

  In order to confirm the effect of the method of the present invention, an experiment was performed using a transient signal generated in a simulated manner. Transient signals are non-Gaussian and non-linear. 8A and 8B show a simulated transient signal Str and a waveform obtained by adding white noise Nwt to the transient signal Str. FIG. 8A shows the generated transient signal, which appears at times of 5, 10, 15, and 20 seconds elapsed from a certain starting point. FIG. 8B shows a waveform obtained by adding white noise Nwt to the generated transient signal Str, and the conventional method and the method of the present invention were applied to this waveform.

FIGS. 9A and 9B are diagrams showing the results of applying the conventional method and the method of the present invention to the generated signal (the signal shown in FIG. 8B). FIG. 9A shows the result of applying the conventional method, and FIG. 9B shows the result of applying the method of the present invention. Calculation parameters were set as follows.
-Number of samples per data block N = 32
-Number of data blocks P = 32
・ Overlap ratio L = 50%
・ Number of frequency pairs Q = 49
-Adjustment amount bias used for calculation of test amount z _γ = 0.0625
・Adjustment amount σ ² = 0.1875 used to calculate the verification amount z _γ
The sampling frequency is 50 kHz,
Arranged on a two-dimensional plane of the pair of frequencies is as listed in FIG. 3, the frequency f _{m (m} = 1~14) shown in FIG. 3 respectively have the following values.
When m = 1, f _m = 0
When 1 <m <14, f _m = f _(m−1) + ΔF
However, ΔF = 1562.5 Hz.

  Comparing the results of the two processes in FIGS. 9A and 9B, the level when only white noise is observed is the same. On the other hand, looking at the response level when observing a transient signal, the method of the present invention is approximately twice as large as the conventional method. Therefore, it can be confirmed that the method of the present invention has an output SNR higher than that of the conventional method and has a high ability to detect non-Gaussian and non-linearity of observation data.

  Although the present invention has been described above as a data analysis apparatus, a data analysis method implemented by the data analysis apparatus also forms part of the present invention. Furthermore, a program for causing a computer to execute the processing of each unit in the data analysis device or each processing in the data analysis method, and a computer-readable recording medium recording such a program also form part of the present invention. .

  The case where the observation data is an audio signal in the audible range has been described above, but the present invention is not limited to this. For example, the present invention can be applied to analysis of ultrasonic signals, analysis of mechanical vibrations, and analysis of brain waves.

2 input signal processing unit, 4 time series data memory, 6 detection processing unit, 12 data block generation unit, 14 Fourier transform unit, 16 frequency component memory, 18 frequency component selection unit, 20 bicoherence calculation unit, 22 bispectral calculation unit , 24 power spectrum calculation unit, 26 product square average calculation unit, 28 normalization unit, 32 test amount calculation unit, 34 comparison unit, 36 parameter memory, 38 adjustment amount calculation unit, 42 white noise generation unit, 44 average variance calculation unit 46 Parameter calculation unit.

Claims (20)

A Fourier transform unit that Fourier-transforms a plurality of data blocks each composed of successive sample values in time-series data representing observation values, and calculates a plurality of frequency components for each data block;
A bicoherence calculator that calculates the bicoherence of each frequency pair consisting of two frequencies of the frequency components for each data block calculated by the Fourier transform unit;
A test amount calculator for calculating a test amount for the time series data based on the bicoherence for each of the frequency pairs;
A comparison unit that determines whether or not the verification amount calculated by the verification amount calculation unit is greater than a predetermined threshold;
The determination result by the comparison unit is output as the analysis result,
The test amount calculation unit determines whether the average value and variance of the test amount are based on the number of sample values constituting each data block, the number of data blocks constituting the time series data, and the overlap ratio of the data blocks. using the adjustment amount of order to be kept constant without, intends line calculation of the test statistic
Data analysis device. A Fourier transform unit that Fourier-transforms a plurality of data blocks each composed of successive sample values in time-series data representing observation values, and calculates a plurality of frequency components for each data block;
A bicoherence calculator that calculates the bicoherence of each frequency pair consisting of two frequencies of the frequency components for each data block calculated by the Fourier transform unit;
A test amount calculator for calculating a test amount for the time series data based on the bicoherence for each of the frequency pairs;
A comparison unit that determines whether or not the verification amount calculated by the verification amount calculation unit is greater than a predetermined threshold;
The determination result by the comparison unit is output as the analysis result,
The test statistic calculation unit, when representing the test amount z _gamma,

(However,
Q is the number of frequency pairs,
γ _q is the bicoherence calculated by the bicoherence calculator for each of the frequency pairs,
bias is an adjustment amount for setting the average of the test amount to 0,
σ ² is an adjustment amount for correcting the variance of the test amount )
Calculate the test statistic z _gamma by
Data analysis device. The adjustment amounts bias and σ ² are

(However,
P is the number of data blocks constituting the time series data,
N is the number of sample values that make up each data block,
α and β are parameters)
The data analysis apparatus according to claim 2 , further comprising an adjustment amount calculation unit that calculates by the following. Using time series data representing white noise instead of time series data representing the observed value,
Performing the Fourier transform and the bicoherence calculation;
Calculate a test verification amount based on the bicoherence obtained by the calculation,
Calculating the mean zwa and variance zwv of the test calibrate for multiple trials;
Based on the average zwa and variance zwv determined by the calculation,

A parameter generating device for generating the parameters α and β,
The data analysis apparatus according to claim 3 , further comprising: a parameter memory that stores the parameters α and β generated by the parameter generation apparatus. The bicoherence calculator is
A bispectrum calculation unit that calculates an average for each of the plurality of data blocks of each of the frequency pairs for each data block and outputs the bispectrum for the time-series data;
An average of the power spectrum of the frequency component equal to the sum of the frequencies constituting each of the frequency pairs for each data block is calculated for the plurality of data blocks and output as a power spectrum for the time series data A power spectrum calculation unit,
The average of the square of the absolute value of the product of the frequency components constituting each of the frequency pairs for each data block is calculated for the plurality of data blocks and output as the product square average for the time series data Product square average calculator,
The bispectrum for the time-series data output from the bispectrum calculation unit, the power spectrum for the time-series data output from the power spectrum calculation unit, and the product square average calculation unit based on the product mean square of the time-series data, the data analyzing apparatus according to any one of claims 1 to 4, characterized in that it has a normalized unit for calculating the bicoherence. When the bispectrum calculation unit represents the bispectrum for each data block as B _i [k ₁ , k ₂ ],

(However,
k ₁ and k ₂ are frequencies constituting the frequency pair,
X _i [k ₁ ], X _i [k ₂ ], and X _i [k ₁ + k ₂ ] are the frequencies k ₁ , k ₂ , k of the i-th data block (i is any one of 1 to P), respectively. ₁ + k ₂ frequency component)
To calculate the bispectrum B _i [k ₁ , k ₂ ] for each data block,
When the bispectrum for the time series data is represented by Ba [k ₁ , k ₂ ],

(Where P is the number of data blocks constituting the time series data)
The data analysis apparatus according to claim 5 , wherein the bispectrum Ba [k ₁ , k ₂ ] for the time series data is calculated by: When the product square average calculation unit represents the product square average of the time series data as Ca [k ₁ , k ₂ ],

(However,
k ₁ and k ₂ are frequencies constituting the frequency pair,
X _i [k ₁ ] and X _i [k ₂ ] are the frequency components of the frequencies k ₁ and k ₂ of the i-th data block (i is any one of 1 to P),
P is the number of data blocks constituting the time series data)
7. The data analysis apparatus according to claim 5 , wherein the product root mean square Ca [k ₁ , k ₂ ] for the time series data is calculated by: When the power spectrum of the time series data is represented by Sa [k ₁ + k ₂ ],

(However,
k ₁ and k ₂ are frequencies constituting the frequency pair,
X _i [k ₁ + k ₂ ] is a frequency component of the frequency k ₁ + k ₂ of the i-th data block (i is any one of 1 to P),
X _i ^* [k ₁ + k ₂ ] is a complex conjugate of X _i [k ₁ + k ₂ ],
P is the number of data blocks. )
The data analysis apparatus according to claim 5 , wherein the power spectrum Sa [k ₁ + k ₂ ] for the time series data is calculated by: When the normalization unit represents the bicoherence as γ [k ₁ , k ₂ ],

The data analysis apparatus according to claim 5 , wherein the bicoherence γ [k ₁ , k ₂ ] is calculated by: A Fourier transform step of Fourier transforming a plurality of data blocks each composed of successive sample values in time-series data representing observation values, and calculating a plurality of frequency components for each data block;
A bicoherence calculation step of calculating the bicoherence of each frequency pair consisting of two frequencies of the frequency components for each data block calculated in the Fourier transform step;
A test amount calculation step of calculating a test amount for the time-series data based on the bicoherence for each of the frequency pairs;
A comparison step for determining whether or not the test amount calculated in the test amount calculation step is larger than a predetermined threshold,
The determination result by the comparison step is output as the analysis result,
In the test amount calculation step, the mean value and variance of the test amount are determined by the number of sample values constituting each data block, the number of data blocks constituting the time series data, and the overlap ratio of the data blocks. using the adjustment amount of order to be kept constant without, intends line calculation of the test statistic
Data analysis method. A Fourier transform step of Fourier transforming a plurality of data blocks each composed of successive sample values in time-series data representing observation values, and calculating a plurality of frequency components for each data block;
A bicoherence calculation step of calculating the bicoherence of each frequency pair consisting of two frequencies of the frequency components for each data block calculated in the Fourier transform step;
A test amount calculation step of calculating a test amount for the time-series data based on the bicoherence for each of the frequency pairs;
A comparison step for determining whether or not the test amount calculated in the test amount calculation step is larger than a predetermined threshold,
The determination result by the comparison step is output as the analysis result,
In the test amount calculation step, when the test amount is represented by z _γ ,

(However,
Q is the number of frequency pairs,
γ _q is the bicoherence calculated in the bicoherence calculation step for each of the frequency pairs,
bias is an adjustment amount for setting the average of the test amount to 0,
σ ² is an adjustment amount for correcting the variance of the test amount )
Calculate the test statistic z _gamma by
Data analysis method. The adjustment amounts bias and σ ² are

(However,
P is the number of data blocks constituting the time series data,
N is the number of sample values that make up each data block,
α and β are parameters)
The data analysis method according to claim 11 , further comprising an adjustment amount calculation step of calculating by the following. Using time series data representing white noise instead of time series data representing the observed value,
Performing the Fourier transform and the bicoherence calculation;
Calculate a test verification amount based on the bicoherence obtained by the calculation,
Calculating the mean zwa and variance zwv of the test calibrate for multiple trials;
Based on the average zwa and variance zwv determined by the calculation,

A parameter generating step for generating the parameters α and β,
The data analysis method according to claim 12 , further comprising a storage step of storing the parameters α and β generated in the parameter generation step in a parameter memory. The bicoherence calculation step includes:
A bispectrum calculation step of calculating an average for each of the plurality of data blocks of each of the frequency pairs for each data block and outputting the bispectrum for the time-series data; and
An average of the power spectrum of the frequency component equal to the sum of the frequencies constituting each of the frequency pairs for each data block is calculated for the plurality of data blocks and output as a power spectrum for the time series data Power spectrum calculation step to perform,
The average of the square of the absolute value of the product of the frequency components constituting each of the frequency pairs for each data block is calculated for the plurality of data blocks and output as the product square average for the time series data Product square average calculation step;
The bispectrum for the time series data output from the bispectrum calculation step, the power spectrum for the time series data output from the power spectrum calculation step, and the product root mean square calculation step. The data analysis method according to claim 10 , further comprising a normalization step of calculating the bicoherence based on the product root mean square of time series data. The bispectrum calculation step represents the bispectrum for each data block as B _i [k ₁ , k ₂ ],

(However,
k ₁ and k ₂ are frequencies constituting the frequency pair,
X _i [k ₁ ], X _i [k ₂ ], and X _i [k ₁ + k ₂ ] are the frequencies k ₁ , k ₂ , k of the i-th data block (i is any one of 1 to P), respectively. ₁ + k ₂ frequency component)
To calculate the bispectrum B _i [k ₁ , k ₂ ] for each data block,
When the bispectrum for the time series data is represented by Ba [k ₁ , k ₂ ],

(Where P is the number of data blocks constituting the time series data)
The data analysis method according to claim 14 , wherein the bispectrum Ba [k ₁ , k ₂ ] for the time series data is calculated by: When the product root mean square calculation step expresses the product root mean square for the time series data as Ca [k ₁ , k ₂ ],

(However,
k ₁ and k ₂ are frequencies constituting the frequency pair,
X _i [k ₁ ] and X _i [k ₂ ] are the frequency components of the frequencies k ₁ and k ₂ of the i-th data block (i is any one of 1 to P),
P is the number of data blocks constituting the time series data)
The data analysis method according to claim 14 or 15 , wherein the product root mean square Ca [k ₁ , k ₂ ] for the time series data is calculated by: In the power spectrum calculation step, when the power spectrum for the time series data is represented by Sa [k ₁ + k ₂ ],

(However,
k ₁ and k ₂ are frequencies constituting the frequency pair,
X _i [k ₁ + k ₂ ] is a frequency component of the frequency k ₁ + k ₂ of the i-th data block (i is any one of 1 to P),
X _i ^* [k ₁ + k ₂ ] is a complex conjugate of X _i [k ₁ + k ₂ ],
P is the number of data blocks. )
The data analysis method according to claim 14, 15, or 16 , wherein the power spectrum Sa [k ₁ + k ₂ ] for the time series data is calculated by: In the normalization step, when the bicoherence is represented by γ [k ₁ , k ₂ ],

The data analysis method according to claim 14 , wherein the bicoherence γ [k ₁ , k ₂ ] is calculated by: The program for making a computer perform each process of the data analysis method in any one of Claims 10 thru | or 18 . A computer-readable recording medium on which the program according to claim 19 is recorded.

Images