CN113820004B - A robust method for initial phase estimation of vibration signals - Google Patents
A robust method for initial phase estimation of vibration signals Download PDFInfo
- Publication number
- CN113820004B CN113820004B CN202111085988.2A CN202111085988A CN113820004B CN 113820004 B CN113820004 B CN 113820004B CN 202111085988 A CN202111085988 A CN 202111085988A CN 113820004 B CN113820004 B CN 113820004B
- Authority
- CN
- China
- Prior art keywords
- signal
- initial phase
- ind
- frequency
- phase
- 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
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01H—MEASUREMENT OF MECHANICAL VIBRATIONS OR ULTRASONIC, SONIC OR INFRASONIC WAVES
- G01H17/00—Measuring mechanical vibrations or ultrasonic, sonic or infrasonic waves, not provided for in the other groups of this subclass
Landscapes
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Measurement Of Mechanical Vibrations Or Ultrasonic Waves (AREA)
Abstract
本发明提供一种鲁棒的振动信号初始相位估计方法，首先根据原始信号进行硬阈值带通滤波，提取信号趋势部分，再利用Hilbert变换提取信号瞬时相位，利用线性插值和信号中间位置估计初始相位，降低信号边界“吉布斯”效应对估计结果的影响。本发明全过程只有一个特征频率参数，提高了相位估计算法鲁棒性；涉及的FFT、IFFT、Hilbert变换等均有快速算法，满足在线处理，可有效应用于故障诊断、设备健康管理等实时工业监测情景；通过滤波提高指定特征频率处信噪比，降低其它非特征频谱对初始相位估计的干扰；通过滤波信号长度一半附近位置估计初始相位，避免因为“吉布斯效应”产生的信号边界震荡带来的计算误差，提高了该方法的适用场景。
The present invention provides a robust initial phase estimation method for a vibration signal. First, a hard threshold bandpass filter is performed on the original signal to extract the trend part of the signal, and then the instantaneous phase of the signal is extracted using the Hilbert transform. The initial phase is estimated using linear interpolation and the middle position of the signal to reduce the influence of the "Gibbs" effect of the signal boundary on the estimation result. The present invention has only one characteristic frequency parameter in the whole process, which improves the robustness of the phase estimation algorithm; the FFT, IFFT, Hilbert transform, etc. involved all have fast algorithms to meet online processing requirements and can be effectively applied to real-time industrial monitoring scenarios such as fault diagnosis and equipment health management; the signal-to-noise ratio at the specified characteristic frequency is improved by filtering, and the interference of other non-characteristic spectra on the initial phase estimation is reduced; the initial phase is estimated by filtering the position near half the length of the signal to avoid the calculation error caused by the signal boundary oscillation caused by the "Gibbs effect", thereby improving the applicable scenarios of the method.
Description
Technical Field
The invention relates to the technical field of measurement and test, in particular to a robust vibration signal initial phase estimation method.
Background
The traditional mechanical fault diagnosis technology mainly relies on the person who has experience of working on site to analyze through the human ear and the touch equipment instrument before, and this mode needs to have a large amount of accumulation of working experience, relies on expert experience's accumulation and can not perfectly promote, and secondly has subjective factor and objective factor of people, and fault diagnosis conclusion human interference is big. Vibration signal based analysis is one of the most widely used fault monitoring methods in equipment diagnostics and predictive maintenance. The main technology is to analyze the characteristic frequency of the equipment through a time domain diagram and a frequency domain diagram to correspond to different fault types.
The time domain waveform reflects the characteristic that the amplitude and the phase of the signal change along with time, the time domain waveform can generally acquire the signal of the equipment through processing modes such as signal filtering and amplifying, the time domain signal can intuitively show the running state and the parameter quantity of the equipment at different moments, the collected vibration waveform can reflect the health state of the equipment, and the change of the vibration signal at different moments can be obtained through intuitively watching the vibration signal waveform. The vibration signal of a typical mechanical device is a non-stationary signal whose waveform of vibration is produced by rotation of multiple parts of the device
Because the time domain waveforms are single in equipment fault analysis, vibration signals of equipment parts are superimposed, and the state of the equipment at the moment cannot be intuitively analyzed, the time domain signals are required to be subjected to Fourier change, dynamic complex signals are decomposed, vibration signals with different frequencies are reflected in a frequency domain diagram, and the fault characteristic frequency of the frequency domain diagram can reflect the running state of the equipment at the moment. Time-frequency converting the signal waveform can convert the acquired waveform in one period into a plurality of frequency components. According to Fourier transform, mapping the time domain signals to a frequency domain, wherein the relation between the amplitude of the signals and the characteristic frequency is reflected in the frequency domain diagram, and the frequency domain signal analysis directly reflects the operation state of the equipment through the frequency signals without considering whether the signals are in a linear relation.
The vibration signal of the starting-up and stopping process is generally called a transient signal, which is the response of the rotor system to the change of the rotating speed and is the reflection of the dynamic characteristics and fault symptoms of the rotor. In the start-stop process of the instrument, the rotor experiences various rotating speeds, and the vibration signal is the response of the rotor system to the change of the rotating speeds, is the external expression of the dynamic characteristics and fault symptoms of the rotor, and contains rich information which is difficult to obtain at ordinary times. Therefore, start-stop process analysis is an important task for rotor detection. By analysing the amplitude and phase shift of the vibration signal of the device, the resonance frequency and critical rotational speed of the rotor system can be found. Based on the amplitude and phase at low speeds, knowledge of the degree of bending of the rotor is facilitated. The response of the rotor system to rotor imbalance surge during the actual start-up and shutdown process can be seen through the phase change. (the amplitude response of the system is a distinct formant when passing through the critical speed, while the phase changes approximately 180 degrees back and forth.)
In order to avoid vibration abnormality caused by resonance and improve the accuracy of fault diagnosis, amplitude and phase estimation is required to be carried out on transient signals (in the condition of start-up and shutdown), and by changing a machine structure, the basic rigidity is increased, so that good balance of a rotor is ensured, and the running rotating speed of a fan is enabled to avoid a resonance region.
The phase calculation method of the vibration signal is as follows:
In a Fast Fourier Transform (FFT) spectrum analysis method, when performing spectrum analysis on a discrete signal, since the discrete signal is time-domain truncated, problems such as energy leakage, a smaller spectrum amplitude, and a reduced measurement accuracy occur. In order to reduce the errors, scholars at home and abroad propose a signal windowing interpolation FFT analysis algorithm based on rectangular windows, hanning windows, blackman windows and Blackman-Harris windows, so that spectrum leakage is restrained to a certain extent, and the accuracy of signal parameter estimation is improved. By adjusting the port parameters, a window function with small side lobe peak level and large side lobe asymptotic attenuation rate is selected. The long-range leakage caused by mutual interference between adjacent frequency spectrum side lobes when the signal time domain is truncated can be effectively reduced by processing the signal.
Sinusoidal approximations were developed based entirely on the development of modern computer technology and digital signal processing technology. The sine approximation method using the least square method as the core can realize the absolute measurement of the phase frequency characteristic of the vibration sensor within the wide frequency band range of 1 Hz-10 kHz, and is particularly suitable for measuring the amplitude-phase characteristic of the vibration sensor under the condition of low frequency and large amplitude with fewer sampling periods.
The full-phase FFT spectrum analysis method can inhibit the problem of spectrum leakage to a certain extent. The analysis method analyzes the frequency spectrums of two sequences with a time delay relationship, corrects the results of the two phase spectrum analysis by using a time shift phase difference method, and can obtain accurate correction of phase and amplitude according to the known sinusoidal signal frequency.
The analog method is to use a phase measuring circuit to measure the phase of an input alternating signal, and the circuit generally relates to a plurality of high-performance components, so that the circuit has complex design, expensive components, limited measuring precision and easy interference from the surrounding environment.
The zero crossing method mainly detects the phase difference between two groups of same-frequency signals when the signal amplitude is 0, the measurement regularity of the method mainly depends on the interval size of sampling time, and the calculation result is easily interfered by external factors. When the correlation method is used for calculating the phase difference, a reference signal with the same frequency is often needed to be used as a reference, the correlation degree is calculated on the input signal and the reference signal, and then the phase difference between the input signal and the reference signal is solved.
The frequency spectrum method is to calculate the phase frequency characteristics of two groups of input signals, then calculate and obtain the phase value of the current input signal according to the peak position in the signal spectrogram to be analyzed, and when the phase difference of the two groups of signals is needed, calculate the phase of a single signal, and then subtract. The latter two detection methods show good anti-interference performance in practical measurement experiments. The spectral method can be used not only to measure the difference between two sets of signals, but also to analyze parameters of each set of input signals independently, such as: amplitude, phase, frequency, etc.
In the prior art, under the condition of mixing different frequencies and noise, the phase information extracted by the full-phase time-shift phase difference method is very accurate, but the accuracy of the amplitude is greatly reduced; the sine approximation method is more accurate in terms of extracting amplitude information, but the extraction of the phase is poorer; the accuracy of FFT spectral analysis is relatively worst. When the frequency resolution of the signal analysis is insufficient to resolve the signal frequency, the accuracy of amplitude and phase estimation is reduced. The traditional method for carrying out FFT spectrum analysis inevitably causes energy leakage, reduces amplitude, reduces precision and can not realize real whole period truncation. Although the accuracy of the measurement of the spectral amplitude is improved after the various window functions are applied to the discrete signal, the maximum error in phase is still large. The sine approximation method is mature in theory and wide in measurement range, but the sine approximation method adopts a least square method to fit the amplitude and the phase of a sine signal, and the error of the estimated quantity determined by the least square method depends on the error of measurement data and the functional relation given by an equation set.
Disclosure of Invention
The invention provides a robust vibration signal initial phase estimation method for on-line fault diagnosis, which aims to solve the problems of low robustness, poor real-time performance and the like of the traditional phase estimation method. According to the method, hard threshold bandpass filtering is carried out according to an original signal, a main trend part of the signal is extracted, then the transient phase of the signal is extracted by Hilbert transformation, the initial phase is estimated (linear interpolation) by utilizing the intermediate position of the signal, and the influence of the Gibbs effect of the signal boundary on an estimation result is reduced.
The invention provides a robust vibration signal initial phase estimation method, which comprises the steps of carrying out FFT on an original signal f (t) and adopting hard threshold band-pass filtering to obtain a filtered signal
Will filter the signalHilbert transform extraction of instantaneous phase
Will instantaneous phaseObtaining initial phase value/>, by using linear interpolation method and signal center position
The invention relates to a robust vibration signal initial phase estimation method, which comprises the following steps as a preferable mode:
S1, fourier transformation: performing FFT on the original signal F (t) to obtain a frequency spectrum F (omega);
s2: hard threshold bandpass filtering: the frequency spectrum F (omega) is subjected to hard threshold bandpass filtering to obtain a transformed frequency spectrum Transformed spectrumPerforming IFFT to obtain filtered signal
S3, hilbert transformation: will filter the signalHilbert transform is carried out to obtain Hilbert transform signalsTransform Hilbert signalsExtraction of instantaneous phase after resolution
S4, linear interpolation: constructing an analytic signal z (t) and using linear interpolation to extract instantaneous phaseAnalyzing, and obtaining an initial phase value/> according to signal intermediate position estimation
In a robust vibration signal initial phase estimation method according to the present invention, as a preferred mode, in step S4,
The analysis signal z (t) is:
wherein the polar coordinates of the resolved signal z (t) are:
A (t) is the envelope of the original signal f (t);
the expression of A (t) is:
in the robust vibration signal initial phase estimation method of the present invention, in step S4, as an optimal mode, the instantaneous phase The method comprises the following steps:
In the robust vibration signal initial phase estimation method, in the step S4, an initial phase value is obtained according to signal intermediate position estimation The specific method of (a) is as follows:
Calculating an index ind about half of the signal length, judging whether the index ind is an integer, if so, then ind epsilon [ t 1,t2,t3...,tL ], recording ind=t k, k=L/2, and initializing the phase value The method comprises the following steps: /(I)
If not, indexing the initial phase value of ind d and ind u, which are the lower neighbor integer and the upper neighbor integer of ind respectivelyThe method comprises the following steps:
the invention relates to a robust vibration signal initial phase estimation method, which is characterized in that as a preferable mode, an index ind calculation formula is as follows:
ind＝ntemp*nt，
Where n temp is the number of whole cycles around half the length of the whole signal, and nt is the number of sampling points in one cycle.
The invention relates to a robust vibration signal initial phase estimation method, which is used as a preferable mode, and the calculation formula of the whole period number n temp near half of the whole signal length is as follows:
wherein L is the number of sampling points;
the calculation formula for one period containing the sampling point number nt is as follows:
Where Fs is the sampling frequency and f 0 is the characteristic frequency.
In the robust vibration signal initial phase estimation method, as a preferred mode, in step S1, the FFT formula is:
wherein x (t) is the signal point of the original signal f (t),
f(t)＝f(t1,t2,t3...,tL)＝[x1,x2,x3...,xL].
In the robust vibration signal initial phase estimation method, in the step S2, as a preferred mode, the hard threshold bandpass filtering is as follows: setting the frequency spectrum F (omega) to 0 outside the range of the [0.75 x F 0,1.25*f0 ] area to obtain a converted frequency spectrumWherein f 0 is a characteristic frequency, i.e., a frequency component to be retained;
The IFFT is:
Wherein, For the filtered signalSignal points of (2)
In the robust vibration signal initial phase estimation method, as an optimal mode, in step S3, hilbert is transformed into:
fourier Transform (FFT) and inverse transform (IFFT)
The Fourier analysis is utilized to realize the mutual conversion between the time domain and the frequency domain of the vibration signal, so that the research of the original time domain signal can be converted into the research of the Fourier coefficient on the frequency domain. In the field of signal processing, fourier transformation plays an important role, has a milestone meaning, and is regarded as a bridge between the time domain and the frequency domain of a signal. For signal x (t), its continuous Fourier transform (FFT) is
Its inverse transform (IFFT) is:
In practical applications, signals tend to be discrete in the time domain and the frequency domain, so that Discrete Fourier Transform (DFT) is commonly used, and in view of the operation speed and the system consumption, fast Fourier Transform (FFT) is widely used in frequency domain analysis of signals.
Hard threshold bandpass filtering
The invention calculates the FFT spectrum of the signal, carries out hard threshold processing on a high frequency spectrum region, removes signal noise by utilizing the FFT inverse transformation IFFT, improves the signal-to-noise ratio of the signal and improves the accuracy of the initial phase estimation.
Given a signal F (t), where t e (- + -infinity, + -infinity) whose fourier transform is F (ω), we assume that we should preserve the frequency component as part F 0, let the spectrum F (ω) be 0 outside the [0.75 x F 0,1.25*f0 ] region (hard thresholding), note that the transformed spectrum isForPerforming an inverse FFT (i.e., IFFT) to compute the filtered signalHilbert transform
Setting a real value signal f (t), wherein t epsilon (- +, +), its Hilbert transform (Hilbert) is defined as:
Phase offset
Recording deviceFourier transformWhere ω is the frequency, according to the theorem F [ F (t) ×g (t) ]=f [ F (t) ] F [ g (t) ],
So that the real value signal f (t) is subjected to Hilbert transform with unchanged amplitude but withTo generate/>, for the positive frequency part of the original signalPhase shift, generating/>, on the negative frequency part of the original signalPhase shift, while the original signal f (t) is shifted from its Hilbert transform (Hilbert)Are orthogonal.
Resolving signals
Definition of resolved signalsWhereinIs the hilbert transform of the real-valued signal f (t),
The z (t) polar coordinate is expressed as: a (t) is the envelope of the signal f (t)/> Is the instantaneous phase (angle) of the signal f (t):
Instantaneous frequency is defined as the derivative of phase:
Linear interpolation
Linear interpolation is a simpler interpolation method, whose interpolation function is a first order polynomial. Let the values of the function y=f (x) at the two points x 0,x1 be y 0,y1, respectively, and calculate the polynomialMake satisfyCalculated to obtain
From the values of the two points, the value of any point between the two points can be estimated by linear interpolation.
Scheme flow
In order to effectively perform initial phase estimation on a vibration signal and improve the subsequent fault diagnosis effect, the invention firstly improves the signal to noise ratio through hard threshold bandpass filtering, then constructs an analysis signal through Hilbert transformation, and calculates the instantaneous phase of an original signal through the analysis signal. In order to overcome the influence of 'spectrum leakage' on boundary part signals, the signal center nearby position is utilized to estimate the initial trigger position of the signals, and the specific scheme flow is as follows: assuming that the signal f (t) has a length of l=2048 and a sampling frequency of Fs＝2.56KHz： f(t)＝f(t1,t2,t3...,tL)＝[x1,x2,x3...,xL].
1. Calculating the Fourier transform F (omega) of the original signal F (x)
2. Setting a hard threshold filtering frequency range according to the rotating speed of equipment: in this example, assuming that the rotation speed of the device is 3000 rpm, the corresponding 1-fold frequency (1X) of the device is 3000/60=50 Hz, 2-fold frequency (2X) is 3000/60×2=100 Hz, … -fold frequency (8X) is 3000/60×8=400 Hz; if we need to estimate the initial phase of 1 multiple (1X: F 0 =50 HZ), let the spectrum F (ω) be 0 outside the [0.75×f 0,1.25*f0 ] region (hard thresholding), note that the transformed spectrum isFor a pair ofPerforming an inverse FFT (i.e., IFFT) to compute the filtered signal
3. Calculating a filtered signalHilbert transform
4. Structure analysis signalThe z (t) polar coordinate is expressed as: /(I)A (t) is the envelope of the signal f (t)Is the instantaneous phase (in angle units) of the signal f (t)
5. Calculating a period includes sampling pointsTaking the number of whole periods around half the length of the whole signal
6. Calculating an index around half the signal length: ind=n temp nt, if ind is an integer
Ind e [ t 1,t2,t3...,tL ], note ind = t k, where k e [1,2,3,..l ]; k=l/2 then the estimated initial phase is: if ind is not an integer, note that the ind lower neighbor integer and the upper neighbor integer are ind d and ind u, respectively, the estimated initial phase is,
The invention has the following advantages:
(1) The whole process of the method only has one parameter of the characteristic frequency f 0, and other intermediate parameters are self-adaptive. The robustness of the phase estimation algorithm is improved.
(2) The method has rapid algorithm related to FFT, IFFT, hilbert transformation and the like in the calculation process, meets the requirement of online processing, and can be effectively applied to real-time industrial monitoring situations such as fault diagnosis, equipment health management and the like.
(3) The method improves the signal-to-noise ratio at the appointed characteristic frequency by filtering, and reduces the interference of other non-characteristic frequency spectrums on the initial phase estimation
(4) The method estimates the initial phase through the position near half of the length of the filtering signal, and avoids calculation errors caused by signal boundary oscillation generated by 'Gibbs effect'.
Drawings
FIG. 1 is a flow chart of a robust vibration signal initial phase estimation method;
FIG. 2 is a diagram of a signal sample A1 single frequency signal of a robust vibration signal initial phase estimation method;
FIG. 3 is a diagram of a signal sample A1 FFT spectrum of a robust vibration signal initial phase estimation method;
FIG. 4 is a schematic diagram of a robust vibration signal initial phase estimation method signal sample A1;
FIG. 5 is a graph of a signal sample A2 single frequency signal of a robust vibration signal initial phase estimation method;
FIG. 6 is a graph of a signal sample A2FFT spectrum of a robust vibration signal initial phase estimation method;
FIG. 7 is a graph of the transient phase of a signal sample A2 of a robust vibration signal initial phase estimation method;
FIG. 8 is a graph of a signal sample A3 single frequency signal of a robust vibration signal initial phase estimation method;
FIG. 9 is a graph of a signal sample A3 FFT spectrum of a robust vibration signal initial phase estimation method;
fig. 10 is a graph of the transient phase of a signal sample A3 of a robust vibration signal initial phase estimation method.
Detailed Description
The following description of the embodiments of the present invention will be made clearly and completely with reference to the accompanying drawings, in which it is apparent that the embodiments described are only some embodiments of the present invention, but not all embodiments.
Example 1
A robust vibration signal initial phase estimation method comprises the steps of performing FFT on an original signal f (t) and then performing hard threshold band-pass filtering to obtain a filtered signal
Will filter the signalHilbert transform extraction of instantaneous phase
Will instantaneous phaseObtaining initial phase value/>, by using linear interpolation method and signal center position
The method comprises the following steps:
S1, fourier transformation: performing FFT on the original signal F (t) to obtain a frequency spectrum F (omega);
the FFT formula is:
wherein x (t) is the signal point of the original signal f (t),
f(t)＝f(t1,t2,t3...,tL)＝[x1,x2,x3...,xL];
S2: hard threshold bandpass filtering: the frequency spectrum F (omega) is subjected to hard threshold bandpass filtering to obtain a transformed frequency spectrumTransformed spectrumPerforming IFFT to obtain filtered signal
The hard threshold bandpass filtering is: setting the frequency spectrum F (omega) to 0 outside the range of the [0.75 x F 0,1.25*f0 ] area to obtain a converted frequency spectrumWherein f 0 is a characteristic frequency, i.e., a frequency component to be retained;
The IFFT is:
Wherein, For the filtered signalSignal points of (2)
S3, hilbert transformation: will filter the signalHilbert transform is carried out to obtain Hilbert transform signalsTransform Hilbert signalsExtraction of instantaneous phase after resolution
The Hilbert transform is:
S4, linear interpolation: constructing an analytic signal z (t) and using linear interpolation to extract instantaneous phase Analyzing, and obtaining an initial phase value/> according to signal intermediate position estimation
The analysis signal z (t) is:
wherein the polar coordinates of the resolved signal z (t) are:
A (t) is the envelope of the original signal f (t);
the expression of A (t) is:
Instantaneous phase The method comprises the following steps:
obtaining initial phase value according to signal intermediate position estimation The specific method of (a) is as follows:
Calculating an index ind about half of the signal length, judging whether the index ind is an integer, if so, then ind epsilon [ t 1,t2,t3...,tL ], recording ind=t k, k=L/2, and initializing the phase value The method comprises the following steps: /(I)
If not, indexing the initial phase value of ind d and ind u, which are the lower neighbor integer and the upper neighbor integer of ind respectivelyThe method comprises the following steps:
the index ind is calculated as follows:
ind＝ntemp*nt，
Wherein n temp is the number of whole periods around half the length of the whole signal, and nt is the number of sampling points in one period;
the calculation formula of the whole period number n temp around half the whole signal length is as follows:
wherein L is the number of sampling points;
the calculation formula for one period containing the sampling point number nt is as follows:
Where Fs is the sampling frequency and f 0 is the characteristic frequency.
Example 2
The method for estimating the initial phase of the robust vibration signal adopts an intelligent operation and maintenance big data cloud platform of a space flight intelligent control (Beijing) monitoring technology limited company to collect real-time data, and the real-time data are respectively used for a single frequency signal A1 (figure 2:25Hz, sampling frequency Fs=2560 Hz, sampling point number L=2048 and initial phase 45 °); single frequency signal plus random noise A2 (fig. 5:25Hz, sampling frequency fs=2560 Hz, sampling point number l=2048, initial phase 45 °, addition of random noise); the random noise A3 is added to the mixture of the plurality of frequency signals (fig. 8: f1=25 Hz, initial phase 45 °; f2=50 Hz, initial phase 60 °; f3=100 Hz, initial phase 30 °; sampling frequency fs=2560 Hz, sampling point number l=2048), and the random noise is added thereto)
1. The signal sample A1 is analyzed online as shown in fig. 2, and the FFT spectrum is calculated as shown in fig. 3. The present signal does not add noise and ignores the hard threshold bandpass filtering step of the scheme. Estimation of its instantaneous phase by Hilbert transform As shown in FIG. 4, the initial phase was estimated to be 44.99 DEG (error 0.01 DEG) from linear interpolation of positions around half the A1 length.
2. On-line analysis of the signal sample A2 is shown in fig. 5, the calculated FFT spectrum is shown in fig. 6, the A2 is relatively large in noise-interference ratio, and the f= [18.75,31.25] partial frequency is maintained through hard threshold band-pass filtering. Estimation of its instantaneous phase by Hilbert transform As shown in FIG. 7, the initial phase was estimated to be 44.67 DEG (error 0.33 DEG) from the position around half the A2 length.
3. On-line analysis of the signal sample A3 is shown in fig. 8, the calculated FFT spectrum is shown in fig. 9, the A3 is relatively noisy, the initial phase of the f2=50 Hz frequency component is specified to be calculated, and the f= [43.75, 56.25] partial frequency is maintained through hard threshold bandpass filtering. Estimation of its instantaneous phase by Hilbert transform As shown in FIG. 10, the initial phase was estimated to be 60.8 DEG (error 0.8 DEG) from the position around half the A3 length.
The foregoing is only a preferred embodiment of the present invention, but the scope of the present invention is not limited thereto, and any person skilled in the art, who is within the scope of the present invention, should make equivalent substitutions or modifications according to the technical scheme of the present invention and the inventive concept thereof, and should be covered by the scope of the present invention.
Claims (4)
1. A robust vibration signal initial phase estimation method is characterized in that: after FFT, the original signal f (t) is subjected to hard threshold band-pass filtering to obtain a filtered signal
The filtered signal is processedHilbert transform extraction of instantaneous phase
-Phase of the instantaneous phaseObtaining initial phase value/>, by using linear interpolation method and signal center position
The method comprises the following steps:
s1, fourier transformation: performing FFT on the original signal F (t) to obtain a frequency spectrum F (omega);
S2: hard threshold bandpass filtering: the frequency spectrum F (omega) is subjected to hard threshold band-pass filtering to obtain a converted frequency spectrum And then transforming the transformed spectrumPerforming IFFT to obtain the filtered signal
S3, hilbert transformation: the filtered signal is processedHilbert transform is carried out to obtain Hilbert transform signalsTransforming the Hilbert transform signalExtracting the instantaneous phase/>, after parsing
S4, linear interpolation: constructing an analytic signal z (t) and using linear interpolation to interpolate the instantaneous phaseAnalyzing, and estimating according to the signal intermediate position to obtain the initial phase value
The analysis signal z (t) is:
Wherein the polar coordinates of the resolved signal z (t) are:
a (t) is the envelope of the original signal f (t);
the expression of A (t) is:
The instantaneous phase The method comprises the following steps:
obtaining the initial phase value according to the signal intermediate position estimation The specific method of (a) is as follows:
Calculating an index ind about half of the signal length and judging whether the index ind is an integer, if so, then ind epsilon [ t 1,t2,t3...,tL ], recording ind=t k, k=L/2, and the initial phase value The method comprises the following steps: /(I)
If not, indexing the initial phase value into ind d and ind u, respectively, which are the lower neighbor integer and the upper neighbor integer of indThe method comprises the following steps:
the index ind has the following calculation formula:
ind＝ntemp*nt，
Wherein n temp is the number of whole periods around half the length of the whole signal, and nt is the number of sampling points in one period;
the calculation formula of the whole period number n temp around half the whole signal length is as follows:
wherein L is the number of sampling points;
the calculation formula for one period containing the sampling point number nt is as follows:
Where Fs is the sampling frequency and f 0 is the characteristic frequency.
2.A method of robust vibration signal initial phase estimation according to claim 1, characterized by:
In step S1, the FFT formula is:
Wherein x (t) is the signal point of the original signal f (t),
f(t)＝f(t1,t2,t3...,tL)＝[x1,x2,x3...,xL].
3. A method of robust vibration signal initial phase estimation according to claim 1, characterized by:
In step S2, the hard threshold bandpass filtering is: setting the frequency spectrum F (omega) to 0 outside the range of the [0.75 x F 0,1.25*f0 ] area to obtain the converted frequency spectrum Wherein f 0 is a characteristic frequency, i.e., a frequency component to be retained;
The IFFT is:
Wherein x to (t) are the filtered signals Signal points of (2)
4. A method of robust vibration signal initial phase estimation according to claim 1, characterized by:
in step S3, hilbert transform is:
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202111085988.2A CN113820004B (en) | 2021-09-16 | 2021-09-16 | A robust method for initial phase estimation of vibration signals |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202111085988.2A CN113820004B (en) | 2021-09-16 | 2021-09-16 | A robust method for initial phase estimation of vibration signals |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN113820004A CN113820004A (en) | 2021-12-21 |
| CN113820004B true CN113820004B (en) | 2024-05-28 |
Family
ID=78914688
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202111085988.2A Active CN113820004B (en) | 2021-09-16 | 2021-09-16 | A robust method for initial phase estimation of vibration signals |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN113820004B (en) |
Families Citing this family (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN114357368A (en) * | 2021-12-31 | 2022-04-15 | 核工业西南物理研究院 | A Nonlinear Mode-Mode Coupling Analysis Method Based on Hilbert Transform |
| CN117168780A (en) * | 2023-08-08 | 2023-12-05 | 国网浙江省电力有限公司电力科学研究院 | GIS circuit breaker spring operating mechanism mechanical performance early warning method and system |
Citations (8)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| KR20100022691A (en) * | 2008-08-20 | 2010-03-03 | 한국표준과학연구원 | Method and apparatus for determining phase sensitivity of an accelerometer based on an analysis of the harmonic components of the interference signal |
| CN101825722A (en) * | 2010-03-25 | 2010-09-08 | 清华大学 | Robust method for estimating instantaneous frequency of seismic signal |
| CN102721462A (en) * | 2012-06-14 | 2012-10-10 | 西安交通大学 | Method for quickly computing Bode plot and Nyquist plot of rotary mechanical vehicle starting and parking processes |
| CN105068571A (en) * | 2015-08-26 | 2015-11-18 | 中国工程物理研究院总体工程研究所 | Multi-dimensional sinusoidal vibration control method and control apparatus |
| CN105738696A (en) * | 2016-04-18 | 2016-07-06 | 天津大学 | Frequency estimation method and device for all-phase time-shift phase difference |
| WO2020183489A1 (en) * | 2019-03-08 | 2020-09-17 | INDIAN INSTITUTE OF TECHNOLOGY MADRAS (IIT Madras) | Method for impedance measurement using multiple phase shifted chirp signals |
| CN111693136A (en) * | 2020-05-20 | 2020-09-22 | 南京航空航天大学 | Acoustic surface wave resonator frequency estimation algorithm adopting echo signal autocorrelation phase spectrum |
| CN113125179A (en) * | 2021-03-11 | 2021-07-16 | 同济大学 | Keyless phase order tracking method for rotating speed fluctuation of rotary machine |
Family Cites Families (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US10263624B2 (en) * | 2017-06-27 | 2019-04-16 | Intel IP Corporation | Phase synchronization between two phase locked loops |
-
2021
- 2021-09-16 CN CN202111085988.2A patent/CN113820004B/en active Active
Patent Citations (8)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| KR20100022691A (en) * | 2008-08-20 | 2010-03-03 | 한국표준과학연구원 | Method and apparatus for determining phase sensitivity of an accelerometer based on an analysis of the harmonic components of the interference signal |
| CN101825722A (en) * | 2010-03-25 | 2010-09-08 | 清华大学 | Robust method for estimating instantaneous frequency of seismic signal |
| CN102721462A (en) * | 2012-06-14 | 2012-10-10 | 西安交通大学 | Method for quickly computing Bode plot and Nyquist plot of rotary mechanical vehicle starting and parking processes |
| CN105068571A (en) * | 2015-08-26 | 2015-11-18 | 中国工程物理研究院总体工程研究所 | Multi-dimensional sinusoidal vibration control method and control apparatus |
| CN105738696A (en) * | 2016-04-18 | 2016-07-06 | 天津大学 | Frequency estimation method and device for all-phase time-shift phase difference |
| WO2020183489A1 (en) * | 2019-03-08 | 2020-09-17 | INDIAN INSTITUTE OF TECHNOLOGY MADRAS (IIT Madras) | Method for impedance measurement using multiple phase shifted chirp signals |
| CN111693136A (en) * | 2020-05-20 | 2020-09-22 | 南京航空航天大学 | Acoustic surface wave resonator frequency estimation algorithm adopting echo signal autocorrelation phase spectrum |
| CN113125179A (en) * | 2021-03-11 | 2021-07-16 | 同济大学 | Keyless phase order tracking method for rotating speed fluctuation of rotary machine |
Non-Patent Citations (1)
| Title |
|---|
| 基于DFT相位的正弦波频率和初相的高精度估计方法;齐国清, 贾欣乐;《电子学报》;20010925(第09期);第1164-1167页 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN113820004A (en) | 2021-12-21 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN105547698B (en) | The method for diagnosing faults and device of rolling bearing | |
| US8121813B2 (en) | System and method for clearance estimation between two objects | |
| CN110763462B (en) | A time-varying vibration signal fault diagnosis method based on synchronous compression operator | |
| CN103576002B (en) | A kind of computing method of capacitive insulator arrangement dielectric loss angle | |
| CN110987438B (en) | Method for detecting periodical vibration impact signals of hydraulic generator in variable rotating speed process | |
| CN102353500B (en) | Extraction method of unbalanced signal for dynamic balance measurement | |
| CN113820004B (en) | A robust method for initial phase estimation of vibration signals | |
| Rodopoulos et al. | A parametric approach for the estimation of the instantaneous speed of rotating machinery | |
| CN112465068A (en) | Rotating equipment fault feature extraction method based on multi-sensor data fusion | |
| US10365297B2 (en) | System and method for generation of a tachometer signal and reduction of jitter | |
| CN115791155B (en) | A method for estimating the instantaneous speed of a gearbox based on short-time Fourier transform spectrum | |
| Lv et al. | Generalized synchroextracting-based stepwise demodulation transform and its application to fault diagnosis of rotating machinery | |
| CN114487589B (en) | Power grid broadband signal adaptive measurement method, device and system | |
| CN108398260B (en) | A Fast Evaluation Method of Gearbox Instantaneous Angular Velocity Based on Mixed Probabilistic Method | |
| CN112556930B (en) | A data quality calculation method for vibration signals of helicopter moving parts | |
| CN117686232A (en) | Method, device and storage medium for extracting vibration fundamental frequency of gas turbine in real time | |
| CN116383629B (en) | Method for diagnosing faults of variable-rotation-speed rolling bearing | |
| CN110376437B (en) | Order analysis method for overcoming non-order frequency component interference | |
| CN120200514A (en) | A fast rotor imbalance adjustment method based on Kalman filtering | |
| CN105467209B (en) | A kind of new metal oxide arrester leakage current analysis method | |
| CN119291294A (en) | Iterative Vold-Kalman filtering method for detecting faults in rotating machinery based on current signals | |
| CN114486252B (en) | Rolling bearing fault diagnosis method of vector mode maximum envelope | |
| CN114383718B (en) | A high-frequency blade passing frequency extraction method based on the vibration signal of the outer casing of the gas turbine | |
| CN114923690A (en) | High-precision bearing fault characteristic frequency estimation diagnosis method | |
| CN116989885A (en) | A time-frequency representation method and system for bearing vibration signals containing harmonics and pulses |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant | ||
| EE01 | Entry into force of recordation of patent licensing contract |
Application publication date: 20211221 Assignee: CHINA TECHNOLOGY EXCHANGE Co.,Ltd. Assignor: Aerospace Intelligent Control (Beijing) Monitoring Technology Co.,Ltd. Contract record no.: X2024110000044 Denomination of invention: A robust initial phase estimation method for vibration signals Granted publication date: 20240528 License type: Exclusive License Record date: 20241107 |
|
| EE01 | Entry into force of recordation of patent licensing contract | ||
| PE01 | Entry into force of the registration of the contract for pledge of patent right |
Denomination of invention: A robust initial phase estimation method for vibration signals Granted publication date: 20240528 Pledgee: CHINA TECHNOLOGY EXCHANGE Co.,Ltd. Pledgor: Aerospace Intelligent Control (Beijing) Monitoring Technology Co.,Ltd. Registration number: Y2024110000392 |
|
| PE01 | Entry into force of the registration of the contract for pledge of patent right | ||
| EC01 | Cancellation of recordation of patent licensing contract |
Contract record no.: X2024110000044 Date of cancellation: 20260114 |
|
| EC01 | Cancellation of recordation of patent licensing contract | ||
| PC01 | Cancellation of the registration of the contract for pledge of patent right |
Granted publication date: 20240528 Pledgee: CHINA TECHNOLOGY EXCHANGE Co.,Ltd. Pledgor: Aerospace Intelligent Control (Beijing) Monitoring Technology Co.,Ltd. Registration number: Y2024110000392 |
|
| PC01 | Cancellation of the registration of the contract for pledge of patent right |