JP3036985B2 - Spectrum calculation device - Google Patents
Spectrum calculation deviceInfo
- Publication number
- JP3036985B2 JP3036985B2 JP4219417A JP21941792A JP3036985B2 JP 3036985 B2 JP3036985 B2 JP 3036985B2 JP 4219417 A JP4219417 A JP 4219417A JP 21941792 A JP21941792 A JP 21941792A JP 3036985 B2 JP3036985 B2 JP 3036985B2
- Authority
- JP
- Japan
- Prior art keywords
- signal
- spectrum
- band
- frequency
- component
- 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.)
- Expired - Fee Related
Links
Landscapes
- Measurement Of Mechanical Vibrations Or Ultrasonic Waves (AREA)
- Complex Calculations (AREA)
Description
【０００１】[0001]
【産業上の利用分野】本発明は、雑音を含む信号を演算
処理することにより、該信号の、雑音の影響が低減化さ
れたスペクトルを求めるスペクトル演算装置に関する。BACKGROUND OF THE INVENTION 1. Field of the Invention The present invention relates to a spectrum calculation device for calculating a signal containing noise by reducing the influence of the noise on the signal by performing arithmetic processing.
【０００２】[0002]
【従来の技術】得られた信号のスペクトルを求め、その
スペクトルに基づいて分析、解析等を行う手法が広範な
分野で採用されている。例えば機械装置が稼動している
場合に、その機械装置の不良等に起因して正常稼動時に
は存在しない振動や音が発生する場合があり、この振動
や音等をピックアップすることによる故障診断が行われ
ている。2. Description of the Related Art Techniques for obtaining a spectrum of an obtained signal and performing analysis and analysis based on the spectrum have been adopted in a wide range of fields. For example, when a mechanical device is operating, vibration or sound that does not exist during normal operation may be generated due to a defect or the like of the mechanical device, and failure diagnosis is performed by picking up the vibration or sound. Have been done.
【０００３】[0003]
【発明が解決しようとする課題】ところがスペクトル分
析のための信号をピックアップする際、往往にして大き
な雑音が混入した信号しか得られない場合がある。例え
ば、機械装置等の不良に起因する振動や音のほか正常な
稼動に起因する振動や音、あるいは対象としている機械
装置以外の周囲にある種々の装置により発生された振動
や音も同時に検出されてしまう場合がある。However, when picking up a signal for spectrum analysis, there is a case where only a signal mixed with large noise may be obtained. For example, vibrations and sounds caused by malfunctions of mechanical devices and the like, vibrations and sounds caused by normal operation, and vibrations and sounds generated by various peripheral devices other than the target mechanical device are simultaneously detected. In some cases.
【０００４】このような場合に、単純にスペクトルを求
めても、あるいは、スペクトルを繰り返し求めて重ね合
せることにより多数の周期に亘る平均的なスペクトルを
求めても、信号に含まれた雑音成分の大きさにより定ま
る所定のＳＮ比のスペクトルしか求めることができない
という問題がある。本発明は、上記問題を解決し、例え
ば本来の信号レベルを上まわる程のレベルの雑音が信号
に混入している場合であっても、従来よりも十分高いＳ
Ｎ比のスペクトルを求めることのできるスペクトル演算
装置を提供することを目的とする。In such a case, even if a spectrum is simply obtained, or an average spectrum over a large number of periods is obtained by repeatedly obtaining and superimposing the spectra, the noise component included in the signal is obtained. There is a problem that only a spectrum having a predetermined SN ratio determined by the size can be obtained. The present invention solves the above-described problem. For example, even when noise of a level higher than the original signal level is mixed in the signal, the S signal is sufficiently higher than the conventional one.
It is an object of the present invention to provide a spectrum calculation device capable of obtaining a spectrum of N ratio.
【０００５】[0005]
【課題を解決するための手段】図１は、本発明のスペク
トル演算装置のブロック図である。本来の信号ｘ0
（ｔ）に雑音ｎ（ｔ）が混入した信号ｘ（ｔ）が、第１
の帯域通過フィルタ１０と第２の帯域通過フィルタ１２
に入力される。第１の帯域通過フィルタ１０は、入力さ
れた信号ｘ（ｔ）のうち所定の第１の中心周波数ｋa の
近傍の周波数成分を抽出することにより帯域制限信号ｘ
a （ｔ）を出力するフィルタであり、第２の帯域通過フ
ィルタ１２は、入力された信号ｘ（ｔ）のうち所定の第
２の中心周波数ｋb 近傍の周波数成分を抽出することに
より第２の帯域制限信号ｘb （ｔ）を出力するフィルタ
である。ここで、第２の帯域通過フィルタ１２は上記第
２の中心周波数ｋb を変更することが自在なフィルタで
あることが好ましい。これら第１及び第２の帯域通過フ
ィルタ１０，１２から出力された第１および第２の帯域
制限信号ｘa （ｔ），ｘb （ｔ）は、それぞれ第１およ
び第２の二乗演算手段１４，１６に入力される。第１の
二乗演算回路１４では、第１の帯域通過フィルタ１０で
抽出された第１の帯域制限信号ｘa （ｔ）の二乗｜ｘ a
（ｔ）｜2 を表わす第１のパワー信号ｚa （ｔ）が求め
られ、第２の二乗演算手段１６では、第２の帯域通過フ
ィルタ１２で抽出された第２の帯域制限信号ｘ b （ｔ）
の二乗｜ｘb （ｔ）｜2 を表わす第２のパワー信号ｚb
（ｔ）が求められる。FIG. 1 is a block diagram of the present invention.
FIG. 2 is a block diagram of a torque calculation device. Original signal x0
The signal x (t) obtained by mixing the noise n (t) into (t) is the first signal x (t).
Band pass filter 10 and second band pass filter 12
Is input to The first band pass filter 10 receives the input signal.
Predetermined first center frequency k of the obtained signal x (t)a of
By extracting nearby frequency components, the band-limited signal x
a (T), and a second band-pass filter.
The filter 12 outputs a predetermined signal out of the input signal x (t).
Center frequency k of 2b To extract nearby frequency components
The second band-limited signal xb Filter that outputs (t)
It is. Here, the second band-pass filter 12 is
Center frequency k of 2b With a filter that can be changed
Preferably, there is. These first and second band-pass filters
First and second bands output from filters 10 and 12
Limit signal xa (T), xb (T) is the first and
And the second square calculation means 14 and 16. First
In the square operation circuit 14, the first band-pass filter 10
The extracted first band-limited signal xa Square of (t) | x a
(T) |Two A first power signal za (T)
In the second square calculating means 16, the second band-pass filter is used.
The second band-limited signal x extracted by the filter 12 b (T)
Square of | xb (T) |Two A second power signal zb
(T) is required.
【０００６】第１および第２の二乗演算手段１４，１６
で求められた第１および第２のパワー信号ｚa （ｔ），
ｚb （ｔ）は、それぞれ第１および第２のスペクトル演
算手段１８，２０に入力される。第１のスペクトル演算
手段１８では、入力された第１のパワー信号ｚa （ｔ）
のスペクトルＺb （ｋ）が求められて伝達関数演算手段
２２に入力され、第２のスペクトル演算手段２０では、
入力された第２のパワー信号ｚb （ｔ）のスペクトルＺ
b （ｋ）が求められて、やはり伝達関数演算手段２２に
入力される。First and second square calculating means 14, 16
The first and second power signals z a (t),
z b (t) is input to the first and second spectrum calculation means 18 and 20, respectively. In the first spectrum calculating means 18, the input first power signal z a (t)
Spectrum Z b (k) is inputted to the transfer function calculating unit 22 sought, in the second spectral computing unit 20,
Spectrum Z of input second power signal z b (t)
b (k) is obtained and also input to the transfer function calculating means 22.
【０００７】伝達関数演算手段２２では、入力された２
つのスペクトルＺa （ｋ），Ｚb （ｋ）を用いて、信号
ｘ（ｔ）の周波数ｋa 成分から周波数ｋb 成分への伝達
関数Ｈab（ｋ）が、式 Ｈab（ｋ）＝ ｛Ｅ［Ｚa *（ｋ）・Ｚb （ｋ）］／Ｅ［Ｚa *（ｋ）・Ｚa （ｋ）］｝1/2 …（１） 但し、Ｅ［……］は集合平均、＊は共役複素数を表わす に従って求められる。In the transfer function calculating means 22, the input 2
Using the two spectra Z a (k) and Z b (k), the transfer function H ab (k) from the frequency k a component to the frequency k b component of the signal x (t) is represented by the equation H ab (k) = {E [Z a * (k ) · Z b (k)] / E [Z a * (k) · Z a (k)]} 1/2 ... (1) where, E [......] is set average , * Are determined according to
【０００８】ここで、上記各フィルタ１０，１２および
各手段１４〜２２は、アナログ信号にフィルタリングあ
るいは各演算を施すものであってもよいが、入力された
信号ｘ（ｔ）をディジタル信号に変換し、このディジタ
ルの信号に基づいてフィルタリングあるいは各演算を施
すものであってもよい。また、これら各フィルタ１０，
１２および各手段１４〜２２は、ハードウェアで構成さ
れたものであってもよくあるいはコンピュータを用いソ
フトウェアで実現したものであってもよい。Here, the filters 10 and 12 and the means 14 to 22 may filter or perform an operation on an analog signal, but convert the input signal x (t) into a digital signal. Then, filtering or each operation may be performed based on this digital signal. Each of these filters 10,
The unit 12 and each of the units 14 to 22 may be configured by hardware, or may be realized by software using a computer.
【０００９】[0009]
【作用】たとえば機械装置等からピックアップされた、
雑音成分を除く本来の信号は所定の周期でほぼ同一の波
形が繰り返される周期的な信号である場合が多い。また
雑音成分は本来の信号よりもその繰り返しサイクル間で
相関が小さい場合が多い。本発明は、信号成分と雑音成
分の、この繰り返しサイクル間の相関の相違に着目し完
成されたものである。即ち、信号成分ｘ0 （ｔ）は相関
が高いため、第１の帯域通過フィルタ１０から出力され
た第１の帯域制限信号ｘa （ｔ）と第２の帯域通過フィ
ルタ１２から出力された第２の帯域制限信号ｘb （ｔ）
との間でも、その本来の信号ｘ0 （ｔ）に起因する成分
は相関が大きく、雑音成分ｎ（ｔ）は相関が低い。した
がってこれら第１の帯域制限信号ｘa （ｔ）と第２の帯
域制限信号ｘb （ｔ）との間のクロススペクトルを適切
な手法により繰り返し求め、この繰り返し求めたクロス
スペクトルの平均を求め、このクロススペクトルを用い
て第１の帯域制限信号ｘa （ｔ）と第２の帯域制限信号
ｘb （ｔ）との間の伝達関数を求めることにより、雑音
成分は徐々に低減化され、高いＳＮ比をもって本来の信
号ｘ0 （ｔ）の周波数ｋa 成分に対する周波数ｋb の成
分の比率が求められることとなる。[Action] For example, picked up from a machine or the like,
The original signal excluding the noise component is often a periodic signal in which substantially the same waveform is repeated at a predetermined cycle. The noise component often has a smaller correlation between its repetition cycles than the original signal. The present invention has been completed by paying attention to the difference in correlation between signal components and noise components between the repetition cycles. That is, since the signal component x 0 (t) has high correlation, the first band-limited signal x a (t) output from the first band-pass filter 10 and the second band-pass signal output from the second band-pass filter 12 are used. 2 band-limited signal x b (t)
Also, the component due to the original signal x 0 (t) has a large correlation, and the noise component n (t) has a low correlation. Therefore, a cross spectrum between the first band limited signal x a (t) and the second band limited signal x b (t) is repeatedly obtained by an appropriate method, and an average of the repeatedly obtained cross spectrum is obtained. By using this cross spectrum to determine the transfer function between the first band-limited signal x a (t) and the second band-limited signal x b (t), the noise component is gradually reduced and the noise component is increased. so that the ratio of the component of the frequency k b is determined for the frequency k a component of the original signal x 0 (t) with the SN ratio.
【００１０】ここで、第１の帯域制限信号ｘa （ｔ）と
第２の帯域制限信号ｘb （ｔ）との間の伝達関数をどの
ようにして求めるかが問題となる。１つには、そのま
ま、第１の帯域制限信号ｘa （ｔ）のスペクトルＸa
（ｋ）と第２の帯域制限信号ｘb （ｔ）のスペクトルＸ
b （ｋ）を求め、式 Ｅ［Ｘa *（ｋ）・Ｘb （ｋ）］／Ｅ［Ｘa *（ｋ）・Ｘa （ｋ）］…（２） 但し、Ｅ［……］は集合平均、＊は共役複素数を表わす に従って信号ｘ（ｔ）の周波数ｋa 成分から周波数ｋb
成分への伝達関数を求めることが考えられる。しかし、
第１の帯域通過フィルタ１０と第２の帯域通過フィルタ
１２は、通常、信号ｘ（ｔ）の、互いに異なる周波数成
分を抽出するものであるため、２つのスペクトルＸa
（ｋ），Ｘb （ｋ）は、通常、重なるスペクトル成分を
有していない。したがってこのような単純な手法では伝
達関数を求めることはできない。Here, how to determine a transfer function between the first band-limited signal x a (t) and the second band-limited signal x b (t) becomes a problem. For one, the spectrum X a of the first band-limited signal x a (t) as it is
(K) and spectrum X of second band-limited signal x b (t)
b (k) is calculated, and the equation E [ Xa * (k) .Xb (k)] / E [ Xa * (k) .Xa (k)] (2) where E [...] ensemble average is * frequency k b from the frequency k a component of the signal x (t) in accordance represents complex conjugate
It is conceivable to determine the transfer function to the component. But,
Since the first band-pass filter 10 and the second band-pass filter 12 usually extract different frequency components of the signal x (t), two spectra X a are used.
(K) and X b (k) usually do not have overlapping spectral components. Therefore, the transfer function cannot be obtained by such a simple method.
【００１１】また他の１つとして、第１の帯域制限信号
ｘa （ｔ），第２の帯域制限信号ｘ b （ｔ）を通常の方
法にしたがって検波することによりそれらの信号ｘa
（ｔ），ｘb （ｔ）の周波数帯域を直流成分近傍に遷移
させ、これらの検波信号のスペクトル間で上記（２）式
に従う伝達関数を演算することが考えられる。この手法
では検波信号はいずれも直流成分近傍の信号であるため
それらの信号のスペクトルどうしは直流成分近傍で重な
ることとなる。本来の信号ｘ0 （ｔ）の繰り返し周期が
正確に常に一定であり、各周期内において発生する信号
ｘ0 （ｔ）の位相も正確に常に一定であるという数学的
モデルにおいては、この手法を用いることにより伝達関
数を求めることができる。しかし、実際の物理現象にお
いては、信号ｘ0 （ｔ）の繰り返し周期がふらつく場合
がほとんどであり、各繰り返し周期内での信号ｘ0
（ｔ）の位置も各繰り返し周期毎にふらつく場合がほと
んどである。このような実際の物理現象においてこの手
法に基づいて伝達関数を求め、ＳＮ比を向上させるため
に多数繰り返し求めた伝達関数の平均を求めると、雑音
ｎ（ｔ）に起因する成分だけでなく本来の信号ｘ0
（ｔ）に起因する成分も零に漸近してしまい、結局これ
ら２つの信号ｘa （ｔ），ｘb （ｔ）間の伝達関数を推
定することはできないこととなる。As another one, a first band-limited signal
xa (T), the second band-limited signal x b (T) is normal
The signals xa
(T), xb Transition the frequency band of (t) to near DC component
And the above equation (2)
It is conceivable to calculate a transfer function according to This technique
Then, the detection signals are all signals near the DC component.
The spectra of these signals overlap near the DC component.
The Rukoto. Original signal x0 The repetition period of (t) is
A signal that is exactly constant at all times and occurs within each cycle
x0 Mathematical that the phase of (t) is also exactly constant at all times
In the model, the transfer function is
You can find the number. However, actual physical phenomena
And the signal x0 When the repetition period of (t) fluctuates
And the signal x within each repetition period0
The position of (t) often fluctuates in each repetition cycle.
It is. This kind of real physical phenomenon
To find the transfer function based on the method and to improve the signal-to-noise ratio
The average of the transfer function obtained many times in
n (t) as well as the original signal x0
The component due to (t) also approaches zero, and eventually this
Two signals xa (T), xb Estimate the transfer function between (t)
Cannot be determined.
【００１２】そこで、本発明では、二乗演算手段１４，
１６を備え、帯域制限信号ｘa （ｔ），ｘb （ｔ）の二
乗信号｜ｘa （ｔ）｜2 ，｜ｘb （ｔ）｜2 を求め、こ
れによりこれらの信号の周波数帯域を直流近傍に遷移さ
せると共に位相情報を捨象し、二乗演算手段１４，１６
により求められたパワー信号Ｚa （ｔ），Ｚb （ｔ）間
の伝達関数を求めることとしたものである。ここで本来
求めるべき伝達関数は、二乗前の信号ｘa （ｔ），ｘb
（ｔ）の間の伝達関数であるため、伝達関数演算手段２
２では、本来の伝達関数の式の平方根が求められる。Therefore, in the present invention, the square operation means 14,
16 to obtain squared signals | x a (t) | 2 , | x b (t) | 2 of the band-limited signals x a (t) and x b (t), whereby the frequency bands of these signals are determined. The signal is shifted to the vicinity of DC and the phase information is discarded.
The transfer function between the power signals Z a (t) and Z b (t) obtained by the following equation is obtained. Here, the transfer function to be originally obtained is the signal x a (t), x b before the square.
(T), the transfer function calculating means 2
In 2, the square root of the original transfer function equation is determined.
【００１３】尚、帯域通過フィルタ１０，１２の通過帯
域幅を狭く設定すると周波数分解能は向上するが分散
（バラつき）が大きくなり、通過帯域幅を広く設定する
と分散（バラつき）は小さくなるが周波数分解能が低下
する。したがって、帯域通過フィルタ１０，１２の通過
帯域幅は、必要とする周波数分解能や、伝達関数を繰り
返し求めて平均するための伝達関数の演算回数、即ち、
許容される演算時間等に応じて適応的に定められる。When the pass band widths of the band-pass filters 10 and 12 are set to be narrow, the frequency resolution is improved but the dispersion (variation) is increased. When the pass band width is set wide, the dispersion (variation) is reduced but the frequency resolution is reduced. Decrease. Therefore, the pass bandwidths of the bandpass filters 10 and 12 are determined by the required frequency resolution and the number of calculations of the transfer function for repeatedly calculating and averaging the transfer function, that is,
It is determined adaptively according to the allowable calculation time and the like.
【００１４】[0014]
【実施例】以下、本発明の実施例について説明する。こ
こでは、信号ｘ（ｔ）＝ｘ0 （ｔ）＋ｎ（ｔ）が所定の
標本化周波数ｆs をもってサンプリングされ、離散的な
ディジタル信号ｘ（ｎ）＝ｘ0 （ｎ）＋ｎ（ｎ）（ｎ
は、ｎ番目の時刻におけるサンプリングであることを表
わす）を得るものとし、ここでは、信号ｘ0 （ｎ）の一
周期を機械系の固有振動モデル ｘ0 （ｎ）＝Ｒe ［ｅｘｐ｛−（α＋ｊω0 ）ｎＴs ｝］ …（３） Ｒe ［……］は実数部を表わす とする。ここでは、標本化周期Ｔs ＝１／ｆs ＝１、減
衰係数α＝０．１、周波数ω0 ＝π／２とし、その固有
振動の波形ｘ0 （ｎ）、スペクトルの振幅｜Ｘ0（ｋ）
｜、及びスペクトルの位相ａｒｇ｛Ｘ0 （ｋ）｝を、そ
れぞれ図２（Ａ）、（Ｂ）、（Ｃ）に示す。ここで、固
有振動ｘ0 （ｎ）の共振角周波数はｋ0 ＝π／２であ
る。Embodiments of the present invention will be described below. Here, the signal x (t) = x 0 ( t) + n (t) is sampled with a predetermined sampling frequency f s, a discrete digital signal x (n) = x 0 ( n) + n (n) ( n
Represents that the sampling is performed at the n-th time). Here, one cycle of the signal x 0 (n) is defined as a natural vibration model x 0 (n) = R e [expＲ− (Α + jω 0 ) nT s }] (3) Let Re [...] represent a real part. Here, the sampling period T s = 1 / f s = 1, the attenuation coefficient alpha = 0.1, and the frequency ω 0 = π / 2, the natural vibration waveform x 0 (n), the amplitude spectrum | X 0 (K)
| And the phase arg {X 0 (k)} of the spectrum are shown in FIGS. 2 (A), (B) and (C), respectively. Here, the resonance angular frequency of the natural vibration x 0 (n) is k 0 = π / 2.
【００１５】ここでは、（３）式に示される固有振動ｘ
0 （ｎ）を周期Ｔ0 で周期的に発生させ、さらに白色雑
音ｎ（ｎ）を加算して観測信号ｘ（ｎ）のモデルとす
る。Here, the natural vibration x expressed by the equation (3)
0 (n) is periodically generated with a period T 0 , and white noise n (n) is added to form a model of the observation signal x (n).
【００１６】[0016]
【数１】 (Equation 1)
【００１７】ここでｍは整数，δ（…）はデルタ関数を
表わす。図３（Ａ）は、（４）式第１項の固有振動周期
波形、図３（Ｂ）は白色雑音ｎ（ｎ）が重畳された
（４）式全体の波形である。ここでは、Ｓ／Ｎ＝−５ｄ
Ｂ（即ち（４）式第１項のパワーよりも第２項のパワー
の方が５ｄＢ大きい）、Ｔ0＝２００Ｔs である。Here, m is an integer, and δ (...) represents a delta function. FIG. 3A shows the natural oscillation period waveform of the first term of the expression (4), and FIG. 3B shows the entire waveform of the expression (4) on which the white noise n (n) is superimposed. Here, S / N = -5d
B (i.e., (4) than the power of the first term is more power in the second term large 5dB formula), is T 0 = 200T s.
【００１８】ここで、固有振動ｘ0 （ｎ）の共振周波数
ｋ0 成分から周波数ｋ1 成分への伝達関数をＨ01（ｋ）
とすると、前述した（１）式より、固有振動スペクトル
の振幅二乗の推定値｜Ｙ0 （ｋ）｜2 は、 ｜Ｙ0 （ｋ）｜2 ＝｜Ｘ0 （ｋ）｜2 ・｜Ｈ01（ｋ）｜2 ＝｜Ｘ0 （ｋ）｜2 ・｜Ｅ［Ｚ0 *（ｋ）・Ｚ1 （ｋ）］／Ｅ［Ｚ0 *（ｋ）・Ｚ0 （ｋ）］｜ …（５） により求められる。ここでＺ0 （ｋ）は、観測信号ｘ
（ｎ）の共振周波数ｋ0 の成分を二乗した信号のスペク
トル、Ｚ1 （ｋ）は、観測信号ｘ（ｎ）の周波数ｋ 1 の
成分を二乗した信号のスペクトル、Ｅ［…］は、集合平
均、＊は共役複素数である。Here, the natural vibration x0 (N) resonance frequency
k0 From component to frequency k1 The transfer function to the component is H01(K)
Then, from the above-mentioned equation (1), the natural vibration spectrum
Estimated value of the amplitude square of | Y0 (K) |Two Is | Y0 (K) |Two = | X0 (K) |Two ・ | H01(K) |Two = | X0 (K) |Two ・ | E [Z0 *(K) · Z1 (K)] / E [Z0 *(K) · Z0 (K)] | (5) Where Z0 (K) is the observation signal x
(N) resonance frequency k0 Of the signal obtained by squaring the components of
Torr, Z1 (K) is the frequency k of the observation signal x (n) 1 of
The spectrum of the signal whose components have been squared, E [...] is the set
And * are conjugate complex numbers.
【００１９】図４，図５，図６は、それぞれ、Ｓ／Ｎ＝
０ｄＢ，−５ｄＢ，−１０ｄＢの場合のシミュレーショ
ン結果を示した図である。ここではいずれも、同期加算
の回数は５００回、帯域通過フィルタの帯域幅は±０．
０１ｆs の場合について示してある。各図中、従来と記
載されたグラフは、通常のパワースペクトル同期加算演
算による推定値｜Ｘ（ｋ）｜2 、本発明と記載されたグ
ラフは、（５）式に基づく推定値｜Ｙ（ｋ）｜2 であ
る。図４〜図６に示すいずれの場合も、本発明の手法を
用いることにより、通常のパワースペクトルによる推定
値｜Ｘ（ｋ）｜2 に比較して、はるかに真値に近い推定
値が得られている。FIGS. 4, 5, and 6 show S / N =
It is a figure showing a simulation result in the case of 0dB, -5dB, and -10dB. In each case, the number of synchronous additions is 500, and the bandwidth of the band-pass filter is ± 0.
It is shown for the case of 01f s. In each figure, the graph described as conventional is an estimated value | X (k) | 2 obtained by ordinary power spectrum synchronous addition operation, and the graph described as the present invention is an estimated value | Y ( k) | 2 . In either case shown in FIGS. 4 to 6, by using the technique of the present invention, the estimated value by normal power spectrum | X (k) |, compared to 2 to give a much estimated value close to the true value Have been.
【００２０】また、Ｓ／Ｎ＝−１０ｄＢの場合に関し
て、同期加算の回数による変化を示すために、図７，図
８に、加算平均回数５００回の図６の場合との比較のた
めのそれぞれ加算平均回数１００回，１０回の場合を示
す。実際の信号波形のスペクトルを推定する際におい
て、加算平均回数をどの程度とすればよいかということ
については、コヒーレンスγ01 2 （ｋ）を、式 γ01 2 （ｋ）＝｛｜Ｅ［Ｚ0 *（ｋ）・Ｚ1 （ｋ）］｜2 ｝／ ｛Ｅ［｜Ｚ0 *（ｋ）｜2 ］・Ｅ［｜Ｚ1 *（ｋ）｜2 ］｝ …（６） に従って求め、このコヒーレンスγ01 2 （ｋ）の値に応
じて加算平均回数を定めることができる。FIGS. 7 and 8 show, for the case of S / N = −10 dB, a comparison with the case of FIG. 6 in which the number of averaging is 500 times, in order to show a change due to the number of synchronous additions. The case where the average number of times is 100 or 10 is shown. In estimating the spectrum of the actual signal waveform, the coherence γ 01 2 (k) is calculated by the equation γ 01 2 (k) = ｛| E [Z 0 * (k) · Z 1 (k)] | 2 ｝ / {E [| Z 0 * (k) | 2 ] · E [| Z 1 * (k) | 2 ]} (6) The number of times of averaging can be determined according to the value of the coherence γ 01 2 (k).
【００２１】以上から、周波数帯域成分間の伝達関数Ｈ
01（ｋ）を用いて、観測信号ｘ（ｎ）から固有振動のス
ペクトルＸ0 （ｋ）の推定を行うことが、特にＳＮ比が
低い場合に有効であることが示される。尚、ここでは、
本発明にいう第１の中心周波数ｋa の一例として共振周
波数ｋ o を用いたが、この第１の中心周波数ｋa は共振
周波数と一致する必要はない。ただしスペクトルの推定
値の精度を向上させるためには、第１の中心周波数ｋa
として、例えば共振周波数やその近傍の周波数等、ＳＮ
比の大きな周波数を選択することが好ましい。From the above, the transfer function H between the frequency band components
01Using (k), the frequency of the natural vibration from the observation signal x (n) is calculated.
Vector X0 Performing the estimation of (k), especially when the SN ratio is
It is shown to be effective when low. Here,
The first center frequency k according to the present inventiona As an example of resonance
Wave number k o But the first center frequency ka Is resonance
It does not need to match the frequency. However, spectrum estimation
To improve the accuracy of the value, the first center frequency ka
For example, SN such as a resonance frequency or a frequency near the resonance frequency
It is preferable to select a frequency having a large ratio.
【００２２】またここでは雑音成分ｎ（ｎ）は白色雑音
であるとしたが、白色雑音でなくてもよく、繰り返し周
期間における相関が本来観察したい信号ｘ0 （ｎ）より
も相関の低い雑音であれば多数回加算平均することによ
り真値に近いスペクトルを推定することができる。この
ときの加算平均回数は、例えば上記（６）式のコヒーレ
ンスを求めることにより定めることができる。Although the noise component n (n) is assumed to be white noise here, the noise component need not be white noise, and the correlation between the repetition periods is a noise having a lower correlation than the signal x 0 (n) to be originally observed. If so, a spectrum close to the true value can be estimated by performing averaging many times. The average number of additions at this time can be determined, for example, by obtaining the coherence of the above equation (6).
【００２３】[0023]
【発明の効果】以上説明したように、本発明のスペクト
ル演算装置は、前述した（１）式に基づいて伝達関係Ｈ
ab（ｋ）を求める伝達関数演算手段を備えたものである
ため、従来よりも雑音の影響が大幅に低減化されたスペ
クトルを求めることができる。As described above, according to the spectrum calculation device of the present invention, the transmission relation H is calculated based on the aforementioned equation (1).
Since a transfer function calculating means for obtaining ab (k) is provided, it is possible to obtain a spectrum in which the influence of noise is significantly reduced as compared with the conventional case.
【図１】本発明のスペクトル演算装置のブロック図であ
る。FIG. 1 is a block diagram of a spectrum calculation device according to the present invention.
【図２】固有振動の波形ｘ0 （ｎ）、そのスペクトルの
振幅｜Ｘ0 （ｋ）｜及びそのスペクトルの位相ａｒｇ
｛Ｘ0 （ｋ）｝を示した図である。FIG. 2 shows the waveform x 0 (n) of the natural vibration, the amplitude | X 0 (k) | of the spectrum, and the phase arg of the spectrum
It is a view showing a {X 0 (k)}.
【図３】固有振動の周期波形、及び白色雑音が重畳され
た波形を示した図である。FIG. 3 is a diagram showing a periodic waveform of a natural vibration and a waveform on which white noise is superimposed.
【図４】Ｓ／Ｎ＝０ｄＢの場合のシミュレーション結果
を示した図である。FIG. 4 is a diagram showing a simulation result when S / N = 0 dB.
【図５】Ｓ／Ｎ＝−５ｄＢの場合のシミュレーション結
果を示した図である。FIG. 5 is a diagram showing a simulation result when S / N = -5 dB.
【図６】Ｓ／Ｎ＝−１０ｄＢの場合のシミュレーション
結果を示した図である。FIG. 6 is a diagram showing a simulation result when S / N = −10 dB;
【図７】加算平均回数１００回の場合のシミュレーショ
ン結果を示した図である。FIG. 7 is a diagram showing a simulation result when the number of averages is 100;
【図８】加算平均回数１０回の場合のシミュレーション
結果を示した図である。FIG. 8 is a diagram showing a simulation result when the number of averages is 10;
１０ 第１の帯域通過フィルタ １２ 第２の帯域通過フィルタ １４ 第１の二乗演算手段 １６ 第２の二乗演算手段 １８ 第１のスペクトル演算手段 ２０ 第２のスペクトル演算手段 ２２ 伝達関数演算手段 DESCRIPTION OF SYMBOLS 10 1st band pass filter 12 2nd band pass filter 14 1st square calculation means 16 2nd square calculation means 18 1st spectrum calculation means 20 2nd spectrum calculation means 22 Transfer function calculation means
───────────────────────────────────────────────────── フロントページの続き (72)発明者 中鉢 憲賢 宮城県仙台市青葉区貝ケ森４丁目３−18 (56)参考文献 特開 平６−207957（ＪＰ，Ａ) 特開 平５−252064（ＪＰ，Ａ) 特開 昭63−101768（ＪＰ，Ａ) 特開 昭59−56170（ＪＰ，Ａ) (58)調査した分野(Int.Cl.7，ＤＢ名) G01R 23/16 G01H 17/00 G06F 17/10 ──────────────────────────────────────────────────続 き Continuation of the front page (72) Inventor Kenken Nakabachi 4-3-1-18 Kaigamori, Aoba-ku, Sendai, Miyagi Prefecture (56) References JP-A-6-207957 (JP, A) JP-A-5-252064 ( JP, A) JP-A-63-101768 (JP, A) JP-A-59-56170 (JP, A) (58) Fields investigated (Int. Cl. 7 , DB name) G01R 23/16 G01H 17/00 G06F 17/10
Claims (2)
し、該信号ｘ（ｔ）の、所定の第１の中心周波数ｋa 近
傍の周波数成分を抽出する第１の帯域通過フィルタと、 前記信号ｘ（ｔ）を入力し、該信号ｘ（ｔ）の、所定の
第２の中心周波数ｋb近傍の周波数成分を抽出する第２
の帯域通過フィルタと、 前記第１の帯域通過フィルタで抽出された第１の帯域制
限信号ｘa （ｔ）の二乗｜ｘa （ｔ）｜2 を表わす第１
のパワー信号ｚa （ｔ）を求める第１の二乗演算手段
と、 前記第２の帯域通過フィルタで抽出された第２の帯域制
限信号ｘb （ｔ）の二乗｜ｘb （ｔ）｜2 を表わす第２
のパワー信号ｚb （ｔ）を求める第２の二乗演算手段
と、 前記第１の二乗演算手段で求められた前記第１のパワー
信号ｚa （ｔ）のスペクトルＺa （ｋ）を求める第１の
スペクトル演算手段と、 前記第２の二乗演算手段で求められた前記第２のパワー
信号ｚb （ｔ）のスペクトルＺb （ｋ）を求める第２の
スペクトル演算手段と、 前記第１のスペクトル演算手段で求められた前記スペク
トルＺa （ｋ）と前記第２のスペクトル演算手段で求め
られた前記スペクトルＺb （ｋ）を用いて、前記信号ｘ
（ｔ）の周波数ｋa 成分から周波数ｋb 成分への伝達関
数Ｈab（ｋ）を、式 Ｈab（ｋ）＝｛Ｅ［Ｚa *（ｋ）・Ｚb （ｋ）］／Ｅ［Ｚ
a *（ｋ）・Ｚa （ｋ）］｝1/2 但し、Ｅ［……］は集合平均、＊は共役複素数を表わす に従って求める伝達関数演算手段とを備えたことを特徴
とするスペクトル演算装置。1. A type noise n a signal x (t) containing (t), the signal x of (t), a first band pass for extracting a predetermined frequency component of the first center frequency k a near A filter, a second component for receiving the signal x (t) and extracting a frequency component of the signal x (t) near a predetermined second center frequency k b
And a first | x a (t) | 2 representing the square of the first band limited signal x a (t) extracted by the first band pass filter.
First square calculating means for obtaining the power signal z a (t) of the second band-limited signal x b (t) extracted by the second band-pass filter | x b (t) | 2 The second representing
A second square calculating means for obtaining the power signal z b (t) of the first power signal; and a second square calculating means for obtaining the spectrum Z a (k) of the first power signal z a (t) obtained by the first square calculating means. (1) spectrum calculating means, second spectrum calculating means for obtaining a spectrum Z b (k) of the second power signal z b (t) obtained by the second square calculating means, Using the spectrum Z a (k) obtained by the spectrum calculation means and the spectrum Z b (k) obtained by the second spectrum calculation means, the signal x
The transfer function H ab (k) from the frequency k a component to the frequency k b component of (t) is expressed by the equation H ab (k) = ｛E [Z a * (k) · Z b (k)] / E [ Z
a * (k) · Z a (k)]｝ 1/2, where E [...] is a set average, and * is a conjugate complex number. apparatus.
２の中心周波数ｋbを変更自在なフィルタであることを
特徴とする請求項１記載のスペクトル演算装置。Wherein said second band-pass filter, the spectrum calculation device according to claim 1, characterized in that the said second center frequency k b freely change filter.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| JP4219417A JP3036985B2 (en) | 1992-08-18 | 1992-08-18 | Spectrum calculation device |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| JP4219417A JP3036985B2 (en) | 1992-08-18 | 1992-08-18 | Spectrum calculation device |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| JPH0666853A JPH0666853A (en) | 1994-03-11 |
| JP3036985B2 true JP3036985B2 (en) | 2000-04-24 |
Family
ID=16735071
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| JP4219417A Expired - Fee Related JP3036985B2 (en) | 1992-08-18 | 1992-08-18 | Spectrum calculation device |
Country Status (1)
| Country | Link |
|---|---|
| JP (1) | JP3036985B2 (en) |
Families Citing this family (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| WO2006093152A1 (en) * | 2005-02-28 | 2006-09-08 | Pioneer Corporation | Characteristics measuring device and characteristics measuring program |
-
1992
- 1992-08-18 JP JP4219417A patent/JP3036985B2/en not_active Expired - Fee Related
Also Published As
| Publication number | Publication date |
|---|---|
| JPH0666853A (en) | 1994-03-11 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| US7099417B2 (en) | Trace video filtering using wavelet de-noising techniques | |
| JP3338370B2 (en) | Frequency analysis method and swept spectrum analyzer using the method | |
| JP2001137243A (en) | Method and apparatus for adaptive wall filtering process in spectrum doppler ultrasonic imaging | |
| EP0451955A2 (en) | Swept signal analysis instrument and method | |
| US7176670B2 (en) | Method and apparatus for zero-mixing spectrum analysis with Hilbert transform | |
| JPH0389174A (en) | Electronic measuring apparatus and estimation of frequency | |
| Matania et al. | Algorithms for spectrum background estimation of non-stationary signals | |
| AU2005289186B2 (en) | Method and device for analysing the spectrum in several frequency ranges having different resolutions | |
| JP3036985B2 (en) | Spectrum calculation device | |
| JPH0252971B2 (en) | ||
| JPH0616782B2 (en) | Method and apparatus for combining continuous estimated signals in real time | |
| JP2002541692A (en) | Digital phase sensitive rectification of signals from AC driven transducers | |
| JP7062302B2 (en) | Filtering device and filtering method | |
| JP4077092B2 (en) | Doppler frequency measurement method and Doppler sonar | |
| JP4350336B2 (en) | Resolution filter for spectrum analyzer | |
| US7424404B2 (en) | Method for determining the envelope curve of a modulated signal in time domain | |
| KR101918559B1 (en) | Method and apparatus for measuring power using maximum instantaneous power and minimum instantaneous power | |
| JP4344356B2 (en) | Detector, method, program, recording medium | |
| JP4572536B2 (en) | Sampling type measuring device | |
| JPH05264335A (en) | Method and apparatus for analyzing rotational order ratio | |
| JP3642834B2 (en) | Ultrasonic Doppler diagnostic device | |
| JP3389127B2 (en) | Radar signal processing apparatus and signal processing method for radar signal processing apparatus | |
| JPS5858625B2 (en) | signal processing device | |
| JP2957572B1 (en) | Earthquake response spectrum calculator | |
| JP2914979B2 (en) | Frequency converter |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| A01 | Written decision to grant a patent or to grant a registration (utility model) |
Free format text: JAPANESE INTERMEDIATE CODE: A01 Effective date: 20000208 |
|
| 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 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20090225 Year of fee payment: 9 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20100225 Year of fee payment: 10 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20110225 Year of fee payment: 11 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20110225 Year of fee payment: 11 |
|
| FPAY | Renewal fee payment (event date is renewal date of database) |
Free format text: PAYMENT UNTIL: 20120225 Year of fee payment: 12 |
|
| LAPS | Cancellation because of no payment of annual fees |