CN113820004A - Robust vibration signal initial phase estimation method - Google Patents
Robust vibration signal initial phase estimation method Download PDFInfo
- Publication number
- CN113820004A CN113820004A CN202111085988.2A CN202111085988A CN113820004A CN 113820004 A CN113820004 A CN 113820004A CN 202111085988 A CN202111085988 A CN 202111085988A CN 113820004 A CN113820004 A CN 113820004A
- Authority
- CN
- China
- Prior art keywords
- signal
- initial phase
- estimation method
- phase estimation
- ind
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Granted
Links
- 238000000034 method Methods 0.000 title claims abstract description 79
- 238000001228 spectrum Methods 0.000 claims abstract description 38
- 238000001914 filtration Methods 0.000 claims abstract description 27
- 238000004364 calculation method Methods 0.000 claims abstract description 15
- 230000009466 transformation Effects 0.000 claims abstract description 6
- 238000005070 sampling Methods 0.000 claims description 22
- 238000004458 analytical method Methods 0.000 claims description 13
- 238000006243 chemical reaction Methods 0.000 claims description 8
- 238000000605 extraction Methods 0.000 claims description 5
- 239000000126 substance Substances 0.000 claims description 3
- 238000003745 diagnosis Methods 0.000 abstract description 8
- 230000008569 process Effects 0.000 abstract description 7
- 238000012545 processing Methods 0.000 abstract description 7
- 230000000694 effects Effects 0.000 abstract description 5
- 238000012544 monitoring process Methods 0.000 abstract description 4
- 230000036541 health Effects 0.000 abstract description 3
- 230000010355 oscillation Effects 0.000 abstract description 2
- 238000010586 diagram Methods 0.000 description 11
- 238000005259 measurement Methods 0.000 description 11
- 230000008859 change Effects 0.000 description 7
- 238000010183 spectrum analysis Methods 0.000 description 6
- 238000005516 engineering process Methods 0.000 description 5
- 230000004044 response Effects 0.000 description 4
- 238000009825 accumulation Methods 0.000 description 2
- 238000001514 detection method Methods 0.000 description 2
- 238000012423 maintenance Methods 0.000 description 2
- 230000010363 phase shift Effects 0.000 description 2
- 238000011160 research Methods 0.000 description 2
- 230000001052 transient effect Effects 0.000 description 2
- 241000764238 Isis Species 0.000 description 1
- 230000005856 abnormality Effects 0.000 description 1
- 230000003044 adaptive effect Effects 0.000 description 1
- 230000003321 amplification Effects 0.000 description 1
- 238000005452 bending Methods 0.000 description 1
- 238000012937 correction Methods 0.000 description 1
- 238000013461 design Methods 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 238000002474 experimental method Methods 0.000 description 1
- 238000005562 fading Methods 0.000 description 1
- 238000003199 nucleic acid amplification method Methods 0.000 description 1
- 230000003595 spectral effect Effects 0.000 description 1
- 238000004611 spectroscopical analysis Methods 0.000 description 1
- 208000024891 symptom Diseases 0.000 description 1
- 238000012360 testing method Methods 0.000 description 1
Images
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 preceding groups
Landscapes
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Measurement Of Mechanical Vibrations Or Ultrasonic Waves (AREA)
Abstract
The invention provides a robust vibration signal initial phase estimation method, which comprises the steps of firstly carrying out hard threshold band-pass filtering according to an original signal, extracting a signal trend part, then extracting a signal instantaneous phase by using Hilbert transformation, estimating an initial phase by using linear interpolation and a signal intermediate position, and reducing the influence of a signal boundary 'Gibbs' effect on an estimation result. The whole process of the invention has only one characteristic frequency parameter, thus improving the robustness of the phase estimation algorithm; the related FFT, IFFT, Hilbert transform and the like have fast algorithms, meet online processing and can be effectively applied to real-time industrial monitoring situations such as fault diagnosis, equipment health management and the like; the signal-to-noise ratio at the specified characteristic frequency is improved through filtering, and the interference of other non-characteristic frequency spectrums on the initial phase estimation is reduced; the initial phase is estimated by filtering the position near half of the signal length, so that the calculation error caused by signal boundary oscillation generated by the Gibbs effect is avoided, and the applicable scene of the method is improved.
Description
Technical Field
The invention relates to the technical field of measurement and testing, in particular to a robust vibration signal initial phase estimation method.
Background
The traditional mechanical fault diagnosis technology mainly relies on a person with work experience on site to analyze through the human ear and a touch device instrument, the method needs accumulation of a large amount of work experience, the accumulation of expert experience cannot be well popularized, and secondly, subjective factors and objective factors of the person exist, and the fault diagnosis conclusion has high man-made interference. In equipment diagnosis and predictive maintenance, vibration signal analysis is one of the most widely used fault monitoring methods. The main technology is to analyze the characteristic frequency of the equipment through a time domain graph and a frequency domain graph to correspond to different fault types.
The time domain waveform reflects the characteristics of signal amplitude and phase change along with time, the time domain waveform can generally acquire signals of equipment through processing modes such as signal filtering, signal amplification and the like, the time domain signal can be visually represented in the running state and 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 by visually watching the vibration signal waveform. The vibration signal of a general mechanical device is a non-stationary signal, and the waveform of the vibration is generated by the rotation of a plurality of parts of the device
The time domain waveform is relatively single in equipment fault analysis, vibration signals of equipment components are superposed, and the current state of the equipment cannot be analyzed visually, so that the time domain signal needs to be subjected to Fourier change, dynamic complex signals are decomposed, the 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 current operation state of the equipment. The time-frequency conversion of the signal waveform can convert the acquired waveform in one period into a plurality of frequency components. According to Fourier transform, time domain signals are mapped to a frequency domain, the relationship between the amplitude of the signals and the characteristic frequency is reflected in a frequency domain diagram, whether the signals are in a linear relationship or not is not required to be considered in the analysis of the frequency domain signals, and the running state of equipment is directly reflected through the frequency signals.
The vibration signal during the starting and stopping process is usually called as 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 signs of the rotor. In the starting and stopping process of the instrument, the rotor experiences various rotating speeds, the vibration signal of the rotor is the response of a rotor system to the change of the rotating speed, the vibration signal is the external expression of the dynamic characteristics and the fault symptoms of the rotor, and the vibration signal contains rich information which is difficult to obtain at ordinary times. Therefore, start-stop process analysis is an important task for rotor detection. By analyzing the amplitude and phase shift of the vibration signal of the device, the resonant frequency and critical rotational speed of the rotor system can be found. From the amplitude and phase at low speed, knowledge of the degree of rotor bending is facilitated. The response of the rotor system to the rotor imbalance sloshing in the real startup and shutdown process can be seen through the phase change. (the system has a pronounced formant in amplitude response and a nearly 180 change in phase before and after it when passing through the critical speed.) the method of the invention
In order to avoid vibration abnormality caused by resonance and improve the accuracy of fault diagnosis, amplitude and phase estimation needs to be carried out on transient signals (starting and stopping conditions), the basic rigidity is increased by changing the structure of a machine, good balance of a rotor is guaranteed, and the running rotating speed of a fan is enabled to avoid a resonance region.
The phase calculation of the vibration signal is commonly performed as follows:
when a Fast Fourier Transform (FFT) spectrum analysis method is used to perform spectrum analysis on a discrete signal, the discrete signal is time-domain truncated, which may cause problems such as energy leakage, small spectrum amplitude, and low measurement accuracy. In order to reduce the errors, scholars at home and abroad propose a signal windowing interpolation FFT analysis algorithm based on a rectangular window, a Hanning window, a Blackman window and a Blackman-Harris window, so that the frequency spectrum leakage is inhibited to a certain extent, and the accuracy of signal parameter estimation is improved. By adjusting the parameters of the port, a window function with small peak level of the side lobe and large fading rate of the side lobe is selected. The signal processing can effectively reduce the long-range leakage caused by the mutual interference between adjacent frequency spectrum sidelobes when the signal time domain is truncated.
The sine approximation method is completely generated based on the development of modern computer technology and digital signal processing technology. The sine approximation method taking the least square method as the core can realize the absolute measurement of the phase-frequency characteristic of the vibration sensor in the broadband range of 1Hz to 10kHz, and is particularly suitable for the measurement of the amplitude-phase characteristic of the vibration sensor under the condition of low frequency and large amplitude with less sampling period number.
The full-phase FFT spectrum analysis method can restrain the spectrum leakage problem to a certain extent. The analysis method analyzes the frequency spectrums of two sequences with a time delay relationship, and corrects the results of two-time phase spectrum analysis by using a time-shift phase difference method, so that accurate phase and amplitude correction can be obtained according to the known sinusoidal signal frequency.
The analog method is to measure the phase of an input alternating signal by using a phase measurement circuit, and the circuit generally involves many high-performance components, and has the disadvantages of complex circuit design, expensive components, limited measurement 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 size degree of the method mainly depends on the interval of sampling time, and the calculation result is easily interfered by external factors. In the correlation method, when phase difference calculation is performed, a reference signal with the same frequency is often required as a reference, correlation calculation is performed on an 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 that the phase-frequency characteristics of two groups of input signals are calculated firstly, then the phase value of the current input signal is calculated and obtained according to the peak position in the spectrogram of the signal to be analyzed, when the phase difference of the two groups of signals needs to be carried out, the phase of a single signal can be calculated, and then the phase difference is subtracted. The latter two detection methods show good anti-interference performance in practical measurement experiments. Not only can spectroscopy be used to measure the difference between two sets of signals, but also to analyze parameters of each set of input signals independently, for example: amplitude, phase, frequency, etc.
In the prior art, under the condition of mixing different frequencies and noises, although phase information extracted by a full-phase time-shift phase difference method is very accurate, the accuracy of amplitude is greatly reduced; the sine approximation method is more accurate in amplitude information extraction, but the phase extraction is poor; the accuracy of FFT spectral analysis is relatively worst. When the frequency resolution of signal analysis is not enough to resolve the signal frequency, the amplitude and phase estimation accuracy is reduced. The traditional method for performing FFT spectrum analysis inevitably has energy leakage, small amplitude and low precision, and can not realize real whole period truncation. Although the accuracy of the measurement of the spectral amplitude is improved after various window functions are applied to the discrete signal, the maximum error of the phase is still large. The sine approximation method is mature in theory and wide in measurement range, but the sine approximation method adopts least square multiplication to fit the amplitude and the phase of a sine signal, and the error of an estimation quantity determined by the least square method depends on the error of measurement data and a functional relation given by an equation set.
Disclosure of Invention
The invention provides a robust vibration signal initial phase estimation method for online 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, firstly, hard threshold band-pass filtering is carried out according to an original signal, a main trend part of the signal is extracted, then, a Hilbert transformation is used for extracting an instantaneous phase of the signal, an initial phase is estimated (linear interpolation) by using a signal middle position, and the influence of a signal boundary 'Gibbs' effect 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 (fast Fourier transform) on an original signal f (t), and then obtaining a filtered signal by adopting hard threshold band-pass filtering
Will filter the signal
Hilbert transform extraction of instantaneous phase
Will instantaneously phase
Obtaining initial phase by linear interpolation and signal center positionValue of
The invention relates to a robust vibration signal initial phase estimation method, which comprises the following steps as a preferred mode:
s1, Fourier transform: carrying out FFT on an original signal F (t) to obtain a frequency spectrum F (omega);
s2: hard threshold bandpass filtering: carrying out hard threshold band-pass filtering on the frequency spectrum F (omega) to obtain a converted frequency spectrum
Then the transformed spectrum is processed
IFFT is carried out to obtain filtered signals
S3, Hilbert transformation: will filter the signal
Hilbert conversion is carried out to obtain a Hilbert conversion signal
Transforming the Hilbert signal
Extracting instantaneous phase after analysis
S4, linear interpolation: constructing an analytic signal z (t) and interpolating the instantaneous phase by linear interpolation
Analyzing, and estimating to obtain initial phase value according to the intermediate position of the signal
In a preferred embodiment of the robust vibration signal initial phase estimation method of the present invention, in step S4,
the analytic signal z (t) is:
wherein the polar coordinates of the analytic signal z (t) are:
a (t) is the envelope of the original signal f (t);
the expression of A (t) is:
in the method for estimating the initial phase of the robust vibration signal according to the present invention, as a preferred mode, in step S4, the instantaneous phase
Comprises the following steps:
in the method for estimating the initial phase of the robust vibration signal according to the present invention, as a preferred mode, in step S4, the initial phase value is obtained according to the signal intermediate position estimation
The specific method comprises the following steps:
calculating an index ind near a half of the signal length and judging whether the index ind is an integer, if so, determining ind belongs to [ t ∈ [ [ t ]1，t2，t3...，tL]Let' ind ═ tkK is L/2, initial phase value
Comprises the following steps:
if not, marking the lower adjacent integer and the upper adjacent integer of the index ind as inddAnd induInitial phase value
Comprises the following steps:
the invention relates to a robust vibration signal initial phase estimation method, which is a preferred mode, and an index ind calculation formula is as follows:
ind＝ntemp*nt，
wherein n istempThe number of whole periods is about half of the length of the whole signal, and nt is the number of sampling points in one period.
The invention discloses a robust vibration signal initial phase estimation method, which is an optimal mode, wherein the number n of the whole period is around half of the length of the whole signaltempThe calculation formula of (a) is as follows:
wherein L is the number of sampling points;
the calculation formula of the number of sampling points nt in one period is as follows:
wherein Fs is the sampling frequency f0Is the characteristic frequency.
In the method for estimating the initial phase of the robust vibration signal, which is 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 method for estimating the initial phase of the robust vibration signal, as a preferred mode, in step S2, the hard threshold band-pass filtering includes: the frequency spectrum F (omega) is set at [ 0.75F0，1.25*f0]Setting the outside of the region as 0 to obtain a transformed frequency spectrum
Wherein f is0Is a characteristic frequency, i.e., a frequency component to be preserved;
the IFFT is as follows:
wherein the content of the first and second substances,
for the filtered signal
Signal point of
In the method for estimating the initial phase of the robust vibration signal according to the present invention, as a preferred embodiment, in step S3, Hilbert transform is:
fourier Transform (FFT) and inverse transform (IFFT)
Mutual conversion between the time domain and the frequency domain of the vibration signal is realized by Fourier analysis, and the research of the original time domain signal can be converted into the research of Fourier coefficients on the frequency domain. In the field of signal processing, Fourier transform plays an important role, has milestone significance, and is regarded as a bridge between a signal time domain and a signal frequency domain. For signal x (t), its successive Fourier transforms (FFT) are
The inverse transform (IFFT) is:
in practical applications, signals tend to be discrete in time domain and frequency domain, so a Discrete Fourier Transform (DFT) is commonly used, and considering the operation speed and the system consumption, a Fast Fourier Transform (FFT) is widely applied to frequency domain analysis of signals.
Hard threshold bandpass filtering
In the invention, the FFT spectrum of the signal is calculated, the hard threshold processing is carried out on the high-frequency spectrum region, the IFFT is used for removing the noise of the signal, the signal-to-noise ratio of the signal is improved, and the accuracy of the subsequent estimation initial phase is improved.
Given a signal F (t), where t ∈ (— infinity, + ∞), its fourier transform is F (ω), assuming we want to preserve the frequency component as F (ω)0In part, let the frequency spectrum F (ω) be [ 0.75F [ ]0，1.25*f0]Out of range is 0 (hard thresholding), and the transformed spectrum is recorded as
For the
Computing a filtered signal by performing an inverse FFT (i.e., IFFT)
Hilbert transform
Let the real-valued signal f (t), where t ∈ (— infinity, + ∞), whose Hilbert transform (Hilbert) is defined as:
Note the book
Fourier transform of
Where ω is the frequency according to the theorem F (t) g (t)]＝F[f(t)]F[g(t)]，
Therefore, after the real-valued signal f (t) is Hilbert transformed, its amplitude is unchanged, but there is
Is generated for the positive frequency part of the original signal
Phase shift, for negative frequency of the original signalPartial generation
Phase shifting while the original signal f (t) is transformed with its Hilbert transform (Hilbert)
Are orthogonal.
Resolving a signal
Definition resolution signal
Wherein
Is a Hilbert transform of the real-valued signal f (t),
z (t) in polar coordinates:
a (t) is the envelope of the signal f (t),
instantaneous phase (angle) of signal f (t):
instantaneous frequency is defined as the derivative of phase:
linear interpolation
Linear interpolation is a simpler interpolation method, and the interpolation function is a first-order polynomial. Let the function y be f (x) at two points x0,x1Are each y0,y1Solving a polynomial equation
To satisfy
Is calculated to
Scheme flow
In order to effectively carry out initial phase estimation on a vibration signal and improve the subsequent fault diagnosis effect, the signal-to-noise ratio of the signal is improved through hard threshold band-pass filtering, an analytic signal is constructed through Hilbert transformation, and the instantaneous phase of the original signal is calculated through the analytic signal. In order to overcome the influence of spectrum leakage on the boundary part signals, the initial trigger position of the signals is estimated by using the position near the center of the signals, and the specific scheme flow is as follows: assuming that the length of the signal f (t) is 2048, and the sampling frequency Fs is 2.56 KHz: f (t) ═ f (t)1,t2,t3...,tL)＝[x1,x2,x3...,xL]。
1. Computing the original signal F (x) Fourier transform F (omega)
2. Setting a hard threshold filtering frequency range according to the rotating speed of the equipment: in this example, assuming that the rotation speed of the equipment is 3000 rpm, the corresponding 1-time multiplier (1X) of the equipment is 3000/60-50 Hz, 2-time multiplier (2X) is 3000/60-2-100 Hz, and … 8-time multiplier (8X) is 3000/60-8-400 Hz; if we need to estimate the 1 octave (1X: f)050HZ), the frequency spectrum F (ω) is set at [0.75 × F)0，1.25*f0]Out of range is 0 (hard thresholding), and the transformed spectrum is recorded as
To pair
Computing a filtered signal by performing an inverse FFT (i.e., IFFT)
3. Calculating a filtered signal
Hilbert transform of
4. Structure analysis signal
z (t) in polar coordinates:
a (t) is the envelope of the signal f (t),
is the instantaneous phase (angle unit) of the signal f (t)
5. Calculating the number of sampling points in one period
Taking the whole period number near half of the whole signal length
6. Calculate the index around half the signal length: ind ═ ntempNt, if ind is an integer
ind∈[t1，t2，t3...，tL]Let' ind ═ tkWherein k ∈ [1, 2, 3](ii) a And when K is L/2, the estimated initial phase is as follows:
if ind is not an integer, let ind be lower neighbor integer and upperAdjacent integers are ind respectivelydAnd induThen the estimated initial phase is,
the invention has the following advantages:
(1) the whole process of the method only has the characteristic frequency f0One parameter, the other intermediate parameters are adaptive. The robustness of the phase estimation algorithm is improved.
(2) In the calculation process of the method, fast algorithms such as FFT, IFFT, Hilbert transform and the like are provided, online processing is met, and the method 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 designated characteristic frequency through filtering, and reduces the interference of other non-characteristic frequency spectrums on the initial phase estimation
(4) The method estimates the initial phase by the position near half of the length of the filtering signal, and avoids the calculation error caused by signal boundary oscillation generated by the 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 robust vibration signal initial phase estimation method signal sample A1 FFT spectrum;
FIG. 4 is a diagram of the instantaneous phase of a signal sample A1 according to a robust vibration signal initial phase estimation method;
FIG. 5 is a diagram of a signal sample A2 single frequency signal of a robust vibration signal initial phase estimation method;
FIG. 6 is a signal sample A2FFT spectrogram of a robust vibration signal initial phase estimation method;
FIG. 7 is a diagram of the instantaneous phase of a signal sample A2 according to a robust vibration signal initial phase estimation method;
FIG. 8 is a diagram of a signal sample A3 single frequency signal for a robust vibration signal initial phase estimation method;
FIG. 9 is a diagram of a robust vibration signal initial phase estimation method signal sample A3 FFT spectrum;
fig. 10 is a diagram of the instantaneous phase of a signal sample a3 according to a robust vibration signal initial phase estimation method.
Detailed Description
The technical solutions in the embodiments of the present invention will be clearly and completely described below with reference to the drawings in the embodiments of the present invention, and it is obvious that the described embodiments are only a part of the embodiments of the present invention, and not all of the embodiments.
Example 1
A robust vibration signal initial phase estimation method is characterized in that an original signal f (t) is subjected to FFT and then filtered by adopting hard threshold band-pass filtering to obtain a filtered signal
Will filter the signal
Hilbert transform extraction of instantaneous phase
Will instantaneously phase
Obtaining initial phase value by linear interpolation and signal center position
The method comprises the following steps:
s1, Fourier transform: carrying out FFT on an 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: carrying out hard threshold band-pass filtering on the frequency spectrum F (omega) to obtain a converted frequency spectrum
Then the transformed spectrum is processed
IFFT is carried out to obtain filtered signals
The hard threshold bandpass filtering is: the frequency spectrum F (omega) is set at [ 0.75F0，1.25*f0]Setting the outside of the region as 0 to obtain a transformed frequency spectrum
Wherein f is0Is a characteristic frequency, i.e., a frequency component to be preserved;
the IFFT is as follows:
wherein the content of the first and second substances,
for the filtered signal
Signal point of
S3, Hilbert transformation: will filter the signal
Hilbert conversion is carried out to obtain a Hilbert conversion signal
Transforming the Hilbert signal
Extracting instantaneous phase after analysis
The Hilbert transform is:
s4, linear interpolation: constructing an analytic signal z (t) and interpolating the instantaneous phase by linear interpolation
Analyzing, and estimating to obtain initial phase value according to the intermediate position of the signal
The analytic signal z (t) is:
wherein the polar coordinates of the analytic signal z (t) are:
a (t) is the envelope of the original signal f (t);
the expression of A (t) is:
instantaneous phase
Comprises the following steps:
estimating and obtaining an initial phase value according to the intermediate position of the signal
The specific method comprises the following steps:
calculating an index ind near a half of the signal length and judging whether the index ind is an integer, if so, determining ind belongs to [ t ∈ [ [ t ]1，t2，t3...，tL]Let' ind ═ tkK is L/2, initial phase value
Comprises the following steps:
if not, marking the lower adjacent integer and the upper adjacent integer of the index ind as inddAnd induInitial phase value
Comprises the following steps:
the index ind is calculated as follows:
ind＝ntemp*nt，
wherein n istempThe number of the whole period is about half of the length of the whole signal, and nt is the number of sampling points in one period;
overall signal lengthNumber n of whole cycles around half degreetempThe calculation formula of (a) is as follows:
wherein L is the number of sampling points;
the calculation formula of the number of sampling points nt in one period is as follows:
wherein Fs is the sampling frequency f0Is the characteristic frequency.
Example 2
A robust vibration signal initial phase estimation method is characterized in that real-time data are collected by an intelligent operation and maintenance big data cloud platform of a space intelligent control (Beijing) monitoring technology Limited company, and single frequency signals A1 (figure 2: 25Hz, sampling frequency Fs is 2560Hz, the number of sampling points L is 2048, and the initial phase is 45 degrees) are respectively subjected to the method; a single frequency signal is added with random noise a2 (fig. 5: 25Hz, sampling frequency Fs is 2560Hz, sampling point number L is 2048, initial phase is 45 °, and random noise is added); a plurality of frequency signals are mixed and added with random noise a3 (fig. 8: f1 is 25Hz, initial phase 45 °, (f 2 is 50Hz, initial phase 60 °, (f 3 is 100Hz, initial phase 30 °, (sampling frequency Fs is 2560Hz, number of sampling points L is 2048, random noise is added)
1. The signal sample a1 is analyzed on-line as shown in fig. 2, and its FFT spectrum is calculated as shown in fig. 3. No noise is added to the signal, and the step of hard threshold band-pass filtering in the scheme is omitted. Estimation of instantaneous phase by Hilbert transform as shown in fig. 4, the initial phase was estimated to be 44.99 ° (error 0.01 °) by linear interpolation from a position around half the length of a 1.
2. The signal sample a2 is analyzed on line as shown in fig. 5, the FFT spectrum is calculated as shown in fig. 6, a2 is relatively interfered by noise, and f ═ 18.75,31.25] partial frequency is maintained by hard threshold band-pass filtering. Estimation of instantaneous phase by Hilbert transform as shown in fig. 7, the initial phase was estimated to be 44.67 ° (error 0.33 °) from a position around half the length of a 2.
3. The signal sample A3 is analyzed on line as shown in fig. 8, the FFT spectrum is calculated as shown in fig. 9, A3 is relatively large in noise interference, the initial phase of the frequency component with f2 equal to 50Hz is specified to be calculated, and the partial frequency with f equal to [43.75, 56.25] is maintained by hard threshold band-pass filtering. Estimation of instantaneous phase by Hilbert transform as shown in fig. 10, the initial phase was estimated to be 60.8 ° (error 0.8 °) from a position around half the length of a 3.
The above description is only for the preferred embodiment of the present invention, but the scope of the present invention is not limited thereto, and any person skilled in the art should be considered to be within the technical scope of the present invention, and the technical solutions and the inventive concepts thereof according to the present invention should be equivalent or changed within the scope of the present invention.
Claims (10)
1. A robust vibration signal initial phase estimation method is characterized in that: FFT is carried out on the original signal f (t), and then the filtered signal is obtained by adopting the band-pass filtering of the hard threshold value
Filtering the filtered signal
Hilbert transform extraction of instantaneous phase
The instantaneous phase is measured
Obtaining initial phase value by linear interpolation and signal center position
2. A robust vibration signal initial phase estimation method as defined in claim 1, wherein: the method comprises the following steps:
s1, Fourier transform: performing FFT on the original signal F (t) to obtain a frequency spectrum F (omega);
s2: hard threshold bandpass filtering: carrying out hard threshold band-pass filtering on the frequency spectrum F (omega) to obtain a converted frequency spectrum
And then converting the converted frequency spectrum
IFFT is carried out to obtain the filtered signal
S3, Hilbert transformation: filtering the filtered signal
Hilbert conversion is carried out to obtain a Hilbert conversion signal
Transforming the Hilbert transform signal
Extracting the instantaneous phase after analysis
S4, linear interpolation: constructing an analytic signal z (t) and interpolating the instantaneous phase by linear interpolation
Analyzing, and estimating according to the signal intermediate position to obtain the signalInitial phase value
3. A robust vibration signal initial phase estimation method as defined in claim 2, wherein:
in the step S4, in the step S,
the analytic signal z (t) is:
wherein the polar coordinates of the analytic signal z (t) are:
a (t) is the envelope of the original signal f (t);
the expression of A (t) is:
4. a robust vibration signal initial phase estimation method as defined in claim 2, wherein:
in step S4, the instantaneous phase
Comprises the following steps:
5. a robust vibration signal initial phase estimation method as defined in claim 2, wherein: in step S4, the signal intermediate position is estimated according to the signal intermediate positionThe initial phase value
The specific method comprises the following steps:
calculating an index ind near a half of the signal length and judging whether the index ind is an integer, if so, determining ind belongs to [ t ∈ [ [ t ]1，t2，t3...，tL]Let' ind ═ tkK is L/2, the initial phase value
Comprises the following steps:
if not, marking the lower adjacent integer and the upper adjacent integer of the index ind as inddAnd induThe initial phase value
Comprises the following steps:
6. a robust vibration signal initial phase estimation method as defined in claim 5, wherein: the index ind calculation formula is as follows:
ind＝ntemp*nt，
wherein n istempThe number of whole periods is about half of the length of the whole signal, and nt is the number of sampling points in one period.
7. A robust vibration signal initial phase estimation method as defined in claim 6, wherein:
number n of whole cycles around half of the total signal lengthtempThe calculation formula of (a) is as follows:
wherein L is the number of sampling points;
the calculation formula of the number of sampling points nt in one period is as follows:
wherein Fs is the sampling frequency f0Is the characteristic frequency.
8. A robust vibration signal initial phase estimation method as defined in claim 2, wherein:
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]。
9. a robust vibration signal initial phase estimation method as defined in claim 2, wherein:
in step S2, the hard threshold band-pass filtering is: (iii) locating said spectrum F (ω) at [ 0.75F [ ]0，1.25*f0]Setting the outside of the region range as 0 to obtain the transformed frequency spectrum
Wherein f is0Is a characteristic frequency, i.e., a frequency component to be preserved;
the IFFT is as follows:
wherein the content of the first and second substances,
for the filtered signal
Signal point of
10. A robust vibration signal initial phase estimation method as defined in claim 2, wherein:
in step S3, Hilbert transforms to:
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202111085988.2A CN113820004B (en) | 2021-09-16 | 2021-09-16 | Robust vibration signal initial phase estimation method |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202111085988.2A CN113820004B (en) | 2021-09-16 | 2021-09-16 | Robust vibration signal initial phase estimation method |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN113820004A true CN113820004A (en) | 2021-12-21 |
| CN113820004B 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 | Robust vibration signal initial phase estimation method |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN113820004B (en) |
Cited By (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN114357368A (en) * | 2021-12-31 | 2022-04-15 | 核工业西南物理研究院 | Nonlinear mode-mode coupling analysis method based on Hilbert transform |
Citations (9)
| 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 |
| US20180375519A1 (en) * | 2017-06-27 | 2018-12-27 | Intel IP Corporation | Phase synchronization between two phase locked loops |
| 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 |
-
2021
- 2021-09-16 CN CN202111085988.2A patent/CN113820004B/en active Active
Patent Citations (9)
| 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 |
| US20180375519A1 (en) * | 2017-06-27 | 2018-12-27 | Intel IP Corporation | Phase synchronization between two phase locked loops |
| 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相位的正弦波频率和初相的高精度估计方法", 《电子学报》, no. 09, 25 September 2001 (2001-09-25), pages 1164 - 1167 * |
Cited By (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN114357368A (en) * | 2021-12-31 | 2022-04-15 | 核工业西南物理研究院 | Nonlinear mode-mode coupling analysis method based on Hilbert transform |
Also Published As
| Publication number | Publication date |
|---|---|
| CN113820004B (en) | 2024-05-28 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| US8121813B2 (en) | System and method for clearance estimation between two objects | |
| CN108761117B (en) | Portable current detection rotating speed tester | |
| CN111307460A (en) | Rolling bearing fault diagnosis method based on calculation order tracking and spectral kurtosis | |
| Wang et al. | Bearing fault diagnosis under time-varying rotational speed via the fault characteristic order (FCO) index based demodulation and the stepwise resampling in the fault phase angle (FPA) domain | |
| 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 | |
| Rodopoulos et al. | A parametric approach for the estimation of the instantaneous speed of rotating machinery | |
| US10365297B2 (en) | System and method for generation of a tachometer signal and reduction of jitter | |
| CN111693775A (en) | Harmonic detection method, device and medium for power transmission network | |
| CN110224394B (en) | Fourier decomposition algorithm suitable for extracting characteristics of non-stationary power oscillation signal | |
| CN113361331B (en) | Power Frequency Interference Elimination Method, System and Medium Based on Windowed Interpolation FFT | |
| CN110108467A (en) | Active sounding speed-measuring method based on portable mobile apparatus | |
| CN111397877A (en) | Rotary machine beat vibration fault detection and diagnosis method | |
| CN112485028B (en) | Feature spectrum extraction method of vibration signal and mechanical fault diagnosis analysis method | |
| CN114061678B (en) | Digital driving method for Coriolis flowmeter | |
| CN113820004B (en) | Robust vibration signal initial phase estimation method | |
| Choudhury et al. | An overview of fault diagnosis of industrial machines operating under variable speeds | |
| CN111856401A (en) | Time delay estimation method based on cross-spectrum phase fitting | |
| Yang et al. | Resampling Technique based Demodulation Analysis for Planet Bearing Cage Fault Diagnosis under Nonstationary Conditions | |
| CN114486252B (en) | Rolling bearing fault diagnosis method of vector mode maximum envelope | |
| CN114487589A (en) | Power grid broadband signal self-adaptive measurement method, device and system | |
| CN115683644A (en) | Double-source beat vibration characteristic identification method for aircraft engine | |
| Desheng et al. | Time-frequency analysis based on BLDC motor fault detection using Hermite S-method | |
| CN114383718A (en) | High-frequency blade passing frequency extraction method based on vibration signals of external casing of gas turbine | |
| CN114184838A (en) | Power system harmonic detection method, system and medium based on SN mutual convolution window |
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 |