CN104820786B - A kind of instantaneous weighting is synchronous to squeeze small echo double-spectrum analysis method - Google Patents

Abstract

An instantaneous weighted synchronous squeezing wavelet bispectral analysis method, comprising six steps: Step 1 to determine the minimum number of segments of the signal, Step 2 to calculate the weight coefficients corresponding to each segment, and Step 3 to calculate each signal segment The weight coefficients corresponding to the segments form the weight matrix of the signal, and the frequency sequence corresponding to each column of weight coefficients is calculated. Step 4 calculates the synchronously squeezed wavelet coefficients of the signal to obtain a time-frequency domain representation with a more compact distribution. Step 5 pairs The wavelet coefficients obtained by weighting the synchrosqueezed wavelet transform coefficients in the frequency range are corrected wavelet coefficients. Step 6 calculates the weighted synchrosqueezed wavelet bispectrum and the instantaneous weighted synchrosqueezed wavelet bispectrum; Wavelet transform can accurately distinguish the frequency components of the signal, and according to the information of the overall frequency distribution of the signal, add weights to the calculated synchronous squeeze wavelet coefficients in time segments, avoiding the shortcomings of low accuracy of traditional methods.

Description

An Instantaneous Weighted Synchronous Squeezing Wavelet Bispectral Analysis Method

technical field

The invention belongs to the technical field of biomedical signal processing, in particular to an instantaneous weighted synchronous squeezing wavelet bispectral analysis method for detecting the instantaneous phase coupling relationship between non-stationary signals.

Background technique

Bispectrum analysis is an effective technique to quantitatively study the non-normal characteristics of nonlinear system response, and is widely used in econometrics, underwater acoustic signal processing, mechanical fault diagnosis and biomedical signal feature extraction and other fields. In the field of biomedical signal processing, since almost all physiological systems exhibit significant nonlinear behavior, the physiological signals generated by them have significant non-stationary and non-normal characteristics. This type of signal is mostly generated by a phase-coupled harmonic process, so bispectral analysis is often used in the analysis and processing of this type of signal. The traditional bispectrum analysis method is based on the time-averaged fast Fourier transform, which has defects in the ability to identify the instantaneous characteristics of the signal. In order to solve this problem, Janez et al. proposed wavelet bispectrum, an analysis technique with the ability to identify instantaneous phase coupling relationships. Combined with the wavelet transform of the time-frequency analysis capability, the concept of instantaneous wavelet bispectrum is proposed, which makes it possible to study the instantaneous nonlinear phase coupling relationship between non-stationary signals in real time, so that more accurate results can be obtained in practical applications. test results.

Due to the band-pass filtering characteristics of wavelet transform, it will cause spectrum leakage in different frequency bands, and the frequency bands may overlap in adjacent scales, so that the wavelet coefficients of the signal will be distorted, and then interference will be introduced when bispectral calculation is performed. Therefore, the accuracy of the algorithm is still flawed. When applied to the phase coupling detection of frequency-time-varying signals, this defect directly causes wavelet bispectrum to wrongly reflect the coupling relationship between current signals. In the actual non-stationary system, the phase coupling relationship is more complicated. In the process of nonlinear phase coupling, not only internal intermodulation products are generated, but also non-phase coupled self-harmonics are generated. Even if these harmonics do not undergo nonlinear phase coupling with other components, as long as their frequencies satisfy the bispectral condition, And the phase condition is not strict, then there will still be peaks at the corresponding frequency coordinates of the calculated bispectrum, and then give wrong results when studying the nonlinear coupling between signals, so a more effective technology is urgently needed to improve detection accuracy.

Contents of the invention

In order to overcome the defects of the above-mentioned prior art, the purpose of the present invention is to propose an instantaneous weighted synchrosqueezing wavelet bispectrum analysis method, which introduces the synchrosqueezing wavelet transform with higher frequency resolution when calculating the wavelet bispectrum, and the signal The information of the overall frequency distribution avoids the disadvantage of low accuracy of traditional methods.

In order to achieve the above object, the technical scheme adopted in the present invention is:

An instantaneous weighted simultaneous squeeze wavelet bispectral analysis method, comprising the following steps:

Step 1: Given a discrete signal sequence g(n) of length N, use the formula (1) to calculate the minimum number of segments nseg required by the signal:

Where: f _s ——signal sampling rate,

f ₀ ——the required minimum frequency resolution;

