JP6627639B2 - Abnormality diagnosis device and abnormality diagnosis method - Google Patents

Description

  The present invention relates to an abnormality diagnosis apparatus and method for performing an abnormality diagnosis by frequency-analyzing waveform data (time-series data) such as vibration and sound during operation of an abnormality diagnosis target device.

  As a method of diagnosing an abnormality of an abnormality diagnosis target device, for example, a method described in Patent Document 1 has been proposed.

  According to the method disclosed in Patent Document 1, time-series waveform data is divided into short time intervals, and Fourier transform is performed to create a large number of samples of multivariate data of frequency components. An eigenvector (also referred to as a rotation matrix) for obtaining a principal component score by performing principal component analysis on a large number of multivariate samples is prepared.

  At the time of diagnosis, the measurement data is Fourier-transformed in the same manner as described above to create a multivariate data sample of the frequency component, which is converted by a prepared rotation matrix to obtain a principal component score, which is calculated in a normal state. Abnormality diagnosis is performed by comparing with the principal component score of. In short, the main point is that, instead of performing the principal component analysis independently for each diagnosis, a result obtained by performing the principal component analysis in advance is used and the comparison is facilitated by performing the conversion based on the same reference at the time of the diagnosis.

  Note that the component calculation method of constant Q conversion used in the present invention is described in, for example, Non-Patent Documents 1 to 4.

Japanese Patent No. 3382240

Judith C. Brown: Calculation of a constant Q spectral transfer, J. Org. Acoustic. Soc. Am. 89 (1): 425-434, 1991 Judith C. Brown and Miller S. Puckette: An effective algorithm for the calculation of a constant Q transform, J. Mol. Acoustic. Soc. Am. 92 (5): 2698-2701, 1992 yukara — 13: [Python] Constant-Q conversion (logarithmic frequency spectrogram), super introduction to music programming (provisional) in Hatena Blog, 2013-12-01, 2013, Internet <URL: http: // yukara-13. hatenablog. com / entry / 2013/12/01/222742>. [2016-02-22 access] yukara_13: [Python] High-speed Constant-Q conversion (with FFT), super introduction to music programming (provisional) in Hatena Blog, 2014-01-05, 2014, Internet <URL: http: // yukara-13. hatenablog. com / entry / 2014/01/05/0624414>. [2016-02-22 access]

  If a long-term diagnosis is to be performed using a rotation matrix created in advance using multivariate samples based on the Fourier transform, a large amount of data is required for the prior analysis. For example, in diagnosis of a rotating machine, a frequency resolution of about 1 Hz is required at least in a low frequency region.

  For this purpose, each sample needs to be measured for at least about 1 second, and even if the time is duplicated and the number of samples is increased, it is essentially a similar sample and lacks diversity, so diagnosis is based on that. Then, in the long term, it is not possible to distinguish between abnormalities and changes over time.

  The case is described below. In this case, the vibration of a certain rotating machine was measured and analyzed over a long period (about 4 months). Once a hour, a measurement is performed for 5 seconds at a sampling frequency of about 100 kHz. From this, 100 samples of slightly more than a second divided data are created from one measurement while overlapping the time. Approximately 80 multivariate samples of 6 octave bands were obtained. The rotation matrix is calculated by principal component analysis of the data for this week (however, about 8000 samples were used because the data was limited to the time of operation), and the degree of abnormality based on the Hotelling T2 / Q statistic was calculated based on this. Is calculated as shown in FIG.

  Although the target rotating machine had no particular problem, both the Indicator 1 based on the Hotelling T2 statistic and the Indicator 2 based on the Q statistic, as can be seen from the “linear” slope in FIG. It can be seen that it is increasing with the increase. If the rate of this increase is 10 times (70 days) the learning period (7 days), the average value will be abnormal, and cannot be used for long-term diagnosis.

  This problem does not depend on the definition of the degree of anomaly. Even if the distribution of the first principal component × the second principal component in FIG. Obviously, the data cannot cover the entire data range indicated by the black circles.

  In FIG. 7, the horizontal axis PC1 indicates the first main component, and the vertical axis PC2 indicates the second main component.

  This problem is considered to be because the diversity of the data cannot be secured because the Fourier transform is used to create the multivariate data from the time series data.

  The present invention has been made to solve the above problems, and an object of the present invention is to provide an abnormality diagnosis apparatus and method capable of ensuring sample diversity and performing highly reliable abnormality diagnosis over a long period of time. .

An abnormality diagnosis apparatus according to claim 1 for solving the above-mentioned problem, generates a multivariate sample by performing constant-Q conversion on time-series data measured in advance from an abnormality diagnosis target and creates a multivariate sample based on the generated multivariate sample. Diagnostic parameter creation means for creating a diagnostic parameter in
Time-series data measured from freshly abnormality diagnostic object by converting the constant Q to create a multivariate sample, and a diagnostic means for diagnosing an abnormality using the diagnostic parameters created by the diagnostic parameter generation means ,
The diagnostic parameter creating means,
A first constant-Q conversion device that performs constant-Q conversion on time-series data to create a multivariate sample;
A principal component analysis unit that calculates a principal component score of each sample while obtaining a transformation matrix by performing principal component analysis on the created multivariate sample;
A first statistic calculation unit that obtains a statistic serving as an index of a sample based on the principal component score;
A first normalization unit that normalizes the statistic to determine the degree of abnormality,
The diagnostic means includes:
A second constant-Q conversion device that performs constant-Q conversion on the time-series data to create a multivariate sample;
From the created multivariate sample, a principal component calculation unit that obtains a principal component score using a transformation matrix obtained by the principal component analysis unit,
A second statistic calculation unit that obtains a statistic based on the principal component score;
A second normalization unit that normalizes the statistic by a normalization coefficient used in the first normalization unit to calculate an abnormality degree, and based on the abnormality degree, Diagnose,
The first and second constant Q conversion devices of the diagnostic parameter creating means and the diagnostic means are
In the set frequency components K _M or lower frequency band, a first amount of calculation of the constant Q transform components determined by the product-sum calculation between the frequency axis coefficients in the discrete Fourier transform and a constant Q transform of the time series data, the K _In a high-frequency band exceeding _M , a total operation amount including the second operation amount of the constant Q conversion component obtained by the product-sum calculation of the time series data and the time axis coefficient in the constant Q conversion is calculated,
Obtains the total calculation amount for each K _M by each setting the K _M to all frequency components K,
A calculation method boundary determining unit that evaluates a total calculation amount for each of the K _M , and determines a K _M with the minimum calculation amount as a frequency component threshold indicating a calculation method boundary;
The low-frequency-band frequency components K is less than the frequency component threshold value K _M which is the determined in the constant Q transform, when selecting the first calculation method is a product-sum computation of the discrete Fourier transform and the frequency axis coefficient series data and, in the high frequency band in which the frequency components K exceeds the frequency component threshold value K _M which is the determined, the calculation method selection means for selecting the second calculation method is a product-sum computation of the time-series data and the time-axis coefficient ,
Constant Q conversion means for performing a first calculation method and a second calculation method selected by the calculation method selection means on the time-series data to obtain a constant Q conversion component;
It is characterized by having.

