CN111259776B - Deterministic signal extraction method based on synchronous average principal component time-frequency analysis - Google Patents
Deterministic signal extraction method based on synchronous average principal component time-frequency analysis Download PDFInfo
- Publication number
- CN111259776B CN111259776B CN202010032827.6A CN202010032827A CN111259776B CN 111259776 B CN111259776 B CN 111259776B CN 202010032827 A CN202010032827 A CN 202010032827A CN 111259776 B CN111259776 B CN 111259776B
- Authority
- CN
- China
- Prior art keywords
- matrix
- tsa
- principal component
- time
- frequency
- 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
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2218/00—Aspects of pattern recognition specially adapted for signal processing
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F18/00—Pattern recognition
- G06F18/20—Analysing
- G06F18/21—Design or setup of recognition systems or techniques; Extraction of features in feature space; Blind source separation
- G06F18/213—Feature extraction, e.g. by transforming the feature space; Summarisation; Mappings, e.g. subspace methods
- G06F18/2135—Feature extraction, e.g. by transforming the feature space; Summarisation; Mappings, e.g. subspace methods based on approximation criteria, e.g. principal component analysis
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2218/00—Aspects of pattern recognition specially adapted for signal processing
- G06F2218/08—Feature extraction
Abstract
The invention discloses a deterministic signal extraction method based on time domain average principal component construction and time-frequency joint analysis, which belongs to the field of mechanical signal processing and comprises the following steps: (1) constructing a factor matrix; (2) constructing a TSA matrix; (3) matrix rearrangement based on energy indicators; (4) reducing the dimension of the matrix based on the difference spectrum and the energy ratio; (5) matrix decomposition based on principal component analysis; (6) reconstructing a matrix based on the principal component energy index; (7) frequency domain analysis of TSA main components of each order; and (8) performing time-frequency analysis on main components of all orders of the TSA. By using the method and the device, the deterministic signal can be effectively extracted under the condition of low signal-to-noise ratio.
Description
Technical Field
The invention belongs to the field of mechanical signal processing, and particularly relates to a deterministic signal extraction method based on time domain average principal component construction and time-frequency joint analysis.
Background
The mechanical equipment can generate sound and vibration time domain signals with rich information content in the operation process, and the current signal analysis methods of the mechanical equipment mainly comprise three types: time domain analysis, frequency domain analysis, and time-frequency analysis.
The time domain analysis is mainly divided into two types, the first type of time domain analysis method is statistical parameter analysis, a target time domain signal is selected, and various corresponding statistical parameters such as an average value, a peak-to-peak value, a standard deviation, a variance, a root mean square and the like are calculated for analysis. The second type of time domain analysis method is overall observation and trend prediction, which focuses on the size of the signal peak, the amplitude of different periods, the time corresponding to the signal peak, the cycle length of the repeated appearance of the waveform with the same shape, and the like, and focuses on the signal change condition, and predicts the signal change trend on the basis.
The main means of frequency domain analysis is fast fourier transform, and the basic idea of fast fourier transform is to decompose a target time domain signal into trigonometric functions with different periods, so that the frequency distribution situation of a stationary signal can be perfectly described, but the method is not suitable for analysis of a non-stationary signal.
The main methods of time-frequency analysis mainly include short-time Fourier transform and wavelet transform. The basic idea of short-time Fourier transform is to decompose a target time domain signal into a plurality of small segments in a signal windowing manner, approximate each small segment of signal to a stable signal, perform Fourier transform on each segment of signal to obtain a spectrum analysis result, and combine all spectrum analysis results to obtain time-frequency joint distribution of the target signal. Since the window function is fixed, the short-time fourier transform has a drawback that the time resolution and the frequency resolution are mutually limited. Different from Fourier transform, which takes an infinite-length trigonometric function as a basis function to decompose a target signal, wavelet transform takes a wavelet function as the basis function and realizes signal analysis of different frequency bands and different time points through the operations of stretching and translating the basis function, thereby ensuring that a wavelet time-frequency analysis result has higher frequency resolution at a low frequency and higher time resolution at a high frequency and realizing the self-adaptive representation of time-frequency analysis.
In the above methods, the time domain method can only obtain time domain information, and cannot obtain frequency domain information, which has great limitation; while the frequency domain analysis can extract components with determined periodicity in the signals, but is not suitable for non-stationary signal analysis, and noise and vibration signals in the actual operation of mechanical equipment are often non-stationary signals; the time-frequency analysis can obtain the representation of two dimensions of time and frequency, but the time resolution and the frequency resolution are often limited, and under the condition that the signal-to-noise ratio of the deterministic signal is low, the time-frequency analysis method is difficult to effectively extract the deterministic signal, and the running state and fault information of the machine are often reflected as specific deterministic signals.
Disclosure of Invention
The invention provides a deterministic signal extraction method based on synchronous average principal component time-frequency analysis, which can effectively extract deterministic signal components in time domain signals under the condition of low signal-to-noise ratio and can be widely applied to the fields of real-time monitoring, state identification, fault diagnosis and the like of various mechanical equipment.
A deterministic signal extraction method based on synchronous average principal component time-frequency analysis comprises the following steps:
(1) Selecting discrete time domain signals to be analyzed, intercepting the total number of signal length points to obtain target time domain signals, and constructing a factor matrix of the total number of signal length points;
(2) Sequentially taking the factors in the factor matrix as the segment number of the time domain synchronous average, carrying out time domain synchronous average processing on the target time domain signal to obtain a time domain synchronous average array corresponding to each factor, and combining all the time domain synchronous average arrays to form an original TSA matrix;
(3) Calculating the root mean square value of each row of the TSA matrix as an energy index, and rearranging the TSA matrix according to the sequence of the energy indexes from large to small;
(4) Calculating a difference spectrum of the obtained energy index, and selecting a difference spectrum order based on a difference spectrum peak value; setting an energy sum threshold, and selecting an energy order based on the energy sum threshold; selecting a larger value from the difference spectrum order and the energy order as a final dimension reduction order, and performing dimension reduction processing on the rearranged TSA matrix according to the value;
(5) Calculating a covariance matrix corresponding to the dimension-reduced TSA matrix, and performing characteristic decomposition on the covariance matrix to obtain a corresponding characteristic value diagonal matrix and a characteristic vector matrix;
(6) Setting a threshold value of sum of eigenvalues, selecting the number of principal components based on the threshold value of sum of eigenvalues, forming a corresponding eigenvector matrix, and multiplying the reduced dimension TSA matrix by the corresponding eigenvector matrix to obtain a TSA matrix of the principal components;
(7) Performing frequency domain analysis on each principal component in the principal component TSA matrix to obtain a frequency domain two-dimensional curve corresponding to each principal component, and performing comprehensive drawing on each curve to obtain a frequency domain analysis result of the principal component TSA matrix;
(8) Performing time-frequency analysis on each principal component in the principal component TSA matrix to obtain a time-frequency matrix corresponding to each principal component, and drawing a time-frequency graph; and adding and averaging the corresponding position elements of all the time frequency matrixes to obtain a uniform time frequency matrix, and drawing a uniform time frequency graph.
The method can overcome the defect that the existing time domain analysis, frequency domain analysis and time frequency analysis can not comprehensively, accurately and effectively extract the deterministic signal, is simple and easy to operate, and can effectively realize the enhancement and extraction of the deterministic signal under the condition of low signal-to-noise ratio so as to realize the corresponding frequency domain characterization and time frequency characterization.
In the step (1), taking the total number L of discrete points contained in the target time domain signal as the target number, and solving all factors F of the L 1 ,F 2 ,…F n And arranged in descending order to form a factor matrix F = { F = 1 ,F 2 ,…F n }。
In the step (2), the specific steps of constructing the original TSA matrix are as follows:
(2-1) for any factor F in the factor matrix i And constructing a corresponding time domain synchronous average array, which specifically comprises the following steps:
(2-1-1) arbitrary factor F in the factor matrix i Dividing the target time domain signal into F as a number of segments i A segment;
(2-1-2) adding F i Adding points at corresponding positions in the segments to obtain a sum signal;
(2-1-3) dividing the sum signal by the number of segments F i And obtaining a time domain synchronous average array, wherein the calculation formula is as follows:
in the formula, TSA (t) m ,F i ) Representing the mth value in the time domain synchronous average array;
(2-2) transversely copying the time domain synchronous average array corresponding to each factor to extend to the length L of the target time domain signal, wherein the calculation formula is as follows:
in the formula, TSA (F) i ) Representing the time domain synchronous average array after the horizontal copying expansion;
(2-3) vertically stacking the time domain synchronous average arrays after the horizontal copying and expansion to form an original TSA matrix, wherein the calculation formula is as follows:
in the formula, TSA (t, F) represents the original TSA matrix.
The specific steps of the step (3) are as follows:
(3-1) calculating the root mean square of each row in the original TSA matrix as an energy index, wherein the calculation formula is as follows:
in the formula, RMS TSA Represents the corresponding root mean square matrix, rms, of the original TSA matrix i Represents the root mean square value of the ith row of the original TSA matrix;
(3-2) rearranging the RMS matrix according to the sequence from large to small to obtain a rearranged RMS matrix, wherein the calculation formula is as follows:
in the formula, RMS r Representing the rearranged root mean square matrix, rmsr i Representing the root mean square value after rearrangement, satisfying rmsr 1 ≥rmsr 2 ≥…≥rmsr n ；
(3-3) rearranging the original TSA matrix according to the sequence of the root mean square values from large to small to obtain a rearranged TSA matrix, wherein the calculation formula is as follows:
wherein F (rmsr) i ) Represents rmsr i A corresponding factor.
In the step (4), the specific step of selecting the order of the difference spectrum is as follows:
and performing difference processing on adjacent elements in the rearranged root-mean-square value matrix to obtain a difference spectrum, wherein the calculation formula is as follows: d i ＝rmsr i -rmsr i+1 (ii) a Finding the maximum value in the difference spectrum, and taking the corresponding subscript as the order of the difference spectrum, wherein the calculation formula is as follows:
in the formula I d Representing the order of the difference spectrum;
the specific steps for selecting the energy order are as follows:
selecting an energy sum threshold s; energy order I is calculated according to the following formula e ：
In the formula, I represents the selected RMS order, and n represents the total RMS order.
In the step (4), the step of performing dimension reduction on the rearranged TSA matrix comprises the following specific steps:
(4-1) taking the larger value of the difference spectrum order and the energy order as a dimensionality reduction order, wherein the calculation formula is as follows:
I＝max(I d ,I e )
in the formula, I is the dimensionality reduction order;
and (4-2) selecting the first I row of the rearranged TSA matrix to form the dimension-reduced TSA matrix based on the dimension-reduced order I.
The specific steps of the step (5) are as follows:
(5-1) calculating a covariance matrix corresponding to the dimension-reduced TSA matrix, wherein the calculation formula is as follows:
TSA cov ＝cov(TSA e (t,F))
in the formula, TSA cov Representing a covariance matrix, TSA e (t, F) is a dimension-reduced TSA matrix, cov (·) represents covariance operation;
(6-2) performing characteristic decomposition on the covariance matrix to obtain a corresponding characteristic value diagonal matrix and a corresponding characteristic vector matrix, wherein the calculation formula is as follows:
[V,U]＝eig(TSA cov )
U＝[u 1 u 2 … u I ]∈R I×I
in the formula, V represents a characteristic value diagonal matrix, U represents a characteristic vector matrix, and eig (·) represents a characteristic decomposition operation.
The specific steps of the step (6) are as follows:
(6-1) selecting a feature value sum threshold a;
(6-2) calculating the number k of the principal components according to the following formula;
in the formula, k p Indicates the number of selected principal components, lambda i Representing the ith eigenvalue in an eigenvalue diagonal matrix;
(6-3) selecting the first k eigenvectors based on the number k of the principal components to form a reconstructed eigenvector matrix, wherein a calculation formula is as follows;
U s ＝[u 1 u 2 … u k ]
in the formula of U s A reconstructed feature vector matrix is obtained;
(6-4) transposing the dimension-reduced TSA matrix, and multiplying the transposed TSA matrix by the reconstructed eigenvector matrix to obtain a TSA principal component matrix consisting of all principal components, wherein the calculation formula is as follows:
TSAF＝TSA e (t,F) T ·U S ＝[TSAP 1 TSAP 2 … TSAP k ]＝[TSA(F(rmsr 1 )) T TSA(F(rmsr 2 )) T … TSA(F(rmsr I )) T ]·[u 1 u 2 … u k ]
wherein TSAF is TSA principal component matrix, TSAP i Is the ith TSA principal component signal.
The specific steps of the step (7) are as follows:
(7-1) carrying out Fourier transform on each principal component signal in the TSA principal component matrix to obtain a corresponding amplitude spectrum, wherein the calculation formula is as follows:
in the formula, P i (f) The amplitude spectrum corresponding to the ith TSA principal component signal;
and (7-2) in an xyz three-dimensional rectangular coordinate system, taking each principal component signal number as an x coordinate, taking the corresponding principal component signal time domain sequence as a y coordinate, taking the corresponding principal component signal amplitude spectrum sequence as a z coordinate, and drawing a frequency domain analysis result of each principal component signal in the TSA principal component matrix.
The specific steps of the step (8) are as follows:
(8-1) performing short-time Fourier transform on each principal component signal in the TSA principal component matrix to obtain a corresponding time-frequency matrix, and drawing a corresponding time-frequency graph, wherein the calculation formula is as follows:
in the formula, STFT i (t, f) is the time-frequency spectrum corresponding to the ith TSA principal component signal, g f,t (τ) is a window function of the short-time Fourier transform, p i (t j ,f k ) Is the element of the jth row and the kth column in the time-frequency matrix;
(8-2) adding and averaging all the position elements corresponding to the time frequency matrix to obtain a uniform time frequency matrix, drawing a uniform time frequency graph, wherein the calculation formula is as follows:
in the formula, STFT a And (t, f) is a normalized time spectrum corresponding to the TSA principal component matrix.
Compared with the prior art, the invention has the following beneficial effects:
(1) The invention provides a deterministic signal extraction method based on time domain average principal component construction and time-frequency joint analysis, which can effectively extract deterministic signal components in signals.
(2) The method is based on time domain synchronous average processing and principal component analysis technology, can greatly improve the signal-to-noise ratio of the deterministic signal component under the condition of low signal-to-noise ratio, and has excellent anti-noise performance.
(3) The method is based on frequency domain analysis and time-frequency analysis methods, can start from a frequency domain and a time-frequency domain, comprehensively represents the frequency distribution state of a target signal, emphasizes and highlights deterministic signal components, and related information can be used for real-time monitoring, state identification and fault diagnosis of mechanical equipment.
Drawings
FIG. 1 is a schematic flow chart of a deterministic signal extraction method based on synchronous mean principal component time-frequency analysis according to the present invention;
fig. 2 shows the time-frequency analysis result of a single-frequency signal (SNR =0 dB) according to an embodiment of the present invention;
fig. 3 is a frequency domain analysis result of a single frequency signal (SNR =0 dB) according to an embodiment of the present invention;
fig. 4 shows the time-frequency analysis result of a multi-frequency signal (SNR =0 dB) according to an embodiment of the present invention;
fig. 5 shows the frequency domain analysis result of the multi-frequency signal (SNR =0 dB) according to the embodiment of the invention;
fig. 6 shows the time-frequency analysis result of a single wideband signal (SNR =0 dB) according to an embodiment of the present invention;
fig. 7 shows the frequency domain analysis result of a single wideband signal (SNR =0 dB) according to an embodiment of the present invention;
fig. 8 shows the time-frequency analysis result of multiple wideband signals (SNR =0 dB) according to an embodiment of the present invention;
fig. 9 shows the frequency domain analysis result of multi-wideband signal (SNR =0 dB) according to the embodiment of the present invention.
Detailed Description
The invention will be described in further detail below with reference to the drawings and examples, which are intended to facilitate the understanding of the invention and are not intended to limit it in any way.
As shown in fig. 1, a deterministic signal extraction method based on synchronous average principal component time-frequency analysis includes the following steps:
s01, factor matrix construction
Selecting discrete time domain signals to be analyzed, intercepting the total number L of signal length points to obtain target time domain signals, and constructing a factor matrix of the total number L of the signal length points.
In the step, the total number L of discrete points contained in the target time domain signal is used as the target number, and all factors F of the L are solved 1 ,F 2 ,…F n And arranged in descending order to form a factor matrix F = { F = { (F) 1 ,F 2 ,…F n }。
S02, TSA matrix construction
In turn by a factor F in a factor matrix i And i ∈ {1,2, … n } is used as the number of segments of the time domain synchronization average to perform time domain synchronization average processing on the target time domain signal to obtain a time domain synchronization average array corresponding to each factor, each time domain synchronization average array is copied and expanded to the total number L of points of the target signal, and each expanded time domain synchronization average array is used as a row to be filled into the original TSA matrix from small to large according to the factors to form the original TSA matrix with n rows and L columns.
The time domain synchronous averaging processing comprises the following specific steps:
(2-1) arbitrary factor F in factor matrix i Dividing the target time domain signal into F as a number of segments i A segment;
(2-2) adding F i Adding points at corresponding positions in the segments to obtain a sum signal;
(2-3) dividing the sum signal by the number of segments F i And obtaining a time domain synchronous average array, wherein the calculation formula is as follows:
in the formula, TSA (t) m ,F i ) Representing the mth value in the time domain synchronous average array.
The specific steps for constructing the original TSA matrix are as follows:
(2-1)' for any factor F in the factor matrix i Constructing a corresponding time domain synchronous average array;
(2-2)' transversely copying and extending the time domain synchronous average array corresponding to each factor to the length L of the target time domain signal, wherein the calculation formula is as follows:
in the formula, TSA (F) i ) Represents the expanded time domain synchronous average array after horizontal replication.
(2-3)' vertically stacking the time domain synchronous average arrays after the horizontal replication and expansion to form an original TSA matrix, wherein the calculation formula is as follows:
in the formula, TSA (t, F) represents the original TSA matrix.
S03, matrix rearrangement based on energy index
And calculating the root mean square value of each row of the TSA matrix as an energy index, and rearranging the TSA matrix according to the sequence of the energy indexes from large to small. The method comprises the following specific steps:
(3-1) calculating the root mean square of each row in the original TSA matrix as an energy index, wherein the calculation formula is as follows:
in the formula, RMS TSA Represents the corresponding root mean square matrix, rms, of the original TSA matrix i Represents the root mean square value of the ith row of the original TSA matrix.
(3-2) rearranging the RMS matrix according to the sequence from large to small to obtain a rearranged RMS matrix, wherein the calculation formula is as follows:
in the formula, RMS r Representing the rearranged root mean square matrix, rmsr i Representing the rearranged root mean square value and satisfying the rmsr 1 ≥rmsr 2 ≥…≥rmsr n 。
(3-3) rearranging the original TSA matrix according to the sequence of the corresponding root mean square values from large to small to obtain a rearranged TSA matrix, wherein the calculation formula is as follows:
wherein F (rmsr) i ) Stands for rmsr i A corresponding factor.
S04, reducing the dimension of the matrix based on the difference spectrum and the energy ratio
Calculating to obtain a difference spectrum of the energy index, and selecting a difference spectrum order based on a difference spectrum peak value; setting an energy sum threshold, and selecting an energy order based on the energy sum threshold; and selecting a larger value from the difference spectrum order and the energy order as a final dimension reduction order, and carrying out dimension reduction on the rearranged TSA matrix according to the value.
The specific steps for determining the order of the difference spectrum are as follows:
(4-1) carrying out difference processing on adjacent elements in the rearranged root mean square value matrix to obtain a difference spectrum, wherein a calculation formula is as follows:
d i ＝rmsr i -rmsr i+1
(4-2) finding the maximum value in the difference spectrum, taking the corresponding subscript as the order of the difference spectrum, and calculating according to the following formula:
in the formula I d Representing the order of the difference spectrum.
The specific steps for determining the energy order are as follows:
(4-1)' selecting an energy sum threshold s;
(4-2)' the energy order I is calculated according to the following formula e ：
The specific steps of performing matrix dimensionality reduction on the rearranged TSA matrix are as follows:
taking the larger value of the difference spectrum order and the energy order as a dimensionality reduction order, and the calculation formula is as follows:
I＝max(I d ,I e )
in the formula, I is the dimensionality reduction order.
And selecting the first I row of the rearranged TSA matrix to form the dimension-reduced TSA matrix based on the dimension-reducing order I.
S05, matrix decomposition based on principal component analysis
And calculating a covariance matrix corresponding to the dimension-reduced TSA matrix, and performing characteristic decomposition on the covariance matrix to obtain a corresponding characteristic value diagonal matrix and a characteristic vector matrix. The method comprises the following specific steps:
(5-1) calculating a covariance matrix corresponding to the dimension-reduced TSA matrix, wherein the calculation formula is as follows:
TSA cov ＝cov(TSA e (t,F))
in the formula, TSA cov Representing a covariance matrix, TSA e (t, F) is the dimension-reduced TSA matrix, cov (·) stands for covariance operation.
(5-2) performing characteristic decomposition on the covariance matrix to obtain a corresponding characteristic value diagonal matrix and a corresponding characteristic vector matrix, wherein the calculation formula is as follows:
[V,U]＝eig(TSA cov )
U＝[u 1 u 2 … u I ]∈R I×I
in the formula, V represents a characteristic value diagonal matrix, U represents a characteristic vector matrix, and eig (·) represents a characteristic decomposition operation.
S06, matrix reconstruction based on principal component energy index
Setting a threshold value of sum of eigenvalues, selecting the number of principal components based on the threshold value of sum of eigenvalues, forming a corresponding eigenvector matrix, and multiplying the dimension-reduced TSA matrix by the corresponding eigenvector matrix to obtain the principal component TSA matrix. The method comprises the following specific steps:
in the step (6), the matrix reconstruction based on the principal component energy index comprises the following specific steps:
(6-1) selecting a feature value sum threshold a;
(6-2) calculating the number k of the principal components according to the following formula;
(6-3) selecting the first k eigenvectors based on the number k of the principal components to form a reconstructed eigenvector matrix, wherein a calculation formula is as follows;
U s ＝[u 1 u 2 … u k ]
in the formula of U s To reconstruct the eigenvector matrix.
(6-4) transposing the dimension-reduced TSA matrix, and multiplying the transposed TSA matrix by the reconstructed eigenvector matrix to obtain a TSA principal component matrix consisting of all principal components, wherein the calculation formula is as follows:
TSAF＝TSA e (t,F) T ·U S ＝[TSAP 1 TSAP 2 … TSAP k ]＝[TSA(F(rmsr 1 )) T TSA(F(rmsr 2 )) T … TSA(F(rmsr I )) T ]·[u 1 u 2 … u k ]
wherein TSAF is TSA principal component matrix, TSAP i Is the ith TSA principal component signal.
S07, frequency domain analysis of TSA principal components of each order
And performing frequency domain analysis on each principal component in the principal component TSA matrix to obtain a frequency domain two-dimensional curve corresponding to each principal component, and performing comprehensive drawing on each curve to obtain a frequency domain analysis result of the principal component TSA matrix. The method comprises the following specific steps:
(7-1) performing Fourier transform on each principal component signal in the TSA principal component matrix to obtain a corresponding amplitude spectrum, wherein the calculation formula is as follows:
in the formula, P i (f) Amplitude spectrum corresponding to ith TSA principal component signal
And (7-2) in an xyz three-dimensional rectangular coordinate system, drawing the frequency domain analysis result of each principal component signal in the TSA principal component matrix by taking each principal component signal number as an x coordinate, the corresponding principal component signal time domain sequence as a y coordinate and the corresponding principal component signal amplitude spectrum sequence as a z coordinate.
S08, time-frequency analysis of TSA main components of each order
Performing time-frequency analysis on each principal component in the principal component TSA matrix to obtain a time-frequency matrix corresponding to each principal component, and drawing a time-frequency graph; and adding and averaging the corresponding position elements of all the time frequency matrixes to obtain a uniform time frequency matrix, and drawing a uniform time frequency graph. The method comprises the following specific steps:
(8-1) performing short-time Fourier transform on each principal component signal in the TSA principal component matrix to obtain a corresponding time-frequency matrix, and drawing a corresponding time-frequency graph, wherein the calculation formula is as follows:
in the formula, STFT i (t, f) is the time-frequency spectrum corresponding to the ith TSA principal component signal, g f,t (τ) is a window function of the short-time Fourier transform, p i (t j ,f k ) Is the j row and k column element in the time frequency matrix.
(8-2) adding and averaging the position elements corresponding to all the time-frequency matrixes to obtain a uniform time-frequency matrix, and drawing a uniform time-frequency graph, wherein the calculation formula is as follows:
in the formula, STFT a And (t, f) is a normalized time spectrum corresponding to the TSA principal component matrix.
In order to verify the effectiveness of the invention, the time-frequency and frequency-domain analysis is respectively carried out on a single-frequency signal, a multi-frequency signal, a single-broadband signal and a multi-broadband signal, which is as follows:
the above method is used to analyze the single-frequency signal x (t) =10cos 2000 pi t + n (t), (n (t) is gaussian white noise, and the single-frequency signal SNR =0 dB), so as to obtain the corresponding time-frequency analysis result and frequency-domain analysis result, which are shown in fig. 2 and fig. 3, respectively. It can be seen that a significant spectral line exists at the position of 1000Hz in the time-frequency diagram, which represents a deterministic single-frequency signal with 1000Hz frequency component; in the frequency domain diagram, the frequency spectrum corresponding to the first principal component shows a line spectrum at the frequency of 1000Hz, which also represents a deterministic single-frequency signal with the frequency component of 1000 Hz.
The multi-frequency signal x (t) =20cos 2400 pi t +15cos 4200 pi t +18cos 7000 pi t + n (t), (n (t) is white gaussian noise, and the multi-frequency signal to noise ratio SNR =0 dB) is analyzed by using the above method, so as to obtain a corresponding time-frequency analysis result and a corresponding frequency-domain analysis result, which are respectively shown in fig. 4 and fig. 5. It can be seen that significant spectral lines exist at the positions of 1200Hz, 2100Hz and 3500Hz of the frequency in the time-frequency diagram respectively, and represent deterministic signal components with frequency components of 1200Hz, 2100Hz and 3500 Hz; in the frequency domain diagram, line spectrums appear at the positions of 1200Hz, 2100Hz and 3500Hz in the frequency spectrums corresponding to the first two main components, and also represent the deterministic signal components of 1200Hz, 2100Hz and 3500Hz in frequency components.
The single wideband signal x (t) = c (t) + n (t) (c (t) is a wideband signal, the frequency range is 1000Hz-2000hz, n (t) is gaussian white noise, and the single wideband signal-to-noise ratio SNR =0 dB) is analyzed by using the above method, and the corresponding time-frequency analysis result and frequency-domain analysis result are obtained, as shown in fig. 6 and fig. 7, respectively. It can be seen that the frequency range of 1000Hz-2000Hz in the time-frequency diagram has a significant dark broadband, representing a broadband signal component with a frequency component of 1000Hz-2000 Hz; in the frequency domain diagram, the corresponding frequency spectrum of the first three main components has prominent peaks in the range of 1000Hz-2000Hz, and the frequency components represent the broadband signal components with the frequency components of 1000Hz-2000 Hz.
The multi-wideband signal x (t) = c (t) + n (t) (c (t) is a wideband signal, the frequency ranges are 1000Hz-2000Hz and 3000Hz-4000hz, n (t) is gaussian white noise, and the multi-wideband signal SNR =0 dB) is analyzed by using the above method, and the corresponding time-frequency analysis result and frequency-domain analysis result are obtained, as shown in fig. 8 and fig. 9, respectively. It can be seen that the frequency ranges of 1000Hz-2000Hz and 3000Hz-4000Hz in the time-frequency diagram have significant dark color broad bands, which represent broad band signal components with frequency components of 1000Hz-2000Hz and 3000Hz-4000 Hz; in the frequency domain diagram, the corresponding frequency spectrums of the first three main components have prominent peaks in the frequency ranges of 1000Hz-2000Hz and 3000Hz-4000Hz, and also represent broadband signal components with the frequency components of 1000Hz-2000Hz and 3000Hz-4000 Hz.
The embodiments described above are intended to illustrate the technical solutions and advantages of the present invention, and it should be understood that the above-mentioned embodiments are only specific embodiments of the present invention, and are not intended to limit the present invention, and any modifications, additions and equivalents made within the scope of the principles of the present invention should be included in the scope of the present invention.
Claims (9)
1. A deterministic signal extraction method based on synchronous average principal component time-frequency analysis is characterized by comprising the following steps:
(1) Selecting discrete time domain signals to be analyzed, intercepting the total number of signal length points to obtain target time domain signals, and constructing a factor matrix of the total number of signal length points;
(2) Sequentially taking factors in the factor matrix as segment numbers of time domain synchronous average, performing time domain synchronous average processing on a target time domain signal to obtain a time domain synchronous average array corresponding to each factor, and combining all the time domain synchronous average arrays to form an original TSA matrix; the specific steps for constructing the original TSA matrix are as follows:
(2-1) for any factor F in the factor matrix i And constructing a corresponding time domain synchronous average array, which specifically comprises the following steps:
(2-1-1) arbitrary factor F in the factor matrix i Dividing the target time domain signal into F as a number of segments i A segment;
(2-1-2) adding F i Adding points at corresponding positions in the segments to obtain a sum signal;
(2-1-3) dividing the sum signal by the number of segments F i And obtaining a time domain synchronous average array, wherein the calculation formula is as follows:
in the formula, TSA (t) m ,F i ) Representing the mth value in the time domain synchronous average array;
(2-2) transversely copying the time domain synchronous average array corresponding to each factor to extend to the length L of the target time domain signal, wherein the calculation formula is as follows:
wherein TSA (F) i ) Representing the time domain synchronous average array after the horizontal copying expansion;
(2-3) vertically stacking the time domain synchronous average arrays after the horizontal copying and expansion to form an original TSA matrix, wherein the calculation formula is as follows:
in the formula, TSA (t, F) represents the original TSA matrix;
(3) Calculating the root mean square value of each row of the TSA matrix as an energy index, and rearranging the TSA matrix according to the sequence of the energy indexes from large to small;
(4) Calculating a difference spectrum of the obtained energy index, and selecting a difference spectrum order based on a difference spectrum peak value; setting an energy sum threshold, and selecting an energy order based on the energy sum threshold; selecting a larger value from the difference spectrum order and the energy order as a final dimension reduction order, and carrying out dimension reduction on the rearranged TSA matrix according to the value;
(5) Calculating a covariance matrix corresponding to the dimension-reduced TSA matrix, and performing characteristic decomposition on the covariance matrix to obtain a corresponding characteristic value diagonal matrix and a characteristic vector matrix;
(6) Setting a threshold value of sum of eigenvalues, selecting the number of principal components based on the threshold value of sum of eigenvalues, forming a corresponding eigenvector matrix, and multiplying the reduced dimension TSA matrix by the corresponding eigenvector matrix to obtain a TSA matrix of the principal components;
(7) Performing frequency domain analysis on each principal component in the principal component TSA matrix to obtain frequency domain two-dimensional curves corresponding to each principal component, and performing comprehensive drawing on each curve to obtain a frequency domain analysis result of the principal component TSA matrix;
(8) Performing time-frequency analysis on each principal component in the principal component TSA matrix to obtain a time-frequency matrix corresponding to each principal component, and drawing a time-frequency graph; and adding and averaging the corresponding position elements of all the time frequency matrixes to obtain a uniform time frequency matrix, and drawing a uniform time frequency graph.
2. The method for extracting deterministic signals based on synchronous mean principal component time-frequency analysis according to claim 1, wherein in step (1), the total number L of discrete points contained in the target time-domain signal is used as the target number, and all factors F of L are solved 1 ,F 2 ,…F n And arranged in descending order to form a factor matrix F = { F = 1 ,F 2 ,…F n }。
3. The method for extracting deterministic signals based on synchronous mean principal component time-frequency analysis according to claim 2, wherein the step (3) comprises the following steps:
(3-1) calculating the root mean square of each row in the original TSA matrix as an energy index, wherein the calculation formula is as follows:
where RMS TSA Represents the corresponding root mean square matrix, rms, of the original TSA matrix i Represents the root mean square value of the ith row of the original TSA matrix;
(3-2) rearranging the RMS matrix according to the sequence from large to small to obtain a rearranged RMS matrix, wherein the calculation formula is as follows:
in the formula, RMS r Representing the rearranged root mean square matrix, rmsr i Representing the rearranged root mean square value and satisfying the rmsr 1 ≥rmsr 2 ≥…≥rmsr a ；
(3-3) rearranging the original TSA matrix according to the sequence of the root mean square values from large to small to obtain a rearranged TSA matrix, wherein the calculation formula is as follows:
wherein F (rmsr) i ) Represents rmsr i A corresponding factor.
4. The method for extracting deterministic signals based on synchronous mean principal component time-frequency analysis according to claim 2, wherein the step (4) of selecting the order of the difference spectrum comprises the following steps:
and performing difference processing on adjacent elements in the rearranged root mean square value matrix to obtain a difference spectrum, wherein the calculation formula is as follows: d is a radical of i ＝rmsr i -rmsr i+1 (ii) a Finding the maximum value in the difference spectrum, taking the corresponding subscript as the order of the difference spectrum, and calculating according to the following formula:
in the formula I d Representing the order of the difference spectrum;
the specific steps for selecting the energy order are as follows:
selecting an energy sum threshold s; energy order I is calculated according to the following formula e ：
In the formula, I represents the selected RMS order, and n represents the RMS total order.
5. The method for extracting deterministic signals based on synchronous mean principal component time-frequency analysis according to claim 4, wherein in the step (4), the step of performing dimension reduction on the rearranged TSA matrix comprises the following specific steps:
(4-1) taking the larger value of the difference spectrum order and the energy order as a dimensionality reduction order, wherein the calculation formula is as follows:
I＝max(I d ，I e )
in the formula, I is the dimensionality reduction order;
and (4-2) selecting the first I row of the rearranged TSA matrix to form the dimension-reduced TSA matrix based on the dimension-reduced order I.
6. The method for extracting deterministic signals based on synchronous mean principal component time-frequency analysis according to claim 1, wherein the step (5) comprises the following steps:
(5-1) calculating a covariance matrix corresponding to the dimension-reduced TSA matrix, wherein the calculation formula is as follows:
TSA cov ＝cov(TSA e (t，F))
in the formula, TSA cov Representing a covariance matrix, TSA e (t, F) is a dimension-reduced TSA matrix, cov (·) represents covariance operation;
(6-2) performing characteristic decomposition on the covariance matrix to obtain a corresponding characteristic value diagonal matrix and a corresponding characteristic vector matrix, wherein the calculation formula is as follows:
『V，U]＝eig(TSA cov )
u＝[u 1 u 2 … u I ]∈R I×I
in the formula, V represents a characteristic value diagonal matrix, U represents a characteristic vector matrix, and eig (·) represents a characteristic decomposition operation.
7. The method for extracting deterministic signals based on synchronous mean principal component time-frequency analysis according to claim 1, wherein the step (6) comprises the following steps:
(6-1) selecting a feature value sum threshold a;
(6-2) calculating the number k of the principal components according to the following formula;
in the formula, k p Indicates the number of selected principal components, lambda i Representing the ith eigenvalue in the eigenvalue diagonal matrix;
(6-3) selecting the first k eigenvectors based on the number k of the principal components to form a reconstructed eigenvector matrix, wherein a calculation formula is as follows;
U S ＝[u 1 u 2 … u k ]
in the formula of U S A reconstructed feature vector matrix is obtained;
(6-4) transposing the dimension-reduced TSA matrix, and multiplying the transposed TSA matrix by the reconstructed eigenvector matrix to obtain a TSA principal component matrix consisting of all principal components, wherein the calculation formula is as follows:
TSAF＝TSA e (t，F) T ·U S ＝[TSAP 1 TSAP 2 … TSAP k ]
＝[TSA(F(rmsr 1 )) T TSA(F(rmsr 2 )) T … TSA(F(rmsr I )) T ]·[u 1 u 2 … u k ]
wherein TSAF is TSA principal component matrix, TSAP i Is the ith TSA principal component signal.
8. The method for extracting deterministic signals based on synchronous mean principal component time-frequency analysis according to claim 1, wherein the step (7) comprises the following steps:
(7-1) carrying out Fourier transform on each principal component signal in the TSA principal component matrix to obtain a corresponding amplitude spectrum, wherein the calculation formula is as follows:
in the formula, P i (f) The amplitude spectrum corresponding to the ith TSA principal component signal;
and (7-2) in an xyz three-dimensional rectangular coordinate system, taking each principal component signal number as an x coordinate, taking the corresponding principal component signal time domain sequence as a y coordinate, taking the corresponding principal component signal amplitude spectrum sequence as a z coordinate, and drawing a frequency domain analysis result of each principal component signal in the TSA principal component matrix.
9. The method for extracting deterministic signals based on synchronized mean principal component time-frequency analysis according to claim 1, wherein the step (8) comprises the following steps:
(8-1) performing short-time Fourier transform on each principal component signal in the TSA principal component matrix to obtain a corresponding time-frequency matrix, and drawing a corresponding time-frequency graph, wherein the calculation formula is as follows:
in the formula, STFT i (t, f) is the time-frequency spectrum corresponding to the ith TSA principal component signal, g f，t (τ) is a window function of the short-time Fourier transform, p i (t j ，f k ) Is the element of the jth row and the kth column in the time frequency matrix;
(8-2) adding and averaging the position elements corresponding to all the time-frequency matrixes to obtain a uniform time-frequency matrix, and drawing a uniform time-frequency graph, wherein the calculation formula is as follows:
in the formula, STFT a And (t, f) is a normalized time spectrum corresponding to the TSA principal component matrix.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202010032827.6A CN111259776B (en) | 2020-01-13 | 2020-01-13 | Deterministic signal extraction method based on synchronous average principal component time-frequency analysis |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202010032827.6A CN111259776B (en) | 2020-01-13 | 2020-01-13 | Deterministic signal extraction method based on synchronous average principal component time-frequency analysis |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN111259776A CN111259776A (en) | 2020-06-09 |
| CN111259776B true CN111259776B (en) | 2023-04-18 |
Family
ID=70954142
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202010032827.6A Active CN111259776B (en) | 2020-01-13 | 2020-01-13 | Deterministic signal extraction method based on synchronous average principal component time-frequency analysis |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN111259776B (en) |
Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| EP2151822A1 (en) * | 2008-08-05 | 2010-02-10 | Fraunhofer-Gesellschaft zur Förderung der angewandten Forschung e.V. | Apparatus and method for processing and audio signal for speech enhancement using a feature extraction |
| CN102662167A (en) * | 2012-04-11 | 2012-09-12 | 西北工业大学 | Feature extraction method of radiated noise signal of underwater target |
| WO2013107076A1 (en) * | 2012-01-19 | 2013-07-25 | 东南大学 | Adaptive window fourier phase extraction method in optical three-dimensional measurement |
| CN105956624A (en) * | 2016-05-06 | 2016-09-21 | 东南大学 | Motor imagery electroencephalogram classification method based on space-time-frequency optimization feature sparse representation |
| CN107145843A (en) * | 2017-04-20 | 2017-09-08 | 浙江大学 | The rotating machinery frequency domain character method for extracting signal counted based on sequential |
-
2020
- 2020-01-13 CN CN202010032827.6A patent/CN111259776B/en active Active
Patent Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| EP2151822A1 (en) * | 2008-08-05 | 2010-02-10 | Fraunhofer-Gesellschaft zur Förderung der angewandten Forschung e.V. | Apparatus and method for processing and audio signal for speech enhancement using a feature extraction |
| WO2013107076A1 (en) * | 2012-01-19 | 2013-07-25 | 东南大学 | Adaptive window fourier phase extraction method in optical three-dimensional measurement |
| CN102662167A (en) * | 2012-04-11 | 2012-09-12 | 西北工业大学 | Feature extraction method of radiated noise signal of underwater target |
| CN105956624A (en) * | 2016-05-06 | 2016-09-21 | 东南大学 | Motor imagery electroencephalogram classification method based on space-time-frequency optimization feature sparse representation |
| CN107145843A (en) * | 2017-04-20 | 2017-09-08 | 浙江大学 | The rotating machinery frequency domain character method for extracting signal counted based on sequential |
Non-Patent Citations (1)
| Title |
|---|
| 车红昆 ; 吕福在 ; 项占琴 ; .基于顺序向前浮动搜索时频优选特征的缺陷识别.浙江大学学报(工学版).2235-2239. * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN111259776A (en) | 2020-06-09 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN107192553B (en) | Gear-box combined failure diagnostic method based on blind source separating | |
| US20080262371A1 (en) | Method for Adaptive Complex Wavelet Based Filtering of Eeg Signals | |
| CN102519725B (en) | Method for processing vibration signal of bearing equipment through nonlinear redundancy lifting wavelet packet | |
| Yi et al. | Mechanical compound faults extraction based on improved frequency domain blind deconvolution algorithm | |
| CN110151175A (en) | Surface electromyogram signal noise-eliminating method based on CEEMD and improvement wavelet threshold | |
| CN110146812A (en) | A kind of Method of Motor Fault Diagnosis based on the study of characteristic node increment type width | |
| CN110096956A (en) | Signal antinoise method and device based on EEMD and arrangement entropy second differnce | |
| CN110109058A (en) | A kind of planar array deconvolution identification of sound source method | |
| CN113537102B (en) | Feature extraction method of microseismic signals | |
| Liang et al. | The optimal ratio time-frequency mask for speech separation in terms of the signal-to-noise ratio | |
| CN109884592A (en) | A kind of auditory localization emulation mode towards low frequency Gaussian noise source | |
| CN116153329A (en) | CWT-LBP-based sound signal time-frequency texture feature extraction method | |
| CN111259776B (en) | Deterministic signal extraction method based on synchronous average principal component time-frequency analysis | |
| CN103915102B (en) | Method for noise abatement of LFM underwater sound multi-path signals | |
| Dang et al. | Fault diagnosis of power transformer by acoustic signals with deep learning | |
| AT510359A1 (en) | METHOD FOR ACOUSTIC SIGNAL TRACKING | |
| CN112180315A (en) | Fault feature extraction method, device and system for optical fiber current transformer | |
| Li et al. | Interference classification and identification of TDCS based on improved convolutional neural network | |
| Tirtom et al. | Enhancement of time-frequency properties of ECG for detecting micropotentials by wavelet transform based method | |
| Liu et al. | A method for blind source separation of multichannel electromagnetic radiation in the field | |
| Joyce | A separable 2-D autoregressive spectral estimation algorithm | |
| CN103198053A (en) | Instantaneous wavelet bicoherence method based on phase randomization | |
| CN106054132B (en) | A kind of ISM methods based on the selection of effective subband and detection statistic weighting | |
| CN112800862B (en) | Non-stationary signal time-frequency matrix reconstruction method and system | |
| CN112995084B (en) | Signal processing method and processing device |
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 |