Step 2: After determining the number of signal segments, use formula (2) to calculate the weight coefficient corresponding to each signal segment:

Among them: m—the serial number of the weight coefficient,

x—the serial number of the signal segment,

G _x (m)——the absolute value of the Fourier transform amplitude of the xth segment of the signal g(n), calculated by formula (3):

G _max - the maximum value of G _x (m);

Step 3: Use the calculated weight coefficients corresponding to each signal segment as a column vector, and arrange them in sequence according to the time period to form a weight matrix of the signal g(n), w=[w ₁ (m),...,w _nseg (m)], each column in the weight matrix corresponds to the elements in the frequency sequence, and the corresponding frequency sequence is obtained by formula (4):

Step 4: Use the formula (5) to calculate the wavelet transform of the discrete sequence g(n), and obtain the time-frequency domain expression form of the discrete signal sequence:

Among them: ψ——the selected mother wavelet function,

u——time factor,

a _j ——the discretization scale parameter of the mother wavelet function ψ,

n——the discretization translation parameter of the mother wavelet function ψ,

g(u)——the signal sequence to be analyzed,

*——indicates to take the conjugate;

Then use the formula (6) to calculate the synchronous squeeze wavelet transform coefficient of the signal g(n):

Among them: Gn—constant that determines the number of discretization scales a _j ,

a _j —— discretization scale, which can be determined by formula (7):

f(a _j ,n)——the frequency plane obtained by deriving the time-frequency plane of wavelet transform,

f _i ——the frequency corresponding to the scale a _i , satisfying the relation f _i =1/a _i ,

f _i ⁺ ,f _i ^- ——the upper and lower bounds of the frequency range determined according to f _i can be determined by formula (8):

n——time factor,

C _ψ —constant coefficient, which can be calculated by formula (9):

Among them: *—indicates to take the conjugate;

Step 5: Use the formula (10) to weight the coefficients of the synchronous squeezed wavelet transform according to the frequency interval to obtain the corrected wavelet coefficients:

WSW _g (f _l ,n ₀ )=SW _g (f _l ,n ₀ )*w _x (k)(10)

Among them: SW _g (f _l ,n ₀ )——synchronous squeeze wavelet transform coefficient of signal g(n),

w _x (k)——the weight coefficient corresponding to the segment where the time factor n ₀ is located,

f _l ——the frequency factor of the expression form in time-frequency domain obtained by synchrosqueezing wavelet transform,

n ₀ ——the time factor of the time-frequency domain expression obtained by synchrosqueezing wavelet transform,

x—the sequence number of the segment where the time factor n ₀ is located, which can be determined by formula (11):

k—the frequency sequence number of the frequency interval where the frequency factor f _l is located, which can be determined by formula (12):

f(k-1)<f _l ≤f(k) (12)

Step 6: Substituting the result calculated by formula (10) into formula (13) to obtain the weighted simultaneous squeezing wavelet bispectrum:

Where: the frequencies f ₁ , f ₂ , and f ₃ must satisfy the relation f ₃ =f ₁ +f ₂ ,

WSW _g (f ₁ ,n)——the value of the weighted synchronous squeeze wavelet coefficient at frequency f ₁ and time n,

WSW _g (f ₂ ,n)——the value of the weighted synchrosqueezing wavelet coefficient at frequency f ₂ and time n,

WSW _g (f ₃ ,n)——weighted synchrosqueezing wavelet coefficient at frequency f ₃ and time n;

Substituting the result calculated by formula (10) into formula (14), the instantaneous weighted synchronous squeezing wavelet bispectrum of the signal sequence g(n) at time n ₀ is obtained:

Since the calculated instantaneous weighted synchrosqueezed wavelet bispectrum is a complex number, it can be expressed in the form of formula (15):

Among them: A(f ₁ ,f,n ₀ )——Instantaneous weighted synchronously squeezed wavelet bispectral amplitude at double frequency (f ₁ ,f ₂ ) at time n ₀ ,

φ(f ₁ , f ₂ , n ₀ )——At time n ₀ , the instantaneous weighted synchronously squeezed wavelet bispectral phase at dual frequencies (f ₁ , f ₂ ).