The abnormality diagnosis method according to claim 5 , wherein the diagnostic parameter creating means creates a multivariate sample by performing constant Q conversion on the time series data measured in advance from the abnormality diagnosis object, and creates the multivariate sample. Diagnostic parameter creation step of creating a diagnostic parameter based on the
A diagnosing step of diagnosing an abnormality using the diagnostic parameter created by the diagnostic parameter creating means, wherein the diagnostic means creates a multivariate sample by performing constant Q conversion on the time series data newly measured from the abnormality diagnostic object; and, with a,
The diagnostic parameter creating step,
A first constant-Q conversion device performing constant-Q conversion on the time-series data to create a multivariate sample;
A principal component analysis unit for performing a principal component analysis on the created multivariate sample to obtain a transformation matrix and calculating a principal component score of each sample;
A step in which a first statistic calculation unit obtains a statistic serving as an index of a sample based on the principal component score;
A first normalizing unit that normalizes the statistic to obtain an abnormal degree,
The diagnosis step includes:
A second constant-Q conversion device performing constant-Q conversion on the time-series data to create a multivariate sample;
A principal component calculation unit, from the created multivariate sample, using a transformation matrix obtained by the principal component analysis unit to obtain a principal component score;
A second statistic calculation unit for obtaining a statistic based on the principal component score;
A second normalization unit for normalizing the statistic by the normalization coefficient used in the first normalization unit to calculate an abnormality degree, and an abnormality diagnosis target based on the abnormality degree Diagnose abnormalities,
The step of the first constant Q conversion device performing constant Q conversion of the time series data to create a multivariate sample,
In the set frequency components K _M or lower frequency band, a first amount of calculation of the constant Q transform components determined by the product-sum calculation between the frequency axis coefficients in the discrete Fourier transform and a constant Q transform of the time series data, the K _In a high-frequency band exceeding _M , a total operation amount including the second operation amount of the constant Q conversion component obtained by the product-sum calculation of the time series data and the time axis coefficient in the constant Q conversion is calculated,
Obtains the total calculation amount for each K _M by each setting the K _M to all frequency components K,
The evaluated total calculation amount for each K _M, comprising the step of determining the K _M of the amount of calculation is minimized, as a frequency component threshold indicating the calculation method boundary,
The step of the second constant-Q conversion device performing constant-Q conversion on the time-series data to create a multivariate sample,
The low-frequency-band frequency components K is less than the frequency component threshold value K _M which is the determined in the constant Q transform, when selecting the first calculation method is a product-sum computation of the discrete Fourier transform and the frequency axis coefficient series data a step, and the frequency component K is for selecting the second calculation method is a product-sum calculation of said time axis factor, the time series data and high frequency band exceeding the frequency component threshold value K _M which is the determined,
Performing a selected first calculation method and a second calculation method on the time-series data to obtain a component of constant Q conversion .

  According to the above configuration, when creating multivariate data of frequency components from time-series data, a constant Q transform is used instead of a Fourier transform. The diversity of the sample can be ensured, and the abnormality can be diagnosed based on the same standard over a long period of time. This improves the reliability of the abnormality diagnosis.

Further, in the abnormality diagnosis device according to claim 2 , in claim 1 , each of the constant Q conversion devices includes a time axis coefficient in the constant Q conversion,

A time axis coefficient calculation unit calculated by:
This is a discrete Fourier transform when the frequency axis coefficient in the constant Q conversion is regarded as time series data of time axis coefficient T _{n, k} with n = 0,..., N−1.

And a frequency axis coefficient calculation unit that calculates as (1 / N) S _{n, k} from
The constant Q conversion means includes:
A low-frequency band calculation unit that performs the first calculation method by performing a product-sum calculation on the frequency axis coefficient and the discrete Fourier transform of the time-series data calculated by the frequency axis coefficient calculation unit;
A high-frequency band calculation unit that performs a product-sum calculation of the time-axis coefficient and time-series data calculated by the time-axis coefficient calculation unit and performs the second calculation method;
A combining unit that combines the outputs of the low-frequency band calculation unit and the high-frequency band calculation unit to obtain a component of constant Q conversion;
It is characterized by having.

  According to the above configuration, the constant Q conversion component can be calculated by a calculation method with a small amount of calculation in both the low frequency band and the high frequency band, and the calculation is speeded up.

Further, the abnormality diagnosis device according to claim 3, in claim 1, wherein the constant Q conversion device,

A plurality of analysis position of the preset time direction _{p m (m = 1, ...} , M) time axis coefficient calculated for each T _{n, k,} and _m determined, the constant said T _{n, k,} and _m A time axis coefficient calculating unit that sets a time axis coefficient in the Q conversion;
This is a discrete Fourier transform when the time axis coefficient T _{n, k} is regarded as time series data of n = 0,..., N−1.

Was obtained from (1 / N) S _n, the _k, determined frequency axis coefficient (1 / N) S n, _k, and _m are calculated for each said analysis position p _m, the (1 / N) S _n, a frequency axis coefficient calculating unit that sets _{k and m} as frequency axis coefficients in the constant Q conversion,
The constant Q conversion means includes:
When asked to Fourier coefficient X _n by performing a discrete Fourier transform for the entire series data x _n, for each of the analyzed position p _m, the frequency axis coefficient calculated by the frequency axis coefficient calculator (1 / N) S _{n , k, m} and the Fourier coefficient X _n to calculate the sum of products and perform the first calculation method, a low frequency band calculation unit,
For each of the analysis positions p _m , a high-frequency wave for performing the second calculation method by performing a product-sum calculation on the time axis coefficient T _{n, k, m} calculated by the time axis coefficient calculation unit and the time series data x _n A bandwidth calculator,
A combining unit that combines the outputs of the low-frequency band calculation unit and the high-frequency band calculation unit to obtain a component of constant Q conversion;
It is characterized by having.

According to the above configuration, when performing the constant Q conversion at each of a plurality of times, the calculation on the frequency axis in the low-frequency band is performed without performing the discrete Fourier transform of the time-series data at each of the plurality of times. _Since the product-sum calculation of the Fourier coefficients _{Xn obtained} by performing the discrete Fourier transform on the entire _n and the frequency axis coefficients at each analysis position is sufficient, the calculation as a whole is speeded up.

Further, in the abnormality diagnosis device according to the fourth aspect , in the first aspect , each of the constant Q conversion devices includes:

A plurality of analysis position of the preset time direction _{p m (m = 1, ...} , M) time axis coefficient T _n is calculated in one analysis position p ₁ of _the, seeking _{k, 1,} the T _{n , k, 1} as a time axis coefficient in the constant Q conversion,
This is a discrete Fourier transform when the time axis coefficient T _{n, k} is regarded as time series data of n = 0,..., N−1.

Was obtained from (1 / N) S _n, the _k, the one frequency axis coefficient calculated in the analysis position p ₁ of (1 / N) S n, obtains a _{k, 1,} the (1 / N) S _{n , k, 1} as frequency axis coefficients in the constant Q conversion,
The constant Q conversion means includes:
A discrete Fourier transform is performed on the entire time series data x _n to obtain a Fourier coefficient X _n, and for each of the analysis positions p _m ,
A frequency axis coefficient (1 / N) S _{n, k, 1} obtained by the frequency axis coefficient calculation unit;
There the analysis frequency axis coefficient for each position _{p m (1 / N) S} n, k, m are in a relationship permitting consequence the frequency axis coefficients in one analysis position _{p 1 (1 / N) S} n, the _{k, 1} Power of complex exponential function derived from

When,
A low-frequency band calculation unit that performs a product-sum calculation of the Fourier coefficient X _n and performs the first calculation method;
By shifting the range for each of the analyzed position p _m, the time-series data x _{n-p1 +} and _pm by product-sum computation the time axis coefficient T _n calculated by the time-axis coefficient calculating _{section, k, 1} and the A high-frequency band calculation unit that performs the calculation method 2;
A combining unit that combines the outputs of the low-frequency band calculation unit and the high-frequency band calculation unit to obtain a component of constant Q conversion;
It is characterized by having.

According to the above configuration, when performing each constant Q transform a plurality of times, so by utilizing the time-axis coefficient and frequency axis coefficients in one analysis position p ₁ is calculated at each time is possible, a plurality of It is not necessary to prepare a time-axis coefficient and a frequency-axis coefficient for each time (for each of a plurality of analysis positions), and the amount of data can be reduced to save memory.

