JP2017198620A - Abnormality diagnosis device and abnormality diagnosis method - Google Patents
Abnormality diagnosis device and abnormality diagnosis method Download PDFInfo
- Publication number
- JP2017198620A JP2017198620A JP2016091695A JP2016091695A JP2017198620A JP 2017198620 A JP2017198620 A JP 2017198620A JP 2016091695 A JP2016091695 A JP 2016091695A JP 2016091695 A JP2016091695 A JP 2016091695A JP 2017198620 A JP2017198620 A JP 2017198620A
- Authority
- JP
- Japan
- Prior art keywords
- constant
- conversion
- calculation
- time
- series data
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Granted
Links
Images
Landscapes
- Testing Of Devices, Machine Parts, Or Other Structures Thereof (AREA)
- Measurement Of Mechanical Vibrations Or Ultrasonic Waves (AREA)
- Testing And Monitoring For Control Systems (AREA)
Abstract
Description
本発明は、異常診断対象機器の動作時の振動・音響などの波形データ（時系列データ）を周波数解析して異常診断を行う異常診断装置および方法に関する。 The present invention relates to an abnormality diagnosis apparatus and method for performing abnormality diagnosis by performing frequency analysis on waveform data (time-series data) such as vibration and sound during operation of an abnormality diagnosis target device.
異常診断対象機器の異常を診断する方法は、例えば特許文献１に記載のものが提案されている。
As a method for diagnosing an abnormality in an abnormality diagnosis target device, for example, a method described in
この特許文献１の方法では、時系列の波形データを、短い時間毎に分割して、フーリエ変換して周波数成分の多変量データのサンプルを多数作り、機器が正常と考えられる期間に予め計測した多数の多変量サンプルを主成分分析して主成分得点を求めるための固有ベクトル（ローテーション行列ともいう）を用意する。
In the method of
診断の際には、計測データを前記同様にフーリエ変換して周波数成分の多変量データ・サンプルを作り、これを用意しておいたローテーション行列で変換して主成分得点を求め、それを正常時の主成分得点と比較することで異常診断を行っている。要するに診断毎に独立に主成分分析するのではなく、予め主成分分析しておいた結果を利用して、診断時には同じ基準で変換することにより比較を容易にするのが要点である。 At the time of diagnosis, the measurement data is Fourier-transformed in the same manner as described above to create multivariate data samples of frequency components, and this is converted with the rotation matrix prepared to obtain the principal component score, which is obtained under normal conditions. Abnormal diagnosis is performed by comparing with the main component score of. In short, instead of performing principal component analysis independently for each diagnosis, it is important to make comparison easier by using the result of principal component analysis in advance and converting according to the same standard at the time of diagnosis.
尚、本発明で利用する定Ｑ変換の成分演算方法は、例えば非特許文献１〜４に記載されている。
The constant Q conversion component calculation method used in the present invention is described in
フーリエ変換を基にした多変量サンプルで予め作成したローテーション行列によって長期にわたって診断しようとすると、事前の解析には多くのデータを必要とする。例えば回転機の診断においては、少なくとも低い周波数領域では１Ｈｚ程度の周波数解像度が必要とされる。 If a long-term diagnosis is made 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 the 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 increased and the number of samples is increased, it becomes essentially a similar sample and lacks diversity. Then, in the long term, it will not be possible to distinguish between abnormalities and changes over time.
以下に、その事例を述べる。この事例は、ある回転機の振動を長期間（約４ヶ月）にわたって測定して分析したものである。１時間に一回、約１００ｋＨｚのサンプリング周波数で５秒間の測定を行い、これから時間を重複させながら１秒強の分割データを１測定から１００サンプル作成して、各サンプルをフーリエ変換して１／６オクターブバンドの約８０系列の多変量サンプルにした。この１週間分（ただし稼動している時刻のデータに絞ったため都合８千サンプルほど）を主成分分析してローテーション行列を計算し、これを基にホテリングＴ２／Ｑ統計量をベースとした異常度を図６のように算出している。 The following are examples. In this case, the vibration of a certain rotating machine is measured and analyzed over a long period (about 4 months). Once every hour, measurement is performed for 5 seconds at a sampling frequency of about 100 kHz, and from this, divided data of just over 1 second is created from 1 measurement to 100 samples, and each sample is subjected to Fourier transform and 1 / There were about 80 series of multivariate samples in 6 octave bands. The degree of abnormality based on the Hotelling T2 / Q statistic is calculated based on the principal component analysis of this one week (however, about 8,000 samples because it was narrowed down to the data at the time of operation). Is calculated as shown in FIG.
対象回転機は特に問題のないものであったが、ホテリングＴ２統計量をベースとしたＩｎｄｉｃａｔｏｒ１、Ｑ統計量をベースとしたＩｎｄｉｃａｔｏｒ２とも、図６の「線形」の傾きから分かるように、時間の経過とともに増加していることが伺える。この増加の割合は、学習期間（７日）の１０倍（７０日）もあれば平均値が異常になるほどであり、長期的な診断にはとても使えない。
The target rotating machine was not particularly problematic. However, both the
この問題は異常度の定義に依存するものではなく、全期間データを主成分分析した図７の第一主成分×第二主成分の散布状況をみても、白色の丸印で示す学習期間のデータが、黒丸印で示す全データの範囲をカバーできていないことは明白である。 This problem does not depend on the definition of the degree of abnormality. Looking at the distribution of the first principal component × second principal component in FIG. It is clear that the data does not 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 thought to be because the diversity of data cannot be ensured because Fourier transform is used to create multivariate data from time series data.
本発明は上記課題を解決するものであり、その目的は、サンプルの多様性を確保し、長期間にわたって信頼性の高い異常診断を実施することができる異常診断装置および方法を提供することにある。 The present invention solves the above-described problems, and an object of the present invention is to provide an abnormality diagnosis apparatus and method capable of ensuring a variety of samples and performing reliable abnormality diagnosis over a long period of time. .
上記課題を解決するための請求項１に記載の異常診断装置は、予め異常診断対象から測定した時系列データを定Ｑ変換して多変量サンプルを作成し、前記作成された多変量サンプルを基に診断用パラメータを作成する診断用パラメータ作成手段と、
新たに異常診断対象から測定した時系列データを定Ｑ変換して多変量サンプルを作成し、前記診断用パラメータ作成手段により作成された診断用パラメータを用いて異常を診断する診断手段と、を備えたことを特徴とする。
The abnormality diagnosis apparatus according to
A diagnostic means for diagnosing an abnormality using the diagnostic parameter created by the diagnostic parameter creation means by creating a multivariate sample by performing constant Q conversion on the time series data newly measured from the abnormality diagnostic target; It is characterized by that.
また、請求項２に記載の異常診断装置は、請求項１において、前記診断用パラメータ作成手段は、
時系列データを定Ｑ変換して多変量サンプルを作成する第１の定Ｑ変換装置と、
前記作成された多変量サンプルを主成分分析して変換行列を得るとともに各サンプルの主成分得点を計算する主成分分析部と、
前記主成分得点を基にサンプルの指標となる統計量を得る第１の統計値計算部と、
前記統計量を正規化して異常度とする第１の正規化部と、を備え、
前記診断手段は、
時系列データを定Ｑ変換して多変量サンプルを作成する第２の定Ｑ変換装置と、
前記作成された多変量サンプルから、前記主成分分析部で得られた変換行列を使って主成分得点を得る主成分計算部と、
前記主成分得点を基に統計量を得る第２の統計値計算部と、
前記統計量を、前記第１の正規化部で使用した正規化係数によって正規化して異常度を計算する第２の正規化部と、を備え、前記異常度に基づいて異常診断対象の異常を診断することを特徴とする。
Further, the abnormality diagnosis apparatus according to claim 2 is characterized in that, in
A first constant Q conversion device for generating multivariate samples by performing constant Q conversion on time series data;
A principal component analysis unit that obtains a transformation matrix by performing principal component analysis on the created multivariate sample and calculates a principal component score of each sample;
A first statistic calculator that obtains a statistic serving as a sample index based on the principal component score;
A first normalization unit that normalizes the statistics and sets the degree of abnormality,
The diagnostic means includes
A second constant Q conversion device for generating multivariate samples by performing constant Q conversion on time series data;
A principal component calculation unit that obtains a principal component score from the created multivariate sample using the transformation matrix obtained by the principal component analysis unit;
A second statistical value calculation unit for obtaining a statistic based on the principal component score;
A second normalization unit that calculates the degree of abnormality by normalizing the statistic with the normalization coefficient used in the first normalization unit, and detects an abnormality of the abnormality diagnosis target based on the degree of abnormality It is characterized by making a diagnosis.
また、請求項７に記載の異常診断方法は、診断用パラメータ作成手段が、予め異常診断対象から測定した時系列データを定Ｑ変換して多変量サンプルを作成し、前記作成された多変量サンプルを基に診断用パラメータを作成する診断用パラメータ作成ステップと、
診断手段が、新たに異常診断対象から測定した時系列データを定Ｑ変換して多変量サンプルを作成し、前記診断用パラメータ作成手段により作成された診断用パラメータを用いて異常を診断する診断ステップと、を備えたことを特徴とする。
The abnormality diagnosis method according to claim 7, wherein the 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 the created multivariate sample A diagnostic parameter creation step for creating a diagnostic parameter based on
A diagnostic step in which a diagnostic means creates a multivariate sample by performing constant Q conversion on time series data newly measured from an abnormality diagnosis target, and diagnoses an abnormality using the diagnostic parameters created by the diagnostic parameter creation means And.
また、請求項８に記載の異常診断方法は、請求項７において、前記診断用パラメータ作成ステップは、
第１の定Ｑ変換装置が、時系列データを定Ｑ変換して多変量サンプルを作成するステップと、
主成分分析部が、前記作成された多変量サンプルを主成分分析して変換行列を得るとともに各サンプルの主成分得点を計算するステップと、
第１の統計値計算部が、前記主成分得点を基にサンプルの指標となる統計量を得るステップと、
第１の正規化部が、前記統計量を正規化して異常度とするステップと、を備え、
前記診断ステップは、
第２の定Ｑ変換装置が、時系列データを定Ｑ変換して多変量サンプルを作成するステップと、
主成分計算部が、前記作成された多変量サンプルから、前記主成分分析部で得られた変換行列を使って主成分得点を得るステップと、
第２の統計値計算部が、前記主成分得点を基に統計量を得るステップと、
第２の正規化部が、前記統計量を、前記第１の正規化部で使用した正規化係数によって正規化して異常度を計算するステップと、を備え、前記異常度に基づいて異常診断対象の異常を診断することを特徴とする。
The abnormality diagnosis method according to
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 obtains a transformation matrix by performing principal component analysis on the created multivariate sample and calculates a principal component score of each sample;
A first statistic calculation unit obtaining a statistic serving as an index of a sample based on the principal component score;
A first normalization unit comprising: normalizing the statistic to obtain an abnormality level;
The diagnostic 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 obtains a principal component score from the created multivariate sample using the transformation matrix obtained by the principal component analysis unit;
A second statistic calculator calculates a statistic based on the principal component score;
A second normalization unit comprising: calculating the degree of abnormality by normalizing the statistic with the normalization coefficient used in the first normalization unit; and an abnormality diagnosis target based on the degree of abnormality It is characterized by diagnosing abnormalities.
上記構成によれば、時系列データから周波数成分の多変量データを作成する際に、フーリエ変換でなく定Ｑ変換を使ったことにより、診断用パラメータ作成手段で使用するデータが比較的少なくてもサンプルの多様性を確保できるようになり、長期間にわたって同じ基準で異常診断が可能になる。これによって異常診断の信頼性が向上する。 According to the above configuration, when creating multivariate data of frequency components from time series data, the constant Q transform is used instead of the Fourier transform, so that relatively little data is used in the diagnostic parameter creation means. The diversity of samples can be secured, and abnormality diagnosis can be performed based on the same standard over a long period of time. This improves the reliability of abnormality diagnosis.
また、請求項３に記載の異常診断装置は、請求項２において、前記診断用パラメータ作成手段および診断手段の第１、第２の定Ｑ変換装置は、ともに、
設定した周波数成分ＫM以下の低周波数帯域において、時系列データの離散フーリエ変換と定Ｑ変換における周波数軸係数との積和計算により求めた定Ｑ変換成分の第１の演算量と、前記ＫMを超える高周波帯域において、時系列データと定Ｑ変換における時間軸係数との積和計算により求めた定Ｑ変換成分の第２の演算量とを含む合計演算量を算出し、
前記ＫMを全周波数成分ｋに各々設定して各ＫM毎の前記合計演算量を求め、
前記各ＫM毎の合計演算量を評価し、該演算量が最小となるＫMを、計算法境界を示す周波数成分閾値として決定する計算法境界決定手段と、
定Ｑ変換における周波数成分ｋが前記決定された周波数成分閾値ＫM以下の低周波帯域では、時系列データの離散フーリエ変換と前記周波数軸係数との積和計算である第１の計算法を選択し、前記周波数成分ｋが前記決定された周波数成分閾値ＫMを超える高周波帯域では、時系列データと前記時間軸係数との積和計算である第２の計算法を選択する計算法選択手段と、
時系列データに対して、前記計算法選択手段により選択された第１の計算法、第２の計算法を実施して定Ｑ変換の成分を求める定Ｑ変換手段と、
を備えたことを特徴とする。
Further, the abnormality diagnosis apparatus according to claim 3 is characterized in that, in claim 2, both the first and second constant Q conversion apparatuses of the diagnostic parameter creation means and the diagnosis 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 calculation amount including the second calculation 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,
To evaluate the total amount of computation of the respective K M, the K M of the calculation amount is minimized, and the calculation method boundary determining means for determining a frequency component threshold indicating the calculation method boundary,
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, in the high frequency band in which the frequency component 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 ,
A 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 provided with.
また、請求項４に記載の異常診断装置は、請求項３において、前記各定Ｑ変換装置は、
前記定Ｑ変換における時間軸係数を、
The abnormality diagnosis device according to claim 4 is the abnormality diagnosis device according to claim 3, wherein each of the constant Q conversion devices is
The time axis coefficient in the constant Q conversion is
によって算出する時間軸係数算出部と、
前記定Ｑ変換における周波数軸係数を、時間軸係数Ｔn,kをｎ＝０，…，Ｎ−１の時系列データとみたときの離散フーリエ変換である
A time axis coefficient calculation unit calculated by
The frequency axis coefficient in the constant Q conversion is a discrete Fourier transform when the time axis coefficient T n, k is regarded as time series data of n = 0,..., N−1.
から（１／Ｎ）Ｓn,kとして算出する周波数軸係数算出部と、を備え、
前記定Ｑ変換手段は、
前記周波数軸係数算出部で算出された周波数軸係数と時系列データの離散フーリエ変換とを積和計算して前記第１の計算法を実施する低周波帯域計算部と、
前記時間軸係数算出部で算出された時間軸係数と時系列データとを積和計算して前記第２の計算法を実施する高周波帯域計算部と、
前記低周波帯域計算部と高周波帯域計算部の各出力を合成して定Ｑ変換の成分を求める合成部と、
を備えていることを特徴とする。
And (1 / N) S n, k , a frequency axis coefficient calculation unit
The constant Q conversion means includes:
A low frequency band calculation unit for performing the first calculation method by multiplying and multiplying the frequency axis coefficient calculated by the frequency axis coefficient calculation unit and the discrete Fourier transform of the time series data;
A high-frequency band calculation unit that performs a product-sum calculation on the time axis coefficient calculated by the time axis coefficient calculation unit and time-series data and implements the second calculation method;
A synthesis unit for synthesizing outputs of the low frequency band calculation unit and the high frequency band calculation unit to obtain a constant Q conversion component;
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 speed can be increased.
また、請求項５に記載の異常診断装置は、請求項３において、前記各定Ｑ変換装置は、
The abnormality diagnosis device according to
を、予め設定した時間方向の複数の解析位置ｐm（ｍ＝１，…，Ｍ）毎に算出して時間軸係数Ｔn,k,mを求め、該Ｔn,k,mを前記定Ｑ変換における時間軸係数とする時間軸係数算出部と、
時間軸係数Ｔn,kをｎ＝０，…，Ｎ−１の時系列データとみたときの離散フーリエ変換である
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 calculation unit as a time axis coefficient in Q conversion;
Discrete Fourier transform when the time axis coefficient T n, k is regarded as time series data of n = 0,..., N−1.
から得た（１／Ｎ）Ｓn,kを、前記解析位置ｐm毎に算出して周波数軸係数（１／Ｎ）Ｓn,k,mを求め、該（１／Ｎ）Ｓn,k,mを前記定Ｑ変換における周波数軸係数とする周波数軸係数算出部と、を備え、
前記定Ｑ変換手段は、
時系列データｘn全体に対して離散フーリエ変換を行ってフーリエ係数Ｘnを求め、前記解析位置ｐm毎に、前記周波数軸係数算出部で算出された周波数軸係数（１／Ｎ）Ｓn,k,mと前記フーリエ係数Ｘnとを積和計算して前記第１の計算法を実施する低周波帯域計算部と、
前記解析位置ｐm毎に、前記時間軸係数算出部で算出された時間軸係数Ｔn,k,mと時系列データｘnとを積和計算して前記第２の計算法を実施する高周波帯域計算部と、
前記低周波帯域計算部と高周波帯域計算部の各出力を合成して定Ｑ変換の成分を求める合成部と、
を備えていることを特徴とする。
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 calculation 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, and a low frequency band calculation unit for performing the first calculation method by calculating the sum of products,
For each of the analyzed position p m, implementing the time-axis coefficient calculating unit time axis coefficients calculated in T n, k, m and the time-series data x n and then product sum calculation of the second calculation method RF A bandwidth calculator;
A synthesis unit for synthesizing outputs of the low frequency band calculation unit and the high frequency band calculation unit to obtain a constant Q conversion component;
It is characterized by having.
上記構成によれば、複数の時刻で各々定Ｑ変換を行う際に、低周波帯域における周波数軸での計算は、時系列データの離散フーリエ変換を、複数の時刻毎に行わなくても、ｘn全体に対して一括して離散フーリエ変換したフーリエ係数Ｘnと各解析位置での周波数軸係数との積和計算でよいため、全体としての計算が高速化される。 According to the above configuration, when performing the constant Q transform at each of a plurality of times, the calculation on the frequency axis in the low frequency band can be performed without performing the discrete Fourier transform of the time series data for each of the plurality of times. Since the product sum calculation of the Fourier coefficient Xn obtained by performing discrete Fourier transform on the whole n and the frequency axis coefficient at each analysis position may be performed, the calculation as a whole is accelerated.
また、請求項６に記載の異常診断装置は、請求項３において、前記各定Ｑ変換装置は、 The abnormality diagnosis device according to claim 6 is the abnormality diagnosis device according to claim 3, wherein each of the constant Q conversion devices is
を、予め設定した時間方向の複数の解析位置ｐm（ｍ＝１，…，Ｍ）のうち１つの解析位置ｐ1において算出して時間軸係数Ｔn,k,1を求め、該Ｔn,k,1を前記定Ｑ変換における時間軸係数とする時間軸係数算出部と、
時間軸係数Ｔn,kをｎ＝０，…，Ｎ−１の時系列データとみたときの離散フーリエ変換である
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 1 of the, seeking k, 1, the T n , k, 1 as a time axis coefficient in the constant Q conversion,
Discrete Fourier transform when the time axis coefficient T n, k is regarded as time series data of n = 0,..., N−1.
から得た（１／Ｎ）Ｓn,kを、前記１つの解析位置ｐ1において算出して周波数軸係数（１／Ｎ）Ｓn,k,1を求め、該（１／Ｎ）Ｓn,k,1を前記定Ｑ変換における周波数軸係数とする周波数軸係数算出部と、を備え、
前記定Ｑ変換手段は、
時系列データｘn全体に対して離散フーリエ変換を行ってフーリエ係数Ｘnを求め、 前記解析位置ｐm毎に、
前記周波数軸係数算出部で求められた周波数軸係数（１／Ｎ）Ｓn,k,1と、
前記解析位置ｐm毎の周波数軸係数（１／Ｎ）Ｓn,k,mは前記１つの解析位置ｐ1における周波数軸係数（１／Ｎ）Ｓn,k,1に帰結できる関係にあることから導かれる、複素指数関数のべき乗
(1 / N) S n, k obtained from the above is calculated at the one analysis position p 1 to obtain the frequency axis coefficient (1 / N) S n, k, 1 , and the (1 / N) S n , k, 1 as a frequency axis coefficient in the constant Q conversion, and a frequency axis coefficient calculation unit,
The constant Q conversion means includes:
Calculated Fourier coefficients X n by performing a discrete Fourier transform for the entire time-series data x n, for each of the analyzed position p m,
The frequency axis coefficient (1 / N) S n, k, 1 obtained by the frequency axis coefficient calculating 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 the complex exponential function
と、
前記フーリエ係数Ｘnとを積和計算して前記第１の計算法を実施する低周波帯域計算部と、
前記解析位置ｐm毎に範囲をずらして、前記時間軸係数算出部で算出された時間軸係数Ｔn,k,1と時系列データｘn-p1+pmとを積和計算して前記第２の計算法を実施する高周波帯域計算部と、
前記低周波帯域計算部と高周波帯域計算部の各出力を合成して定Ｑ変換の成分を求める合成部と、
を備えていることを特徴とする。
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 for performing the calculation method of 2,
A synthesis unit for synthesizing outputs of the low frequency band calculation unit and the high frequency band calculation unit to obtain a constant Q conversion component;
It is characterized by having.
上記構成によれば、複数の時刻で各々定Ｑ変換を行う際に、１つの解析位置ｐ1における時間軸係数および周波数軸係数を利用して各時刻での計算が可能となるので、複数の時刻毎（複数の解析位置毎）に時間軸係数および周波数軸係数を用意する必要はなく、データ量が少なくなってメモリを節約することができる。 According to the above configuration, when performing constant Q conversion at a plurality of times, calculation at each time can be performed using the time axis coefficient and the frequency axis coefficient at one analysis position p 1 . 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.
（１）請求項１〜８に記載の発明によれば、時系列データから周波数成分の多変量データを作成する際に、フーリエ変換でなく定Ｑ変換を使ったことにより、診断用パラメータ作成手段で使用するデータが比較的少なくてもサンプルの多様性を確保できるようになり、長期間にわたって同じ基準で異常診断が可能になる。これによって異常診断の信頼性が向上する。
（２）請求項３、４に記載の発明によれば、低周波帯域、高周波帯域のどちらの場合でも演算量の少ない計算法で定Ｑ変換の成分を演算することができ、演算の高速化が実現される。
（３）請求項５に記載の発明によれば、複数の時刻で各々定Ｑ変換を行う際に、低周波帯域における周波数軸での計算は、時系列データの離散フーリエ変換を、複数の時刻毎に行わなくても、ｘn全体に対して一括して離散フーリエ変換したフーリエ係数Ｘnと各解析位置での周波数軸係数との積和計算でよいため、全体としての計算が高速化される。
（４）請求項６に記載の発明によれば、複数の時刻で各々定Ｑ変換を行う際に、１つの解析位置ｐ1における時間軸係数および周波数軸係数を利用して各時刻での計算が可能となるので、複数の時刻毎（複数の解析位置毎）に時間軸係数および周波数軸係数を用意する必要はなく、データ量が少なくなってメモリを節約することができる。
(1) According to the first to eighth aspects of the present invention, when creating multivariate data of frequency components from time series data, a constant Q transform is used instead of a Fourier transform. The sample diversity can be ensured even if the data used in the system is relatively small, and the abnormality can be diagnosed with the same standard over a long period of time. This improves the reliability of abnormality diagnosis.
(2) According to the third and fourth aspects of the 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 speed is increased. Is realized.
(3) According to the invention described in
(4) according According to the invention described in claim 6, 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 1 Therefore, it is not necessary to prepare a time axis coefficient and a frequency axis coefficient for each of a plurality of times (for each of a plurality of analysis positions), and the data amount can be reduced and the memory can be saved.
以下、図面を参照しながら本発明の実施の形態を説明するが、本発明は下記の実施形態例に限定されるものではない。本発明では、時系列データ（波形データ）から周波数成分の多変量サンプルを生成する際に、従来のフーリエ変換の代わりに定Ｑ変換を使う。 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 does not analyze data in the same period in all frequency bands as in Fourier transform, but refers to data for each frequency so that it has the same number of cycles in all frequency bands. Analyze by changing the number. In this method, since data is analyzed with a shorter period as the high frequency band, when a long period (1 second or more) of measurement is required to increase the frequency resolution in the low frequency band (for example, about 1 Hz), the period is set. Even if a large number of samples are overlapped, the periods do not overlap in the high frequency band, so that diversity is ensured as a multivariate sample.
具体的な診断手法の構成を図１に示す。図１において、本実施例１の手法は、予め測定して集めた時系列データから診断のためのパラメータを計算する準備フェーズと、新たに測定された時系列データを診断する診断フェーズからなる。 A specific configuration of the diagnostic technique is shown in FIG. In FIG. 1, the method of the first embodiment includes a preparation phase for calculating parameters for diagnosis from time series data collected in advance and a diagnosis phase for diagnosing newly measured time series data.
準備フェーズは、予め異常診断対象から測定して集めた時系列データを定Ｑ変換して多変量サンプルを作成し、前記作成された多変量サンプルを基に診断用パラメータを作成する診断用パラメータ作成手段１００として構成している。 The preparation phase creates a multivariate sample by performing constant Q conversion on time-series data collected and measured in advance from an abnormality diagnosis target, and creates a diagnostic parameter based on the created multivariate sample The means 100 is configured.
診断フェーズは、新たに異常診断対象から測定した時系列データを定Ｑ変換して多変量サンプルを作成し、前記診断用パラメータ作成手段により作成された診断用パラメータを用いて異常を診断する診断手段２００として構成している。
The diagnostic phase is a diagnostic means for diagnosing an abnormality using the diagnostic parameter created by the diagnostic parameter creation means by creating a multivariate sample by performing constant Q conversion on the time series data newly measured from the abnormality
診断用パラメータ作成手段１００は、予め異常診断対象から測定して集めた時系列データを定Ｑ変換して多変量サンプルを作成する第１の定Ｑ変換装置１０１と、前記作成された多変量サンプルを主成分分析して変換行列を得るとともに各サンプルの主成分得点を計算する主成分分析部１０２と、前記主成分得点を基にサンプルの指標となる統計量を得る第１の統計値計算部１０３と、前記統計量を正規化して異常度とする第１の正規化部１０４と、を備えている。
The diagnostic parameter creation means 100 includes a first constant
前記診断手段２００は、新たに測定された時系列データを定Ｑ変換して多変量サンプルを作成する第２の定Ｑ変換装置２０１と、前記作成された多変量サンプルから、前記主成分分析部１０２で得られた変換行列を使って主成分得点を得る主成分計算部２０２と、前記主成分得点を基に統計量を得る第２の統計値計算部２０３と、前記統計量を、前記第１の正規化部１０４で使用した正規化係数によって正規化して異常度を計算する第２の正規化部２０４と、を備え、前記異常度に基づいて異常診断対象の異常を診断するものである。
The
図１の異常診断装置は、例えばコンピュータにより構成され、通常のコンピュータのハードウェアリソース、例えばＲＯＭ、ＲＡＭ、ＣＰＵ、入力装置、出力装置、通信インターフェース、ハードディスク、記録媒体およびその駆動装置を備えている。 The abnormality diagnosis apparatus of FIG. 1 is configured by a computer, for example, and includes hardware resources of a normal computer, such as ROM, RAM, CPU, input device, output device, communication interface, hard disk, recording medium, and driving device thereof. .
このハードウェアリソースとソフトウェアリソース（ＯＳ、アプリケーションなど）との協働の結果、異常診断装置は、図１のように第１の定Ｑ変換装置１０１、主成分分析部１０２、第１の統計値計算部１０３、第１の正規化部１０４、第２の定Ｑ変換装置２０１、主成分計算部２０２、第２の統計値計算部２０３、第２の正規化部２０４を実装する。
As a result of the cooperation between the hardware resource and the software resource (OS, application, etc.), the abnormality diagnosis apparatus has a first constant
前記多変量サンプル、主成分得点、統計量、異常度、変換行列、正規化係数などのデータはハードディスクあるいはＲＡＭなどの保存手段・記憶手段に格納されるものとする。 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 a storage unit / storage unit such as a hard disk or a RAM.
次に、図１の装置の各部の動作を説明する。準備フェーズ（診断用パラメータ作成手段１００）では、予め集めた時系列データを、第１の定Ｑ変換装置１０１において順次定Ｑ変換して多数の多変量サンプルを作成する。定Ｑ変換装置１０１は、例えば非特許文献１〜４に開示されている方式を実行するものでもよく、また後述する実施例２〜４の高速演算が可能な定Ｑ変換装置を採用しても良い。尚、実施例２〜４に示す方法で高速に求める場合には、時系列データのデータ長はそろっていることが望ましい。
Next, the operation of each part of the apparatus of FIG. 1 will be described. In the preparation phase (diagnostic parameter creation means 100), time series data collected in advance is sequentially subjected to constant Q conversion in the first constant
例えば、１０２．４ｋＨｚサンプリングで５秒間のデータが１００件など。定Ｑ変換の高周波領域での参照データ数の少なさにより、１回の計測から回転周期との位相や電源周期との位相がそれぞれ異なる１００〜４００程度の多様なサンプルを得ることができる。 For example, there are 100 data for 5 seconds with 102.4 kHz sampling. Due to the small number of reference data in the high frequency region of the constant Q conversion, various samples of about 100 to 400 having different phases from the rotation cycle and the power cycle can be obtained from one measurement.
次に、主成分分析部１０２において、多変量サンプルを主成分分析して変換行列を得るとともに、各サンプルの主成分得点を計算する。そして統計値計算部１０３が、主成分得点を基にホテリングＴ２／Ｑ統計量のようなサンプルの指標となる統計量を得る。
Next, the principal
最後に、算出した統計量が準備フェーズの時系列データで平均０になるように、第１の正規化部１０４において統計量を正規化して異常度とする。この際に、正規化に使用した正規化係数を保存して診断フェーズ（診断手段２００）で使う。尚、異常度の取り得る数値範囲を限定するために、異常度を統計量の対数値にすることが望ましい。
Finally, the
診断フェーズ（診断手段２００）では、新たに測定した時系列データを第２の定Ｑ変換装置２０１において定Ｑ変換して多変量サンプルを作成する。
In the diagnosis phase (diagnostic means 200), the newly measured time-series data is subjected to constant Q conversion in the second constant
定Ｑ変換装置２０１は、定Ｑ変換装置１０１と同様に、例えば非特許文献１〜４に開示されている方式を実行するものでもよく、また後述する実施例２〜４の高速演算が可能な定Ｑ変換装置を採用しても良い。
The constant
次に、多変量サンプルに準備フェーズの主成分分析部１０２で作成した変換行列を使って、主成分計算部２０２で主成分得点を得る。そしてこれを基に、第２の統計値計算部２０３が準備フェーズと同様の計算で統計量を得る。最後に準備フェーズの第１の正規化部１０４で作成した正規化係数を使って、第２の正規化部２０４が異常度を計算し、その値により機器が正常か異常かを診断する。
Next, a principal component score is obtained by the principal
異常診断に使う異常値の閾値を決めるのは一般的に困難であるが、異常値を統計値の対数で取る場合には、簡易の診断としては閾値をｌｏｇ３、ｌｏｇ５のような値としても良い。これは振動実効値での異常判断基準が通常時の３倍、５倍などとするのと同様である。 Although it is generally difficult to determine a threshold value of an abnormal value used for abnormality diagnosis, when the abnormal value is taken as a logarithm of a statistical value, the threshold value may be a value such as log 3 or log 5 for simple diagnosis. . This is the same as the case where the abnormality judgment standard based on the vibration effective value is 3 times, 5 times, or the like in the normal state.
本発明の効果を示すため、図７の第一主成分×第二主成分散布図と同じデータについて、定Ｑ変換による多変量化データを主成分分析した場合の第一主成分×第二主成分散布図を図２に示す。図２によれば、図７の場合と異なり、白色の丸印で示す学習期間のデータが、黒丸印で示す全データの分布範囲を概ねカバーできていることが分かる。 In order to show the effect of the present invention, for the same data as the first principal component × second principal component scatter diagram of FIG. The component scatter diagram is shown in FIG. According to FIG. 2, it can be seen that, unlike the case of FIG. 7, the data of the learning period indicated by white circles can generally cover the distribution range of all data indicated by black circles.
また、図６に示したものと同様の異常度を提案手法により計算した結果を図３に示す。図３によれば、図６の場合と異なり、時間がたっても異常度が上昇傾向にないことが分かる。線形回帰による傾きはほとんどゼロであり、正の傾きを取るＩｎｄｉｃａｔｏｒ１でも、異常度の平均が閾値に達するのに１０年ほど掛かる。
FIG. 3 shows the result of calculating the degree of abnormality similar to that shown in FIG. 6 by the proposed method. According to FIG. 3, it can be seen that, unlike the case of FIG. 6, the degree of abnormality does not tend to increase over time. The slope by linear regression is almost zero, and even with
次に、提案手法による異常診断例を図４に示す。図４は、図６や図３とは別の回転機の計算例である。こちらも同様に１時間に一回、約１００ｋＨｚのサンプリング周波数で５秒間の測定を行ったが、測定開始から２５日目に故障により緊急停止した。２５日目は稼動直後から異常度が非常に大きくなっているが、良く見ると、１３日目にＩｎｄｉｃａｔｏｒ１が注意ラインを越えていることが分かる。
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 those in FIGS. 6 and 3. Similarly, the measurement was performed once every hour at a sampling frequency of about 100 kHz for 5 seconds, but was stopped due to a failure on the 25th day from the start of the measurement. On the 25th day, the degree of abnormality is very large immediately after the operation, but if you look closely, you can see that on the 13th day,
尚、この診断結果はフーリエ変換による手法でもほぼ同様の結果が得られる。 This diagnosis result can be obtained by using a method using Fourier transform.
本実施例２は、図１の第１、第２の定Ｑ変換装置１０１、２０１を、より精度の高い定Ｑ変換を、少ない演算量で高速に行うことができる図５の定Ｑ変換装置で構成したものである。
In the second embodiment, the first and second constant
まず定Ｑ変換の成分について説明する。 First, the constant Q conversion component will be described.
非特許文献１によると、定Ｑ変換は次で定式化される（表記は少し変更している）。
According to
０以外の値は各ｋに対して、ｎ＝（Ｎ−Ｎk）／２，…，（Ｎ＋Ｎk）／２−１の範囲でのみ取り、この範囲の窓関数の形状として具体的にはハミング窓、ハニング窓、矩形窓など、多くの種類が提案されている
Ｎk：窓幅，Ｎk＝Ｑ（ｆs／ｆk）となるように取る
Ｎ≧Ｎ1≧…≧Ｎk…≧ＮK≧２Ｑとなる
Ｑ：窓の周期数，定Ｑ変換ではこれを全ての周波数で共通にする
非特許文献１、２ではＱ＝（２1/B−１）-1程度に取るようにしている
Ｂはビン数で、１オクターブを分割する数のこと
ｆs：サンプリング周波数
ｆk：第ｋ成分の周波数，ｆk＝（２1/B）k-1ｆmin
ｆmin＝ｆ1は定Ｑ変換で解析する最小の周波数で、解析する範囲から決める
ｆkは解析する最大の周波数で、ｆk≦ｆs／２となるようにする。
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: Number of window periods, common to all frequencies in constant Q conversion In
f min = f 1 is the minimum frequency to be analyzed by the constant Q conversion, and f k determined from the analysis range is the maximum frequency to be analyzed, so that f k ≦ f s / 2.
定Ｑ変換は特に低周波帯域の成分で窓幅Ｎkが非常に大きくなることにより演算量が多い。 The constant Q conversion has a large amount of calculation due to the extremely large window width N k especially in the low frequency band component.
これに対して非特許文献２では、高速フーリエ変換（ＦＦＴ；Ｆａｓｔ Ｆｏｕｒｉｅｒ ｔｒａｎｓｆｏｒｍ）を利用した演算の高速化方法が提案されている。 On the other hand, Non-Patent Document 2 proposes a method for speeding up the operation using Fast Fourier Transform (FFT).
この方法ではパーセバルの公式により、定Ｑ変換を次の数式（２）のように時間軸での計算から周波数軸での計算に置き換える。 In this method, the constant Q transformation is replaced with the calculation on the frequency axis from the calculation on the time axis as shown in the following formula (2) by the Parseval formula.
尚、明細書中の以下の文章では、前記数式（１）、（２）の各左辺の定Ｑ変換の成分を、「Ｘk cq」と表記することもある。 In the following sentences in the specification, the constant Q conversion component on each left side of the mathematical expressions (1) and (2) may be expressed as “X k cq ”.
ここでＳn,kは（小さなｋに対して特に）ほとんどの要素が０に近いので、一定値以下の項を０として疎行列表現すると（実質非ゼロ項の積和計算だけになるため）非常に高速に計算できるようになる。 Here, since most elements of Sn , k are close to 0 (especially for small k), expressing a sparse matrix with a term of a certain value or less as 0 (because it is only a product-sum calculation of non-zero terms) It becomes possible to calculate very fast.
さらに（１／Ｎ）Ｓn,kは、計算するデータ長と解析する周波数範囲が決まっていれば、入力データ列ｘnに寄らず事前に計算できるためこれを予め用意しておけば、入力データの離散フーリエ変換はＦＦＴで実行するとして、定Ｑ変換を入力データに対する１回のＦＦＴと、ＦＦＴ結果ベクトルと予め用意した疎行列との積演算により計算できることになる。 Furthermore, (1 / N) S n, k can be calculated in advance without depending on the input data string 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 executed by FFT, the constant Q transform can be calculated by a product operation of one FFT on the input data, an FFT result vector, and a sparse matrix prepared in advance.
離散でないフーリエ変換でも同等の公式が成り立ち、これらはパーセバルの定理、プランシュレルの定理など様々な呼び方がされる。 Equivalent formulas also hold for non-discrete Fourier transforms, and these are called variously, such as Perseval's theorem and Planschler's theorem.
非特許文献２による高速化のアイデアは周波数軸での係数Ｓn,kが疎行列とみなせることに依拠しているが、実際にはＳn,kが疎行列というのは粗い評価である。確かにその値は一部を除いて小さいが、精度を高めると必ずしも疎にならない。特に高周波帯域でこの傾向が顕著である。実際、非特許文献２、４で係数Ｓn,kが十分小さいとしている閾値はそれぞれ０．１５、０．００５４であり、通常の評価には問題ないが有効数字７桁以上の精度を得るには十分とは言えない。このため、要求精度によってはこの高速化手法がかえって遅くなる可能性がある。 The idea of speeding up by 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 reality it is a rough evaluation that S n, k is a sparse matrix. Certainly the value is small except for a part, but it does not necessarily become sparse when the accuracy is increased. This tendency is particularly remarkable in the high frequency band. Actually, the thresholds for which the coefficient Sn , k is sufficiently small in Non-Patent Documents 2 and 4 are 0.15 and 0.0054, respectively, and there is no problem in normal evaluation, but an accuracy of 7 significant digits or more is obtained. Is not enough. For this reason, depending on the required accuracy, there is a possibility that this speed-up method may become slower.
実際に周波数軸での計算でどの程度の演算量を必要とするかを調べるため、非特許文献４に準じた設定により時間軸での計算と周波数軸での計算でそれぞれ計算対象となる非ゼロの項数を評価する。 In order to investigate how much computation is actually required in the calculation on the frequency axis, the setting according to Non-Patent Document 4 is used, and the non-zero that is the calculation target for the calculation on the time axis and the calculation on the frequency axis, respectively. Evaluate the number of terms.
設定パラメータは、本明細書の記号法で表記してＢ＝２４、Ｑ＝２８、ｆs＝１６０００、ｆmin＝６０、Ｋ＝１６０とし、定Ｑ変換をそれぞれ矩形窓、ハミング窓、ハニング窓で計算する場合において、時間軸での計算に必要な非ゼロ項数と、周波数軸での計算で所定の精度の係数まで計算する際に必要な項数を表１に示す。 The setting parameters are expressed by the symbolic method of this specification, and B = 24, Q = 28, f s = 16000, f min = 60, K = 160, and the constant Q conversion is a rectangular window, a Hamming window, and a Hanning window, respectively. Table 1 shows the number of non-zero terms necessary for the calculation on the time axis and the number of terms necessary for calculating the coefficient with a predetermined accuracy by the calculation on the frequency axis.
なお、非特許文献１、２に合わせるとＢ＝２４のときＱ＝３４となるが、非特許文献３、４では追加の係数ｆｒａｔｉｏを導入してＱの値を変えているので、その設定値Ｑ＝２８を使っている。また解析する最大周波数６０００ＨｚからＫ＝１６０を導いた。ｆsの値は非特許文献３，４に明記されていなかったので、最大周波数６０００Ｈｚまで解析できるサンプリング周波数１２０００Ｈｚ以上のＷＡＶ音声データとしてｆs＝１６０００Ｈｚを選択した。
Note that Q = 34 when B = 24 when combined with
表１中の周波数軸右の数値は係数をゼロと評価する係数の大きさの閾値であり、各数値右のパーセンテージは矩形窓の時間軸計算の項数を基準としている。 The numerical value on the right side of the frequency axis in Table 1 is a threshold value of the coefficient size for evaluating 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.
表１から、窓形状により程度が違うが、精度が高いほど（係数をゼロと評価する係数Ｓn,kの大きさの閾値が小さいほど）周波数軸での計算に必要な項数が増えることが分かる。特に矩形窓やハミング窓では精度が高くなると時間軸での計算より周波数軸での計算の方が積和すべき項数が大きくなり、周波数軸での計算で高速にならなくなる。 From Table 1, although the degree varies depending on the window shape, the higher the accuracy (the smaller the threshold value of the coefficient Sn, k that evaluates the coefficient to zero), the more terms required for the calculation on the frequency axis. I understand. In particular, when the accuracy is high in a rectangular window or a Hamming window, the number of terms to be summed up in the frequency axis is larger than that in the time axis, and the calculation in the frequency axis does not become faster.
尚、周波数軸での計算にはＦＦＴ計算と疎行列の非ゼロ係数抽出のオーバーヘッドもあるので、必要項数に差がなければ時間軸計算が有利である。 Note that the calculation on the frequency axis also has the overhead of FFT calculation and non-zero coefficient extraction of the sparse matrix, so the time axis calculation is advantageous if there is no difference in the number of required terms.
より精度の高い定Ｑ変換を高速に行うには工夫が必要になる。 Ingenuity is required to perform constant Q conversion with higher accuracy at high speed.
そこで本実施例２では、定Ｑ変換を高精度でも演算量を少なく計算するために、定Ｑ変換の成分Ｘk cqを数式（２）の時間軸での計算と周波数軸での計算のどちらの方法で計算するか、予め周波数成分ｋ毎に演算量の小さい方を選択しておくことで実現する。 Therefore, in the second embodiment, in order to calculate the constant Q conversion with high accuracy even with a small amount of calculation, the constant X conversion component X k cq is calculated on the time axis or the frequency axis in equation (2). This can be realized by the above method or by selecting the smaller calculation amount for each frequency component k in advance.
定Ｑ変換の成分Ｘk cqの計算はおおよそ低周波帯域では周波数軸での計算の演算量が小さく、高周波帯域では時間軸での計算の演算量が小さいから、ｋの境界値ＫMを決めて、ｋ≦ＫMでは周波数軸での計算、ｋ＞ＫMでは時間軸での計算、というようにすれば良い。 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, 10 is a preparation unit for determining parameters for performing constant Q conversion in advance and preparing for constant Q conversion calculation, and 20 is an individual calculation unit for calculating constant Q conversion for individual time series data. is there.
準備部１０は、
設定した周波数成分ＫM以下の低周波数帯域において、時系列データの離散フーリエ変換と定Ｑ変換における周波数軸係数との積和計算により求めた定Ｑ変換成分の第１の演算量と、前記ＫMを超える高周波帯域において、時系列データと定Ｑ変換における時間軸係数との積和計算により求めた定Ｑ変換成分の第２の演算量とを含む合計演算量を算出し、前記ＫMを全周波数成分ｋに各々設定して各ＫM毎の前記合計演算量を求め、前記各ＫM毎の合計演算量を評価し、該演算量が最小となるＫMを、計算法境界を示す周波数成分閾値として決定する計算法境界決定手段の機能と、
定Ｑ変換における周波数成分ｋが前記決定された周波数成分閾値ＫM以下の低周波帯域では、時系列データの離散フーリエ変換と前記周波数軸係数との積和計算である第１の計算法を選択し、前記周波数成分ｋが前記決定された周波数成分閾値ＫMを超える高周波帯域では、時系列データと前記時間軸係数との積和計算である第２の計算法を選択する計算法選択手段の機能とを具備している。
The
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 calculation amount including the second calculation 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, and the above K M is calculated. 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 The function of the calculation method boundary determination means for determining the 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.
前記計算法境界決定手段および計算法選択手段の各機能を実行する準備部１０は、具体的には、定Ｑ変換のパラメータを決定するパラメータ決定部１１、時間軸係数Ｔn,kを計算する時間軸係数算出部１２、周波数軸係数（１／Ｎ）Ｓn,kを計算する周波数軸係数算出部１３、時間軸での計算と周波数軸での計算との境界のｋ値ＫMを決定する計算法境界決定部１４からなる。
Specifically, the
個別計算部２０は、時系列データに対して、前記計算法選択手段により選択された第１の計算法、第２の計算法を実施して定Ｑ変換の成分を求める定Ｑ変換手段の機能を具備している。 The individual calculation unit 20 performs the first calculation method and the second calculation method selected by the calculation method selection unit on the time series data, and functions as a constant Q conversion unit for obtaining a constant Q conversion component. It has.
前記定Ｑ変換手段の機能を実行する個別計算部２０は、具体的には、時系列データの長さを整える前処理部２１、高周波帯域の定Ｑ変換を計算する高周波帯域計算部２２、低周波帯域の定Ｑ変換を計算する低周波帯域計算部２３、高周波帯域での結果と低周波帯域での結果をまとめる定Ｑ変換合成部２４からなる。
Specifically, the individual calculation unit 20 that executes the function of the constant Q conversion means includes a
図５の定Ｑ変換の成分演算装置は、例えばコンピュータにより構成され、通常のコンピュータのハードウェアリソース、例えばＲＯＭ、ＲＡＭ、ＣＰＵ、入力装置、出力装置、通信インターフェース、ハードディスク、記録媒体およびその駆動装置を備えている。 5 is configured 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 driving device thereof It has.
このハードウェアリソースとソフトウェアリソース（ＯＳ、アプリケーションなど）との協働の結果、定Ｑ変換の成分演算装置は、図５に示すように、パラメータ決定部１１、時間軸係数算出部１２、周波数軸係数算出部１３、計算法境界決定部１４、前処理部２１、高周波帯域計算部２２、低周波帯域計算部２３、定Ｑ変換合成部２４を実装する。
As a result of the cooperation between the hardware resource and the software resource (OS, application, etc.), as shown in FIG. 5, the constant Q conversion component calculation device includes a
また各部での処理結果は、ハードディスク或いはＲＡＭなどの保存手段・記憶手段に格納され、随時利用されるものとする。 In addition, the processing results in each unit are stored in a storage unit / storage unit such as a hard disk or a RAM and used as needed.
まず、準備部１０のパラメータ決定部１１では、予め定Ｑ変換で解析する時系列データのサンプリング周波数ｆs、解析する周波数範囲（ｆmin，ｆmax）と１オクターブの周波数を区切るビン数Ｂ、解析対象データの周期数Ｑ、窓関数を決定する。このときｆmax≦ｆs／２が必要で、Ｋ＝ｆｌｏｏｒ（Ｂ＊ｌｏｇ2（ｆmax／ｆmin））＋１となる。これを基に計算する入力データの長さをＮ＝Ｎ1＝Ｑ（ｆs／ｆmin）とする。
First, in the
また周波数軸係数をゼロとみなすべき係数の大きさの閾値Ｔｈもここで決定しておく。 In addition, a threshold value Th for the magnitude of the coefficient that should be regarded as zero as the frequency axis coefficient is determined here.
時間軸係数算出部１２では、パラメータ決定部１１で決定したパラメータを基に、時間軸係数Ｔn,kを次の定義式に従って計算する。
The time axis
周波数軸係数算出部１３では、パラメータ決定部１１で決定したパラメータを基に、周波数軸係数（１／Ｎ）Ｓn,kを、
The frequency axis
の定義式に従って計算する。この際、｜Ｓn,k｜≦Ｔｈとなる項は周波数軸係数（１／Ｎ）Ｓn,k＝０とする。 Calculate according to the definition formula. In this case, a term satisfying | S n, k | ≦ Th is a frequency axis coefficient (1 / N) S n, k = 0.
計算法境界決定部１４では、周波数軸での計算を行う最大のｋ値ＫM（周波数成分閾値）を求める。ＫMは０，…，Ｋの何れかの値を持ち、ＫM＝０は周波数軸での計算を行わずに非特許文献１、２と同様に全て時間軸での計算で求める場合を意味する。
The calculation method
周波数成分閾値ＫMの決定方法としては、ＫM＝０，…，Ｋの各場合について、定Ｑ変換の成分Ｘk cqを求めるのに必要な演算量ＣX（ＫM）を評価してそれが最小となるようにＫMを決定する。 As a method of determining the frequency component threshold value K M , for each case of K M = 0,..., K, the amount of calculation C X (K M ) necessary for obtaining the constant Q conversion component X k cq is evaluated. K M is determined so that it is minimized.
定Ｑ変換の成分Ｘk cqの演算量ＣX（ＫM）（合計演算量）は、長さＮのＦＦＴの演算量ＣFFT（Ｎ）、非ゼロの時間軸係数毎の演算量ＣT、非ゼロの時間軸係数の数ＮT（ｋ）、非ゼロの周波数軸係数毎の演算量ＣF、ｋで非ゼロとみなす周波数軸係数の数ＮF（ｋ）を使っておおよそ The computation amount C X (K M ) (total computation amount) of the constant X conversion component X k cq is the FFT computation amount C FFT (N) of length N, and the computation amount C T for each non-zero time axis coefficient. , Using the number N T (k) of non-zero time axis coefficients, the calculation amount 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
と表せる。 It can be expressed.
ここでＣFFTは計算を実行する計算機やそれぞれの算法の詳細により変わり得るが、実測からの一例としては、ＣFFT（Ｎ）＝（１／４）Ｎｌｏｇ2（Ｎ）、ＣT＝１、ＣF＝１としても良い。 Here, C FFT may vary depending on the computer that performs the calculation and the details of each calculation method. As an example from actual measurement, C FFT (N) = (1/4) Nlog 2 (N), C T = 1, C F = 1 may be set.
上記演算量ＣX（ＫM）を示す数式（３）において、右辺第１項は、ＫMが非ゼロのときＣFFTとし、そうでないときは０とすることを表し、右辺第２項は、周波数成分ｋが周波数成分閾値ＫM以下の低周波帯域における周波数軸係数毎の演算量を表し、右辺第３項は、周波数成分ｋが周波数成分閾値ＫMを超える高周波帯域における時間軸係数毎の演算量を表している。 In the mathematical expression (3) indicating the calculation amount C X (K M ), the first term on the right side represents C FFT when K M is non-zero, and represents 0 when not, and the second term on the right side is 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 computation.
このためＫMが非ゼロのときは、数式（３）の右辺第１項のＦＦＴの演算量および右辺第２項の演算量（周波数軸係数による演算量）からなる第１の演算量と、数式（３）の右辺第３項の時間軸係数による第２の演算量との合計が演算量ＣX（ＫM）となる。 For this reason, when KM is non-zero, the first amount of computation composed of the amount of computation of the first term on the right side of Equation (3) and the amount of computation of the second term on the right side (the amount of computation based on the frequency axis coefficient); The sum of the second calculation amount based on the time axis coefficient in the third term on the right side of Equation (3) is the calculation amount C X (K M ).
またＫMがゼロのときは、数式（３）の右辺第１項および右辺第２項がともに０となって右辺第３項の時間軸係数による第２の演算量のみが演算量ＣX（ＫM）となる。 When KM is zero, the first term and the second term on the right side of Equation (3) are both 0, and only the second computation amount based on the time axis coefficient of the third term on the right side is the computation amount C X ( K M ).
この演算量ＣX（ＫM）は、非ゼロの時間軸係数の数ＮT（ｋ）がｋに関して指数的に減少し、非ゼロの周波数軸係数の数ＮF（ｋ）がｋに関しておおよそ指数的に増加するため、ＫM＝１，…，Ｋの範囲ではおおよそ下に凸の特性カーブとなる。ただし、ＫM＝０ではＣFFTがなくなるためＣX（０）が小さくなる。結果として演算量ＣX（ＫM）はＫM＝０あるいはＫM＝１，…，Ｋでの最小点の何れかで最小値を取ることになる。 This amount of computation C X (K M ) is exponentially decreased with respect to k as the number of non-zero time-axis coefficients N T (k), and the number N F (k) of non-zero frequency-axis coefficients as approximate with respect to k. Since it increases exponentially, it becomes a characteristic curve that is approximately downward in the range of K M = 1,. However, C X (0) becomes small because C FFT disappears when K M = 0. As a result, the calculation amount C X (K M ) takes a minimum value at any of the minimum points at K M = 0 or K M = 1,.
個別計算部２０では、入力された時系列データの定Ｑ変換を次のようにして求める。 The individual calculation unit 20 obtains the constant Q conversion of the input time series data as follows.
前処理部２１では、入力された時系列データの長さがＮより短ければ前後に０を補完して、Ｎより長ければ前後の値を除去して長さがＮの時系列データｘnになるようにする。このとき入力された時系列データの中央部分が時系列データｘnの中央部分に一致するようにする。
In the
次に定Ｑ変換の成分Ｘk cqを、前記計算法境界決定部１４で決定した周波数成分閾値ＫMにより、ｋ≦ＫMについては低周波帯域計算部２３で計算し、ｋ＞ＫMについては高周波帯域計算部２２で計算する。
Next, the constant Q conversion component X k cq is calculated by the low frequency
高周波帯域計算部２２は、前記時間軸係数算出部１２で算出された時間軸係数Ｔn,kと時系列データｘnを積和計算して定Ｑ変換の成分Ｘk cqを求める。各ｋについて非ゼロの時間軸係数Ｔn,kは高々Ｎk個なので、積和計算を非ゼロの範囲に限定することで効率的に計算できる。
The high frequency
低周波帯域計算部２３は、まず時系列データｘnをＦＦＴしてフーリエ係数Ｘnを求める。次に前記周波数軸係数算出部１３で算出された周波数軸係数（１／Ｎ）Ｓn,kとフーリエ係数Ｘnを積和計算して定Ｑ変換の成分Ｘk cqを求める。各ｋについて非ゼロの周波数軸係数（１／Ｎ）Ｓn,kは少数なので、積和計算を非ゼロの範囲に限定することで効率的に計算できる。
The low frequency
尚、周波数軸係数は多くの場合、両端（周波数ゼロ近辺とサンプリング周波数近辺）部分では非ゼロ項とゼロ項が入り混じるが中間部分ではゼロ項が続くので、非ゼロ項の範囲を少し広げて、両端部分は全て非ゼロ項として計算しても良い。この場合、非ゼロ項が少し増えるがその数はそれほど多くはなく、また非ゼロ判定ロジックが簡略化できるので演算時間は大差なくなる。 In many cases, the frequency axis coefficient is mixed with non-zero terms and zero terms at both ends (near the frequency zero and near the sampling frequency), but the zero term continues in the middle, so the range of the non-zero terms is slightly expanded. Both end portions may be calculated as non-zero terms. In this case, the number of non-zero terms increases slightly, but the number thereof is not so large, and the non-zero determination logic can be simplified, so that the calculation time does not differ greatly.
最後に定Ｑ変換合成部２４で、高周波帯域計算部２２と低周波帯域計算部２３とで計算した定Ｑ変換の成分Ｘk cqをまとめて、定Ｑ変換を完了する。
Finally, the constant Q
上記のように本実施例２によれば、従来の高速化技術の問題点である、高い精度の定Ｑ変換を計算しようとすると演算量が膨れ上り通常の時間軸での計算よりもむしろ遅くなるという問題が改善される。 As described above, according to the second embodiment, when calculating a high-accuracy constant Q conversion, which is a problem of the conventional high-speed technology, the amount of calculation increases and is slower than the calculation on the normal time axis. The problem of becoming is improved.
その効果を示すため、定Ｑ変換の計算方法（時間軸計算、周波数軸計算、提案手法）毎に、必要な演算量を表２のように試算した。試算は、時系列データのパラメータを表１と同じにして、演算量パラメータをＣFFT（Ｎ）＝（１／４）Ｎｌｏｇ2（Ｎ）、ＣT＝１、ＣF＝１として計算した定Ｑ変換の成分Ｘk cqの演算量ＣX（ＫM）を表２に示す。 In order to show the effect, the necessary amount of calculation was calculated as shown in Table 2 for each constant Q conversion calculation method (time axis calculation, frequency axis calculation, proposed method). In the trial calculation, the parameters of the time series data are the same as in Table 1, and the calculation amount parameters are C FFT (N) = (1/4) Nlog 2 (N), C T = 1, 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 coefficient magnitude for evaluating 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 appears slightly due to the amount of FFT calculation, the Hanning window has an amount of calculation of approximately 1/5 (22.5%) even at the threshold of 1.0e-5. Yes. Even if the calculation amount of the FFT operation is heavy and the effect on the frequency axis is not effective, the calculation on the time axis is automatically performed, so that the problem of slowing down does not occur.
実施例２は、ある一時刻における定Ｑ変換の各成分値を求める際には効率的であるが、実際の定Ｑ変換では、時系列データの周波数成分の時間変化を捉えるために少しずつ時刻をずらしながら繰り返し定Ｑ成分を計算する。このときずらす時刻は、最短ではサンプリング周波数で、時系列データの情報を余すことなく計算するには最長でも最高周波数での窓幅にしたい。もちろん、時系列データの特徴サンプル抽出などでは、これよりも長い間隔で時刻をずらして定Ｑ変換する場合もあるが、いずれにせよ時系列データ全体に対しては多くの繰り返し計算をすることになり、まだまだ負荷が高い。 The second embodiment is efficient in obtaining each component value of the constant Q conversion at a certain time, but in the actual constant Q conversion, the time is gradually measured in order to capture the time change of the frequency component of the time series data. The constant Q component is calculated repeatedly while shifting. The time to be shifted at this time is the sampling frequency at the shortest, and it is desired to set the window width at the maximum frequency at the longest to calculate the time series data without leaving any information. Of course, when extracting feature samples of time-series data, there is a case where constant Q conversion is performed by shifting the time at intervals longer than this, but in any case, many repetitive calculations are performed on the entire time-series data. The load is still high.
多数の時刻で同時に定Ｑ変換を計算する際には更なる効率化が必要である。 When calculating the constant Q conversion at a large number of times at the same time, further efficiency is required.
そこで本実施例３では、数式（２）の周波数軸での計算で定Ｑ変換を行う場合に、定Ｑ変換を計算する時刻ごとに計算に使用する時系列データの範囲でＤＦＴ（あるいはＦＦＴ）する代わりに、複数時刻での定Ｑ変換に必要な範囲全体で一括してＤＦＴ（あるいはＦＦＴ）を行うことで高速化を図った。 Therefore, in the third embodiment, when constant Q conversion is performed in the calculation on the frequency axis of Equation (2), DFT (or FFT) is performed within the range of time series data used for calculation at each time when constant Q conversion is calculated. Instead, the DFT (or FFT) is performed collectively over the entire range necessary for the constant Q conversion at a plurality of times to increase the speed.
数式（１）で定義する定Ｑ変換はデータ列ｘn，ｎ＝０，…，Ｎ−１の中央ｘN/2での成分Ｘk cqを計算するように定式化しているが、窓関数の非ゼロ範囲を時刻方向に平行移動することで他の時刻での成分を計算するように調整できる。すなわちＴn,kとｘnとの積和のインデックスをずらすことで実現できる。 The constant Q transformation defined by Equation (1) is formulated so as to calculate the component X k cq at the center x N / 2 of the data string x n , n = 0 ,. It is possible to adjust so as to calculate the components at other times by translating the non-zero range in the time direction. That is , it can be realized by shifting the product sum index of T n, k and x n .
数式（２）による周波数軸での計算に関しては、初めに定Ｑ変換する時刻全体で必要な全データに対して一括してＤＦＴ（あるいはＦＦＴ）を行い、予め各時刻用に用意した（１／Ｎ）Ｓn,kと順次積和することで実現する。 With respect to the calculation on the frequency axis according to the equation (2), DFT (or FFT) is performed on all the data necessary for the entire time for the first constant Q conversion, and prepared for each time in advance (1 / N) Realized by sequentially multiplying and summing with Sn , k .
本実施例３による定Ｑ変換の成分演算装置の具体的な計算の構成は、図５と同様であるが、各部で実行される機能が以下に示すように異なる。 The specific calculation configuration of the constant Q conversion component arithmetic apparatus according to the third embodiment is the same as that shown in FIG. 5, but the functions executed in each unit are different as shown below.
まず、準備部１０のパラメータ決定部１１では、予め定Ｑ変換で解析する時系列データのサンプリング周波数ｆｓ、解析する周波数範囲（ｆmin，ｆmax）と１オクターブの周波数を区切るビン数Ｂ、解析対象データの周期数Ｑ、窓関数を決定する。このときｆmax≦ｆs／２が必要で、Ｋ＝ｆｌｏｏｒ（Ｂ＊ｌｏｇ2（ｆmax／ｆmin））＋１となる。
First, the
さらに入力データｘnの全長Ｎと、解析位置ｐm，ｍ＝１，…，Ｍ（時間方向の解析位置；時刻）を決定しておく。このとき、各位置ｐmでの定Ｑ変換はｘpmを中心とする長さＮ1＝Ｑ（ｆs／ｆmin）のデータ列をもとに計算し、（Ｎ1／２）≦ｐ1≦…≦ｐM≦Ｎ−（Ｎ1／２）であるものとする。すなわち、全ての解析位置ｐm，ｍ＝１，…，Ｍで定Ｑ変換の計算に必要なデータがｘ0，…，ｘN-1に含まれるものとする。尚、入力データの全点で定Ｑ変換を計算したいような場合は予め入力データを前後にゼロ拡張して上記条件を満たすように変換しておけばよい。
Further, the total length N of the input data x n and the analysis positions p m , m = 1,..., M (analysis position in the time direction; time) are determined. At this time, the constant Q transformation at each position p m is computed based on the data sequence of length N 1 = Q centered at x pm (f s / f min ), (
また、周波数軸係数をゼロとみなすべき係数の大きさの閾値Ｔｈもここで決定しておく。 In addition, a threshold value Th of the coefficient magnitude that should be regarded as a frequency axis coefficient is also determined here.
時間軸係数算出部１２ではパラメータ決定部１１で決定したパラメータを基に、時間軸係数Ｔn,kを各解析位置ｐmについて次の定義式に従って計算する。
Parameter based on the determined time axis
このときｐmに対するＴn,kをＴn,k,mと書くと、非ゼロ項に関してＴn,k,m＝Ｔn-pm+p1,k,1である。 T n for the time p m, k and T n, k, when written as m, T n with respect to non-zero terms, k, a m = T n-pm + p1 , k, 1.
周波数軸係数算出部１３ではパラメータ決定部１１で決定したパラメータを基に、各解析位置ｐmについて周波数軸係数（１／Ｎ）Ｓn,k,mを定義式に従って計算する。この際、｜Ｓn,k,m｜≦Ｔｈとなる項は周波数軸係数（１／Ｎ）Ｓn,k,m＝０とする。
Based on the parameters determined in the frequency axis
計算法境界決定部１４では、周波数軸での計算を行う最大のｋ値ＫM（周波数成分閾値）を求める。ＫMは０，…，Ｋの何れかの値を持ち、ＫM＝０は周波数軸での計算を行わずに非特許文献１、２と同様に全て時間軸での計算で求める場合を意味する。
The calculation method
周波数成分閾値ＫMの決定方法としては、ＫM＝０，…，Ｋの各場合について、定Ｑ変換の成分Ｘk cqを求めるのに必要な演算量ＣX（ＫM）を評価してそれが最小となるようにＫMを決定する。尚、この条件は複数のｐmに寄らないので、ｐ1（１つの解析位置）について求めておけば十分である。 As a method of determining the frequency component threshold value K M , for each case of K M = 0,..., K, the amount of calculation C X (K M ) necessary for obtaining the constant Q conversion component X k cq is evaluated. K M is determined 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 (1 single analysis position).
定Ｑ変換の成分Ｘk cqの演算量ＣX（ＫM）（合計演算量）は、長さＮのＦＦＴの演算量ＣFFT（Ｎ）、非ゼロの時間軸係数毎の演算量ＣT、非ゼロの時間軸係数の数ＮT（ｋ）、非ゼロの周波数軸係数毎の演算量ＣF、ｋで非ゼロとみなす周波数軸係数の数ＮF（ｋ）を使っておおよそ、前記実施例２で述べた数式（３） The computation amount C X (K M ) (total computation amount) of the constant X conversion component X k cq is the FFT computation amount C FFT (N) of length N, and the computation amount C T for each non-zero time axis coefficient. , Using the number N T (k) of non-zero time axis coefficients, the calculation amount 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, Formula (3) described in the second embodiment
と表せる。 It can be expressed.
ここでＣFFTは計算を実行する計算機やそれぞれの算法の詳細により変わり得るが、実測からの一例としては、ＣFFT（Ｎ）＝（１／４）Ｎｌｏｇ2（Ｎ）、ＣT＝１、ＣF＝１としても良い。 Here, C FFT may vary depending on the computer that performs the calculation and the details of each calculation method. As an example from actual measurement, C FFT (N) = (1/4) Nlog 2 (N), C T = 1, C F = 1 may be set.
この演算量ＣX（ＫM）は、時間軸係数および周波数軸係数の、ハニング窓での項数と周波数（周波数成分ｋ；対数周波数）の関係が、非ゼロの時間軸係数の数ＮT（ｋ）がｋに関して指数的に減少し、非ゼロの周波数軸係数の数ＮF（ｋ）がｋに関しておおよそ指数的に増加するため、ＫM＝１，…，Ｋの範囲ではおおよそ下に凸の特性カーブとなる。ただし、ＫM＝０ではＣFFTがなくなるためＣX（０）が小さくなる。結果として演算量ＣX（ＫM）はＫM＝０あるいはＫM＝１，…，Ｋでの最小点の何れかで最小値を取ることになる。 This calculation amount C X (K M ) is the number N T of time axis coefficients whose 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 N F (k) of non-zero frequency axis coefficients increases approximately exponentially with respect to k, it is approximately lower in the range of K M = 1,. It becomes a convex characteristic curve. However, C X (0) becomes small because C FFT disappears when K M = 0. As a result, the calculation amount C X (K M ) takes a minimum value at any of the minimum points at K M = 0 or K M = 1,.
個別計算部２０では入力された時系列の定Ｑ変換を次のようにして求める。 The individual calculation unit 20 obtains the input time-series constant Q conversion as follows.
前処理部２１では、入力された時系列データの長さがＮより短ければ前後にゼロを補完して、Ｎより長ければ前後の値を除去して長さがＮの時系列データｘnになるようにする。ゼロを補完する場合には、入力された時系列データの中央部分が時系列データｘnの中央部分に一致するようにする。データを除去する場合は、解析したい範囲に寄るが、単純に後ろのデータを除去するのでも良い。
In the
次に定Ｑ変換の成分Ｘk cqを、計算法境界決定部１４で決定したＫMにより、ｋ≦ＫMについては低周波帯域計算部２３で、ｋ＞ＫMについては高周波帯域計算部２２で計算する。
Then component X k cq constant Q transform, the K M determined by calculation method
尚、前記定Ｑ変換の成分は以下、「Ｘk,m cq」と表記することもある。 The constant Q conversion component may also be expressed as “X k, m cq ” below.
各ｋについて、前記時間軸係数算出部１２で算出された非ゼロの時間軸係数Ｔn,k,mは高々Ｎk個なので、積和計算を非ゼロの範囲に限定することで効率的に計算できる。
For each k, there are at most N k non-zero time-axis coefficients T n, k, m calculated by the time-axis
低周波帯域計算部２３は、まず時系列データｘnを全体に対して一度ＦＦＴしてフーリエ係数Ｘnを求める。次にｐm毎に、前記周波数軸係数算出部１３で算出された周波数軸係数（１／Ｎ）Ｓn,k,mとフーリエ係数Ｘnを積和計算して定Ｑ変換の成分Ｘk,m cqを求める。
Low frequency
各ｋについて非ゼロの周波数軸係数（１／Ｎ）Ｓn,k,mは少数なので、積和計算を非ゼロの範囲に限定することで効率的に計算できる。 Since there are a small number of non-zero frequency axis coefficients (1 / N) S n, k, m for each k, it can be efficiently calculated by limiting the product-sum calculation to a non-zero range.
尚、周波数軸係数は多くの場合、両端（周波数ゼロ近辺とサンプリング周波数近辺）部分では非ゼロ項とゼロ項が入り混じるが中間部分ではゼロ項が続くので、非ゼロ項の範囲を少し広げて、両端部分は全て非ゼロ項として計算しても良い。この場合、非ゼロ項が少し増えるがその数はそれほど多くはなく、また非ゼロ判定ロジックが簡略化できるので演算時間は大差なくなる。 In many cases, the frequency axis coefficient is mixed with non-zero terms and zero terms at both ends (near the frequency zero and near the sampling frequency), but the zero term continues in the middle, so the range of the non-zero terms is slightly expanded. Both end portions may be calculated as non-zero terms. In this case, the number of non-zero terms increases slightly, but the number thereof is not so large, and the non-zero determination logic can be simplified, so that the calculation time does not differ greatly.
最後に定Ｑ変換合成部２４で、高周波帯域計算部２２と低周波帯域計算部２３とで計算した定Ｑ変換の成分Ｘk,m cqをまとめて、定Ｑ変換を完了する。
Finally, the constant Q
前記実施例３では、時間軸係数Ｔn,k,m、周波数軸係数（１／Ｎ）Ｓn,k,mは解析位置ｐm毎に用意するので、事前に用意して保管しておくデータ量が大きくなる。本実施例４では、係数データのメモリ使用量を大幅に減らすように以下のような対策を講じた。 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 amount of data increases. In the fourth embodiment, the following measures are taken so as to greatly reduce the memory usage of coefficient data.
尚、数式（４）中の In addition, in Formula (4)
は、複素指数関数の時刻および周波数毎に異なるべき乗を表している。 Represents a different power for each time and frequency of the complex exponential function.
本実施例４による定Ｑ変換の成分演算装置の具体的な構成は図５と同様であり、以下に図５と異なる部分のみを説明する。 The specific configuration of the constant Q conversion component arithmetic apparatus according to the fourth embodiment is the same as that shown in FIG. 5, and only the parts different from FIG. 5 will be described below.
時間軸係数算出部１２では、パラメータ決定部１１で決定したパラメータを基に、時間軸係数Ｔn,k,1を定義式に従って求める。すなわち、
The time axis
を、予め設定した解析位置ｐmのうち１つの解析位置（例えば１番目の解析位置）ｐ1において算出し、時間軸係数Ｔn,k,1を求める。 And calculated in one analysis position (e.g. the first analysis position) p 1 of the preset analysis position p m, the time axis coefficient T n, obtains the k, 1.
周波数軸係数算出部１３では、パラメータ決定部１１で決定したパラメータを基に、周波数軸係数（１／Ｎ）Ｓn,k,1を定義式に従って計算する。
The frequency axis
すなわち、 That is,
から得た（１／Ｎ）Ｓn,kを、前記１つの解析位置ｐ1において算出して周波数軸係数（１／Ｎ）Ｓn,k,1を求める。 (1 / N) S n, k obtained from (1) is calculated at the one analysis position p 1 to obtain the frequency axis coefficient (1 / N) S n, k, 1 .
この際、｜Ｓn,k,1｜≦Ｔｈとなる項は周波数軸係数（１／Ｎ）Ｓn,k,1＝０とする。 In this case, a term satisfying | S n, k, 1 | ≦ Th is a frequency axis coefficient (1 / N) S n, k, 1 = 0.
計算法境界決定部１４では、周波数軸計算の際に複素べき乗の積和が追加されるために非ゼロの周波数軸係数毎の演算量ＣFを実施例３の場合より大きく、例えば５（複素べき乗の演算に＋３、複素べき乗の積和に＋１）とする。
The calculation method
高周波帯域計算部２２は、複数の解析位置ｐm毎に範囲をずらして、前記時間軸係数算出部１２で算出された時間軸係数Ｔn,k,1と時系列データｘn-p1+pmを積和計算して定Ｑ変換の成分Ｘk,m cqを求める。
High frequency
低周波帯域計算部２３は、まず時系列データｘnを全体に対して一度ＦＦＴしてフーリエ係数Ｘnを求める。
Low frequency
前記実施例３、４の手法によれば、定Ｑ変換を高速に実行するにあたって、先行の高速化技術（実施例２）が時間軸での計算でも周波数軸での計算でもそれらを併用する計算方法でも、一時刻毎に定Ｑ変換を計算するのに対して、複数時刻での定Ｑ変換を同時に計算する際に、ＤＦＴ（あるいはＦＦＴ）を一時刻の計算範囲毎ではなく全体に対して一度だけ実施することで、全体としての計算が高速化される。 According to the methods of the third and fourth embodiments, when the constant Q conversion is executed at high speed, the preceding high-speed technique (second embodiment) uses both the time-axis calculation and the frequency-axis calculation. In the method, the constant Q conversion is calculated every hour, but when calculating the constant Q conversion at a plurality of times at the same time, the DFT (or FFT) is not calculated for every calculation range of one time. Performing only once speeds up the overall calculation.
実施例４では、時間軸と周波数軸との係数データを計算時刻毎に持つ代わりに、それら複素数のべき乗計算に置き換えて（数式（４））、演算量は少し増える代わりに係数データのメモリ使用量を大幅に減らすことができる。 In the fourth embodiment, instead of having the coefficient data of the time axis and the frequency axis for each calculation time, the complex data is replaced with a power of the complex number (Equation (4)), and the amount of calculation is increased slightly. The amount can be greatly reduced.
１０…準備部
１１…パラメータ決定部
１２…時間軸係数算出部
１３…周波数軸係数算出部
１４…計算法境界決定部
２０…個別計算部
２１…前処理部
２２…高周波帯域計算部
２３…低周波帯域計算部
２４…定Ｑ変換合成部
１００…診断用パラメータ作成手段
１０１…第１の定Ｑ変換装置
１０２…主成分分析部
１０３…第１の統計値計算部
１０４…第１の正規化部
２００…診断手段
２０１…第２の定Ｑ変換装置
２０２…主成分計算部
２０３…第２の統計値計算部
２０４…第２の正規化部
DESCRIPTION OF
Claims (8)
新たに異常診断対象から測定した時系列データを定Ｑ変換して多変量サンプルを作成し、前記診断用パラメータ作成手段により作成された診断用パラメータを用いて異常を診断する診断手段と、を備えたことを特徴とする異常診断装置。 A diagnostic parameter creating 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;
A diagnostic means for diagnosing an abnormality using the diagnostic parameter created by the diagnostic parameter creation means by creating a multivariate sample by performing constant Q conversion on the time series data newly measured from the abnormality diagnostic target; An abnormality diagnosis apparatus characterized by that.
時系列データを定Ｑ変換して多変量サンプルを作成する第１の定Ｑ変換装置と、
前記作成された多変量サンプルを主成分分析して変換行列を得るとともに各サンプルの主成分得点を計算する主成分分析部と、
前記主成分得点を基にサンプルの指標となる統計量を得る第１の統計値計算部と、
前記統計量を正規化して異常度とする第１の正規化部と、を備え、
前記診断手段は、
時系列データを定Ｑ変換して多変量サンプルを作成する第２の定Ｑ変換装置と、
前記作成された多変量サンプルから、前記主成分分析部で得られた変換行列を使って主成分得点を得る主成分計算部と、
前記主成分得点を基に統計量を得る第２の統計値計算部と、
前記統計量を、前記第１の正規化部で使用した正規化係数によって正規化して異常度を計算する第２の正規化部と、を備え、前記異常度に基づいて異常診断対象の異常を診断することを特徴とする請求項１に記載の異常診断装置。 The diagnostic parameter creation means includes
A first constant Q conversion device for generating multivariate samples by performing constant Q conversion on time series data;
A principal component analysis unit that obtains a transformation matrix by performing principal component analysis on the created multivariate sample and calculates a principal component score of each sample;
A first statistic calculator that obtains a statistic serving as a sample index based on the principal component score;
A first normalization unit that normalizes the statistics and sets the degree of abnormality,
The diagnostic means includes
A second constant Q conversion device for generating multivariate samples by performing constant Q conversion on time series data;
A principal component calculation unit that obtains a principal component score from the created multivariate sample using the transformation matrix obtained by the principal component analysis unit;
A second statistical value calculation unit for obtaining a statistic based on the principal component score;
A second normalization unit that calculates the degree of abnormality by normalizing the statistic with the normalization coefficient used in the first normalization unit, and detects an abnormality of the abnormality diagnosis target based on the degree of abnormality The abnormality diagnosis apparatus according to claim 1, wherein diagnosis is performed.
設定した周波数成分ＫM以下の低周波数帯域において、時系列データの離散フーリエ変換と定Ｑ変換における周波数軸係数との積和計算により求めた定Ｑ変換成分の第１の演算量と、前記ＫMを超える高周波帯域において、時系列データと定Ｑ変換における時間軸係数との積和計算により求めた定Ｑ変換成分の第２の演算量とを含む合計演算量を算出し、
前記ＫMを全周波数成分Ｋに各々設定して各ＫM毎の前記合計演算量を求め、
前記各ＫM毎の合計演算量を評価し、該演算量が最小となるＫMを、計算法境界を示す周波数成分閾値として決定する計算法境界決定手段と、
定Ｑ変換における周波数成分Ｋが前記決定された周波数成分閾値ＫM以下の低周波帯域では、時系列データの離散フーリエ変換と前記周波数軸係数との積和計算である第１の計算法を選択し、前記周波数成分Ｋが前記決定された周波数成分閾値ＫMを超える高周波帯域では、時系列データと前記時間軸係数との積和計算である第２の計算法を選択する計算法選択手段と、
時系列データに対して、前記計算法選択手段により選択された第１の計算法、第２の計算法を実施して定Ｑ変換の成分を求める定Ｑ変換手段と、
を備えたことを特徴とする請求項２に記載の異常診断装置。 The diagnostic parameter creation means and the first and second constant Q conversion devices of the diagnostic means 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 calculation amount including the second calculation 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,
To evaluate the total amount of computation of the respective K M, the K M of the calculation amount is minimized, and the calculation method boundary determining means for determining a frequency component threshold indicating the 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 ,
A 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;
The abnormality diagnosis device according to claim 2, further comprising:
前記定Ｑ変換における時間軸係数を、
前記定Ｑ変換における周波数軸係数を、時間軸係数Ｔn,kをｎ＝０，…，Ｎ−１の時系列データとみたときの離散フーリエ変換である
前記定Ｑ変換手段は、
前記周波数軸係数算出部で算出された周波数軸係数と時系列データの離散フーリエ変換とを積和計算して前記第１の計算法を実施する低周波帯域計算部と、
前記時間軸係数算出部で算出された時間軸係数と時系列データとを積和計算して前記第２の計算法を実施する高周波帯域計算部と、
前記低周波帯域計算部と高周波帯域計算部の各出力を合成して定Ｑ変換の成分を求める合成部と、
を備えていることを特徴とする請求項３に記載の異常診断装置。 Each of the constant Q conversion devices is
The time axis coefficient in the constant Q conversion is
The frequency axis coefficient in the constant Q conversion is a discrete Fourier transform when the time axis coefficient T n, k is regarded as time series data of n = 0,..., N−1.
The constant Q conversion means includes:
A low frequency band calculation unit for performing the first calculation method by multiplying and multiplying the frequency axis coefficient calculated by the frequency axis coefficient calculation unit and the discrete Fourier transform of the time series data;
A high-frequency band calculation unit that performs a product-sum calculation on the time axis coefficient calculated by the time axis coefficient calculation unit and time-series data and implements the second calculation method;
A synthesis unit for synthesizing outputs of the low frequency band calculation unit and the high frequency band calculation unit to obtain a constant Q conversion component;
The abnormality diagnosis device according to claim 3, comprising:
時間軸係数Ｔn,kをｎ＝０，…，Ｎ−１の時系列データとみたときの離散フーリエ変換である
前記定Ｑ変換手段は、
時系列データｘn全体に対して離散フーリエ変換を行ってフーリエ係数Ｘnを求め、前記解析位置ｐm毎に、前記周波数軸係数算出部で算出された周波数軸係数（１／Ｎ）Ｓn,k,mと前記フーリエ係数Ｘnとを積和計算して前記第１の計算法を実施する低周波帯域計算部と、
前記解析位置ｐm毎に、前記時間軸係数算出部で算出された時間軸係数Ｔn,k,mと時系列データｘnとを積和計算して前記第２の計算法を実施する高周波帯域計算部と、
前記低周波帯域計算部と高周波帯域計算部の各出力を合成して定Ｑ変換の成分を求める合成部と、
を備えていることを特徴とする請求項３に記載の異常診断装置。 Each of the constant Q conversion devices is
Discrete Fourier transform when the time axis coefficient T n, k is regarded as time series data of n = 0,..., N−1.
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, and a low frequency band calculation unit for performing the first calculation method by calculating the sum of products,
For each of the analyzed position p m, implementing the time-axis coefficient calculating unit time axis coefficients calculated in T n, k, m and the time-series data x n and then product sum calculation of the second calculation method RF A bandwidth calculator;
A synthesis unit for synthesizing outputs of the low frequency band calculation unit and the high frequency band calculation unit to obtain a constant Q conversion component;
The abnormality diagnosis device according to claim 3, comprising:
時間軸係数Ｔn,kをｎ＝０，…，Ｎ−１の時系列データとみたときの離散フーリエ変換である
前記定Ｑ変換手段は、
時系列データｘn全体に対して離散フーリエ変換を行ってフーリエ係数Ｘnを求め、 前記解析位置ｐm毎に、
前記周波数軸係数算出部で求められた周波数軸係数（１／Ｎ）Ｓn,k,1と、
前記解析位置ｐm毎の周波数軸係数（１／Ｎ）Ｓn,k,mは前記１つの解析位置ｐ1における周波数軸係数（１／Ｎ）Ｓn,k,1に帰結できる関係にあることから導かれる、複素指数関数のべき乗
前記フーリエ係数Ｘnとを積和計算して前記第１の計算法を実施する低周波帯域計算部と、
前記解析位置ｐm毎に範囲をずらして、前記時間軸係数算出部で算出された時間軸係数Ｔn,k,1と時系列データｘn-p1+pmとを積和計算して前記第２の計算法を実施する高周波帯域計算部と、
前記低周波帯域計算部と高周波帯域計算部の各出力を合成して定Ｑ変換の成分を求める合成部と、
を備えていることを特徴とする請求項３に記載の異常診断装置。 Each of the constant Q conversion devices is
Discrete Fourier transform when the time axis coefficient T n, k is regarded as time series data of n = 0,..., N−1.
The constant Q conversion means includes:
Calculated Fourier coefficients X n by performing a discrete Fourier transform for the entire time-series data x n, for each of the analyzed position p m,
The frequency axis coefficient (1 / N) S n, k, 1 obtained by the frequency axis coefficient calculating 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 the complex exponential function
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 for performing the calculation method of 2,
A synthesis unit for synthesizing outputs of the low frequency band calculation unit and the high frequency band calculation unit to obtain a constant Q conversion component;
The abnormality diagnosis device according to claim 3, 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,
A diagnostic step in which a diagnostic means creates a multivariate sample by performing constant Q conversion on time series data newly measured from an abnormality diagnosis target, and diagnoses an abnormality using the diagnostic parameters created by the diagnostic parameter creation means And an abnormality diagnosis method.
第１の定Ｑ変換装置が、時系列データを定Ｑ変換して多変量サンプルを作成するステップと、
主成分分析部が、前記作成された多変量サンプルを主成分分析して変換行列を得るとともに各サンプルの主成分得点を計算するステップと、
第１の統計値計算部が、前記主成分得点を基にサンプルの指標となる統計量を得るステップと、
第１の正規化部が、前記統計量を正規化して異常度とするステップと、を備え、
前記診断ステップは、
第２の定Ｑ変換装置が、時系列データを定Ｑ変換して多変量サンプルを作成するステップと、
主成分計算部が、前記作成された多変量サンプルから、前記主成分分析部で得られた変換行列を使って主成分得点を得るステップと、
第２の統計値計算部が、前記主成分得点を基に統計量を得るステップと、
第２の正規化部が、前記統計量を、前記第１の正規化部で使用した正規化係数によって正規化して異常度を計算するステップと、を備え、前記異常度に基づいて異常診断対象の異常を診断することを特徴とする請求項７に記載の異常診断方法。 The diagnostic parameter creation step includes:
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 obtains a transformation matrix by performing principal component analysis on the created multivariate sample and calculates a principal component score of each sample;
A first statistic calculation unit obtaining a statistic serving as an index of a sample based on the principal component score;
A first normalization unit comprising: normalizing the statistic to obtain an abnormality level;
The diagnostic 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 obtains a principal component score from the created multivariate sample using the transformation matrix obtained by the principal component analysis unit;
A second statistic calculator calculates a statistic based on the principal component score;
A second normalization unit comprising: calculating the degree of abnormality by normalizing the statistic with the normalization coefficient used in the first normalization unit; and an abnormality diagnosis target based on the degree of abnormality The abnormality diagnosis method according to claim 7, wherein the abnormality is diagnosed.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| JP2016091695A JP6627639B2 (en) | 2016-04-28 | 2016-04-28 | Abnormality diagnosis device and abnormality diagnosis method |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| JP2016091695A JP6627639B2 (en) | 2016-04-28 | 2016-04-28 | Abnormality diagnosis device and abnormality diagnosis method |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| JP2017198620A true JP2017198620A (en) | 2017-11-02 |
| JP6627639B2 JP6627639B2 (en) | 2020-01-08 |
Family
ID=60237779
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| JP2016091695A Active JP6627639B2 (en) | 2016-04-28 | 2016-04-28 | Abnormality diagnosis device and abnormality diagnosis method |
Country Status (1)
| Country | Link |
|---|---|
| JP (1) | JP6627639B2 (en) |
Cited By (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JP2019101777A (en) * | 2017-12-04 | 2019-06-24 | 株式会社明電舎 | Abnormality diagnosing apparatus and abnormality diagnosing method |
| JP2020079718A (en) * | 2018-11-12 | 2020-05-28 | 東ソー株式会社 | Joint failure detection method of metal-resin composite |
| JP2020181443A (en) * | 2019-04-26 | 2020-11-05 | 株式会社豊田中央研究所 | Abnormality detection apparatus, abnormality detection method, and computer program |
| CN113557414A (en) * | 2019-03-22 | 2021-10-26 | Abb瑞士股份有限公司 | Device for monitoring equipment |
Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JPH07219581A (en) * | 1994-01-31 | 1995-08-18 | Babcock Hitachi Kk | Discrimination method for acoustic signal and device therefor |
| JP2004020193A (en) * | 2002-06-12 | 2004-01-22 | Takayoshi Yamamoto | Method for diagnosing object facility, computer program, and device for diagnosing subject facility |
| JP2008502927A (en) * | 2004-06-14 | 2008-01-31 | フラウンホッファー−ゲゼルシャフト ツァ フェルダールング デァ アンゲヴァンテン フォアシュンク エー．ファオ | Apparatus and method for converting an information signal into a spectral representation with variable resolution |
-
2016
- 2016-04-28 JP JP2016091695A patent/JP6627639B2/en active Active
Patent Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JPH07219581A (en) * | 1994-01-31 | 1995-08-18 | Babcock Hitachi Kk | Discrimination method for acoustic signal and device therefor |
| JP2004020193A (en) * | 2002-06-12 | 2004-01-22 | Takayoshi Yamamoto | Method for diagnosing object facility, computer program, and device for diagnosing subject facility |
| JP2008502927A (en) * | 2004-06-14 | 2008-01-31 | フラウンホッファー−ゲゼルシャフト ツァ フェルダールング デァ アンゲヴァンテン フォアシュンク エー．ファオ | Apparatus and method for converting an information signal into a spectral representation with variable resolution |
Cited By (9)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JP2019101777A (en) * | 2017-12-04 | 2019-06-24 | 株式会社明電舎 | Abnormality diagnosing apparatus and abnormality diagnosing method |
| JP7009961B2 (en) | 2017-12-04 | 2022-01-26 | 株式会社明電舎 | Abnormality diagnosis device and abnormality diagnosis method |
| JP2020079718A (en) * | 2018-11-12 | 2020-05-28 | 東ソー株式会社 | Joint failure detection method of metal-resin composite |
| JP7218546B2 (en) | 2018-11-12 | 2023-02-07 | 東ソー株式会社 | METHOD FOR DETECTING POOR JOINING OF METAL MEMBER-RESIN MEMBER COMPOSITE |
| CN113557414A (en) * | 2019-03-22 | 2021-10-26 | Abb瑞士股份有限公司 | Device for monitoring equipment |
| JP2022525963A (en) * | 2019-03-22 | 2022-05-20 | アーベーベー・シュバイツ・アーゲー | Equipment for equipment monitoring |
| US11835429B2 (en) | 2019-03-22 | 2023-12-05 | Abb Schweiz Ag | Apparatus for equipment monitoring |
| JP7485693B2 (en) | 2019-03-22 | 2024-05-16 | アーベーベー・シュバイツ・アーゲー | Equipment for monitoring equipment |
| JP2020181443A (en) * | 2019-04-26 | 2020-11-05 | 株式会社豊田中央研究所 | Abnormality detection apparatus, abnormality detection method, and computer program |
Also Published As
| Publication number | Publication date |
|---|---|
| JP6627639B2 (en) | 2020-01-08 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| US20210089917A1 (en) | Heuristic Inference of Topological Representation of Metric Relationships | |
| EP2499504B1 (en) | A precision measurement of waveforms using deconvolution and windowing | |
| JP6627639B2 (en) | Abnormality diagnosis device and abnormality diagnosis method | |
| US5938594A (en) | Method and apparatus for detecting nonlinearity and chaos in a dynamical system | |
| Uğuz et al. | A biomedical system based on hidden Markov model for diagnosis of the heart valve diseases | |
| CN108009122B (en) | Improved HHT method | |
| Leise | Wavelet-based analysis of circadian behavioral rhythms | |
| EP2880578A2 (en) | Estimating remaining useful life from prognostic features discovered using genetic programming | |
| US11055382B2 (en) | Methods and systems that estimate a degree of abnormality of a complex system | |
| JP2018501561A (en) | Quality control engine for complex physical systems | |
| Pang et al. | Fault state recognition of wind turbine gearbox based on generalized multi-scale dynamic time warping | |
| JP7140426B2 (en) | Time-varying structure instantaneous frequency determination method, system, device and storage medium | |
| CN113657244B (en) | Fan gearbox fault diagnosis method and system based on improved EEMD and speech spectrum analysis | |
| Liu et al. | Bearing failure diagnosis at time-varying speed based on adaptive clustered fractional Gabor transform | |
| JP6677069B2 (en) | Constant Q conversion component operation device and constant Q conversion component operation method | |
| US20220350987A1 (en) | Feature amount extraction device, time-sequential inference apparatus, time-sequential learning system, time-sequential feature amount extraction method, time-sequential inference method, and time-sequential learning method | |
| CN110675858A (en) | Terminal control method and device based on emotion recognition | |
| TWI743842B (en) | Method and system for processing electroencephalogram signal | |
| CN108629441A (en) | Prediction technique and device based on clustering and the improved fan noise of small echo | |
| Marchi et al. | Adaptive synchrosqueezing wavelet transform for real-time applications | |
| Bopaiah et al. | Precision/Recall trade-off analysis in Abnormal/Normal heart sound classification | |
| JP2021047100A (en) | Diagnostic device and diagnostic method | |
| JP7009961B2 (en) | Abnormality diagnosis device and abnormality diagnosis method | |
| CN117349661B (en) | Method, device, equipment and storage medium for extracting vibration signal characteristics of plunger pump | |
| JP6385113B2 (en) | An apparatus, method, and computer program for analyzing the relaxation elastic modulus of a polymer model. |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| A621 | Written request for application examination |
Free format text: JAPANESE INTERMEDIATE CODE: A621 Effective date: 20190219 |
|
| A977 | Report on retrieval |
Free format text: JAPANESE INTERMEDIATE CODE: A971007 Effective date: 20190917 |
|
| A131 | Notification of reasons for refusal |
Free format text: JAPANESE INTERMEDIATE CODE: A131 Effective date: 20191001 |
|
| A521 | Request for written amendment filed |
Free format text: JAPANESE INTERMEDIATE CODE: A523 Effective date: 20191021 |
|
| TRDD | Decision of grant or rejection written | ||
| A01 | Written decision to grant a patent or to grant a registration (utility model) |
Free format text: JAPANESE INTERMEDIATE CODE: A01 Effective date: 20191105 |
|
| A61 | First payment of annual fees (during grant procedure) |
Free format text: JAPANESE INTERMEDIATE CODE: A61 Effective date: 20191118 |
|
| R150 | Certificate of patent or registration of utility model |
Ref document number: 6627639 Country of ref document: JP Free format text: JAPANESE INTERMEDIATE CODE: R150 |