The advantages of the present invention are: in order to solve the problem of overlapping wavelet decomposition coefficients existing in wavelet bispectrum, a synchronous squeezing wavelet transform with higher frequency resolution is introduced so as to accurately distinguish the frequency components of signals. Secondly, in order to suppress the interference caused by a large number of harmonic components produced in the phase coupling process of non-stationary signals, weights are added to the calculated synchrosqueezed wavelet coefficients by time segments. Therefore, the defect that the traditional method cannot identify the phase coupling relationship between non-stationary signals due to low accuracy is avoided. The invention can not only be used in the technical field of biomedical signal processing, but also can be used in the detection of instantaneous phase coupling relationship in other nonlinear systems.

Description of drawings

Figure 1 is the time-domain waveform of the simulated signal.

Figure 2 shows the weight coefficients of each segmented signal.

Fig. 3 is the wavelet transform time-frequency diagram of the simulated signal.

Fig. 4 is the time-frequency diagram of the synchrosqueezing wavelet transform of the simulated signal.

Fig. 5 is a three-dimensional image of instantaneous bispectrum based on traditional wavelet transform.

Fig. 6 is a three-dimensional image of instantaneous bispectrum based on synchrosqueezing wavelet transform.

Fig. 7 is the time-varying curve of the peak phase and amplitude of the instantaneous wavelet bispectrum of the simulated signal.

Fig. 8 is the time-varying curve of the peak phase and amplitude of the simulated signal instantaneously weighted synchronously squeezed wavelet bispectrum.

Fig. 9 is a flowchart of the present invention.

specific embodiment

The present invention will be described in detail below in conjunction with accompanying drawings and examples.

Taking the simulation signal as an example, the expression of the simulation signal is shown in formula (16)

X(n)＝cos(2πf ₁ n)+η(cos(2πf ₁ n)-cos(2πf ₂ n)) ² +ξ(n) (16)

Among them: ξ(n)——Gaussian white noise with a mean value of zero;

f ₁ ——time-varying coupling frequency, its value is shown in formula (17):

f ₂ ——time-varying coupling frequency, its value is shown in formula (18):

η—coupling strength coefficient;

The duration of the simulation signal X(n) is 300 seconds, the sampling rate is 10Hz, and the time-varying coupling frequencies f ₁ and f ₂ are respectively modulated by sinusoidal signals of different frequencies, which increase from 0.85 to 0.95 and 0.25 to 0.35. The value of the coupling strength coefficient η changes alternately, η = 0 (no coupling exists), η = 0.6 (weak coupling exists), and changes alternately every 20 seconds. The time-domain waveform and its frequency spectrum of the simulated signal are shown in Figure 1. In order to identify the secondary phase coupling alternately appearing between simulated frequency and time-varying signals, the invention is used to analyze the data.

An instantaneous weighted simultaneous squeeze wavelet bispectral analysis method, comprising the following steps:

Step 1: Given a simulation signal sequence X(n), take the minimum frequency resolution f ₀ =0.1, f _s =10, N=3000, and use formula (1) to calculate the minimum number of segments nseg=30 required for the signal.

Step 2: After determining the number of signal segments, formula (2) is used to calculate the weight coefficients corresponding to each signal segment, and the calculated weight matrix is shown in FIG. 2 .

Step 3: Use the calculated weight coefficients corresponding to each signal segment as a column vector, and arrange them according to the order of the time period to form a weight matrix of the signal X(n), w=[w ₁ (m),...,w ₃₀ (m)], and then calculate the frequency sequence corresponding to each column of weight coefficients.

Step 4: Use formula (5) to calculate the wavelet transform of the simulated discrete sequence X(n), and obtain the time-frequency domain expression form of the discrete signal sequence. The time-frequency diagram of the signal is shown in Figure 3, as can be seen from Figure 3, Due to the band-pass characteristic of wavelet transform, the wavelet coefficients are scattered in adjacent frequency intervals, which will cause the overlap of wavelet coefficients between frequency intervals. At the same time, it can be seen that when the secondary phase coupling occurs (η=0.6), Due to the time-varying frequency of the signal, a large number of high-order harmonics are generated during the coupling process, which are reflected in the time-frequency diagram as wavelet coefficients with a wide range of frequency distribution.

Then use formula (6) to calculate the synchronous squeeze wavelet transform of signal X(n), and the obtained time-frequency diagram is shown in Figure 4. It can be seen from Figure 4 that the synchrosqueezing wavelet transform obtains a very compact time-frequency surface, and the distribution of wavelet coefficients is compressed, which avoids the overlapping problem of wavelet coefficients between different frequency intervals; During the time period when the coupling occurs, high-order harmonics distributed on the time-frequency plane are observed. These harmonics will mislead the subsequent bispectral calculation and need to be suppressed by adding weights to the wavelet coefficients distributed in different frequency intervals.