(1) According to the first to fifth aspects of the present invention, when generating multivariate data of frequency components from time-series data, a constant Q-transform instead of a Fourier transform is used, so that diagnostic parameter generating means It is possible to ensure the diversity of the sample even if the data used in the method is relatively small, and the abnormality can be diagnosed based on the same standard over a long period of time. This improves the reliability of the abnormality diagnosis.
(2) According to the first and second aspects of the present invention, the constant Q conversion component can be calculated by a calculation method with a small amount of calculation in both the low frequency band and the high frequency band, and the calculation is speeded up. Is realized.
(3) According to the third aspect of the invention, when performing the constant Q conversion at each of a plurality of times, the calculation on the frequency axis in the low frequency band is performed by performing the discrete Fourier transform of the time series data by the plurality of times. even without every, for good in product-sum calculation between the frequency axis coefficient on the Fourier coefficients X _n and each analysis position discrete Fourier transform collectively on the entire x _n, is calculated as a whole faster You.
(4) according According to the invention described in claim 4, when performing each constant Q transform a plurality of times, the computation of a single analysis position p each time by using the time axis coefficient and frequency axis coefficients in ₁ Therefore, it is not necessary to prepare a time-axis coefficient and a frequency-axis coefficient for each of a plurality of times (each of a plurality of analysis positions), and the amount of data can be reduced to save memory.

FIG. 1 is a block diagram showing a configuration of an abnormality diagnosis device according to an embodiment of the present invention. FIG. 4 is a scatter diagram of first principal component × second principal component of all period data according to the embodiment of the present invention. Explanatory drawing which shows the period change of the abnormal degree by embodiment of this invention. Explanatory drawing which shows the abnormality diagnosis example by embodiment of this invention. FIG. 1 is a configuration diagram illustrating an example of a constant Q conversion device according to an embodiment of the present invention. Explanatory drawing which shows the period change of the abnormal degree based on the conventional method. FIG. 4 is a scatter diagram of a first principal component × a second principal component of all period data by a conventional Fourier transform.

  Hereinafter, embodiments of the present invention will be described with reference to the drawings, but the present invention is not limited to the following embodiments. In the present invention, when generating multivariate samples of frequency components from time-series data (waveform data), a constant Q transform is used instead of the conventional Fourier transform.

  Constant Q transform (Constant Q transform) does not analyze data of the same period in all frequency bands unlike Fourier transform, but refers to data referred to for each frequency so that the number of periods is the same in all frequency bands. Analyze with different numbers. In this method, since analysis is performed using data of a shorter period in a higher frequency band, a longer period (1 second or longer) is required to increase the frequency resolution in a low frequency band (for example, about 1 Hz). Even if a large number of samples are overlapped, the periods do not overlap in the high frequency band, so that diversity is secured as a multivariate sample.

  FIG. 1 shows the configuration of a specific diagnosis method. In FIG. 1, the method of the first embodiment includes a preparation phase for calculating a parameter for diagnosis from time series data collected by measuring in advance, and a diagnosis phase for diagnosing newly measured time series data.

  In the preparation phase, a time-series data measured and collected from an abnormality diagnosis target is subjected to constant Q conversion to create a multivariate sample, and a diagnostic parameter is created based on the created multivariate sample. It is configured as means 100.

  The diagnostic phase is a diagnostic means for generating a multivariate sample by performing constant-Q conversion on time series data newly measured from an abnormality diagnostic object and diagnosing an abnormality using the diagnostic parameters created by the diagnostic parameter creating means. 200.

  The diagnostic parameter creating means 100 includes a first constant Q conversion device 101 for performing constant Q conversion on time series data collected by measuring from an abnormality diagnosis target in advance to create a multivariate sample, and the created multivariate sample And a principal component analysis unit 102 that obtains a transformation matrix by calculating the principal component and calculates a principal component score of each sample, and a first statistical value calculation unit that obtains a statistic serving as an index of the sample based on the principal component score. 103 and a first normalization unit 104 that normalizes the statistic to determine the degree of abnormality.

  The diagnostic means 200 comprises: a second constant-Q conversion device 201 for performing constant-Q conversion on newly measured time-series data to generate a multivariate sample; and the principal component analysis section from the generated multivariate sample. A principal component calculation unit 202 that obtains a principal component score using the transformation matrix obtained in step 102; a second statistical value calculation unit 203 that obtains a statistic based on the principal component score; A second normalization unit 204 that normalizes with the normalization coefficient used in the first normalization unit 104 to calculate the degree of abnormality, and diagnoses the abnormality to be diagnosed based on the degree of abnormality. .

  The abnormality diagnosis apparatus in FIG. 1 is configured by, for example, a computer, and includes hardware resources of a normal computer, such as a ROM, a RAM, a CPU, an input device, an output device, a communication interface, a hard disk, a recording medium, and a driving device thereof. .

  As a result of the cooperation between the hardware resources and the software resources (OS, application, etc.), the abnormality diagnosis device includes a first constant Q conversion device 101, a principal component analysis unit 102, a first statistical value as shown in FIG. A calculation unit 103, a first normalization unit 104, a second constant Q conversion device 201, a principal component calculation unit 202, a second statistic value calculation unit 203, and a second normalization unit 204 are mounted.

  The data such as the multivariate sample, the principal component score, the statistic, the degree of abnormality, the transformation matrix, and the normalization coefficient are stored in storage means such as a hard disk or RAM.

  Next, the operation of each unit of the apparatus of FIG. 1 will be described. In the preparation phase (diagnosis parameter creation means 100), the first constant Q conversion device 101 sequentially converts the time series data collected in advance into a large number of multivariate samples. The constant Q conversion device 101 may execute a method disclosed in Non-Patent Documents 1 to 4, for example, or may adopt a constant Q conversion device capable of performing high-speed calculations in Examples 2 to 4 described later. good. In the case of obtaining the data at high speed by the method described in the second to fourth embodiments, it is desirable that the data lengths of the time-series data are uniform.

  For example, there are 100 data for 5 seconds at 102.4 kHz sampling. Due to the small number of reference data in the high frequency region of the constant Q conversion, it is possible to obtain various samples of about 100 to 400 having different phases from the rotation cycle and the power supply cycle from one measurement.

  Next, the principal component analysis unit 102 performs a principal component analysis on the multivariate sample to obtain a transformation matrix, and calculates a principal component score of each sample. Then, the statistic calculation unit 103 obtains a statistic serving as an index of the sample such as the Hotelling T2 / Q statistic based on the principal component score.

  Finally, the first normalization unit 104 normalizes the statistic so that the calculated statistic has an average of 0 in the time-series data in the preparation phase, and determines the degree of abnormality. At this time, the normalization coefficient used for normalization is stored and used in the diagnosis phase (diagnosis means 200). Note that, in order to limit the numerical range in which the degree of abnormality can be taken, it is desirable that the degree of abnormality be a logarithmic value of a statistic.

  In the diagnostic phase (diagnosing means 200), the second constant-Q converter 201 performs constant-Q conversion on the newly measured time-series data to create a multivariate sample.

  The constant-Q converter 201 may execute, for example, the methods disclosed in Non-Patent Documents 1 to 4, similarly to the constant-Q converter 101, and can perform high-speed calculations in Examples 2 to 4 described later. A constant Q conversion device may be employed.

  Next, a principal component score is obtained by the principal component calculation unit 202 using the transformation matrix created by the principal component analysis unit 102 in the preparation phase for the multivariate sample. Then, based on this, the second statistic calculation unit 203 obtains a statistic by the same calculation as in the preparation phase. Finally, the second normalization unit 204 calculates the degree of abnormality using the normalization coefficient created by the first normalization unit 104 in the preparation phase, and diagnoses whether the device is normal or abnormal based on the value.

  It is generally difficult to determine a threshold value of an abnormal value used for abnormality diagnosis. However, when an abnormal value is calculated as a logarithm of a statistical value, the threshold value may be log3 or log5 as a simple diagnosis. . This is the same as when the abnormality determination criterion based on the vibration effective value is three times, five times, or the like as compared with the normal case.

  In order to show the effect of the present invention, the same data as the first principal component × second principal component scatter diagram in FIG. FIG. 2 shows the dispersion diagram. According to FIG. 2, unlike the case of FIG. 7, it can be seen that the data in the learning period indicated by the white circles can substantially cover the distribution range of all data indicated by the black circles.

  FIG. 3 shows a result of calculating the degree of abnormality similar to that shown in FIG. 6 by the proposed method. According to FIG. 3, unlike the case of FIG. 6, it can be seen that the degree of abnormality does not tend to increase over time. The slope based on the linear regression is almost zero, and it takes about 10 years for Indicator 1 having a positive slope to reach the threshold value of the average of the abnormalities.

  Next, an example of abnormality diagnosis by the proposed method is shown in FIG. FIG. 4 is a calculation example of a rotating machine different from FIGS. 6 and 3. In this case as well, measurement was performed once an hour at a sampling frequency of about 100 kHz for 5 seconds, but the emergency stop due to a failure occurred on the 25th day from the start of the measurement. On the 25th day, the degree of abnormality has become extremely large immediately after the operation, but if you look closely, you can see that Indicator 1 has crossed the caution line on the 13th day.

  It should be noted that, as a result of this diagnosis, almost the same result can be obtained even by a method using Fourier transform.

  In the second embodiment, the first and second constant-Q converters 101 and 201 of FIG. 1 can perform high-precision constant-Q conversion at a high speed with a small amount of calculation. It consists of.

  First, the components of constant Q conversion will be described.

  According to Non-Patent Document 1, the constant Q conversion is formulated as follows (the notation is slightly changed).

