JP5854551B2 - Real-time frequency analysis method - Google Patents
Real-time frequency analysis method Download PDFInfo
- Publication number
- JP5854551B2 JP5854551B2 JP2011191146A JP2011191146A JP5854551B2 JP 5854551 B2 JP5854551 B2 JP 5854551B2 JP 2011191146 A JP2011191146 A JP 2011191146A JP 2011191146 A JP2011191146 A JP 2011191146A JP 5854551 B2 JP5854551 B2 JP 5854551B2
- Authority
- JP
- Japan
- Prior art keywords
- frequency
- output
- time
- data
- value
- 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.)
- Active
Links
- 238000004458 analytical method Methods 0.000 title claims description 89
- 238000004364 calculation method Methods 0.000 claims description 62
- 238000000034 method Methods 0.000 claims description 58
- 238000012545 processing Methods 0.000 claims description 53
- 238000000605 extraction Methods 0.000 claims description 37
- 230000008569 process Effects 0.000 claims description 23
- 238000012546 transfer Methods 0.000 claims description 23
- 230000004044 response Effects 0.000 claims description 19
- 238000003672 processing method Methods 0.000 claims description 16
- 238000005070 sampling Methods 0.000 claims description 13
- 230000001133 acceleration Effects 0.000 claims description 3
- 230000006870 function Effects 0.000 description 19
- 230000008859 change Effects 0.000 description 14
- 238000001228 spectrum Methods 0.000 description 8
- 238000004422 calculation algorithm Methods 0.000 description 7
- 230000036541 health Effects 0.000 description 6
- 238000005259 measurement Methods 0.000 description 6
- 238000012544 monitoring process Methods 0.000 description 5
- 238000010586 diagram Methods 0.000 description 3
- 239000000284 extract Substances 0.000 description 3
- 238000012935 Averaging Methods 0.000 description 2
- 230000005540 biological transmission Effects 0.000 description 2
- 238000004891 communication Methods 0.000 description 2
- 238000001514 detection method Methods 0.000 description 2
- 238000011156 evaluation Methods 0.000 description 2
- 230000004888 barrier function Effects 0.000 description 1
- 230000008901 benefit Effects 0.000 description 1
- 238000006243 chemical reaction Methods 0.000 description 1
- 238000007796 conventional method Methods 0.000 description 1
- 230000006866 deterioration Effects 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 238000010606 normalization Methods 0.000 description 1
- 230000009467 reduction Effects 0.000 description 1
- 238000009991 scouring Methods 0.000 description 1
- 239000004065 semiconductor Substances 0.000 description 1
- 230000002123 temporal effect Effects 0.000 description 1
Images
Landscapes
- Measurement Of Mechanical Vibrations Or Ultrasonic Waves (AREA)
Description
本発明は、リアルタイム周波数解析方法に関し、特に、計算のアルゴリズムを簡素化して処理装置や記憶装置の負担を少なくしたリアルタイム周波数解析方法に関する。 The present invention relates to a real-time frequency analysis method, and more particularly to a real-time frequency analysis method that simplifies a calculation algorithm and reduces the burden on a processing device and a storage device.
従来、時々刻々に周波数成分が変化する信号に対する周波数解析として使われてきた技術には、主に２つの技術が有ることが電気情報通信学会誌８／’９６（非特許文献１）に記載されている。以下では、この２つの技術について説明する。
（１） 短時間フーリエ変換を用いた解析方法
この方法は、基本的には、定常信号を解析するためのフーリエ変換［下記（１６）式参照］の処理を、図１０に示すように、０〜T秒間の信号をＮ分割したＴ／Ｎの長さ毎（ｉ＝１〜Ｎ）の信号に適用し、その結果を、該時間区間の周波数成分とする解析法である。
Conventionally, it has been described in the Journal of the Institute of Electrical, Information and Communication Engineers 8 / '96 (Non-patent Document 1) that there are mainly two techniques that have been used as frequency analysis for signals whose frequency components change from moment to moment. ing. Hereinafter, these two techniques will be described.
(1) Analysis method using short-time Fourier transform In this method, basically, processing of Fourier transform [see the following equation (16)] for analyzing a stationary signal is performed as shown in FIG. This is an analysis method in which a signal for ˜T seconds is applied to a signal for each T / N length (i = 1 to N) obtained by dividing the signal into N, and the result is used as a frequency component of the time interval.
この方法では、フーリエ変換が積分変換となるため、演算に多大の時間を要する。また、演算を早くする高速フーリエ変換（FFT）のアルゴリズムを使ったとしても、時々刻々の演算は複雑となり、また、ＦＦＴの制約により、周波数分解能ΔＦは、計算に用いる信号の長さの逆数、即ちΔｆ＝Ｎ／Ｔにならざるを得ないため、図１０において、周波数の時間変化を細かく見るために分割数Ｎを大きくする（個々の計算する信号の時間長さを短くする）と、周波数分解能は荒くならざるを得ず、これとは逆に、周波数分解能を細かくしようとする（個々の計算に用いる信号の長さを長くする）と、周波数成分の時間変化の分解能を荒くせざるを得ないといった欠点を有している。
（２） ウェーブレット変換［下記（１７）式参照］を用いた解析方法
In this method, since the Fourier transform becomes an integral transform, a long time is required for the calculation. Even if a fast Fourier transform (FFT) algorithm that speeds up the computation is used, the computation is complicated every moment, and due to the constraints of FFT, the frequency resolution ΔF is the reciprocal of the length of the signal used for computation, That is, Δf = N / T, and in FIG. 10, if the number of divisions N is increased in order to see the time change of the frequency in detail (the time length of each signal to be calculated is shortened), the frequency The resolution must be rough, and conversely, if you try to make the frequency resolution finer (increase the length of the signal used for each calculation), the resolution of the time change of the frequency component must be rough. It has the disadvantage that it cannot be obtained.
(2) Analysis method using wavelet transform [refer to equation (17) below]
この方法も、前述のフーリエ変換と同様で、（１６）式でのｅの乗数である−ｊωｔの代わりに、マザーウェーブレットと呼ばれる関数を、元の信号に掛けた積分変換であるため、計算のアルゴリズムが複雑となり、計算にも多大の時間を要することになる。
このマザーウェーブレットと呼ばれる関数には、幾つかの種類の関数が提案されているが、例えば、スケールと呼ばれている数ａの逆数である１／ａが周波数に対応し、数ｂが時間のシフトに相当するような、ａ，ｂ対応の関数も提案されている。
しかしながら、この解析方法も、前述の計算時間の問題だけでは済まず、周波数分解能が周波数によって倍々に変化することになるため、周波数が高くなるに従って周波数分解能が荒くならざるを得ないことになり、また時間分解能は、低周波になるに従って荒くならざるを得ないといった欠点を有している。
この分野の特許文献としては、例えば特許文献１（特開２０００−１８６９８４号公報）に、簡単な構成および処理によって、洗掘等に起因する橋脚の支持力の変化を適切に且つ定常的に監視し、運行規制の的確な発令／解除を可能とする技術が開示されている。具体的には、ＦＦＴ処理部が、橋脚の天端に取着された加速度計による検出データに高速フーリエ変換を施してフーリエスペクトルを得る。平均化処理部は、フーリエスペクトルを所定の時間についての加算平均により時間的な平均を求め、所定の周波数幅について周波数軸上で移動平均を求めて、フーリエスペクトル波形を平滑化する。重み付け処理部は、健全時の橋脚の伝達関数に基づく１自由度系の応答倍率曲線に基づく重み付け関数で平滑化フーリエスペクトルを重み付けする。面積算定正規化処理部は、重み付けフーリエスペクトルの曲線とベースラインとで囲まれる領域の面積を求めて正規化し、正規化面積減小評価部による閾値との比較評価に供するものとしている。
This method is also the same as the above-described Fourier transform, and instead of -jωt, which is a multiplier of e in equation (16), is an integral transform obtained by multiplying the original signal by a function called a mother wavelet. The algorithm becomes complicated and the calculation takes a lot of time.
Several types of functions have been proposed for this function called mother wavelet. For example, 1 / a, which is the inverse of the number a called scale, corresponds to the frequency, and the number b is the time. A function corresponding to a and b corresponding to a shift has also been proposed.
However, this analysis method is not limited to the above-described problem of calculation time, and the frequency resolution changes twice according to the frequency. Therefore, the frequency resolution must be roughened as the frequency increases. In addition, the time resolution has a drawback that it must be roughened as the frequency becomes low.
As a patent document in this field, for example, in Patent Document 1 (Japanese Patent Laid-Open No. 2000-186984), a change in the supporting force of a bridge pier due to scouring or the like is appropriately and regularly monitored with a simple configuration and processing. However, a technology that enables accurate issuance / cancellation of operation regulations is disclosed. Specifically, the FFT processing unit performs a fast Fourier transform on the detection data from the accelerometer attached to the top of the pier to obtain a Fourier spectrum. The averaging processing unit obtains a temporal average of the Fourier spectrum by addition averaging for a predetermined time, obtains a moving average on the frequency axis for a predetermined frequency width, and smoothes the Fourier spectrum waveform. The weighting processing unit weights the smoothed Fourier spectrum with a weighting function based on a response magnification curve of a one-degree-of-freedom system based on a transfer function of a bridge pier in a healthy state. The area calculation normalization processing unit obtains and normalizes the area of the region surrounded by the curve of the weighted Fourier spectrum and the baseline, and uses it for comparison evaluation with the threshold value by the normalized area reduction evaluation unit.
ところで、リアルタイム周波数解析方法は、構造物の振動を利用したヘルスモニタリングに使用されることが多いが、その普及には構造物に取り付ける振動センサ内にてある程度振動特性を評価できるところまで演算処理した結果が出力できるような安価なスマートセンサの実現が必要であると考えられる。
このためのスマートセンサ内での演算処理は、高度な処理能力や大きな記憶容量は望めないであろうことから、負荷が少なく、また、記憶容量を多く必要としない処理法が必要となる。
但し、構造物の振動データから、この振動特性の変化を見ようとした場合、少なくとも振動データの周波数成分の経時変化を求めることができる時間周波数解析的な手法が必要になる。
しかしながら、上記背景技術で述べた従来のリアルタイム周波数解析方法にあっては、前述のとおり、時間周波数解析としては短時間フーリエ変換によるランニングスペクトル、ウェーブレット解析等が通常使われているが、これらの解析手法は、演算ステップ数がリアルタイム処理を行うには多めであり、かつ、ある程度の時間データが必要なため、処理装置の処理能力や記憶容量の点で、安価なスマートセンサに組込むのは難しいと考えられる。
By the way, the real-time frequency analysis method is often used for health monitoring using the vibration of the structure, but for its widespread use, the vibration sensor attached to the structure has been processed up to a point where vibration characteristics can be evaluated to some extent. It is considered necessary to realize an inexpensive smart sensor that can output the result.
For this purpose, the processing within the smart sensor cannot be expected to have a high processing capacity and a large storage capacity. Therefore, a processing method that requires a small load and does not require a large storage capacity is required.
However, when it is attempted to see the change in the vibration characteristics from the vibration data of the structure, a time frequency analysis method capable of obtaining at least a change with time of the frequency component of the vibration data is required.
However, in the conventional real-time frequency analysis method described in the background art, as described above, a running spectrum by short-time Fourier transform, wavelet analysis, etc. are usually used as the time-frequency analysis. The method has a large number of computation steps to perform real-time processing, and requires a certain amount of time data, so it is difficult to incorporate it into an inexpensive smart sensor in terms of processing capacity and storage capacity of the processing device. Conceivable.
また、周波数分解能に関して、必要な周波数帯域内で必要な周波数分解能を確保することが難しい。
これらの点から、従来のリアルタイム周波数解析方法は、構造物のリアルタイムヘルスモニタリングに用いるには不向きであるという問題点があった。
本発明は、上記従来の問題点に鑑みてなされたものであって、振動の物理現象と直感的に結びつく１自由度振動系を基にした複数の狭帯域バンドパスフィルタを用いる解析法を使用することを可能にして、必要な周波数分解能を確保しつつ、処理装置に対して負荷が少なく、また、記憶容量を多く必要としない、処理能力の高くない処理装置に適用可能で、特に構造物の振動特性評価のために好適なリアルタイム周波数解析方法を提供することを目的としている。
In addition, regarding the frequency resolution, it is difficult to ensure the necessary frequency resolution within the necessary frequency band.
From these points, the conventional real-time frequency analysis method has a problem that it is not suitable for real-time health monitoring of structures.
The present invention has been made in view of the above-described conventional problems, and uses an analysis method using a plurality of narrow-band bandpass filters based on a one-degree-of-freedom vibration system that is intuitively linked to a physical phenomenon of vibration. making it possible to, while maintaining the frequency resolution needed, less load to the processing unit, also not much storage capacity required, possible, in particular structure applicable to not higher processor processing capability An object of the present invention is to provide a real-time frequency analysis method suitable for evaluating vibration characteristics of
請求項１に記載した本発明に係るリアルタイム周波数解析方法は、上述した目的を達成するために、加速度、荷重等の物理量を測定するセンサの出力を一定の時間間隔でサンプリングして取り出すＡ／Ｄ変換器の出力時系列データ、または該出力時系列データを一旦記憶装置に集録したデータに対し、任意の周波数成分の振動特性を抽出するリアルタイム周波数解析法であって、
任意の１つの中心周波数近傍の狭帯域の周波数成分の波形を抽出する狭帯域周波数成分波形抽出ステップと、
前記狭帯域周波数成分波形抽出ステップで抽出された波形の絶対値を求める絶対値回路の機能を有する絶対値算出ステップと、
前記絶対値算出ステップの出力絶対値波形の、任意時間毎の移動平均を求める移動平均算出ステップと、
を備え、
前記狭帯域周波数成分波形抽出ステップは、
用いる時系列データの現時点のデータおよび２サンプル前のデータと、現時点より１サンプル前および２サンプル前の当該ステップの出力から、現時点の前記当該ステップの出力を算出する機能を有し、
用いる時系列データのサンプリング時間間隔と、任意の中心周波数と、中心周波数に対する１より小さい任意の通過帯域比率とを用いて該中心周波数に対する通過帯域比率分の周波数だけ中心周波数から低い周波数側と高い周波数側にずれたカットオフ周波数と、前記任意の通過帯域比率に応じた遮断特性を有し、
前記絶対値算出ステップは、
入力された値が正値の場合はそのままの値を出力し、また、入力された値が負値の場合は負の符号を外した値を出力する機能を有し、
前記移動平均算出ステップは、
前記狭帯域周波数成分波形抽出ステップにおいて設定した任意の中心周波数に相当する周期の整数倍の時間を、前記狭帯域周波数成分振幅抽出ステップに適用するデータのサンプリング時間間隔で除したデータ数分の逆数を当該ステップに用いる時系列データの現時点のデータに乗じ、この値と、１から移動平均を取るデータ数の逆数を差し引いた値を現時点の１サンプル前の当該移動平均ステップの算出結果に乗じた値と、を足し合わせた値を当該移動平均ステップの現時点の算出結果とすることで、順次１サンプリング間隔毎の平均値を算出する機能を有し、
前記移動平均算出ステップの出力結果が前記狭帯域周波数成分波形抽出ステップにおいて設定した任意の１つの中心周波数近傍の周波数成分の時々刻々の振幅の実行値に相当する値を任意の１つの周波数成分の振動特性として抽出するに際し、掛け算の回数および一連のステップで記憶を必要とするデータの数を少なくして、演算負荷と必要記憶容量を小さく抑えるように構成したことを特徴としている。
In order to achieve the above-described object, the real-time frequency analysis method according to the first aspect of the present invention is an A / D that samples and extracts the output of a sensor that measures a physical quantity such as acceleration and load at regular time intervals. input time series data of the transducer or to data acquire time-series data output to one Dan storage device, a real-time frequency analysis method for extracting the vibration characteristics of any frequency component,
A narrowband frequency component waveform extraction step for extracting a waveform of a narrowband frequency component near any one center frequency;
An absolute value calculating step having a function of an absolute value circuit for obtaining an absolute value of the waveform extracted in the narrowband frequency component waveform extracting step;
A moving average calculating step for obtaining a moving average for each arbitrary time of the output absolute value waveform of the absolute value calculating step;
Equipped with a,
The narrowband frequency component waveform extraction step includes:
A function for calculating the current output of the step from the current data of the time-series data to be used and the data of two samples before and the output of the step one sample before and two samples before the current time;
Using the sampling time interval of the time-series data to be used, an arbitrary center frequency, and an arbitrary passband ratio smaller than 1 with respect to the center frequency, the frequency corresponding to the passband ratio with respect to the center frequency is higher on the lower frequency side than the center frequency. A cutoff frequency shifted to the frequency side, and a cutoff characteristic according to the arbitrary passband ratio,
The absolute value calculating step includes:
When the input value is a positive value, the value is output as it is, and when the input value is a negative value, it has a function of outputting a value with a negative sign removed,
The moving average calculation step includes:
The reciprocal of the number of data obtained by dividing the time that is an integral multiple of the period corresponding to the arbitrary center frequency set in the narrowband frequency component waveform extraction step by the sampling time interval of data applied to the narrowband frequency component amplitude extraction step Is multiplied by the current data of the time-series data used in this step, and this value is multiplied by the result of subtracting the reciprocal of the number of data to be taken from 1 from the calculation result of the moving average step one sample before the current time. By adding the value and the value as the current calculation result of the moving average step, it has a function of calculating an average value for each sampling interval sequentially,
The output result of the moving average calculation step is a value corresponding to the execution value of the momentary amplitude of the frequency component in the vicinity of any one center frequency set in the narrowband frequency component waveform extraction step. When extracting as vibration characteristics, the number of multiplications and the number of data that need to be stored in a series of steps are reduced to reduce the calculation load and the required storage capacity .
また、請求項２に記載した本発明に係るリアルタイム周波数解析方法は、請求項１のリアルタイム周波数解析方法であって、請求項１のリアルタイム周波数解析方法の全ステップを、所定の複数個の中間周波数毎に並列に処理できるステップを有し、前記複数並列の処理方法の各々の処理方法で用いる狭帯域周波数成分波形抽出ステップの中心周波数と、中心周波数に対する１より小さい任意の通過帯域比率とは、各々異なる値を設定することができ、設定した任意の複数の中心周波数近傍の周波数成分の時々刻々の振幅の実行値に相当する値を任意の複数の周波数成分の振動特性として抽出することを特徴としている。
Also, real-time frequency analysis method according to the present invention described in
また、請求項３に記載した本発明に係るリアルタイム周波数解析方法は、請求項１のリアルタイム周波数解析方法であって、前記解析対象への入力信号を検出する第１のセンサの出力データと前記解析対象からの出力信号を検出する前記第２のセンサの出力データに対して、請求項１に記載の前記狭帯域周波数成分波形抽出ステップと、前記絶対値算出ステップと、前記移動平均算出ステップとを適用する入出力狭帯域周波数成分波形抽出ステップと、
前記入出力狭帯域周波数成分波形抽出ステップの、前記第１のセンサと前記第２のセンサとに対応した出力の出力比を求める除算ステップと、
を備え、
前記第１のセンサと前記第２のセンサとの出力に適用する、請求項１記載の前記狭帯域周波数成分波形抽出ステップと、前記絶対値算出ステップと、前記移動平均算出ステップにおける各処理方法では、同じ中心周波数と、同じ通過帯域比率と、同じ整数周期分の移動平均データ数とを参照することで、前記解析対象に設定した任意の中心周波数に対する伝達特性応答倍率を、任意の周波数成分の振動特性として抽出することを特徴としている。
The real-time frequency analysis method according to the present invention described in
A division step of obtaining an output ratio of outputs corresponding to the first sensor and the second sensor in the input / output narrowband frequency component waveform extraction step;
Equipped with a,
The processing methods in the narrowband frequency component waveform extraction step, the absolute value calculation step, and the moving average calculation step according to
また、請求項４に記載した本発明に係るリアルタイム周波数解析方法は、請求項１のリアルタイム周波数解析方法であって、前記解析対象への入力信号を検出する前記第１のセンサの出力データと前記第２のセンサの出力データに対して、請求項１に記載の前記狭帯域周波数成分波形抽出ステップと、前記絶対値算出ステップと、前記移動平均算出ステップとを適用する入出力狭帯域周波数成分波形抽出ステップと、
前記入出力狭帯域周波数成分波形抽出ステップの、前記第１のセンサと前記第２のセンサとに対応した出力の出力比を求める除算ステップと、
前記除算ステップを前記第１のセンサに対応するデータに適用した結果の出力と、前記除算ステップを前記第２のセンサに対応するデータに適用した結果の出力とを乗算する伝達特性応答倍率算出ステップと、
請求項１に記載の前記狭帯域周波数成分波形抽出ステップを前記第１のセンサの出力データに適用した結果の出力と、前記狭帯域周波数成分波形抽出ステップを前記第２のセンサの出力データに適用した結果の出力とを乗算する波形乗算ステップと、
前記波形乗算ステップの出力波形に対し任意時間毎の移動平均を求める移動平均算出ステップと、
前記移動平均算出ステップの算出結果の２倍を、前記伝達特性応答倍率算出ステップの算出結果で除した結果の逆余弦を求める位相算出ステップと、
を備えたことを特徴としている。
Also, real-time frequency analysis method according to the present invention described in
A division step of obtaining an output ratio of outputs corresponding to the first sensor and the second sensor in the input / output narrowband frequency component waveform extraction step;
Transfer characteristic response magnification calculation step of multiplying the output of the result obtained by applying the division step to the data corresponding to the first sensor and the output of the result of applying the division step to the data corresponding to the second sensor When,
Applying the narrowband frequency component waveform extracting step of
A moving average calculating step for obtaining a moving average for each arbitrary time with respect to the output waveform of the waveform multiplying step;
A phase calculation step for obtaining an inverse cosine of a result obtained by dividing twice the calculation result of the moving average calculation step by the calculation result of the transfer characteristic response magnification calculation step;
It is characterized by having.
また、請求項５に記載した本発明に係るリアルタイム周波数解析方法は、請求項３に記載のリアルタイム周波数解析方法の全ステップを、所定の複数個の中間周波数毎に並列に処理するステップを有することを特徴としている。
また、請求項６に記載した本発明に係るリアルタイム周波数解析方法は、請求項４に記載のリアルタイム周波数解析方法の全ステップを、所定の複数個の中間周波数毎に並列に処理するステップを有することを特徴としている。
Also, real-time frequency analysis method according to the present invention described in
Also, real-time frequency analysis method according to the present invention described in
以上述べたように、本発明によれば、周波数分解能と時間分解能とを、時々刻々に変化する測定結果の現在値を示す関数に含まれる係数内の中心周波数と、通過帯域比率を設定することにより、任意に定めることが可能となり、計算アルゴリズムについても、１つの周波数に対して、前記関数を基本とした数式を時々刻々計算するだけであり、
特に、前記移動平均算出ステップでは、
前記狭帯域周波数成分波形抽出ステップにおいて設定した任意の中心周波数に相当する周期の整数倍の時間を、前記狭帯域周波数成分振幅抽出ステップに適用するデータのサンプリング時間間隔で除したデータ数分の逆数を当該ステップに用いる時系列データの現時点のデータに乗じ、この値と、１から移動平均を取るデータ数の逆数を差し引いた値を現時点の１サンプル前の当該移動平均ステップの算出結果に乗じた値と、を足し合わせた値を当該移動平均ステップの現時点の算出結果とすることで、順次１サンプリング間隔毎の平均値を算出する機能を有し、
前記移動平均算出ステップの出力結果が前記狭帯域周波数成分波形抽出ステップにおいて設定した任意の１つの中心周波数近傍の周波数成分の時々刻々の振幅の実行値に相当する値を任意の１つの周波数成分の振動特性として抽出するに際し、掛け算の回数および一連のステップで記憶を必要とするデータの数を少なくして、演算負荷と必要記憶容量を小さく抑えるように構成したことにより、従来よりも単純な計算処理となるので、演算負荷を少なくすることができ、処理能力の高くない処理装置に適用可能で、負荷が少なく、記憶容量を多く必要としない構造物の振動特性評価のための時間周波数解析法、伝達関数算出法を構築することができる。また、本発明によれば、構造物に取り付けるヘルスモニタリング用スマートセンサ内に構造物の振動データを用いた健全性評価のための時々刻々の処理アルゴリズムを組み込むことが容易になる。
また、請求項２〜請求項６の記載によれば、いずれも請求項１に規定するように、掛け算の回数および一連のステップで記憶を必要とするデータの数を少なくして、演算負荷と必要記憶量を小さく抑えるように構成することにより、処理能力の高くない処理装置に適用可能で、負荷が少なく、記憶容量を多く必要としない構造物の振動特性評価のために好適なリアルタイム周波数解析方法を提供することができ、延いては、構造物に取り付けるヘルスモニタリング用スマートセンサ内に構造物の振動データを用いた健全性評価のための時々刻々の処理アルゴリズムを組み込むことが容易になる。
As described above, according to the present invention, the frequency resolution and the time resolution are set with the center frequency in the coefficient included in the function indicating the current value of the measurement result that changes from moment to moment, and the passband ratio. Thus, it is possible to arbitrarily determine the calculation algorithm, and for each frequency, only a mathematical expression based on the function is calculated momentarily.
In particular, in the moving average calculation step,
The reciprocal of the number of data obtained by dividing the time that is an integral multiple of the period corresponding to the arbitrary center frequency set in the narrowband frequency component waveform extraction step by the sampling time interval of data applied to the narrowband frequency component amplitude extraction step Is multiplied by the current data of the time-series data used in this step, and this value is multiplied by the result of subtracting the reciprocal of the number of data to be taken from 1 from the calculation result of the moving average step one sample before the current time. By adding the value and the value as the current calculation result of the moving average step, it has a function of calculating an average value for each sampling interval sequentially,
The output result of the moving average calculation step is a value corresponding to the execution value of the momentary amplitude of the frequency component in the vicinity of any one center frequency set in the narrowband frequency component waveform extraction step. When extracting as vibration characteristics, the number of multiplications and the number of data that need to be stored in a series of steps are reduced, and the calculation load and the required storage capacity are reduced, making calculation simpler than before. Time-frequency analysis method for evaluating vibration characteristics of structures that can reduce computation load, can be applied to processing devices with low processing capacity, and do not require much storage capacity A transfer function calculation method can be constructed. In addition, according to the present invention, it is easy to incorporate a processing algorithm for evaluating the soundness using the vibration data of the structure into the smart sensor for health monitoring attached to the structure.
Further, according to
本発明のリアルタイム周波数解析方法は、加速度、荷重等の物理量を測定するセンサの出力を一定の時間間隔でサンプリングして取り出すＡ／Ｄ変換器（図示は省略）を使用し、該Ａ／Ｄ変換器の出力時系列データ、または該出力時系列データを、一旦記憶装置に集録したデータに対し、任意の周波数成分の振動特性を抽出する方法である。
よって、本発明のリアルタイム周波数解析方法は、構造物のヘルスモニタリングのためのリアルタイム周波数解析に適した方法である。
以下、本発明のリアルタイム周波数解析方法の実施の形態について、〔実施の形態１〕〜〔実施の形態７〕の順に、図面を参照して詳細に説明する。
The real-time frequency analysis method of the present invention uses an A / D converter (not shown) that samples and extracts the output of a sensor that measures physical quantities such as acceleration and load at regular time intervals. vessels of input time series data or the time-series data output, with respect to data acquisition to one Dan storage device, a method of extracting the vibration characteristics of the arbitrary frequency component.
Therefore, the real-time frequency analysis method of the present invention is a method suitable for real-time frequency analysis for health monitoring of structures.
Hereinafter, embodiments of the real-time frequency analysis method of the present invention will be described in detail in the order of [Embodiment 1] to [Embodiment 7] with reference to the drawings.
〔実施の形態１〕
本実施の形態１に係るリアルタイム周波数解析方法は、設定した任意の１つの中心周波数近傍の周波数成分の時々刻々の振幅の実行値に相当する値を任意の１つの周波数成分の振動特性として抽出する方法である。
即ち、本方法は、特定の周波数成分の時間変化を計算する方法であり、例えば、回転機器周辺の振動において、回転機器の回転数と同じ周波数成分の変化を得ようとする場合、中心周波数を、時々刻々の回転機器の回転数に相当する値に変えることで、上記の目的を達成することができる。
具体的な処理ステップとしては、任意の１つの中心周波数近傍の狭帯域の周波数成分の波形を抽出する狭帯域周波数成分波形抽出ステップと、前記狭帯域周波数成分波形抽出ステップで抽出された波形の絶対値を求める絶対値回路の機能を有する絶対値算出のステップと、前記絶対値算出ステップの出力絶対値波形に対し任意時間毎の移動平均を求める移動平均算出ステップと、を有する。
[Embodiment 1]
The real-time frequency analysis method according to the first embodiment extracts a value corresponding to the execution value of the amplitude of the frequency component in the vicinity of any one set center frequency as the vibration characteristic of any one frequency component. Is the method.
In other words, this method is a method for calculating a time change of a specific frequency component. For example, in the case of vibration around a rotating device, when trying to obtain the same frequency component change as the rotational speed of the rotating device, the center frequency is calculated. By changing the value to the value corresponding to the rotational speed of the rotating device from moment to moment, the above object can be achieved.
As specific processing steps, a narrowband frequency component waveform extracting step for extracting a waveform of a narrowband frequency component near any one center frequency, and an absolute waveform extracted in the narrowband frequency component waveform extracting step are described. that Yusuke a step of absolute value calculating with the function of the absolute value circuit for obtaining a value, a moving average calculation step of calculating the moving average of each arbitrary time with respect to the output absolute value waveform of the absolute value calculation step.
前記狭帯域周波数成分波形抽出ステップは、時系列データの現時点のデータおよび２サンプル前のデータと、現時点より１サンプル前および２サンプル前の当該ステップの出力とから、現時点の出力を算出し、その狭帯域周波数成分波形抽出の周波数特性は、時系列データのサンプリング時間間隔と、任意の中心周波数と、中心周波数に対する１より小さい任意の通過帯域比率と、を用いることにより、該中心周波数に対する通過帯域比率分の周波数だけ中心周波数から低い周波数側と高い周波数側にずれたカットオフ周波数と、前記任意の通過帯域比率に応じた遮断特性とを有する。
また、絶対値算出ステップは、入力された値が正の値の場合はそのままの値を出力し、入力された値が負の値の場合は該負を示す符号を外した値を出力する。
さらに、移動平均算出ステップは、前記狭帯域周波数成分波形抽出ステップにおいて設定した任意の中心周波数に相当する周期の整数倍の時間を、前記狭帯域周波数成分振幅抽出ステップに適用するデータのサンプリング時間間隔で除したデータ数分の逆数を当該ステップに用いる時系列データの現時点のデータに乗じ、この値と、１から移動平均を取るデータ数の逆数を差し引いた値を現時点の１サンプル前の当該移動平均ステップの算出結果に乗じた値と、を足し合わせた値を当該移動平均ステップの現時点の算出結果とすることで、順次１サンプリング間隔毎の平均値を算出する機能を有し、
前記移動平均算出ステップの出力結果が前記狭帯域周波数成分波形抽出ステップにおいて設定した任意の１つの中心周波数近傍の周波数成分の時々刻々の振幅の実行値に相当する値を任意の１つの周波数成分の振動特性として抽出するに際し、掛け算の回数および一連のステップで記憶を必要とするデータの数を少なくして、演算負荷と必要記憶容量を小さく抑えるように構成したものである。
この移動平均ステップは、信号の移動平均により、振幅実行値や平均値を求めるための移動平均計算を簡略化する方法である。この方法によれば、計算アルゴリズムの簡略化と計算速度の向上を達成することができる。
この（３）式で示す指数移動平均ステップは、移動平均（Ｚｎ）を取るデータ数の逆数を当該方法に用いる時系列データの現時点のデータに乗じ、この値と、１から移動平均を取るデータ数の逆数を差し引いた値を現時点の１サンプル前の前記指数移動平均ステップの算出結果に乗じた値とを足し合わせた値を、現時点の算出結果とするものである。
但し、（３）式の符号の意味を、下記のとおりとする。
ｚ：ｖの移動平均値；本ステップの出力
ｋ：移動平均を取るデータ数の逆数
The narrowband frequency component waveform extraction step calculates the current output from the current data of the time series data and the data of 2 samples before and the output of the
The absolute value calculation step, and outputs the raw value if the entered value is a positive value, when the input value is negative value you outputs a value obtained by removing a code indicating a negative .
Further, the moving average calculating step, the narrow band frequency component waveform extracted arbitrary integer multiple of the period corresponding to the center frequency time set in step, sampling time of data to be applied to the narrow-band frequency component amplitude extraction steps The reciprocal of the number of data divided by the interval is multiplied by the current data of the time-series data used in the step, and this value and the value obtained by subtracting the reciprocal of the number of data taking the moving average from 1 A value obtained by adding the value obtained by multiplying the calculation result of the moving average step and the current calculation result of the moving average step has a function of sequentially calculating an average value for each sampling interval,
The output result of the moving average calculation step is a value corresponding to the execution value of the momentary amplitude of the frequency component in the vicinity of any one center frequency set in the narrowband frequency component waveform extraction step. At the time of extraction as vibration characteristics, the number of multiplications and the number of data that need to be stored in a series of steps are reduced to reduce the calculation load and the required storage capacity.
This moving average step is a method of simplifying the moving average calculation for obtaining the amplitude execution value and the average value by the moving average of the signal. According to this method, the calculation algorithm can be simplified and the calculation speed can be improved.
The exponential moving average step represented by the equation (3) is obtained by multiplying the current data of the time series data used in the method by the reciprocal of the number of data for which the moving average (Zn) is taken, and taking this value and the moving average from 1 A value obtained by adding the value obtained by subtracting the reciprocal of the number to the calculation result of the exponential moving average step one sample before the current time is used as the current calculation result.
However, the meanings of the symbols in the expression (3) are as follows.
z: moving average of v; output of this step
k: Reciprocal of the number of data to take the moving average
図１は、本発明の実施の形態１に係るリアルタイム周波数解析方法の処理手順を示すフローチャート図である。
但し、ここで使用される符号の意味は、下記のとおりとする。
ｘ：使用する時系列データ
ｙ：周波数成分の波形出力
入出力の添え字は、ｎ：現在の値、ｎ−ｊ：現在よりｊ：サンプル前の値
ｖ：ｙの絶対値
ｚ：ｖの移動平均値：本ステップの出力
ｋ：移動平均を取るデータ数の逆数
ａ、ｂ、ｃ：係数
（処理手順）
まず、ステップＳ１の数式〔（１）式〕を計算し、次にステップＳ２の数式〔（２）式〕を計算し、最後にステップＳ３の数式〔（３）式〕を計算する。
FIG. 1 is a flowchart showing a processing procedure of a real-time frequency analysis method according to
However, the meanings of the symbols used here are as follows .
x : Time series data to be used y: Waveform output of frequency component Input / output subscripts are: n: current value, nj: current value j: value before sample v: absolute value of y z: movement of v Average value : Output of this step
k: Reciprocal number a, b, c of coefficient of moving average taking data (coefficient) (processing procedure)
First, the mathematical expression [Equation (1)] in Step S1 is calculated, then the mathematical expression [Equation (2)] in Step S2 is calculated, and finally the mathematical expression [Equation (3)] in Step S3 is calculated.
〔実施の形態２〕
本発明の実施の形態２に係るリアルタイム周波数解析方法は、設定した任意の複数の中心周波数近傍の周波数成分の時々刻々の振幅の実行値に相当する値を任意の複数の周波数成分の振動特性として抽出する方法、即ち、前記中心周波数についての周波数成分の時間変化を得るための方法であるが、従来法に比べ、出力関数の係数ａ、ｂ、ｃ内の中心周波数と通過帯域比率の設定を任意に行うことが可能なため、周波数毎に所望の周波数分解能と時間分解能を設定できる利点を有する。このことは、例えば構造物の損傷により構造物の固有振動数が変化するような現象においては、比較的長い時間間隔における固有振動数の微小な変化を検知する必要があるため、時間分解能は荒くし、周波数分解能を細かくするといった操作を可能にするものである。
具体的な処理ステップとしては、前述の実施の形態１に係るリアルタイム周波数解析方法の全ステップを複数並列に処理するように配置した処理方法を有し、前記複数並列に処理するように配置された前記の処理方法で用いる狭帯域周波数成分波形抽出ステップの中心周波数と、該中心周波数に対する１より小さい任意の通過帯域比率とは、各々、異なる値を設定することができるようにし、これにより、設定した任意の複数の中心周波数近傍の周波数成分の時々刻々の振幅の実行値に相当する値を任意の複数の周波数成分の振動特性として抽出できるようにしている。
[Embodiment 2]
In the real-time frequency analysis method according to the second embodiment of the present invention, a value corresponding to the execution value of the amplitude of each frequency component in the vicinity of any set plurality of center frequencies is set as the vibration characteristic of any set of frequency components. method of extracting, i.e., wherein at the center is a method for obtaining the time variation of the frequency component of the frequency, compared with the conventional method, the coefficient of the output functions a, b, the center frequency and pass band ratios in c Since setting can be arbitrarily performed, there is an advantage that desired frequency resolution and time resolution can be set for each frequency. This is because, for example, in a phenomenon in which the natural frequency of the structure changes due to damage to the structure, it is necessary to detect minute changes in the natural frequency in a relatively long time interval, so the time resolution is rough. In addition, operations such as finer frequency resolution are possible.
As a specific processing step, there is a processing method arranged so as to process all the steps of the real-time frequency analysis method according to the first embodiment described above in parallel, and the processing steps are arranged so as to process in parallel. The center frequency of the narrowband frequency component waveform extraction step used in the processing method and an arbitrary passband ratio smaller than 1 with respect to the center frequency can be set to different values, respectively. The value corresponding to the effective value of the amplitude of the frequency components in the vicinity of the plurality of arbitrary center frequencies can be extracted as the vibration characteristics of the plurality of frequency components.
図２は、本発明の実施の形態２に係るリアルタイム周波数解析方法処理手順を示すフローチャート図である。
但し、ここでは、測定環境を示す一般的な符号、及び図２に記載の符号の意味を下記のとおりとする。
ｘ：使用する時系列データ
ｙ：周波数成分の波形出力
入出力の添え字は、ｎ：現在の値、ｎ−ｊ：現在よりｊ：サンプル前の値
ｖ：ｙの絶対値
ｚ：ｖの移動平均値（ｐ個を出力）
ｋ：移動平均を取るデータ数の逆数
（処理手順）
図２におけるステップＳ１〜Ｓ３の処理を、ｉ＝１〜ｐ個並べて並列に処理する。
FIG. 2 is a flowchart showing a real-time frequency analysis method processing procedure according to the second embodiment of the present invention.
However, here, the meanings of general symbols indicating the measurement environment and the symbols described in FIG. 2 are as follows .
x : Time series data to be used y: Waveform output of frequency component Input / output subscripts are: n: current value, nj: current value j: value before sample v: absolute value of y z: movement of v Average value (p is output)
k: Reciprocal of the number of data for which the moving average is taken (processing procedure)
The processes in steps S1 to S3 in FIG.
〔実施の形態３〕
本発明の実施の形態３に係るリアルタイム周波数解析方法は、任意の中心周波数に対する伝達特性応答倍率を任意の周波数成分の振動特性として抽出できるようにしたものであり、特定の周波数の入力信号に対する出力信号の応答倍率を計算できることから、例えば、回転機器取付け部と周辺の振動において、回転数の変化に対し、それに相当する周波数成分の応答倍率の変化を求めることなどの処理が可能となり、振動対策が必要な場合などの解析に有用となる。
具体的な処理ステップとしては、前記実施の形態１に係るリアルタイム周波数解析方法の全処理ステップと同じ処理ステップを、所定の複数個の中間周波数毎に配して並列に処理できるステップを有し、かつ前記複数並列の処理方法の各々の処理方法で用いる狭帯域周波数成分波形抽出ステップの中心周波数と、中心周波数に対する１より小さい任意の通過帯域比率とは、各々異なる値を設定することを可能とする。これにより、設定した任意の複数の中心周波数近傍の周波数成分の時々刻々の振幅の実行値に相当する値を任意の複数の周波数成分の振動特性として抽出することができる。
図３は、本発明の実施の形態３に係るリアルタイム周波数解析方法の処理手順を示すフローチャート図である。
但し、ここでは、測定環境を示す一般的な符号、及び図３に記載の符号の意味を下記のとおりとする。
[Embodiment 3]
The real-time frequency analysis method according to
As specific processing steps, the same processing steps as all the processing steps of the real-time frequency analysis method according to
FIG. 3 is a flowchart showing a processing procedure of the real-time frequency analysis method according to
However, here, the general symbols indicating the measurement environment and the meanings of the symbols described in FIG. 3 are as follows.
ｘ：使用する時系列データ
ｙ：周波数成分の波形出力
入出力の添え字は、ｎ：現在の値、ｎ−ｊ：現在よりｊ：サンプル前の値、ｉｎ：構造物への入力、ｏｔ：構造物からの出力
ｖ：ｙの絶対値
ｚ：ｖの移動平均値
ｋ：移動平均を取るデータ数の逆数
Ａａｍｐ＿ｉ ：任意の中心周波数に対する伝達特性応答倍率とする。
（処理手順）
まず、ステップＳ４の数式〔（４）式および（５）式〕を計算し、次にステップＳ５の数式〔（６）式および（７）式〕を計算し、次に、ステップＳ６の数式〔（８）式および（９）式〕を計算し、最後にステップＳ７の数式〔（１０）式〕を計算する。
x : Time series data to be used y: Waveform output of frequency component Input / output subscripts are: n: current value, nj: current value j: value before sample, in: input to structure, ot: Output from structure v: absolute value of y z: moving average value of v
k: Reciprocal number of data that takes moving average Aamp_i: Transfer characteristic response magnification with respect to an arbitrary center frequency.
(Processing procedure)
First, the equations [Equation (4) and (5)] in Step S4 are calculated, then the equations [Equation (6) and (7)] in Step S5 are calculated, and then the equations [Step 6] Equations (8) and (9)] are calculated, and finally the equation [Equation (10)] in step S7 is calculated.
〔実施の形態４〕
本発明の実施の形態４に係るリアルタイム周波数解析方法は、設定した任意の中心周波数に対する伝達特性位相差を任意の周波数成分の振動特性として抽出する方法であり、即ち、特定の周波数の入力信号に対する出力信号の位相差を計算するものであり、前述の実施の形態３と同様の用途への適用が可能であり、複数の出力点と入力点の位相関係から振動の伝達経路の推定等に役立たせることができる。
具体的な処理ステップとしては、解析対象への入力信号を検出する第１のセンサの出力データと前記第２のセンサの出力データに対して、前記実施の形態１と同じ、狭帯域周波数成分波形抽出ステップと、絶対値算出ステップと、移動平均算出ステップとを適用する入出力狭帯域周波数成分波形抽出ステップと、前記入出力狭帯域周波数成分波形抽出ステップの、前記第１のセンサと前記第２のセンサとに対応した出力の出力比を求める除算ステップと、前記除算ステップを前記第１のセンサに対応するデータに適用した結果の出力と、前記除算ステップを前記第２のセンサに対応するデータに適用した結果の出力とを乗算する伝達特性応答倍率算出ステップと、前記実施の形態１と同じ狭帯域周波数成分波形抽出ステップを前記第１のセンサの出力データに適用した結果の出力と、前記狭帯域周波数成分波形抽出ステップを前記第２のセンサの出力データに適用した結果の出力とを乗算する波形乗算ステップと、前記波形乗算ステップの出力波形に対し任意時間毎の移動平均を求める移動平均算出ステップと、前記移動平均算出ステップの算出結果の２倍を、前記伝達特性応答倍率算出ステップの算出結果で除した結果の逆余弦を求める位相算出ステップと、を有する。
[Embodiment 4]
The real-time frequency analysis method according to the fourth embodiment of the present invention is a method for extracting a transfer characteristic phase difference for a set arbitrary center frequency as a vibration characteristic of an arbitrary frequency component, that is, for an input signal having a specific frequency. This is to calculate the phase difference of the output signal, and can be applied to the same use as in the third embodiment, and is useful for estimating the vibration transmission path from the phase relationship between a plurality of output points and input points. Can be made.
As specific processing steps, the same narrowband frequency component waveform as that in the first embodiment is applied to the output data of the first sensor for detecting the input signal to the analysis target and the output data of the second sensor. The first sensor and the second of the input / output narrowband frequency component waveform extraction step, the input / output narrowband frequency component waveform extraction step applying the extraction step, the absolute value calculation step, and the moving average calculation step A division step for obtaining an output ratio of the output corresponding to the sensor of the first output, an output resulting from applying the division step to the data corresponding to the first sensor, and data corresponding to the second sensor of the division step The transfer characteristic response magnification calculating step of multiplying the output of the result applied to the above and the same narrowband frequency component waveform extracting step as in the first embodiment are performed in the first sensor. A waveform multiplying step for multiplying the output obtained by applying the output data to the output data obtained by applying the narrowband frequency component waveform extracting step to the output data of the second sensor, and an output waveform of the waveform multiplying step A moving average calculation step for obtaining a moving average for each arbitrary time, and a phase calculation for obtaining an inverse cosine of a result obtained by dividing twice the calculation result of the moving average calculation step by the calculation result of the transfer characteristic response magnification calculation step Steps.
このように構成して、前記各ステップには、前記実施の形態１と同じ中心周波数と、同じ通過帯域比率と、を用い、前記移動平均算出ステップには前記実施の形態１と同じ整数周期分の移動平均データ数を用いることで、前記構造物の設定した任意の中心周波数に対する伝達特性位相差を任意の周波数成分の振動特性として抽出することができるようにしている。
図４は、本発明の実施の形態４に係るリアルタイム周波数解析方法の処理手順を示すフローチャート図である。
但し、ここでは、測定環境を示す一般的な符号、及び図４に記載の符号の意味を下記のとおりとする。
With this configuration, the same center frequency and the same passband ratio as those in the first embodiment are used for each step, and the same average period as that in the first embodiment is used for the moving average calculation step. By using this moving average data number, the transfer characteristic phase difference with respect to an arbitrary center frequency set by the structure can be extracted as a vibration characteristic of an arbitrary frequency component.
FIG. 4 is a flowchart showing the processing procedure of the real-time frequency analysis method according to
However, here, the meanings of general symbols indicating the measurement environment and the symbols described in FIG. 4 are as follows.
ｘ：使用する時系列データ
ｙ：周波数成分の波形出力
入出力の添え字は、ｎ：現在の値、ｎ−ｊ：現在よりｊ：サンプル前の値、ｉｎ：構造物への入力、ｏｔ：構造物からの出力
ｖ：ｙの絶対値
ｚ：ｖの移動平均値
ｋ：移動平均を取るデータ数の逆数
Aamp_inot：任意の中心周波数に対する構造物の入出力の振幅成分の積
Ainot ：任意の中心周波数に対する構造物の入出力の波形の積
cos φ：中心周波数に対する伝達特性位相差の余弦
φ：中心周波数に対する伝達特性位相差；
（処理手順）
まず、ステップＳ４の数式〔（４）式と（５）式〕を計算し、次に、ステップＳ５，６の数式〔（６）式、（７）式、（８）式、（９）式〕とステップＳ８〔（１１）式〕の数式を計算し、また、この計算と並列にステップＳ９〔（１２）式〕の数式とＳ１０〔（１３）式〕の数式を計算し、最後にステップＳ１１〔（１４）式〕の数式を計算する。
x : Time series data to be used y: Waveform output of frequency component Input / output subscripts are: n: current value, nj: current value j: value before sample, in: input to structure, ot: Output from structure v: absolute value of y z: moving average value of v
k: Reciprocal of the number of data to take the moving average
Aamp_inot: The product of the amplitude components of the input and output of the structure for an arbitrary center frequency
Ainot: product of the input and output waveform of the structure for an arbitrary center frequency
cos φ: cosine of transfer characteristic phase difference with respect to center frequency φ: transfer characteristic phase difference with respect to center frequency;
(Processing procedure)
First, the equations [Equation (4) and (5)] in step S4 are calculated, and then the equations [Equations (6), (7), (8), and (9) in steps S5 and S6]. ] And the formula of step S8 [formula (11)], and the formula of step S9 [formula (12)] and the formula of S10 [formula (13)] are calculated in parallel with this calculation. S11 [Expression (14)] is calculated.
〔実施の形態５〕
本発明の実施の形態５に係るリアルタイム周波数解析方法は、狭帯域伝達特性応答倍率スペクトル抽出方法であり、即ち、複数の中心周波数についての入力信号に対する出力信号の応答倍率を計算する方法である。
この方法により、解析対象の時々刻々に伝達関数の内の応答倍率の変化を見ることが可能となるので、この方法は、例えば構造物の損傷により構造物の固有振動数の変化を検知することに適用することができる。
具体的な処理ステップとしては、前記実施の形態３に係るリアルタイム周波数解析方法の全処理ステップと同じ処理ステップを複数配して並列処理できるステップを有し、前記複数並列に処理するように配置された処理方法の各々が用いる狭帯域周波数成分波形抽出ステップの中心周波数と、該中心周波数に対する１より小さい任意の通過帯域比率とは、各々異なる値を設定することができるようにしている。これにより、設定した任意の複数の任意の中心周波数に対する伝達特性応答倍率を任意の周波数成分の振動特性として抽出することができる。
[Embodiment 5]
The real-time frequency analysis method according to the fifth embodiment of the present invention is a narrowband transfer characteristic response magnification spectrum extraction method, that is, a method of calculating the response magnification of the output signal with respect to the input signal for a plurality of center frequencies.
This method makes it possible to see the change in the response magnification of the transfer function from time to time of the analysis target, so this method can detect changes in the natural frequency of the structure due to damage to the structure, for example. Can be applied to.
As specific processing steps, there are steps that can be processed in parallel by arranging a plurality of the same processing steps as all the processing steps of the real-time frequency analysis method according to the third embodiment, and arranged so as to perform the processing in parallel. The center frequency of the narrowband frequency component waveform extraction step used by each of the processing methods and an arbitrary passband ratio smaller than 1 with respect to the center frequency can be set to different values. Thereby, it is possible to extract the transfer characteristic response magnification with respect to any set arbitrary center frequencies as vibration characteristics of an arbitrary frequency component.
図５は、本発明の実施の形態５に係るリアルタイム周波数解析方法の処理手順を示すフローチャート図である。
但し、ここでは、測定環境を示す一般的な符号、及び図５に記載の符号の意味を下記のとおりとする。
ｘ：使用する時系列データ
ｙ：周波数成分の波形出力
入出力の添え字は、ｎ：現在の値、ｎ−ｊ：現在よりｊ：サンプル前の値、ｉｎ：構造物への入力、ｏｔ：構造物からの出力
ｖ：ｙの絶対値
ｚ：ｖの移動平均値
Aamp＿i ：任意の中心周波数に対する伝達特性応答倍率
（処理手順）
図５のステップＳ４〜Ｓ７の処理を、ｉ＝１〜ｐ個並べて並列に処理する。
FIG. 5 is a flowchart showing the processing procedure of the real-time frequency analysis method according to the fifth embodiment of the present invention.
However, here, the meanings of general symbols indicating the measurement environment and the symbols described in FIG. 5 are as follows .
x : Time series data to be used y: Waveform output of frequency component Input / output subscripts are: n: current value, nj: current value j: value before sample, in: input to structure, ot: Output from structure v: absolute value of y z: moving average value of v
A amp_i: transmission characteristic response magnification to an arbitrary center frequency
( Processing procedure)
The processes of steps S4 to S7 in FIG. 5 are performed in parallel by arranging i = 1 to p.
〔実施の形態６〕
本発明の実施の形態６に係るリアルタイム周波数解析方法は、狭帯域周波数成分位相差スペクトル抽出方法であり、即ち、複数の中心周波数についての入力信号に対する出力信号の位相特性を計算するものである。
この方法を使用すると、解析対象の時々刻々に伝達関数の内の応答倍率の変化を見ることが可能となり、よって、例えば構造物の損傷により構造物の固有振動数の変化を位相特性の変化により検知することに適用することができる。
具体的な処理ステップとしては、前記実施の形態４の全ステップと同じ処理ステッブを複数配して並列処理できるようにしたステップを有し、かつ前記複数並列の処理方法の各々の処理方法で用いる狭帯域周波数成分波形抽出ステップの中心周波数と、中心周波数に対する１より小さい任意の通過帯域比率とは、各々異なる値を設定することを可能とする。これにより、設定した任意の複数の任意の中心周波数に対する伝達特性位相差を任意の周波数成分の振動特性として抽出することができる。
図６は、本発明の実施の形態６に係るリアルタイム周波数解析方法の処理手順を示すフローチャート図である。
但し、ここでは、測定環境を示す一般的な符号および図６に記載の符号の意味を下記のとおりとする。
[Embodiment 6]
The real-time frequency analysis method according to
By using this method, it is possible to see the change in the response magnification of the transfer function from time to time in the analysis target, so that the change in the natural frequency of the structure due to the damage of the structure is caused by the change in the phase characteristics, for example. It can be applied to detection.
As specific processing steps, a plurality of processing steps that are the same as all the steps of the fourth embodiment are arranged so that they can be processed in parallel, and are used in each processing method of the plurality of parallel processing methods. It is possible to set different values for the center frequency of the narrowband frequency component waveform extraction step and the arbitrary passband ratio smaller than 1 with respect to the center frequency. Thereby, the transfer characteristic phase difference with respect to the set arbitrary plural arbitrary center frequencies can be extracted as the vibration characteristic of the arbitrary frequency component.
FIG. 6 is a flowchart showing a processing procedure of the real-time frequency analysis method according to the sixth embodiment of the present invention.
However, here, the meanings of general symbols indicating the measurement environment and the symbols described in FIG. 6 are as follows.
ｘ：使用する時系列データ
ｙ：周波数成分の波形出力
入出力の添え字は、ｎ：現在の値、ｎ−ｊ：現在よりｊ：サンプル前の値、ｉｎ：構造物への入力、ｏｔ：構造物からの出力
ｖ：ｙの絶対値
ｚ：ｖの移動平均値
ｋ：移動平均を取るデータ数の逆数
Aamp＿inot：任意の中心周波数に対する構造物の入出力の振幅成分の積
Ainot ：任意の中心周波数に対する構造物の入出力の波形の積
cos φ：中心周波数に対する伝達特性位相差の余弦
φ：中心周波数に対する伝達特性位相差（ｐ個出力）
（処理手順）
図６におけるステップＳ４〜Ｓ６およびステップＳ８〜Ｓ１１の処理を、ｉ＝１〜ｐ個並べて並列に処理する。
x : Time series data to be used y: Waveform output of frequency component Input / output subscripts are: n: current value, nj: current value j: value before sample, in: input to structure, ot: Output from structure v: absolute value of y z: moving average value of v
k: Reciprocal of the number of data to take the moving average
A amp_inot: The product of the amplitude components of the input and output of the structure for an arbitrary center frequency
Ainot: product of the input and output waveform of the structure for an arbitrary center frequency
cos φ: cosine of transfer characteristic phase difference with respect to center frequency φ: transfer characteristic phase difference with respect to center frequency (p outputs )
( Processing procedure)
The processes of steps S4 to S6 and steps S8 to S11 in FIG. 6 are arranged in parallel by arranging i = 1 to p.
図７は、本発明の実施の形態２に係るリアルタイム周波数解析方法での処理を示す説明図である。
同図において、本発明の実施の形態１に係るリアルタイム周波数解析方法での処理は、一点鎖線内のブロックに纏めて示している。
図７に示す処理は、前述のとおり、中心周波数の周波数成分の時々刻々の振幅値を求める処理である。
図８は、本発明の実施の形態５に係るリアルタイム周波数解析方法での処理を示す説明図である。
同図において、本発明の実施の形態３に係るリアルタイム周波数解析方法での処理は、一点鎖線内のブロックに纏めて示している。
図８に示す処理は、前述のとおり、各々異なる中心周波数での本発明の実施の形態１に係るリアルタイム周波数解析方法の処理を纏めたものであり、各々の中心周波数の周波数成分の時々刻々の振幅値を求める処理である。
FIG. 7 is an explanatory diagram showing processing in the real-time frequency analysis method according to
In the figure, the processes in the real-time frequency analysis method according to the first embodiment of the present invention are collectively shown in blocks within a one-dot chain line.
The process shown in FIG. 7 is a process for obtaining the amplitude value of the frequency component of the center frequency every moment as described above.
FIG. 8 is an explanatory diagram showing processing in the real-time frequency analysis method according to the fifth embodiment of the present invention.
In the figure, the processes in the real-time frequency analysis method according to the third embodiment of the present invention are collectively shown in blocks within a one-dot chain line.
The process shown in FIG. 8 is a summary of the processes of the real-time frequency analysis method according to the first embodiment of the present invention at different center frequencies as described above. This is a process for obtaining an amplitude value.
図９は、本発明の実施の形態６に係るリアルタイム周波数解析方法での処理を示す説明図である。
図９において、本発明の実施の形態４に係るリアルタイム周波数解析方法での処理は一点鎖線（但し、画像変換時の画質の劣化により一点鎖線でなく、破線として見える場合もある）内のブロックに纏めて示している。
図９に示す処理は、前述のとおり、各々異なる中心周波数での本発明の実施の形態４に係るリアルタイム周波数解析方法の処理を纏めたものであり、各周波数毎の入出力の位相差を求める処理である。
なお、本発明に係るリアルタイム周波数解析方法を実現する処理の少なくとも一部をコンピュータ制御により実行するものとし、かつ、上記処理を、図１〜６のフローチャートで示した手順および（１５）式によりコンピュータに実行せしめるプログラムは、半導体メモリを始め、ＣＤ−ＲＯＭや磁気テープなどのコンピュータ読み取り可能な記録媒体に格納して配付してもよい。そして、少なくともマイクロコンピュータ，パーソナルコンピュータ，汎用コンピュータを範疇に含むコンピュータが、上記の記録媒体から上記プログラムを読み出して、実行するものとしてもよい。
FIG. 9 is an explanatory diagram showing processing in the real-time frequency analysis method according to the sixth embodiment of the present invention.
In FIG. 9, the processing in the real-time frequency analysis method according to the fourth embodiment of the present invention is performed on a block in a dashed line (however, it may appear as a broken line instead of a dashed line due to image quality deterioration during image conversion). It shows together.
The process shown in FIG. 9 summarizes the processes of the real-time frequency analysis method according to the fourth embodiment of the present invention at different center frequencies as described above, and obtains the input / output phase difference for each frequency. It is processing.
It should be noted that at least a part of the processing for realizing the real-time frequency analysis method according to the present invention is executed by computer control, and the above processing is performed by the procedure shown in the flowcharts of FIGS. The program to be executed may be distributed after being stored in a computer-readable recording medium such as a semiconductor memory, a CD-ROM, or a magnetic tape. A computer including at least a microcomputer, a personal computer, and a general-purpose computer may read the program from the recording medium and execute the program.
Ｔ：サンプリング時間間隔；データ集録時設定
ωｉ：中心周波数；任意設定
ζｉ（ツェータ）；通過帯域比率；任意設定
ｘ：使用する時系列データ
ｙ：周波数成分の波形出力
入出力の添え字は、ｎ：現在の値、ｎ−ｊ：現在よりｊ：サンプル前の値
ｖ：ｙの絶対値
ｚ：ｖの移動平均値
T: Sampling time interval; Data acquisition setting ωi: Center frequency; Arbitrary setting ζi (Zeta); Passband ratio; Arbitrary setting x: Time series data used y: Waveform output of frequency component Input / output subscript is n : Current value, nj: j from current value: value before sample v: absolute value of y z: moving average value of v
Claims (6)
任意の１つの中心周波数近傍の狭帯域の周波数成分の波形を抽出する狭帯域周波数成分波形抽出ステップと、
前記狭帯域周波数成分波形抽出ステップで抽出された波形の絶対値を求める絶対値回路の機能を有する絶対値算出ステップと、
前記絶対値算出ステップの出力絶対値波形の、任意時間毎の移動平均を求める移動平均算出ステップと、
を備え、
前記狭帯域周波数成分波形抽出ステップは、
用いる時系列データの現時点のデータおよび２サンプル前のデータと、現時点より１サンプル前および２サンプル前の当該ステップの出力から、現時点の前記当該ステップの出力を算出する機能を有し、
用いる時系列データのサンプリング時間間隔と、任意の中心周波数と、中心周波数に対する１より小さい任意の通過帯域比率とを用いて該中心周波数に対する通過帯域比率分の周波数だけ中心周波数から低い周波数側と高い周波数側にずれたカットオフ周波数と、前記任意の通過帯域比率に応じた遮断特性を有し、
前記絶対値算出ステップは、
入力された値が正値の場合はそのままの値を出力し、また、入力された値が負値の場合は負の符号を外した値を出力する機能を有し、
前記移動平均算出ステップは、
前記狭帯域周波数成分波形抽出ステップにおいて設定した任意の中心周波数に相当する周期の整数倍の時間を、前記狭帯域周波数成分振幅抽出ステップに適用するデータのサンプリング時間間隔で除したデータ数分の逆数を当該ステップに用いる時系列データの現時点のデータに乗じ、この値と、１から移動平均を取るデータ数の逆数を差し引いた値を現時点の１サンプル前の当該移動平均ステップの算出結果に乗じた値と、を足し合わせた値を当該移動平均ステップの現時点の算出結果とすることで、順次１サンプリング間隔毎の平均値を算出する機能を有し、
前記移動平均算出ステップの出力結果が前記狭帯域周波数成分波形抽出ステップにおいて設定した任意の１つの中心周波数近傍の周波数成分の時々刻々の振幅の実行値に相当する値を任意の１つの周波数成分の振動特性として抽出するに際し、掛け算の回数および一連のステップで記憶を必要とするデータの数を少なくして、演算負荷と必要記憶容量を小さく抑えるように構成したことを特徴とするリアルタイム周波数解析方法。 Acceleration, to data acquisition input time series data of the A / D converter is taken out by sampling the output of a sensor for measuring the physical quantity of the load, etc., at regular time intervals, or time-series data output to one Dan storage device A real-time frequency analysis method for extracting vibration characteristics of an arbitrary frequency component,
A narrowband frequency component waveform extraction step for extracting a waveform of a narrowband frequency component near any one center frequency;
An absolute value calculating step having a function of an absolute value circuit for obtaining an absolute value of the waveform extracted in the narrowband frequency component waveform extracting step;
A moving average calculating step for obtaining a moving average for each arbitrary time of the output absolute value waveform of the absolute value calculating step;
With
The narrowband frequency component waveform extraction step includes:
A function for calculating the current output of the step from the current data of the time-series data to be used and the data of two samples before and the output of the step one sample before and two samples before the current time;
Using the sampling time interval of the time-series data to be used, an arbitrary center frequency, and an arbitrary passband ratio smaller than 1 with respect to the center frequency, the frequency corresponding to the passband ratio with respect to the center frequency is higher on the lower frequency side than the center frequency. A cutoff frequency shifted to the frequency side, and a cutoff characteristic according to the arbitrary passband ratio,
The absolute value calculating step includes:
When the input value is a positive value, the value is output as it is, and when the input value is a negative value, it has a function of outputting a value with a negative sign removed,
The moving average calculation step includes:
The reciprocal of the number of data obtained by dividing the time that is an integral multiple of the period corresponding to the arbitrary center frequency set in the narrowband frequency component waveform extraction step by the sampling time interval of data applied to the narrowband frequency component amplitude extraction step Is multiplied by the current data of the time-series data used in this step, and this value is multiplied by the result of subtracting the reciprocal of the number of data to be taken from 1 from the calculation result of the moving average step one sample before the current time. By adding the value and the value as the current calculation result of the moving average step, it has a function of calculating an average value for each sampling interval sequentially,
The output result of the moving average calculation step is a value corresponding to the execution value of the momentary amplitude of the frequency component in the vicinity of any one center frequency set in the narrowband frequency component waveform extraction step. Real-time frequency analysis method characterized in that the number of multiplications and the number of data that need to be stored in a series of steps are reduced to reduce the calculation load and the required storage capacity when extracting as vibration characteristics .
前記入出力狭帯域周波数成分波形抽出ステップの、前記第１のセンサと前記第２のセンサとに対応した出力の出力比を求める除算ステップと、
を備え、
前記第１のセンサと前記第２のセンサとの出力に適用する、請求項１記載の前記狭帯域周波数成分波形抽出ステップと、前記絶対値算出ステップと、前記移動平均算出ステップにおける各処理方法では、同じ中心周波数と、同じ通過帯域比率と、同じ整数周期分の移動平均データ数とを参照することで、前記解析対象に設定した任意の中心周波数に対する伝達特性応答倍率を、任意の周波数成分の振動特性として抽出することを特徴とする請求項１に記載のリアルタイム周波数解析方法。 The output data of said second sensor for detecting an output signal from the first output data and the analyzed sensor for detecting an input signal to the analyzed, the narrowband frequency according to claim 1 a component waveform extraction step, the absolute value calculation step, a narrowband frequency component waveform extracting step input of applying and the moving average calculating step,
A division step of obtaining an output ratio of outputs corresponding to the first sensor and the second sensor in the input / output narrowband frequency component waveform extraction step;
Equipped with a,
Applied to the output of the previous SL first sensor and the second sensor, and wherein the narrowband frequency component waveform extracting step of claim 1, wherein said absolute value calculating step, each processing method in the moving average calculating step Then, by referring to the same center frequency, the same passband ratio, and the number of moving average data for the same integer period, the transfer characteristic response magnification with respect to the arbitrary center frequency set as the analysis target is set to the arbitrary frequency component. Real-time frequency analysis method according to claim 1, characterized that you extracted as the vibration characteristics of the.
前記入出力狭帯域周波数成分波形抽出ステップの、前記第１のセンサと前記第２のセンサとに対応した出力の出力比を求める除算ステップと、
前記除算ステップを前記第１のセンサに対応するデータに適用した結果の出力と、前記除算ステップを前記第２のセンサに対応するデータに適用した結果の出力とを乗算する伝達特性応答倍率算出ステップと、
請求項１に記載の前記狭帯域周波数成分波形抽出ステップを前記第１のセンサの出力データに適用した結果の出力と、前記狭帯域周波数成分波形抽出ステップを前記第２のセンサの出力データに適用した結果の出力とを乗算する波形乗算ステップと、
前記波形乗算ステップの出力波形に対し任意時間毎の移動平均を求める移動平均算出ステップと、
前記移動平均算出ステップの算出結果の２倍を、前記伝達特性応答倍率算出ステップの算出結果で除した結果の逆余弦を求める位相算出ステップと、
を備えたことを特徴とする請求項１に記載のリアルタイム周波数解析方法。 The output data of the first output data and the second sensor of the sensor for detecting an input signal to the analyzed, and the narrow-band frequency component waveform extracting step according to claim 1, wherein the absolute value a calculation step, a narrowband frequency component waveform extracting step input of applying and the moving average calculating step,
A division step of obtaining an output ratio of outputs corresponding to the first sensor and the second sensor in the input / output narrowband frequency component waveform extraction step;
Transfer characteristic response magnification calculation step of multiplying the output of the result obtained by applying the division step to the data corresponding to the first sensor and the output of the result of applying the division step to the data corresponding to the second sensor When,
Applying the narrowband frequency component waveform extracting step of claim 1 and an output of the result of applying the output data of the first sensor, the narrowband frequency component waveform extracting step to the output data of the second sensor A waveform multiplication step for multiplying the resulting output by
A moving average calculating step for obtaining a moving average for each arbitrary time with respect to the output waveform of the waveform multiplying step;
A phase calculation step for obtaining an inverse cosine of a result obtained by dividing twice the calculation result of the moving average calculation step by the calculation result of the transfer characteristic response magnification calculation step;
The real-time frequency analysis method according to claim 1, further comprising :
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| JP2011191146A JP5854551B2 (en) | 2011-09-01 | 2011-09-01 | Real-time frequency analysis method |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| JP2011191146A JP5854551B2 (en) | 2011-09-01 | 2011-09-01 | Real-time frequency analysis method |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| JP2013053883A JP2013053883A (en) | 2013-03-21 |
| JP5854551B2 true JP5854551B2 (en) | 2016-02-09 |
Family
ID=48131010
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| JP2011191146A Active JP5854551B2 (en) | 2011-09-01 | 2011-09-01 | Real-time frequency analysis method |
Country Status (1)
| Country | Link |
|---|---|
| JP (1) | JP5854551B2 (en) |
Families Citing this family (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| KR101674757B1 (en) | 2014-12-16 | 2016-11-10 | 전남대학교산학협력단 | Apparatus and method for estimating oscilation mode |
| JP2021067598A (en) * | 2019-10-25 | 2021-04-30 | アンリツ株式会社 | Signal analyzer and signal analysis method |
Family Cites Families (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JPH08136409A (en) * | 1994-10-28 | 1996-05-31 | Motorola Inc | Knock detection system |
| GB2434649A (en) * | 2006-01-28 | 2007-08-01 | Martin Paul Simpson | Signal analyser |
| JP4842669B2 (en) * | 2006-02-27 | 2011-12-21 | 株式会社デジタル | Anomaly detection device |
| JP4472769B2 (en) * | 2007-12-28 | 2010-06-02 | 株式会社シグネット | Real-time seismometer and prediction method of seismic intensity using it |
-
2011
- 2011-09-01 JP JP2011191146A patent/JP5854551B2/en active Active
Also Published As
| Publication number | Publication date |
|---|---|
| JP2013053883A (en) | 2013-03-21 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Gryllias et al. | Estimation of the instantaneous rotation speed using complex shifted Morlet wavelets | |
| CN107505135B (en) | Rolling bearing composite fault extraction method and system | |
| Yan et al. | Improved Hilbert–Huang transform based weak signal detection methodology and its application on incipient fault diagnosis and ECG signal analysis | |
| Chen et al. | A performance enhanced time-varying morphological filtering method for bearing fault diagnosis | |
| EP2237002B1 (en) | Method for analysing vibration in rotor blades | |
| Wang et al. | Bearing fault diagnosis under time-varying rotational speed via the fault characteristic order (FCO) index based demodulation and the stepwise resampling in the fault phase angle (FPA) domain | |
| GB2543521A (en) | Method and data processing device for severity assessment of bearing defects using vibration energy | |
| CN108287073B (en) | Resonance band selection method based on singular value component frequency domain spectrum | |
| US11927501B2 (en) | Method and device for monitoring a gear system | |
| US20130006150A1 (en) | Bruxism detection device and bruxism detection method | |
| Zhang et al. | Improved local cepstrum and its applications for gearbox and rolling bearing fault detection | |
| CN111487318B (en) | Time-varying structure instantaneous frequency extraction method | |
| Zhang et al. | Operation conditions monitoring of flood discharge structure based on variance dedication rate and permutation entropy | |
| EP2618269A1 (en) | A method for processing of measurements from several sensors | |
| JP5854551B2 (en) | Real-time frequency analysis method | |
| CN114061947B (en) | Sparse time-frequency analysis-based variable-rotation-speed fault diagnosis method and system for gearbox | |
| CN107941324A (en) | A kind of method of estimation of consumer level inertia sensing unit environments noise | |
| CN114325740A (en) | Speed measurement method and system applied to laser Doppler velocimeter | |
| CN113074941A (en) | Variable-speed gear fault signal extraction method, diagnosis method and system based on adaptive time-varying comb filtering and storage medium | |
| CN112446329A (en) | Method, system, device and storage medium for determining instantaneous frequency of time-varying structure | |
| US8995230B2 (en) | Method of extracting zero crossing data from full spectrum signals | |
| CN116701840A (en) | Method and system for optimizing and calculating cepstrum of mechanical vibration signal | |
| CN116509365A (en) | Vital sign monitoring method and device | |
| CN115266066A (en) | Method, device, equipment and storage medium for identifying fault of rotating body | |
| CN111812404B (en) | Signal processing method and processing device |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| A621 | Written request for application examination |
Free format text: JAPANESE INTERMEDIATE CODE: A621 Effective date: 20140818 |
|
| A977 | Report on retrieval |
Free format text: JAPANESE INTERMEDIATE CODE: A971007 Effective date: 20150409 |
|
| A131 | Notification of reasons for refusal |
Free format text: JAPANESE INTERMEDIATE CODE: A131 Effective date: 20150417 |
|
| A521 | Request for written amendment filed |
Free format text: JAPANESE INTERMEDIATE CODE: A523 Effective date: 20150616 |
|
| 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: 20151204 |
|
| A61 | First payment of annual fees (during grant procedure) |
Free format text: JAPANESE INTERMEDIATE CODE: A61 Effective date: 20151207 |
|
| R150 | Certificate of patent or registration of utility model |
Ref document number: 5854551 Country of ref document: JP Free format text: JAPANESE INTERMEDIATE CODE: R150 |
|
| R250 | Receipt of annual fees |
Free format text: JAPANESE INTERMEDIATE CODE: R250 |
|
| R250 | Receipt of annual fees |
Free format text: JAPANESE INTERMEDIATE CODE: R250 |
|
| R250 | Receipt of annual fees |
Free format text: JAPANESE INTERMEDIATE CODE: R250 |
|
| R250 | Receipt of annual fees |
Free format text: JAPANESE INTERMEDIATE CODE: R250 |