Step 5: Use the formula (10) to weight the coefficients of the synchro-squeezing wavelet transform according to the frequency range to obtain the corrected wavelet coefficients.

Step 6: Substituting the result calculated by formula (10) into formula (13), and calculating the weighted synchrosqueezed wavelet bispectrum representing the overall coupling characteristics of the signal. Substituting the result calculated by formula (10) into formula (14), the instantaneous weighted synchronous squeeze wavelet bispectrum of the simulated signal at time n ₀ is calculated.

Figure 5 is a three-dimensional diagram of the traditional instantaneous wavelet bispectrum of the simulated signal when secondary phase coupling occurs (η=0.6). It can be found that due to the overlapping of wavelet coefficients in adjacent frequency intervals, in the bispectrum calculation results, the peak distribution Also overlap each other, making it difficult to analyze dual frequencies where phase coupling occurs in real time. Figure 6 is a three-dimensional image of the instantaneous weighted synchronously squeezed wavelet bispectrum at the same time as Figure 5. Its peaks are independently distributed, and the three-dimensional image is very compact, which can analyze the dual frequencies where phase coupling occurs in real time.

Applying the calculated instantaneous weighted synchrosqueezing wavelet bispectrum, the coupling characteristics of simulated signals changing with time can be accurately detected. Figure 7 shows the curves of the amplitude and phase of the maximum peak (the main frequency component that produces the secondary phase coupling) of the simulated signal's traditional instantaneous wavelet bispectrum as a function of time, and Figure 8 shows the instantaneous weighted synchronously squeezed wavelet bispectrum A plot of the amplitude and phase of the maximum peak value as a function of time. Theoretically, when there is a secondary phase coupling, the peak value of the bispectrum is not zero, and the phase is a constant. By comparing the instantaneous weighted synchronous squeezing wavelet bispectral analysis proposed by the present invention and the traditional instantaneous wavelet bispectral analysis results show that the instantaneous bispectral analysis based on the traditional wavelet transform cannot distinguish the hidden frequency-time variation in the simulation signal X(n) Intermittent phase coupling between signals, and the instantaneous weighted synchronous squeezing wavelet bispectrum analysis algorithm proposed by the present invention is in the non-coupling time period: 0-20, 40-60, 80-100, 120-140, 160-180 , 200-220, 240-260, 280-300 (in seconds), the phase shows a random change trend, and the amplitude is almost zero. In the time period with coupling, the phase presents a linear change trend, which is constant, and the corresponding amplitude is also greater than zero. It can accurately distinguish the time period of secondary phase coupling in the simulation signal X(n) with time-varying frequency.

The traditional calculation formula of wavelet bispectrum is:

Where: frequency f ₁ , f ₂ , f satisfies the relation f=f ₁ +f ₂ ,

W _g (f,n)——Wavelet transform.

It has two major defects: First, when the wavelet transform it uses is used to process frequency-varying signals, due to the band-pass filtering characteristics of wavelet transform, wavelet coefficients overlap in adjacent scale intervals, so three times of cumulative calculation Sometimes errors will be introduced, resulting in wrong conclusions when studying the phase coupling relationship of non-stationary signals. Secondly, according to the Fourier bispectrum theory, if a harmonic process as shown in formula (20):

Where: ω—angular frequency;

φ——phase;

Satisfying the condition ω ₃ ＝ω ₁ +ω ₂ ; φ ₃ ＝φ ₁ +φ ₂ +φ ₀ , where φ ₀ is the phase delay, then x(n) can be called a second-order phase-coupled harmonic process, and its bispectrum Impulse values will be generated at the frequency coordinates (ω ₁ ,ω ₂ ), (ω ₂ ,ω ₁ ). In the actual non-stationary signal, the phase coupling relationship is more complicated. In the nonlinear phase coupling process, not only the internal Intermodulation products will also produce a large number of their own harmonics. Even if these harmonics do not have nonlinear phase coupling with other components, as long as the above conditions are met between their frequencies, then the calculated frequency coordinates corresponding to the bispectrum are still There will be peaks, which will disturb the essential results of bispectrum in the study of nonlinear coupling between signals, making it easy to make wrong judgments when wavelet bispectrum is analyzing non-stationary signals.