Values other than 0 are taken for each k only in the range of n = (N−N _k ) / 2,..., (N + N _k ) / 2-1. Hamming window, Hanning window, such as a rectangular window, N _k are many types have been proposed: the window _{width, N k = Q (f s} / f k) and so as to take _{N ≧ N 1 ≧ ... ≧ N} k ... ≧ N _K ≧ 2Q Q: The number of periods of the window, which is common to all frequencies in constant Q conversion In Non-Patent Documents 1 and 2, it is assumed that Q = (2 ^{1 / B} −1) ⁻¹ B is the number of bins and is a number that divides one octave. F _s : sampling frequency f _k : frequency of the k-th component, f _k = (2 ^{1 / B} ) ^k−1 f _min
f _min = f ₁ is the smallest frequency to be analyzed at a constant Q transform, the f _k determined from the scope of analyzing the maximum frequency to be analyzed, so that the f _k ≦ f _s / 2.

The constant Q conversion requires a large amount of calculation because the window width _Nk becomes extremely large particularly in a low frequency band component.

  On the other hand, Non-Patent Literature 2 proposes a method for increasing the speed of an arithmetic operation using a fast Fourier transform (FFT; Fast Fourier transform).

  In this method, according to Parseval's formula, the constant Q conversion is replaced from the calculation on the time axis to the calculation on the frequency axis as in the following equation (2).

In the following text of the specification, the component of the constant Q conversion on each of the left sides of the equations (1) and (2) may be described as “X _k ^cq ”.

Here, most of the elements of S _{n, k} (especially for small k) are close to 0, so if terms less than a certain value are expressed as sparse matrices as 0 (since only the product-sum calculation of substantially non-zero terms) You can calculate very fast.

Further, (1 / N) S _{n, k} can be calculated in advance regardless of the input data sequence x _n if the data length to be calculated and the frequency range to be analyzed are determined. Assuming that the discrete Fourier transform of the data is performed by FFT, the constant Q transform can be calculated by one FFT on the input data and a product operation of the FFT result vector and a sparse matrix prepared in advance.

  Equivalent formulas hold for non-discrete Fourier transforms, which are referred to in various ways as Parseval's theorem and Planschler's theorem.

The idea of speeding up according to Non-Patent Document 2 relies on the fact that the coefficient S _{n, k on} the frequency axis can be regarded as a sparse matrix, but in practice, the fact that S _{n, k} is a sparse matrix is a rough evaluation. Certainly, the value is small except for a part, but it is not necessarily sparse when the accuracy is increased. This tendency is particularly remarkable in a high frequency band. In fact, in Non-Patent Documents 2 and 4, the thresholds at which the coefficient S _{n, k} is considered to be sufficiently small are 0.15 and 0.0054, respectively. Is not enough. For this reason, depending on the required accuracy, there is a possibility that the speed-up method may be rather slowed down.

  In order to find out how much calculation amount is actually required in the calculation on the frequency axis, non-zero values to be calculated in the calculation on the time axis and the calculation on the frequency axis are set according to Non-Patent Document 4. Evaluate the number of terms.

Setting parameters, denoted by the symbol method of the present specification and B = 24, Q = 28, f s = 16000, f min = 60, K = 160, respectively rectangular window constant Q transform, Hamming window, Hanning window Table 1 shows the number of non-zero terms required for calculation on the time axis and the number of terms required for calculation up to a coefficient of predetermined accuracy in calculation on the frequency axis.

It should be noted that when combined with Non-Patent Documents 1 and 2, Q = 34 when B = 24. However, in Non-Patent Documents 3 and 4, the value of Q is changed by introducing an additional coefficient “fraction”. Q = 28 is used. K = 160 was derived from the maximum frequency 6000 Hz to be analyzed. Since the value of f _s was not specified in Non-Patent Documents 3 and 4, f _s = 16000 Hz was selected as WAV audio data having a sampling frequency of 12000 Hz or more that can be analyzed up to a maximum frequency of 6000 Hz.

  The numerical value on the right side of the frequency axis in Table 1 is a threshold value of the magnitude of the coefficient that evaluates the coefficient as zero, and the percentage on the right side of each numerical value is based on the number of terms in the time axis calculation of the rectangular window.

From Table 1, the degree depends on the window shape, but the higher the accuracy (the smaller the threshold value of the coefficient S _{n, k} for evaluating the coefficient to be zero), the more the number of terms required for the calculation on the frequency axis is increased. I understand. In particular, in a rectangular window or a Hamming window, when the accuracy is high, the number of terms to be summed up in the calculation on the frequency axis becomes larger than that in the calculation on the time axis, and the calculation on the frequency axis does not become faster.

  Note that the calculation on the frequency axis has an overhead of FFT calculation and non-zero coefficient extraction of a sparse matrix. Therefore, if there is no difference in the number of required terms, the time axis calculation is advantageous.

  To perform high-precision constant Q conversion at high speed, a device is required.

Therefore, in the second embodiment, in order to calculate the constant Q conversion with high accuracy and a small amount of calculation, the constant X conversion component X _k ^cq is calculated using either the time axis calculation or the frequency axis calculation of the equation (2). Or by selecting in advance the smaller calculation amount for each frequency component k.

The calculation of the component X _k ^cq constant Q transform small calculation amount of calculations in the frequency axis at approximately the low frequency band, because in the high frequency band is small calculation amount calculation in the time axis, determines the boundary value K _M for k Te, calculated at k ≦ K _M in the frequency axis, calculations in k> K in _M time axis, may be so called.

  In FIG. 5, reference numeral 10 denotes a preparation unit that determines parameters for performing constant Q conversion in advance to prepare for constant Q conversion calculation, and 20 denotes an individual calculation unit that calculates constant Q conversion for individual time-series data. is there.

