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】[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】[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】[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.
【0004】このような場合に、単純にスペクトルを求
めても、あるいは、スペクトルを繰り返し求めて重ね合
せることにより多数の周期に亘る平均的なスペクトルを
求めても、信号に含まれた雑音成分の大きさにより定ま
る所定のSN比のスペクトルしか求めることができない
という問題がある。本発明は、上記問題を解決し、例え
ば本来の信号レベルを上まわる程のレベルの雑音が信号
に混入している場合であっても、従来よりも十分高いS
N比のスペクトルを求めることのできるスペクトル演算
装置を提供することを目的とする。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】[0005]
【課題を解決するための手段】図1は、本発明のスペク
トル演算装置のブロック図である。本来の信号x0
(t)に雑音n(t)が混入した信号x(t)が、第1
の帯域通過フィルタ10と第2の帯域通過フィルタ12
に入力される。第1の帯域通過フィルタ10は、入力さ
れた信号x(t)のうち所定の第1の中心周波数ka の
近傍の周波数成分を抽出することにより帯域制限信号x
a (t)を出力するフィルタであり、第2の帯域通過フ
ィルタ12は、入力された信号x(t)のうち所定の第
2の中心周波数kb 近傍の周波数成分を抽出することに
より第2の帯域制限信号xb (t)を出力するフィルタ
である。ここで、第2の帯域通過フィルタ12は上記第
2の中心周波数kb を変更することが自在なフィルタで
あることが好ましい。これら第1及び第2の帯域通過フ
ィルタ10,12から出力された第1および第2の帯域
制限信号xa (t),xb (t)は、それぞれ第1およ
び第2の二乗演算手段14,16に入力される。第1の
二乗演算回路14では、第1の帯域通過フィルタ10で
抽出された第1の帯域制限信号xa (t)の二乗|x a
(t)|2 を表わす第1のパワー信号za (t)が求め
られ、第2の二乗演算手段16では、第2の帯域通過フ
ィルタ12で抽出された第2の帯域制限信号x b (t)
の二乗|xb (t)|2 を表わす第2のパワー信号zb
(t)が求められる。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.
【0006】第1および第2の二乗演算手段14,16
で求められた第1および第2のパワー信号za (t),
zb (t)は、それぞれ第1および第2のスペクトル演
算手段18,20に入力される。第1のスペクトル演算
手段18では、入力された第1のパワー信号za (t)
のスペクトルZb (k)が求められて伝達関数演算手段
22に入力され、第2のスペクトル演算手段20では、
入力された第2のパワー信号zb (t)のスペクトルZ
b (k)が求められて、やはり伝達関数演算手段22に
入力される。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.
【0007】伝達関数演算手段22では、入力された2
つのスペクトルZa (k),Zb (k)を用いて、信号
x(t)の周波数ka 成分から周波数kb 成分への伝達
関数Hab(k)が、式 Hab(k)= {E[Za *(k)・Zb (k)]/E[Za *(k)・Za (k)]}1/2 …(1) 但し、E[……]は集合平均、*は共役複素数を表わす に従って求められる。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
【0008】ここで、上記各フィルタ10,12および
各手段14〜22は、アナログ信号にフィルタリングあ
るいは各演算を施すものであってもよいが、入力された
信号x(t)をディジタル信号に変換し、このディジタ
ルの信号に基づいてフィルタリングあるいは各演算を施
すものであってもよい。また、これら各フィルタ10,
12および各手段14〜22は、ハードウェアで構成さ
れたものであってもよくあるいはコンピュータを用いソ
フトウェアで実現したものであってもよい。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】[0009]
【作用】たとえば機械装置等からピックアップされた、
雑音成分を除く本来の信号は所定の周期でほぼ同一の波
形が繰り返される周期的な信号である場合が多い。また
雑音成分は本来の信号よりもその繰り返しサイクル間で
相関が小さい場合が多い。本発明は、信号成分と雑音成
分の、この繰り返しサイクル間の相関の相違に着目し完
成されたものである。即ち、信号成分x0 (t)は相関
が高いため、第1の帯域通過フィルタ10から出力され
た第1の帯域制限信号xa (t)と第2の帯域通過フィ
ルタ12から出力された第2の帯域制限信号xb (t)
との間でも、その本来の信号x0 (t)に起因する成分
は相関が大きく、雑音成分n(t)は相関が低い。した
がってこれら第1の帯域制限信号xa (t)と第2の帯
域制限信号xb (t)との間のクロススペクトルを適切
な手法により繰り返し求め、この繰り返し求めたクロス
スペクトルの平均を求め、このクロススペクトルを用い
て第1の帯域制限信号xa (t)と第2の帯域制限信号
xb (t)との間の伝達関数を求めることにより、雑音
成分は徐々に低減化され、高いSN比をもって本来の信
号x0 (t)の周波数ka 成分に対する周波数kb の成
分の比率が求められることとなる。[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.
【0010】ここで、第1の帯域制限信号xa (t)と
第2の帯域制限信号xb (t)との間の伝達関数をどの
ようにして求めるかが問題となる。1つには、そのま
ま、第1の帯域制限信号xa (t)のスペクトルXa
(k)と第2の帯域制限信号xb (t)のスペクトルX
b (k)を求め、式 E[Xa *(k)・Xb (k)]/E[Xa *(k)・Xa (k)]…(2) 但し、E[……]は集合平均、*は共役複素数を表わす に従って信号x(t)の周波数ka 成分から周波数kb
成分への伝達関数を求めることが考えられる。しかし、
第1の帯域通過フィルタ10と第2の帯域通過フィルタ
12は、通常、信号x(t)の、互いに異なる周波数成
分を抽出するものであるため、2つのスペクトルXa
(k),Xb (k)は、通常、重なるスペクトル成分を
有していない。したがってこのような単純な手法では伝
達関数を求めることはできない。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.
【0011】また他の1つとして、第1の帯域制限信号
xa (t),第2の帯域制限信号x b (t)を通常の方
法にしたがって検波することによりそれらの信号xa
(t),xb (t)の周波数帯域を直流成分近傍に遷移
させ、これらの検波信号のスペクトル間で上記(2)式
に従う伝達関数を演算することが考えられる。この手法
では検波信号はいずれも直流成分近傍の信号であるため
それらの信号のスペクトルどうしは直流成分近傍で重な
ることとなる。本来の信号x0 (t)の繰り返し周期が
正確に常に一定であり、各周期内において発生する信号
x0 (t)の位相も正確に常に一定であるという数学的
モデルにおいては、この手法を用いることにより伝達関
数を求めることができる。しかし、実際の物理現象にお
いては、信号x0 (t)の繰り返し周期がふらつく場合
がほとんどであり、各繰り返し周期内での信号x0
(t)の位置も各繰り返し周期毎にふらつく場合がほと
んどである。このような実際の物理現象においてこの手
法に基づいて伝達関数を求め、SN比を向上させるため
に多数繰り返し求めた伝達関数の平均を求めると、雑音
n(t)に起因する成分だけでなく本来の信号x0
(t)に起因する成分も零に漸近してしまい、結局これ
ら2つの信号xa (t),xb (t)間の伝達関数を推
定することはできないこととなる。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.
【0012】そこで、本発明では、二乗演算手段14,
16を備え、帯域制限信号xa (t),xb (t)の二
乗信号|xa (t)|2 ,|xb (t)|2 を求め、こ
れによりこれらの信号の周波数帯域を直流近傍に遷移さ
せると共に位相情報を捨象し、二乗演算手段14,16
により求められたパワー信号Za (t),Zb (t)間
の伝達関数を求めることとしたものである。ここで本来
求めるべき伝達関数は、二乗前の信号xa (t),xb
(t)の間の伝達関数であるため、伝達関数演算手段2
2では、本来の伝達関数の式の平方根が求められる。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.
【0013】尚、帯域通過フィルタ10,12の通過帯
域幅を狭く設定すると周波数分解能は向上するが分散
(バラつき)が大きくなり、通過帯域幅を広く設定する
と分散(バラつき)は小さくなるが周波数分解能が低下
する。したがって、帯域通過フィルタ10,12の通過
帯域幅は、必要とする周波数分解能や、伝達関数を繰り
返し求めて平均するための伝達関数の演算回数、即ち、
許容される演算時間等に応じて適応的に定められる。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】[0014]
【実施例】以下、本発明の実施例について説明する。こ
こでは、信号x(t)=x0 (t)+n(t)が所定の
標本化周波数fs をもってサンプリングされ、離散的な
ディジタル信号x(n)=x0 (n)+n(n)(n
は、n番目の時刻におけるサンプリングであることを表
わす)を得るものとし、ここでは、信号x0 (n)の一
周期を機械系の固有振動モデル x0 (n)=Re [exp{−(α+jω0 )nTs }] …(3) Re [……]は実数部を表わす とする。ここでは、標本化周期Ts =1/fs =1、減
衰係数α=0.1、周波数ω0 =π/2とし、その固有
振動の波形x0 (n)、スペクトルの振幅|X0(k)
|、及びスペクトルの位相arg{X0 (k)}を、そ
れぞれ図2(A)、(B)、(C)に示す。ここで、固
有振動x0 (n)の共振角周波数はk0 =π/2であ
る。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 [expR− (Α + 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.
【0015】ここでは、(3)式に示される固有振動x
0 (n)を周期T0 で周期的に発生させ、さらに白色雑
音n(n)を加算して観測信号x(n)のモデルとす
る。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】[0016]
【数1】 (Equation 1)
【0017】ここでmは整数,δ(…)はデルタ関数を
表わす。図3(A)は、(4)式第1項の固有振動周期
波形、図3(B)は白色雑音n(n)が重畳された
(4)式全体の波形である。ここでは、S/N=−5d
B(即ち(4)式第1項のパワーよりも第2項のパワー
の方が5dB大きい)、T0=200Ts である。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.
【0018】ここで、固有振動x0 (n)の共振周波数
k0 成分から周波数k1 成分への伝達関数をH01(k)
とすると、前述した(1)式より、固有振動スペクトル
の振幅二乗の推定値|Y0 (k)|2 は、 |Y0 (k)|2 =|X0 (k)|2 ・|H01(k)|2 =|X0 (k)|2 ・|E[Z0 *(k)・Z1 (k)]/E[Z0 *(k)・Z0 (k)]| …(5) により求められる。ここでZ0 (k)は、観測信号x
(n)の共振周波数k0 の成分を二乗した信号のスペク
トル、Z1 (k)は、観測信号x(n)の周波数k 1 の
成分を二乗した信号のスペクトル、E[…]は、集合平
均、*は共役複素数である。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.
【0019】図4,図5,図6は、それぞれ、S/N=
0dB,−5dB,−10dBの場合のシミュレーショ
ン結果を示した図である。ここではいずれも、同期加算
の回数は500回、帯域通過フィルタの帯域幅は±0.
01fs の場合について示してある。各図中、従来と記
載されたグラフは、通常のパワースペクトル同期加算演
算による推定値|X(k)|2 、本発明と記載されたグ
ラフは、(5)式に基づく推定値|Y(k)|2 であ
る。図4〜図6に示すいずれの場合も、本発明の手法を
用いることにより、通常のパワースペクトルによる推定
値|X(k)|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.
【0020】また、S/N=−10dBの場合に関し
て、同期加算の回数による変化を示すために、図7,図
8に、加算平均回数500回の図6の場合との比較のた
めのそれぞれ加算平均回数100回,10回の場合を示
す。実際の信号波形のスペクトルを推定する際におい
て、加算平均回数をどの程度とすればよいかということ
については、コヒーレンスγ01 2 (k)を、式 γ01 2 (k)={|E[Z0 *(k)・Z1 (k)]|2 }/ {E[|Z0 *(k)|2 ]・E[|Z1 *(k)|2 ]} …(6) に従って求め、このコヒーレンスγ01 2 (k)の値に応
じて加算平均回数を定めることができる。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).
【0021】以上から、周波数帯域成分間の伝達関数H
01(k)を用いて、観測信号x(n)から固有振動のス
ペクトルX0 (k)の推定を行うことが、特にSN比が
低い場合に有効であることが示される。尚、ここでは、
本発明にいう第1の中心周波数ka の一例として共振周
波数k o を用いたが、この第1の中心周波数ka は共振
周波数と一致する必要はない。ただしスペクトルの推定
値の精度を向上させるためには、第1の中心周波数ka
として、例えば共振周波数やその近傍の周波数等、SN
比の大きな周波数を選択することが好ましい。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.
【0022】またここでは雑音成分n(n)は白色雑音
であるとしたが、白色雑音でなくてもよく、繰り返し周
期間における相関が本来観察したい信号x0 (n)より
も相関の低い雑音であれば多数回加算平均することによ
り真値に近いスペクトルを推定することができる。この
ときの加算平均回数は、例えば上記(6)式のコヒーレ
ンスを求めることにより定めることができる。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】[0023]
【発明の効果】以上説明したように、本発明のスペクト
ル演算装置は、前述した(1)式に基づいて伝達関係H
ab(k)を求める伝達関数演算手段を備えたものである
ため、従来よりも雑音の影響が大幅に低減化されたスペ
クトルを求めることができる。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.
【図1】本発明のスペクトル演算装置のブロック図であ
る。FIG. 1 is a block diagram of a spectrum calculation device according to the present invention.
【図2】固有振動の波形x0 (n)、そのスペクトルの
振幅|X0 (k)|及びそのスペクトルの位相arg
{X0 (k)}を示した図である。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)}.
【図3】固有振動の周期波形、及び白色雑音が重畳され
た波形を示した図である。FIG. 3 is a diagram showing a periodic waveform of a natural vibration and a waveform on which white noise is superimposed.
【図4】S/N=0dBの場合のシミュレーション結果
を示した図である。FIG. 4 is a diagram showing a simulation result when S / N = 0 dB.
【図5】S/N=−5dBの場合のシミュレーション結
果を示した図である。FIG. 5 is a diagram showing a simulation result when S / N = -5 dB.
【図6】S/N=−10dBの場合のシミュレーション
結果を示した図である。FIG. 6 is a diagram showing a simulation result when S / N = −10 dB;
【図7】加算平均回数100回の場合のシミュレーショ
ン結果を示した図である。FIG. 7 is a diagram showing a simulation result when the number of averages is 100;
【図8】加算平均回数10回の場合のシミュレーション
結果を示した図である。FIG. 8 is a diagram showing a simulation result when the number of averages is 10;
10 第1の帯域通過フィルタ 12 第2の帯域通過フィルタ 14 第1の二乗演算手段 16 第2の二乗演算手段 18 第1のスペクトル演算手段 20 第2のスペクトル演算手段 22 伝達関数演算手段 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)発明者 中鉢 憲賢 宮城県仙台市青葉区貝ケ森4丁目3−18 (56)参考文献 特開 平6−207957(JP,A) 特開 平5−252064(JP,A) 特開 昭63−101768(JP,A) 特開 昭59−56170(JP,A) (58)調査した分野(Int.Cl.7,DB名) 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)
し、該信号x(t)の、所定の第1の中心周波数ka 近
傍の周波数成分を抽出する第1の帯域通過フィルタと、 前記信号x(t)を入力し、該信号x(t)の、所定の
第2の中心周波数kb近傍の周波数成分を抽出する第2
の帯域通過フィルタと、 前記第1の帯域通過フィルタで抽出された第1の帯域制
限信号xa (t)の二乗|xa (t)|2 を表わす第1
のパワー信号za (t)を求める第1の二乗演算手段
と、 前記第2の帯域通過フィルタで抽出された第2の帯域制
限信号xb (t)の二乗|xb (t)|2 を表わす第2
のパワー信号zb (t)を求める第2の二乗演算手段
と、 前記第1の二乗演算手段で求められた前記第1のパワー
信号za (t)のスペクトルZa (k)を求める第1の
スペクトル演算手段と、 前記第2の二乗演算手段で求められた前記第2のパワー
信号zb (t)のスペクトルZb (k)を求める第2の
スペクトル演算手段と、 前記第1のスペクトル演算手段で求められた前記スペク
トルZa (k)と前記第2のスペクトル演算手段で求め
られた前記スペクトルZb (k)を用いて、前記信号x
(t)の周波数ka 成分から周波数kb 成分への伝達関
数Hab(k)を、式 Hab(k)={E[Za *(k)・Zb (k)]/E[Z
a *(k)・Za (k)]}1/2 但し、E[……]は集合平均、*は共役複素数を表わす に従って求める伝達関数演算手段とを備えたことを特徴
とするスペクトル演算装置。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.
2の中心周波数kbを変更自在なフィルタであることを
特徴とする請求項1記載のスペクトル演算装置。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 |