Claims (1)

1. an instantaneous weighted synchronous squeeze wavelet bispectrum analysis method, is characterized in that, comprises the following steps: Step 1: Given a discrete signal sequence g(i) of length N, use the formula (1) to calculate the minimum number of segments nseg required by the signal: Where: f _s ——signal sampling rate, f ₀ ——the required minimum frequency resolution; Step 2: After determining the number of signal segments, use formula (2) to calculate the weight coefficient corresponding to each signal segment: Among them: m—the serial number of the weight coefficient, x——Signal segment number, G _x (m)——the absolute value of the Fourier transform amplitude of the xth segment of the signal sequence g(i), calculated by formula (3): Among them, g _x (q)——the xth segment data of the signal sequence g(i), G _max - the maximum value of G _x (m); Step 3: Use the calculated weight coefficients corresponding to each signal segment as a column vector, and arrange them in sequence according to the time period to form a weight matrix of the signal sequence g(i), w=[w ₁ (m),...,w _nseg (m)], each column in the weight matrix has a one-to-one correspondence with the elements in the frequency sequence, and the corresponding frequency sequence is obtained by formula (4): Step 4: Use the formula (5) to calculate the wavelet transform of the discrete signal sequence g(i), and obtain the time-frequency domain expression of the discrete signal sequence: Where: ψ—the selected mother wavelet function, a _j ——the discretization scale parameter of the mother wavelet function ψ, n——the discretization translation parameter of the mother wavelet function ψ, *——indicates to take the conjugate; Then use the formula (6) to calculate the synchrosqueezing wavelet transform coefficient of the signal sequence g(i): Among them: Gn—constant that determines the number of discretization scales a _j , a _j is determined by formula (7): f(a _j ,n)——the frequency plane obtained by deriving the time-frequency plane of wavelet transform, f _j ——the frequency corresponding to the scale a _j , satisfying the relation f _j =1/a _j , f _j ⁺ , f _j ^- ——the upper and lower bounds of the frequency interval determined according to f _j , determined by the formula (8): C _ψ —constant coefficient, calculated by formula (9): Where: ε—integral variable; Step 5: Use the formula (10) to weight the coefficients of the synchronous squeezed wavelet transform according to the frequency interval to obtain the corrected wavelet coefficients: Among them: SW _g (f _l ,n ₀ )——synchronous squeeze wavelet transform coefficient of signal sequence g(i), ——the weight coefficient corresponding to the segment where the time factor n ₀ is located, f _l ——the frequency factor of the expression form in time-frequency domain obtained by synchrosqueezing wavelet transform, n ₀ ——the time factor of the time-frequency domain expression obtained by synchrosqueezing wavelet transform, x ₀ ——the signal segment sequence number at the moment corresponding to the time factor n ₀ , determined by formula (11): k——the frequency serial number of the frequency interval where the frequency factor f _l is located, satisfying the formula (12): f(k-1)<f _l ≤f(k) (12) Step 6: Substituting the result calculated by formula (10) into formula (13) to obtain the weighted simultaneous squeezing wavelet bispectrum: Where: the frequencies f ₁ , f ₂ , and f ₃ must satisfy the relation f ₃ =f ₁ +f ₂ , WSW _g (f ₁ ,n ₀ )——the value of the weighted synchronous squeeze wavelet coefficient at frequency f ₁ and time factor n ₀ , WSW _g (f ₂ ,n ₀ )——the value of the weighted synchronous squeeze wavelet coefficient at frequency f ₂ and time factor n ₀ , WSW _g (f ₃ ,n ₀ )——the value of the weighted synchrosqueezing wavelet coefficient at frequency f ₃ and time factor n ₀ ; Substituting the result calculated by formula (10) into formula (14), the instantaneous weighted synchronous squeeze wavelet bispectrum of discrete signal sequence g(i) at time factor n ₀ is obtained: Since the calculated instantaneous weighted synchrosqueezed wavelet bispectrum is a complex number, it can be expressed in the form of formula (15): Among them: A(f ₁ ,f,n ₀ )——at the time factor n ₀ , the instantaneous weighted synchronously squeezed wavelet bispectral amplitude at dual frequencies (f ₁ , f ₂ ), φ(f ₁ , f ₂ , n ₀ )——At time factor n ₀ , the instantaneous weighted synchronously squeezed wavelet bispectral phase at dual frequencies (f ₁ , f ₂ ).