The preparation unit 10
In the set frequency components K _M or lower frequency band, a first amount of calculation of the constant Q transform components determined by the product-sum calculation between the frequency axis coefficients in the discrete Fourier transform and a constant Q transform of the time series data, the K in a high frequency band exceeding the _M, when calculating the total amount of calculation and a second calculation of the constant Q transform components determined by the product-sum calculation of the time axis coefficient in the sequence data and the constant Q transform, the K _M respectively set to all frequency components k seek the total calculation amount for each K _M, evaluates the total computation amount of the respective K _M, the K _M of the calculation amount becomes minimum, indicating the calculation method boundary A function of a calculation method boundary determining means for determining as a frequency component threshold,
The low-frequency-band frequency component k is less frequency component threshold value K _M which is the determined in the constant Q transform, when selecting the first calculation method is a product-sum computation of the discrete Fourier transform and the frequency axis coefficient series data and, wherein in a high frequency band in which the frequency component k exceeds the frequency component threshold value K _M which is the determination of the calculation method selection means for selecting the second calculation method is a product-sum computation of the time-series data and the time-axis coefficient Function.

The preparation unit 10 that executes the functions of the calculation method boundary determination unit and the calculation method selection unit specifically includes a parameter determination unit 11 that determines parameters for constant Q conversion _, and calculates a time axis coefficient T _{n, k} . time axis coefficient calculating section 12, the frequency axis coefficient (1 / n) S n, the frequency axis coefficient calculation unit 13 for calculating the _k, determining a k value K _M of the boundary between the calculation of the calculation and the frequency axis in the time axis And a calculation method boundary determining unit 14.

  The individual calculation unit 20 performs a first calculation method or a second calculation method selected by the calculation method selection unit on the time-series data to obtain a constant Q conversion component. Is provided.

  Specifically, the individual calculation unit 20 that performs the function of the constant Q conversion unit includes a preprocessing unit 21 that adjusts the length of the time series data, a high frequency band calculation unit 22 that calculates a constant Q conversion of a high frequency band, It comprises a low-frequency band calculation unit 23 for calculating constant-Q conversion of a frequency band, and a constant-Q conversion synthesis unit 24 for compiling results in a high-frequency band and results in a low-frequency band.

  The constant Q conversion component operation device of FIG. 5 is constituted by, for example, a computer, and hardware resources of a normal computer, for example, ROM, RAM, CPU, input device, output device, communication interface, hard disk, recording medium, and its driving device It has.

  As a result of the cooperation between the hardware resources and the software resources (OS, application, etc.), the constant-Q conversion component operation device includes a parameter determination unit 11, a time axis coefficient calculation unit 12, a frequency axis A coefficient calculation unit 13, a calculation method boundary determination unit 14, a preprocessing unit 21, a high frequency band calculation unit 22, a low frequency band calculation unit 23, and a constant Q conversion synthesis unit 24 are mounted.

  Further, the processing results of the respective units are stored in a storage unit or storage unit such as a hard disk or a RAM, and are used as needed.

First, the parameter determination unit 11 of the preparation unit 10 preliminarily determines the sampling frequency f _s of time-series data to be analyzed by constant Q conversion, the number of bins B that separates the frequency range to be analyzed (f _min , f _max ) and the frequency of one octave, The number of periods Q of the data to be analyzed and the window function are determined. At this time, f _max ≤ f _s / 2 is required, and K = floor (B * log ₂ (f _max / f _min )) + 1. The length of the input data to be calculated based on this and _{N = N 1 = Q (f} s / f min).

  The threshold value Th of the magnitude of the coefficient for which the frequency axis coefficient should be regarded as zero is also determined here.

The time axis coefficient calculation unit 12 calculates the time axis coefficient T _{n, k} based on the parameters determined by the parameter determination unit 11 according to the following definitional expression.

The frequency axis coefficient calculator 13 calculates a frequency axis coefficient (1 / N) S _{n, k} based on the parameters determined by the parameter determiner 11.

Calculate according to the definition formula. At this time, the term of | S _{n, k} | ≦ Th is set to the frequency axis coefficient (1 / N) S _{n, k} = 0.

The calculation method boundary determination unit 14 obtains a maximum k value K _M (frequency component threshold) for performing calculation on the frequency axis. K _M is 0, ..., have one of the values of K, K _M = 0 is means if determined by calculation in the same manner all the time axis and Non-Patent Documents 1 and 2 calculated without the frequency axis I do.

As a method of determining frequency components threshold K _M is, K _M = 0, ..., for each case of K, by evaluating the constant Q transform component X _k calculation amount necessary for obtaining the ^cq C _X (K _M) Determine K _M so that it is minimized.

The operation amount C _X (K _M ) (total operation amount) of the constant X transform component X _k ^cq is the operation amount C _FFT (N) of the FFT of length N, and the operation amount C _{T for} each non-zero time axis coefficient. , The number of non-zero time axis coefficients N _T (k), the amount of computation C _{F for} each non-zero frequency axis coefficient, and the number N _F (k) of frequency axis coefficients regarded as non-zero by k

  Can be expressed as

Here, the C _FFT can vary depending on the computer that performs the calculation and the details of each algorithm, but as an example from actual measurement, C _FFT (N) = (１ ／) Nlog ₂ (N), C _T = 1, C _F = 1 may be set.

In Equation (3) representing the operation amount C _X (K _M), the first term is to K _M is a C _FFT when nonzero, indicates that a 0 and if not, the second term represents the calculation amount of each frequency axis coefficient frequency component k is in the following low-frequency band frequency component threshold value K _M, the third term on the right side, every time axis coefficient in a high frequency band frequency component k exceeds the frequency component threshold value K _M Represents the amount of calculation of.

For this reason, when _KM is non-zero, a first operation amount including the operation amount of the FFT of the first term on the right side and the operation amount of the second term on the right side (operation amount based on the frequency axis coefficient) of Expression (3); The sum of the calculation amount C _X (K _M ) and the second calculation amount based on the time axis coefficient of the third term on the right side of Expression (3) becomes the calculation amount C _X (K _M ).

When _KM is zero, the first term on the right side and the second term on the right side of equation (3) are both 0, and only the second operation amount based on the time axis coefficient of the third term on the right side is the operation amount C _X ( K _M ).

The amount of computation C _X (K _M ) is such that the number of non-zero time axis coefficients N _T (k) decreases exponentially with k, and the number of non-zero frequency axis coefficients N _F (k) is approximately with k. Since it increases exponentially, a characteristic curve which is approximately convex downward in the range of K _M = 1,..., K. However, when K _M = 0, C _FFT disappears and C _X (0) becomes smaller. As a result, the operation amount C _X (K _M ) takes the minimum value at any one of the minimum points at K _M = 0 or K _M = 1,..., K.

  The individual calculation unit 20 calculates the constant Q conversion of the input time-series data as follows.

The pre-processing unit 21 complements the leading and trailing zeros if the length of the input time-series data is shorter than N, and removes the preceding and following values if the length is longer than N to form the time-series data x _n of length N. To be. At this time, the central part of the input time-series data is made to coincide with the central part of the time-series data _xn .

Then component X _k ^cq constant Q transform, by the calculation method boundary determining unit 14 frequency components threshold K _M determined in, for k ≦ K _M was calculated in the low frequency bandwidth calculation unit 23, for k> K _M Is calculated by the high frequency band calculator 22.

The high-frequency band calculation unit 22 calculates the product-sum of the time-axis coefficient T _{n, k} calculated by the time-axis coefficient calculation unit 12 and the time-series data x _n to obtain a constant Q-transform component X _k ^cq . Since there are at most N _k non-zero time-axis coefficients T _n, k for each k, the product-sum calculation can be performed efficiently by limiting it to a non-zero range.

The low-frequency band calculation unit 23 first obtains a Fourier coefficient X _n by performing FFT on the time-series data x _n . Next, the product of the frequency axis coefficient (1 / N) S _{n, k} calculated by the frequency axis coefficient calculation unit 13 and the Fourier coefficient X _n is calculated to obtain a constant Q conversion component X _k ^cq . Since the number of non-zero frequency axis coefficients (1 / N) S _{n, k for each k} is small, efficient calculation can be performed by limiting the product-sum calculation to the non-zero range.

  In many cases, the non-zero term and the zero term are mixed at both ends (near the frequency zero and the vicinity of the sampling frequency), but the zero term continues in the middle part. , Both ends may be calculated as non-zero terms. In this case, the number of non-zero terms slightly increases, but the number is not so large, and since the non-zero determination logic can be simplified, the calculation time is not much different.

Finally, the constant Q conversion synthesis unit 24 puts together the constant X conversion components X _k ^cq calculated by the high frequency band calculation unit 22 and the low frequency band calculation unit 23, and completes the constant Q conversion.

  As described above, according to the second embodiment, when trying to calculate a high-precision constant Q conversion, which is a problem of the conventional high-speed technology, the amount of calculation expands and the calculation is slower than the calculation on the normal time axis. Problem is improved.

In order to show the effect, the required amount of calculation was calculated as shown in Table 2 for each calculation method of constant Q conversion (time axis calculation, frequency axis calculation, proposed method). In the trial calculation, the parameters of the time-series data were set to be the same as those in Table 1, and the calculation parameters were calculated as C _FFT (N) = (1/4) Nlog ₂ (N), C _T = 1, and C _F = 1. Table 2 shows the calculation amount C _X (K _M ) of the component X _k ^cq of the Q conversion.

  The numerical value on the right side of the frequency axis in Table 2 is a threshold value of the magnitude of the coefficient that evaluates the coefficient as zero, and the percentage on the right side of each numerical value is based on the amount of calculation of the time axis calculation of the rectangular window.

  In Table 2, although the effect of speeding-up is small due to the calculation amount of the FFT calculation, the calculation amount of the Hanning window is approximately 1/5 (22.5%) even with the threshold value of 1.0e-5. I have. Further, even when the calculation amount of the FFT calculation is heavy and the calculation on the frequency axis does not produce an effect, the calculation is automatically performed on the time axis, so that there is no problem that the calculation is slowed down.

  The second embodiment is efficient when obtaining each component value of the constant Q conversion at a certain time, but in the actual constant Q conversion, the time is gradually increased in order to capture the time change of the frequency component of the time series data. The constant Q component is repeatedly calculated while shifting. The time to be shifted at this time is the sampling frequency at the shortest, and the window width at the longest is the longest at the maximum in order to calculate the information of the time-series data without leaving it. Of course, when extracting feature samples of time-series data, the time may be shifted at longer intervals to perform constant Q conversion, but in any case, many repetitive calculations are performed on the entire time-series data. No, the load is still high.

  When the constant Q conversion is calculated at many times at the same time, further efficiency is required.

  Therefore, in the third embodiment, when performing the constant Q conversion by the calculation on the frequency axis of Expression (2), the DFT (or FFT) is performed within the range of the time-series data used for the calculation at each time when the constant Q conversion is calculated. Instead, the DFT (or FFT) is performed collectively over the entire range necessary for constant Q conversion at a plurality of times to increase the speed.

The constant Q transformation defined by the equation (1) is formulated to calculate the component X _k ^cq at the center x _{N / 2} of the data sequence x _n , n = 0 ^,. Can be adjusted to calculate components at other times by translating the non-zero range in the time direction. That is _, it can be realized by shifting the index of the product sum of T _{n, k} and x _n .

With respect to the calculation on the frequency axis according to the equation (2), first, DFT (or FFT) is performed on all necessary data over the entire time for constant Q conversion, and is prepared in advance for each time (1/1). N) This is realized by sequentially summing with S _{n, k} .

  The specific calculation configuration of the constant-Q conversion component operation device according to the third embodiment is the same as that of FIG. 5, but the functions executed by the respective units are different as described below.

First, the parameter determination unit 11 of the preparation unit 10 determines in advance the sampling frequency fs of the time-series data to be analyzed by constant Q conversion, the number of bins B separating the frequency range to be analyzed (f _min , f _max ) and the frequency of one octave, The number of periods Q of the target data and the window function are determined. At this time, f _max ≤ f _s / 2 is required, and K = floor (B * log ₂ (f _max / f _min )) + 1.

Further, the total length N of the input data _xn and the analysis positions p _m , m = 1,..., M (analysis position in the time direction; time) are determined in advance. At this time, the constant Q transformation at each position p _m is computed based on the data sequence of length N ₁ = Q centered at _{x pm (f s / f min} ), (N 1/2) ≦ p it is assumed to be _{1 ≦ ... ≦ p M ≦ N-} (N 1/2). That is, all of the analysis position _{p m, m = 1, ...} , data necessary for the calculation of the constant Q transform by M is x _0, ..., are intended to be included in x _N-1. When it is desired to calculate the constant Q conversion at all points of the input data, the input data may be extended to zero before and after and converted so as to satisfy the above condition.

  Further, the threshold value Th of the magnitude of the coefficient for which the frequency axis coefficient should be regarded as zero is also determined here.

Parameter based on the determined time axis coefficient calculating section 12, the parameter determination unit 11, the time axis coefficient T _n, the _k is calculated according to the following defining equation for each analysis position p _m.

At this time, if T _{n, k} for p _m is written as T _{n, k, m} , T _{n, k, m} = T _{n−pm + p1, k, 1} for the non-zero terms.

Based on the parameters determined in the frequency axis coefficient calculation unit 13, the parameter determination unit 11 calculates the frequency axis coefficient (1 / N) S n, k, a _m according defining equation for each analysis position p _m. At this time, the term satisfying | S _{n, k, m} | ≦ Th is set to the frequency axis coefficient (1 / N) S _{n, k, m} = 0.

The calculation method boundary determination unit 14 obtains a maximum k value K _M (frequency component threshold) for performing calculation on the frequency axis. K _M is 0, ..., have one of the values of K, K _M = 0 is means if determined by calculation in the same manner all the time axis and Non-Patent Documents 1 and 2 calculated without the frequency axis I do.

As a method of determining frequency components threshold K _M is, K _M = 0, ..., for each case of K, by evaluating the constant Q transform component X _k calculation amount necessary for obtaining the ^cq C _X (K _M) Determine K _M so that it is minimized. Note that this condition does not depend on a plurality of p _m, it is sufficient if seeking for p ₁ (1 single analysis position).

The operation amount C _X (K _M ) (total operation amount) of the constant X transform component X _k ^cq is the operation amount C _FFT (N) of the FFT of length N, and the operation amount C _{T for} each non-zero time axis coefficient. , The number of non-zero time axis coefficients N _T (k), the amount of computation C _{F for} each non-zero frequency axis coefficient, and the number of frequency axis coefficients N _F (k) regarded as non-zero by k Equation (3) described in the second embodiment

  Can be expressed as

Here, the C _FFT can vary depending on the computer that performs the calculation and the details of each algorithm, but as an example from actual measurement, C _FFT (N) = (１ ／) Nlog ₂ (N), C _T = 1, C _F = 1 may be set.

The amount of operation C _X (K _M ) is based on the number _{NT of} non-zero time axis coefficients where the relationship between the number of terms in the Hanning window and the frequency (frequency component k; logarithmic frequency) of the time axis coefficient and the frequency axis coefficient is non-zero. Since (k) decreases exponentially with respect to k and the number of non-zero frequency axis coefficients N _F (k) increases approximately exponentially with respect to k, it is approximately below in the range of K _M = 1,. It becomes a convex characteristic curve. However, when K _M = 0, C _FFT disappears and C _X (0) becomes smaller. As a result, the operation amount C _X (K _M ) takes the minimum value at any one of the minimum points at K _M = 0 or K _M = 1,..., K.

  The individual calculating unit 20 obtains the input time-series constant Q conversion as follows.

The pre-processing unit 21 complements zeros before and after the input time-series data if the length is shorter than N, and removes values before and after the input time-series data if the length is longer than N to obtain the time-series data x _n having the length N. To be. When complementing zero, the central part of the input time-series data is made to coincide with the central part of the time-series data _xn . When removing the data, it depends on the range to be analyzed, but it is also possible to simply remove the subsequent data.

Then component X _k ^cq constant Q transform, the K _M determined by calculation method boundary determining unit 14, for k ≦ K _M in the low frequency bandwidth calculation unit 23, k> for K _M is the high-frequency bandwidth calculation unit 22 Calculate with

Note that the component of the constant Q conversion may be referred to as “X _{k, m} ^cq ” below.

For each k, the number of non-zero time-axis coefficients T _{n, k, m} calculated by the time-axis coefficient calculation unit 12 is at most N _k , so that the product-sum calculation is efficiently limited to a non-zero range. Can calculate.

Low frequency band calculation unit 23 calculates the Fourier coefficients X _n and FFT once for the entire time-series data x _n first. Then for each p _m, the frequency axis coefficient frequency axis coefficient calculated by the calculating section 13 (1 / N) S n , k, m and component X _k of the Fourier coefficients X _n to the product-sum calculation constant Q transform _{, m} ^cq .

Since the number of non-zero frequency axis coefficients (1 / N) S _{n, k, m} is small for each k, the product-sum calculation can be efficiently performed by limiting the product-sum calculation to a non-zero range.

  In many cases, the non-zero term and the zero term are mixed at both ends (near the frequency zero and the vicinity of the sampling frequency), but the zero term continues in the middle part. , Both ends may be calculated as non-zero terms. In this case, the number of non-zero terms slightly increases, but the number is not so large, and since the non-zero determination logic can be simplified, the calculation time is not much different.

Finally, the constant-Q conversion synthesizing unit 24 puts together the constant-Q conversion components X _{k and} ^mcq calculated by the high-frequency band calculation unit 22 and the low-frequency band calculation unit 23, and completes the constant-Q conversion.

In Example 3, the time axis coefficient T _{n, k, m,} the frequency axis coefficient (1 / N) S n, k, since _m is prepared for each analysis position p _m, keep a prepared beforehand The data volume increases. In the fourth embodiment, the following countermeasures were taken to significantly reduce the memory usage of coefficient data.

  Note that, in equation (4),

  Represents a power which differs for each time and frequency of the complex exponential function.

  The specific configuration of the constant Q conversion component calculation device according to the fourth embodiment is the same as that of FIG. 5, and only the portions different from FIG. 5 will be described below.

The time axis coefficient calculation unit 12 obtains a time axis coefficient T _{n, k, 1} based on the parameters determined by the parameter determination unit 11 according to the definition formula. That is,

And calculated in one analysis position (e.g. the first analysis position) p ₁ of the preset analysis position p _m, the time axis coefficient T _n, obtains the _{k, 1.}

The frequency axis coefficient calculation unit 13 calculates the frequency axis coefficient (1 / N) S _{n, k, 1} based on the parameters determined by the parameter determination unit 11 according to the definition formula.

  That is,

(1 / N) S _{n, k obtained} from the above is calculated at the _one analysis position p ₁ to obtain a frequency axis coefficient (1 / N) S _{n, k, 1} .

At this time, the term of | S _{n, k, 1} | ≦ Th is set to the frequency axis coefficient (1 / N) S _{n, k, 1} = 0.

In the calculation method boundary determination unit 14, since the product sum of the complex exponentiation is added in the frequency axis calculation, the calculation amount C _F for each non-zero frequency axis coefficient is larger than that in the third embodiment, for example, 5 (complex). The calculation of exponentiation is +3, and the product sum of complex exponentiation is +1).

High frequency band calculation unit 22, by shifting the range for each of a plurality of analysis locations p _m, coefficient time axis calculated by the time-axis coefficient calculating section 12 T _{n, k, 1} and the time-series data x _{n-p1 + pm} To obtain the constant X-transform components X _{k, m} ^cq .

Low frequency band calculation unit 23 calculates the Fourier coefficients X _n and FFT once for the entire time-series data x _n first.

  According to the methods of the third and fourth embodiments, in performing the constant Q conversion at a high speed, the preceding speed-up technique (the second embodiment) uses both the calculation on the time axis and the calculation on the frequency axis. In the method, the constant Q conversion is calculated at each time, but when the constant Q conversion at a plurality of times is calculated at the same time, the DFT (or FFT) is applied not to the calculation range at one time but to the whole. Performing it only once speeds up the calculation as a whole.

  In the fourth embodiment, instead of having the coefficient data of the time axis and the frequency axis at each calculation time, the coefficient data is replaced by exponentiation of these complex numbers (Equation (4)), and the amount of calculation is slightly increased and the memory of the coefficient data is used. The amount can be greatly reduced.

DESCRIPTION OF SYMBOLS 10 ... Preparation part 11 ... Parameter determination part 12 ... Time axis coefficient calculation part 13 ... Frequency axis coefficient calculation part 14 ... Calculation method boundary determination part 20 ... Individual calculation part 21 ... Preprocessing part 22 ... High frequency band calculation part 23 ... Low frequency Band calculator 24 Constant Q conversion synthesizer 100 Diagnosis parameter creation means 101 First constant Q converter 102 Principal component analyzer 103 First statistical value calculator 104 First normalizer 200 Diagnostic means 201 Second constant Q converter 202 Principal component calculating unit 203 Second statistical value calculating unit 204 Second normalizing unit

Claims (5)

Diagnostic parameter creation means for creating a multivariate sample by performing constant Q conversion on time series data measured in advance from an abnormality diagnosis target and creating a diagnostic parameter based on the created multivariate sample;
Time-series data measured from freshly abnormality diagnostic object by converting the constant Q to create a multivariate sample, and a diagnostic means for diagnosing an abnormality using the diagnostic parameters created by the diagnostic parameter generation means ,
The diagnostic parameter creating means,
A first constant-Q conversion device that performs constant-Q conversion on time-series data to create a multivariate sample;
A principal component analysis unit that calculates a principal component score of each sample while obtaining a transformation matrix by performing principal component analysis on the created multivariate sample;
A first statistic calculation unit that obtains a statistic serving as an index of a sample based on the principal component score;
A first normalization unit that normalizes the statistic to determine the degree of abnormality,
The diagnostic means includes:
A second constant-Q conversion device that performs constant-Q conversion on the time-series data to create a multivariate sample;
From the created multivariate sample, a principal component calculation unit that obtains a principal component score using a transformation matrix obtained by the principal component analysis unit,
A second statistic calculation unit that obtains a statistic based on the principal component score;
A second normalization unit that normalizes the statistic by a normalization coefficient used in the first normalization unit to calculate an abnormality degree, and based on the abnormality degree, Diagnose,
The first of said diagnostic parameter generation means and diagnostic means, the second constant Q converter are both
In the set frequency components K _M or lower frequency band, a first amount of calculation of the constant Q transform components determined by the product-sum calculation between the frequency axis coefficients in the discrete Fourier transform and a constant Q transform of the time series data, the K _In a high-frequency band exceeding _M , a total operation amount including the second operation amount of the constant Q conversion component obtained by the product-sum calculation of the time series data and the time axis coefficient in the constant Q conversion is calculated,
Obtains the total calculation amount for each K _M by each setting the K _M to all frequency components K,
A calculation method boundary determining unit that evaluates a total calculation amount for each of the K _M , and determines a K _M with the minimum calculation amount as a frequency component threshold indicating a calculation method boundary;
The low-frequency-band frequency components K is less than the frequency component threshold value K _M which is the determined in the constant Q transform, when selecting the first calculation method is a product-sum computation of the discrete Fourier transform and the frequency axis coefficient series data and, in the high frequency band in which the frequency components K exceeds the frequency component threshold value K _M which is the determined, the calculation method selection means for selecting the second calculation method is a product-sum computation of the time-series data and the time-axis coefficient ,
Constant Q conversion means for performing a first calculation method and a second calculation method selected by the calculation method selection means on the time-series data to obtain a constant Q conversion component;
An abnormality diagnosis device comprising: Each said constant Q conversion device,
The time axis coefficient in the constant Q conversion is

A time axis coefficient calculation unit calculated by:
This is a discrete Fourier transform when the frequency axis coefficient in the constant Q conversion is regarded as time series data of time axis coefficient T _{n, k} with n = 0,..., N−1.

And a frequency axis coefficient calculation unit that calculates as (1 / N) S _{n, k} from
The constant Q conversion means includes:
A low-frequency band calculation unit that performs the first calculation method by performing a product-sum calculation on the frequency axis coefficient and the discrete Fourier transform of the time-series data calculated by the frequency axis coefficient calculation unit;
A high-frequency band calculation unit that performs a product-sum calculation of the time-axis coefficient and time-series data calculated by the time-axis coefficient calculation unit and performs the second calculation method;
A combining unit that combines the outputs of the low-frequency band calculation unit and the high-frequency band calculation unit to obtain a component of constant Q conversion;
The abnormality diagnosis device according to claim 1 , further comprising: Each said constant Q conversion device,

A plurality of analysis position of the preset time direction _{p m (m = 1, ...} , M) time axis coefficient calculated for each T _{n, k,} and _m determined, the constant said T _{n, k,} and _m A time axis coefficient calculating unit that sets a time axis coefficient in the Q conversion;
This is a discrete Fourier transform when the time axis coefficient T _{n, k} is regarded as time series data of n = 0,..., N−1.

Was obtained from (1 / N) S _n, the _k, determined frequency axis coefficient (1 / N) S n, _k, and _m are calculated for each said analysis position p _m, the (1 / N) S _n, a frequency axis coefficient calculating unit that sets _{k and m} as frequency axis coefficients in the constant Q conversion,
The constant Q conversion means includes:
When asked to Fourier coefficient X _n by performing a discrete Fourier transform for the entire series data x _n, for each of the analyzed position p _m, the frequency axis coefficient calculated by the frequency axis coefficient calculator (1 / N) S _{n , k, m} and the Fourier coefficient X _n to calculate the sum of products and perform the first calculation method, a low frequency band calculation unit,
For each of the analysis positions p _m , a high-frequency wave for performing the second calculation method by performing a product-sum calculation on the time axis coefficient T _{n, k, m} calculated by the time axis coefficient calculation unit and the time series data x _n A bandwidth calculator,
A combining unit that combines the outputs of the low-frequency band calculation unit and the high-frequency band calculation unit to obtain a component of constant Q conversion;
The abnormality diagnosis device according to claim 1 , further comprising: Each said constant Q conversion device,

A plurality of analysis position of the preset time direction _{p m (m = 1, ...} , M) time axis coefficient T _n is calculated in one analysis position p ₁ of _the, seeking _{k, 1,} the T _{n , k, 1} as a time axis coefficient in the constant Q conversion,
This is a discrete Fourier transform when the time axis coefficient T _{n, k} is regarded as time series data of n = 0,..., N−1.

Was obtained from (1 / N) S _n, the _k, the one frequency axis coefficient calculated in the analysis position p ₁ of (1 / N) S n, obtains a _{k, 1,} the (1 / N) S _{n , k, 1} as frequency axis coefficients in the constant Q conversion,
The constant Q conversion means includes:
A discrete Fourier transform is performed on the entire time series data x _n to obtain a Fourier coefficient X _n, and for each of the analysis positions p _m ,
A frequency axis coefficient (1 / N) S _{n, k, 1} obtained by the frequency axis coefficient calculation unit;
There the analysis frequency axis coefficient for each position _{p m (1 / N) S} n, k, m are in a relationship permitting consequence the frequency axis coefficients in one analysis position _{p 1 (1 / N) S} n, the _{k, 1} Power of complex exponential function derived from

When,
A low-frequency band calculation unit that performs a product-sum calculation of the Fourier coefficient X _n and performs the first calculation method;
By shifting the range for each of the analyzed position p _m, the time-series data x _{n-p1 +} and _pm by product-sum computation the time axis coefficient T _n calculated by the time-axis coefficient calculating _{section, k, 1} and the A high-frequency band calculation unit that performs the calculation method 2;
A combining unit that combines the outputs of the low-frequency band calculation unit and the high-frequency band calculation unit to obtain a component of constant Q conversion;
The abnormality diagnosis device according to claim 1 , further comprising: Diagnostic parameter creation means creates a multivariate sample by performing constant Q conversion on time series data measured in advance from an abnormality diagnosis target, and creates a diagnostic parameter based on the created multivariate sample. Steps and
A diagnosing step of diagnosing an abnormality using the diagnostic parameter created by the diagnostic parameter creating means, wherein the diagnostic means creates a multivariate sample by performing constant Q conversion on the time series data newly measured from the abnormality diagnostic object; and, with a,
The diagnostic parameter creating step,
A first constant-Q conversion device performing constant-Q conversion on the time-series data to create a multivariate sample;
A principal component analysis unit for performing a principal component analysis on the created multivariate sample to obtain a transformation matrix and calculating a principal component score of each sample;
A step in which a first statistic calculation unit obtains a statistic serving as an index of a sample based on the principal component score;
A first normalizing unit that normalizes the statistic to obtain an abnormal degree,
The diagnosis step includes:
A second constant-Q conversion device performing constant-Q conversion on the time-series data to create a multivariate sample;
A principal component calculation unit, from the created multivariate sample, using a transformation matrix obtained by the principal component analysis unit to obtain a principal component score;
A second statistic calculation unit for obtaining a statistic based on the principal component score;
A second normalization unit for normalizing the statistic by the normalization coefficient used in the first normalization unit to calculate an abnormality degree, and an abnormality diagnosis target based on the abnormality degree Diagnose abnormalities ,
The step of the first constant Q conversion device performing constant Q conversion of the time series data to create a multivariate sample,
In the set frequency components K _M or lower frequency band, a first amount of calculation of the constant Q transform components determined by the product-sum calculation between the frequency axis coefficients in the discrete Fourier transform and a constant Q transform of the time series data, the K _In a high-frequency band exceeding _M , a total operation amount including the second operation amount of the constant Q conversion component obtained by the product-sum calculation of the time series data and the time axis coefficient in the constant Q conversion is calculated,
Obtains the total calculation amount for each K _M by each setting the K _M to all frequency components K,
The evaluated total calculation amount for each K _M, comprising the step of determining the K _M of the amount of calculation is minimized, as a frequency component threshold indicating the calculation method boundary,
The step of the second constant-Q conversion device performing constant-Q conversion on the time-series data to create a multivariate sample,
The low-frequency-band frequency components K is less than the frequency component threshold value K _M which is the determined in the constant Q transform, when selecting the first calculation method is a product-sum computation of the discrete Fourier transform and the frequency axis coefficient series data a step, and the frequency component K is for selecting the second calculation method is a product-sum calculation of said time axis factor, the time series data and high frequency band exceeding the frequency component threshold value K _M which is the determined,
Performing a selected first calculation method and a second calculation method on the time-series data to obtain a component of constant Q conversion .

Images