阅读说明：本技术 一种时频分析方法和系统 (Time-frequency analysis method and system ) 是由 轩建平 李锐 唐律 张青 于 2020-11-25 设计创作，主要内容包括：本发明公开了一种时频分析方法和系统。该方法包括步骤：对输入信号进行第一次信号分解；对第一次信号分解得到的信号进行第二次信号分解,获得三个长度相等的信号；对第二次信号分解得到的三个信号分别进行加窗,得到三个加窗信号；对三个加窗信号分别进行对称离散傅里叶变换,得到三个频域信号；分别对三个频域信号的虚部进行求和,分别对三个频域信号的实部进行求和；根据虚部求和结果计算信号瞬时频率；根据计算的瞬时频率及实部求和结果计算信号的直流分量、原点相位和幅值。本发明实现了瞬时频率、直流分量、瞬时幅值和瞬时相位的计算,适用于平稳信号和非平稳信号,计算精度高,对信号长度几乎无要求,抗直流干扰能力强,计算时间短。(The invention discloses a time-frequency analysis method and a time-frequency analysis system. The method comprises the following steps: performing first signal decomposition on an input signal; performing second signal decomposition on the signal obtained by the first signal decomposition to obtain three signals with equal length; windowing three signals obtained by the second signal decomposition respectively to obtain three windowed signals; symmetric discrete Fourier transform is respectively carried out on the three windowed signals to obtain three frequency domain signals; summing the imaginary parts of the three frequency domain signals respectively, and summing the real parts of the three frequency domain signals respectively; calculating the instantaneous frequency of the signal according to the result of the summation of the imaginary parts; and calculating the direct current component, the origin phase and the amplitude of the signal according to the calculated instantaneous frequency and the real part summation result. The invention realizes the calculation of instantaneous frequency, direct current component, instantaneous amplitude and instantaneous phase, is suitable for stable signals and non-stable signals, has high calculation precision, almost no requirement on signal length, strong direct current interference resistance and short calculation time.)

1. A time-frequency analysis method is characterized by comprising the following steps:

acquiring an input signal, identifying the signal length of the input signal, and performing first signal decomposition on the input signal;

performing second signal decomposition on the signal obtained by the first signal decomposition to obtain three signals with equal length;

windowing three signals obtained by the second signal decomposition respectively to obtain three windowed signals;

symmetric discrete Fourier transform is respectively carried out on the three windowed signals to obtain three frequency domain signals;

summing the imaginary parts of the three frequency domain signals respectively, and summing the real parts of the three frequency domain signals respectively;

calculating the instantaneous frequency of the signal according to the summation result of the three imaginary parts of the three frequency domain signals;

and calculating the direct current component, the origin phase and the amplitude of the signal according to the calculated instantaneous frequency and the real parts of the three frequency domain signals.

2. The time-frequency analysis method according to claim 1, wherein the first signal decomposition comprises the steps of:

the input signal before the first signal decomposition is denoted as S (i), the length of the input signal S (i) is denoted as K, and the signal S (i) is decomposed into K-L +1 short signals s with the length of L_jAs an output signal of the first signal decomposition, wherein j takes the value of 1,2,3, …, K-L + 1; the method comprises the following specific steps: selecting S (1) to S (L) as the 1 st short signal S₁Selecting S (2) to S (L +1) as the 2 nd short signal S₂Selecting S (3) to S (L +2) as the 3 rd short signal S₃And so on until a K-L +1 short signal s is obtained_K-L+1(ii) a Wherein L is an odd number.

3. The time-frequency analysis method according to claim 1, wherein the second signal decomposition comprises the steps of:

the signal obtained by the first signal decomposition is recorded as s_jSignal s_jIs denoted as L, signal s_jDecomposed into three signals of length N, denoted x respectively₁，x₂And x₃Wherein the signal x₁For the 1 st to L-2k th samples of the input sequence, where the signal x₂For the 1+ k to L-k samples of the input sequence, where the signal x₃Is the 1+2k to lth sample of the input sequence, where k is a positive integer and N is an odd number.

4. A time-frequency analysis method according to claim 1, wherein the discrete fourier transform comprises the steps of:

let three windowed signals be respectively denoted as y₁，y₂And y₃Let the three frequency domain signals be respectively denoted as Y₁，Y₂And Y₃The length of the windowed signal is recorded as N, and the frequency domain signal Y₁The calculation formula of (2) is as follows:

frequency domain signal Y₂The calculation formula of (2) is as follows:

frequency domain signal Y₃The calculation formula of (2) is as follows:

wherein, the value of m is { m belongs to Z | - (N-1)/2 ≦ m ≦ N-1)/2 }.

5. A time-frequency analysis method according to claim 4, characterized in that the calculation formula for summing the imaginary parts of the three frequency domain signals is:

wherein Im () represents an imaginary component operation; im is₁As a frequency domain signal Y₁Result of the summation of imaginary parts, Im₂As a frequency domain signal Y₂Result of the summation of imaginary parts, Im₃As a frequency domain signal Y₃The result of the imaginary summation of (a);

the calculation formula for summing the real parts of the three frequency domain signals is as follows:

Re₁as a frequency domain signal Y₁Sum of real parts of (Re)₂As a frequency domain signal Y₂Sum of real parts of (Re)₃As a frequency domain signal Y₃The real part of (d) is summed.

6. A time-frequency analysis method according to claim 5, wherein said calculated instantaneous signal frequency is obtained according to the following equation:

is a signal x₁Sum signal x₂The phase difference of (a) is greater than (b),is a signal s_jA is the signal s_jGamma is a predetermined parameter, the parameter havingAnd (4) preparing translation invariance.

7. The time-frequency analysis method according to claim 6, wherein the DC component, the phase of the origin and the amplitude of the calculated signal are obtained according to the following equation:

wherein DC represents the signal s_jOf the dc component.

8. The time-frequency analysis method according to claim 1, further comprising the steps of:

after calculating the phase of the origin, phase unwrapping is carried out on the phase of the origin, and the instantaneous frequency of the signal is recalculated according to the phase after phase unwrapping.

9. The time-frequency analysis method according to claim 8, wherein the phase unwrapping is obtained by calculating discrete partial derivatives of the phase.

10. A time-frequency analysis system, comprising:

the first signal decomposition module is used for acquiring an input signal, identifying the signal length of the input signal and performing first signal decomposition on the input signal;

the second signal decomposition module is used for carrying out second signal decomposition on the signal obtained by the first signal decomposition to obtain three signals with equal length;

the windowing module is used for respectively windowing the three signals obtained by the second signal decomposition to obtain three windowed signals;

the symmetrical discrete Fourier transform module is used for respectively carrying out symmetrical discrete Fourier transform on the three windowed signals to obtain three frequency domain signals;

the summation module is used for respectively summing the imaginary parts of the three frequency domain signals and respectively summing the real parts of the three frequency domain signals;

the first calculation module is used for calculating the instantaneous frequency of the signal according to the summation result of the three imaginary parts of the three frequency domain signals;

and the second calculation module is used for calculating the direct current component, the origin phase and the amplitude of the signal according to the summation result of the three real parts of the three frequency domain signals and the instantaneous frequency obtained by calculation.

Technical Field

The invention belongs to the technical field of signal analysis, and particularly relates to a time-frequency analysis method and a time-frequency analysis system.

Background

Signals are the carriers of information, and it is recognized that everything in the world and communication are based on an understanding of the information. Acquiring information in a signal relies on signal analysis. Usually, the actual signals are time-varying signals, and many scientists use non-stationary signals to refer to time-varying signals. Time-frequency analysis is widely used for analyzing non-stationary signals, such as gravitational waves, animal languages, brain waves and other signals. With the understanding of the world, the requirement for the time-frequency analysis method is higher and higher.

Various time-frequency analysis methods exist in the prior art. For example, Linear Time-frequency representation, Bilinear Time-frequency distribution, Time-varying high-order spectroscopy, Adaptive parameter Time-frequency analysis, Time-frequency ARMA model, Adaptive non-parametric Time-frequency analysis, and the like.

The existing time frequency analysis method has the following problems: there is a large error in the calculation of the frequency parameter for non-stationary signals. Constrained by the heisenberg uncertainty criterion (the resolution in time and frequency domains cannot reach the highest level at the same time), the existing method cannot obtain accurate local frequency parameters (amplitude, frequency and phase) by using a few samples. When the signal is a stable signal, the existing time frequency analysis method can give a more accurate result. When the signal is a non-stationary signal, an analysis result given by the existing time-frequency analysis method has an error. Generally, the faster the frequency parameter of the signal changes, the greater the error. The theoretical basis of time-frequency analysis of non-stationary signals is that any non-stationary signal can be regarded as a short-time stationary signal and the head and the tail are connected, and more samples are selected to shake the theoretical basis, so that the estimation of frequency parameters of the non-stationary signals has larger errors.

The existing time frequency analysis method has the following problems: amplitude and phase change information is lost. The amplitude reflects the energy intensity of the signal and the phase reflects the position of the signal. Taking short-time fourier transform (STFT) as an example, the time-frequency diagram can only approximately reflect the relation between time and frequency, and the amplitude and phase information is lost.

Disclosure of Invention

Aiming at least one defect or improvement requirement in the prior art, the invention provides a time-frequency analysis method and a time-frequency analysis system, wherein the time-frequency analysis is carried out by adopting symmetric discrete Fourier transform, the frequency parameter calculation is realized by utilizing the integral characteristic of the symmetric discrete Fourier transform, the instantaneous amplitude and the instantaneous phase can also be calculated, the method is suitable for stable signals and non-stable signals, the calculation precision is high, the signal length is not required, the direct-current interference resistance is strong, and the calculation time is short.

To achieve the above object, according to a first aspect of the present invention, there is provided a time-frequency analysis method, including the steps of:

acquiring an input signal, identifying the signal length of the input signal, and performing first signal decomposition on the input signal according to the signal length;

performing second signal decomposition on the signal obtained by the first signal decomposition to obtain three signals with equal length;

windowing three signals obtained by the second signal decomposition respectively to obtain three windowed signals;

symmetric discrete Fourier transform is respectively carried out on the three windowed signals to obtain three frequency domain signals;

summing the imaginary parts of the three frequency domain signals respectively, and summing the real parts of the three frequency domain signals respectively;

calculating the instantaneous frequency of the signal according to the summation result of the three imaginary parts of the three frequency domain signals;

and calculating the direct current component, the origin phase and the amplitude of the signal according to the summation result of the three real parts of the three frequency domain signals and the calculated instantaneous frequency.

Preferably, the first signal decomposition includes the steps of:

the input signal before the first signal decomposition is denoted as S (i), the length of the input signal S (i) is denoted as K, and the signal S (i) is decomposed into K-L +1 short signals s with the length of L_jAs an output signal of the first signal decomposition, wherein j takes the value of 1,2,3, …, K-L + 1; the method comprises the following specific steps: selecting S (1) to S (L) as the 1 st short signal S₁Selecting S (2) to S (L +1) as the 2 nd short signal S₂Selecting S (3) to S (L +2) as the 3 rd short signal S₃Thus, according toGo on until K-L +1 short signal s is obtained_K-L+1(ii) a Wherein L is an odd number.

Preferably, the second signal decomposition comprises the steps of:

the signal obtained by the first signal decomposition is recorded as s_jSignal s_jIs denoted as L, signal s_jDecomposed into three signals of equal length, denoted x respectively₁，x₂And x₃Wherein the signal x₁For the 1 st to L-2k th samples of the input sequence, where the signal x₂For the 1+ k to L-k samples of the input sequence, where the signal x₃Is the 1+2k to lth sample of the input sequence, where k is a positive integer and N is an odd number.

Preferably, the discrete fourier transform comprises the steps of:

let three windowed signals be respectively denoted as y₁，y₂And y₃Let the three frequency domain signals be respectively denoted as Y₁，Y₂And Y₃The length of the windowed signal is recorded as N, and the frequency domain signal Y₁The calculation formula of (2) is as follows:

frequency domain signal Y₂The calculation formula of (2) is as follows:

frequency domain signal Y₃The calculation formula of (2) is as follows:

wherein, the value of m is { m belongs to Z | - (N-1)/2 ≦ m ≦ N-1)/2 }.

Preferably, the calculation formula for summing the imaginary parts of the three frequency domain signals is as follows:

wherein Im () represents an imaginary component operation; im is₁As a frequency domain signal Y₁Result of the summation of imaginary parts, Im₂As a frequency domain signal Y₂Result of the summation of imaginary parts, Im₃As a frequency domain signal Y₃The result of the imaginary summation of (a);

the calculation formula for summing the real parts of the three frequency domain signals is as follows:

Re₁as a frequency domain signal Y₁Sum of real parts of (Re)₂As a frequency domain signal Y₂Sum of real parts of (Re)₃As a frequency domain signal Y₃The real part of (d) is summed.

Preferably, the calculated signal phase difference and instantaneous frequency are obtained according to the following equation system:

is a signal x₁Sum signal x₂The phase difference of (a) is greater than (b),is a signal s_jA is the signal s_jIs a predetermined parameter, which has a translational invariance.

Preferably, the dc component, the origin phase and the amplitude of the calculation signal are obtained according to the following equation:

where DC represents the direct current component in the signal.

Preferably, the time-frequency analysis method further comprises the steps of:

after calculating the phase of the origin, phase unwrapping is carried out on the phase of the origin, and the instantaneous frequency of the signal is recalculated according to the phase after phase unwrapping.

Preferably, the phase unwrapping is obtained by calculating discrete partial derivatives of the phase.

According to a second aspect of the present invention, there is provided a time-frequency analysis system, comprising:

the first signal decomposition module is used for acquiring an input signal, identifying the signal length of the input signal and performing first signal decomposition on the input signal;

the second signal decomposition module is used for carrying out second signal decomposition on the signal obtained by the first signal decomposition to obtain three signals with equal length;

the windowing module is used for respectively windowing the three signals obtained by the second signal decomposition to obtain three windowed signals;

the symmetrical discrete Fourier transform module is used for respectively carrying out symmetrical discrete Fourier transform on the three windowed signals to obtain three frequency domain signals;

the summation module is used for respectively summing the imaginary parts of the three frequency domain signals and respectively summing the real parts of the three frequency domain signals;

the first calculation module is used for calculating the instantaneous frequency of the signal according to the summation result of the three imaginary parts of the three frequency domain signals;

and the second calculation module is used for calculating the direct current component, the origin phase and the amplitude of the signal according to the summation result of the three real parts of the three frequency domain signals and the instantaneous frequency obtained by calculation.

In general, compared with the prior art, the invention has the following beneficial effects:

(1) the invention is an analytic algorithm, when the signal is a steady signal, an analytic solution with extremely high precision can be obtained, and the error of frequency is less than 1 x 10-12 Hz;

(2) when the signal is a non-stationary signal, although the instantaneous frequency obtained by direct calculation has a large error, the instantaneous phase has extremely high precision, and the instantaneous frequency is recalculated by differentiating the instantaneous phase, namely, the instantaneous frequency is calculated indirectly, which shows that the recalculated instantaneous frequency has higher precision.

(3) The method can calculate the instantaneous amplitude and the instantaneous phase while calculating the instantaneous frequency, and is superior to the traditional time-frequency signal analysis method based on the time-frequency diagram.

(4) The invention has almost no requirement on the length of the signal, and can calculate and obtain a local frequency value by adopting five samples and also can obtain a local amplitude, a local phase and a direct current component.

(5) The invention is an analytic algorithm, which uses the frequency of the signal calculated by the imaginary part of SDFT without the DC component information to eliminate the influence of the DC component on the frequency estimation, thus the algorithm has strong DC interference resistance.

(6) The calculation time of the invention is very short, and the total time of 500 independent simulation tests is 19.6 seconds.

(7) The method can be used for correcting the sensor, and although the zero drift problem of the sensor can be corrected by adopting a method of subtracting the mean value of the signal, the error still exists after the correction. The invention provides a method for calculating the analytic solution of a direct current component in a signal, which can completely eliminate the problem of zero drift of a sensor. Assuming that the amplitude of the simple harmonic signal is 1 unit, the calculation error of the dc component is less than 1 × 10-10 units.

Drawings

FIG. 1 is a flow chart of a time-frequency analysis method according to an embodiment of the present invention;

fig. 2 is a schematic diagram of a second signal decomposition according to an embodiment of the present invention.

Detailed Description

In order to make the objects, technical solutions and advantages of the present invention more apparent, the present invention is described in further detail below with reference to the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are merely illustrative of the invention and are not intended to limit the invention. In addition, the technical features involved in the embodiments of the present invention described below may be combined with each other as long as they do not conflict with each other.

The embodiment of the invention provides a time frequency analysis method, which is used for solving the problems that the existing time frequency analysis method cannot adopt a few samples to calculate instantaneous frequency parameters, the direct current component in a signal is calculated and the like.

Referring to fig. 1, an embodiment of the present invention provides a time-frequency analysis method, which includes a signal first decomposition step, a signal second decomposition step, a signal windowing step, a symmetric discrete fourier transform step, a frequency domain integration and summation step, a frequency calculation step, a step of calculating other frequency parameters, a phase unwrapping step, a frequency recalculation step, and nine steps in total. The last two steps are not necessary, when the stationary signal analysis is carried out, the frequency, the phase and the amplitude of the signal can be calculated by adopting the first seven steps, and the frequency domain parameters obtained by adopting the algorithm have extremely high precision. When non-stationary signal analysis is carried out, the instantaneous frequency, the instantaneous phase and the instantaneous amplitude of the signal can be calculated by adopting the first seven steps, and the instantaneous frequency can be recalculated by adopting the seventh step and the eighth step, so that the calculation accuracy of the instantaneous frequency is improved. The preferred implementation of the various steps is as follows:

(1) signal first decomposition step:

obtaining input signal, assuming that the length of input signal S (i) is K, decomposing signal S (i) into K-L +1 short signals s with length L_j(ii) a The method comprises the following specific steps: selecting S (1) to S (L) as the 1 st short signal S₁Selecting S (2) to S (L +1) as the 2 nd short signal S₂Selecting S (3) to S (L +2) as the 3 rd short signal S₃And so on until a K-L +1 short signal s is obtained_K-L+1(ii) a Algorithm of the invention for signal s_jThe length (L) of (2) is almost not required, and only L is required to be an odd number of 5 or more.

(2) And signal second decomposition:

suppose the input signal is s_jSignal s_jIs L; a second signal decomposition step of decomposing the input signal s_jDecomposed into three signals (x) of equal length₁，x₂And x₃) As shown in fig. 2; wherein the signal x₁For the 1 st to L-2k th samples of the input sequence, where the signal x₂For the 1+ k to L-k samples of the input sequence, where the signal x₃Is the 1+2k th to Lth samples of the input sequence; where k is an integer of 1 or more, it is preferable that the value of k be 1 in actual calculation.

(3) Signal windowing:

the window function: in signal processing, the window function is a real-valued function that is all 0 outside a given interval; the length of a given interval is called the length of the window function, which is usually equal to the length of the signal; in the discrete signal analysis, the length of the signal is the number of samples; by inquiring the signal analysis teaching material, window functions in various forms can be obtained; common window functions include rectangular, hanning, gaussian, hamming, and triangular window.

Three signals (x) obtained by the decomposition₁，x₂And x₃) Sequentially windowing to sequentially obtain y₁，y₂And y₃The signal windowing is to perform a Hadamard product on the signal and a window function, wherein the Hadamard product is a vector multiplication operation, namely, corresponding elements of two signals with equal length are multiplied to obtain a signal with the length equal to that of the two signals; the method is not limited to a certain window function, and preferably a rectangular window is used as a default window function of the algorithm in the embodiment of the invention;

(4) a symmetrical discrete Fourier transform step:

the Continuous Fourier Transform (CFT) is a mathematical transform whose infinite integration results in poor practical applicability. In the digital age, the Fourier transforms used in various industries are Discrete Fourier Transforms (DFT). There are currently two main forms of DFT, the first being the conventional Discrete Fourier Transform (ODFT) and the second being the Symmetric Discrete Fourier Transform (SDFT). ODFT is widely used in Fast Fourier Transform (FFT), SDFT is less used, and existing time-frequency analysis methods are all based on ODFT. Compared to ODFT, SDFT is much closer to CFT, with symmetry and integration characteristics similar to CFT, whereas ODFT does not have these properties.

This step combines the three windowed signals (y)₁，y₂And y₃) Sequentially performing Symmetric Discrete Fourier Transform (SDFT) to obtain Y₁，Y₂And Y₃(ii) a With signal y₁For example, signal y₁The symmetric discrete Fourier transform is shown in a formula (1), wherein the value of m is { m belongs to Z | - (N-1)/2 ≦ m ≦ N-1)/2 };

existing fast algorithms for Discrete Fourier Transform (DFT) are fast discrete fourier transform (FFT), which is a fast algorithm for conventional discrete fourier transform (ODFT), and is not suitable for SDFT; the calculation time can be reduced by adopting a fast algorithm of SDFT; theoretically, the fast algorithms of the SDFT are various and cannot be listed one by one; all fast algorithms used to calculate equation (1) should be equivalent to this step.

(5) And a frequency domain integral summation step:

in sequence to Y₁，Y₂And Y₃Is subjected to integral summation to obtain Im₁，Im₂And Im₃(ii) a In sequence to Y₁，Y₂And Y₃Is subjected to integral summation to obtain Re₁，Re₂And Re₃(ii) a With Y₁For example, the imaginary part integral calculation formula is shown as formula (2); wherein Im () represents an imaginary component operation;

with Y₁For example, the partial integral calculation formula is shown in formula (3),

the imaginary parts of the frequency domain signals (Y1, Y2 and Y3) are odd symmetric, so that the imaginary parts of Y1, Y2 and Y3 and the calculation formula are respectively three to calculate Im₁For example, the other two calculation formulas are shown in formulas (4) and (5); im is₂And Im₃Similar calculation formulas are provided, which are not listed; wherein Im () represents fetchPerforming imaginary part operation;

(6) and a frequency calculating step:

according to the SDFT frequency domain imaginary part integral characteristic, the parameters of the simple harmonic signal have the following three equation relations:

the real parts of the frequency domain signals (Y1, Y2 and Y3) are even-symmetric, so that the real parts of Y1, Y2 and Y3 and the calculation formulas are respectively three in three to calculate Re₁For example, the other two calculation formulas are shown in formulas (9) and (10); re₂And Re₃Similar calculation formulas are provided, which are not listed; wherein Re () represents the real part operation;

wherein N is the length of the equal length signal of three,is a signal x₁Sum signal x₂The phase difference of (a) is greater than (b),is a signal x₂Is the phase of the origin of (A) is the signal x₂γ is a parameter with translational invariance; the method is a phase difference method, and the frequency of a signal can be calculated by utilizing the phase difference of the signal in a certain time period; signal x₁Sum signal x₂Is equal to the signal x₂Sum signal x₃The time interval T between the two signals is k/fs, wherein fs is the sampling frequency; signal x₁Sum signal x₂Is not equal toThe calculation formula is shown as formula (11);

according to the principle of the phase difference method, the frequency calculation formula of the signal is shown as formula (12);

equations (6) - (8) can be simplified into a nonlinear system of equations with and without three unknown parameters, and thus the phase differenceThere are many calculation methods, and the formula (11) is adopted to calculate the phase differenceBut one of them is a simpler form; when the equations (6) - (8) are adopted to construct an equality relationship to obtain the phase difference, the algorithm can be regarded as being equivalent to the algorithm of the invention;

(7) calculating other frequency parameters:

according to the SDFT frequency domain real part integral characteristic, the parameters of the simple harmonic signal have the following three equation relations:

wherein N is the length of the signal with three equal lengths, DC represents the direct current component in the signal, and the phase difference obtained by the calculation of the formula (7) isSubstituting equations (13) - (15) into, the signal s can be obtained_jDC, phase of originAnd the calculation formula of the amplitude A are shown in formulas (16) to (18);

calculating DC component DC, origin phase using equations (13) - (15)And amplitude A, according toDifference in solving order, DC component, phase of originThe calculation formula of the sum amplitude A has various forms; in this embodiment, the DC component DC is calculated first, and in practice, the phase or amplitude may also be calculated first; the method adopts equations (13) - (15) to construct an equality relationship to obtain signal parameters, and the algorithm can be regarded as equivalent to the algorithm of the invention;

(8) phase unwrapping step:

phase unwrapping: with a simple harmonic signalFor example, where A is the amplitude of the simple harmonic signal, f is the frequency of the simple harmonic signal, π is the circumference,in order to be the initial phase of the signal,for the phase at time t, the range of the phase is (— ∞, ∞) theoretically; the phase is constrained to its principal value in the actual calculation (whether it is the interval (-pi, pi)]Or [0,2 π ])), such a phase is called the phase of the wrap; the process of reducing the winding phase to the range (— infinity, infinity) is phase unwrapping. Discrete partial derivatives of phase are widely used for phase unwrapping. When the phase of the simple harmonic signal is limited to its main value, in the range (-pi, pi)]For example, the phase changes periodically. The process from-pi to pi is a gradual change process, and pi is mutated into-pi and then gradually changed into pi, and the process appears repeatedly and is periodic; the change profile is obvious, the layers are uniform, and the catastrophe points are phase period demarcation points. Therefore, under ideal conditions, the discrete partial derivative of the phase can be extracted, and the partial derivative is integrated, so that the purpose of phase unwrapping can be achieved.

Preferably, the phase unwrapping is performed by calculating discrete partial derivatives of the phase, the calculation being as shown in equation (19);

wherein the value range of i is (1: K-N); then, the size of delta phi (i) is judged, if the delta phi (i) is smaller than-pi, the value of delta phi (i) is judgedTo2 pi is increased;

(9) and a frequency recalculation step:

a calculation formula for calculating the instantaneous frequency by differentiating the instantaneous phase (recalculating the instantaneous frequency) is shown in formula (20),

wherein the value range of i is (2: K-N);

the time-frequency analysis system of the embodiment of the invention comprises:

the first signal decomposition module is used for acquiring an input signal, identifying the signal length of the input signal and performing first signal decomposition on the input signal;

the second signal decomposition module is used for carrying out second signal decomposition on the signal obtained by the first signal decomposition to obtain three signals with equal length;

the windowing module is used for respectively windowing the three signals obtained by the second signal decomposition to obtain three windowed signals;

the symmetrical discrete Fourier transform module is used for respectively carrying out symmetrical discrete Fourier transform on the three windowed signals to obtain three frequency domain signals;

the summation module is used for respectively summing the imaginary parts of the three frequency domain signals and respectively summing the real parts of the three frequency domain signals;

the first calculation module is used for calculating the instantaneous frequency of the signal according to the summation result of the three imaginary parts of the three frequency domain signals;

and the second calculation module is used for calculating the direct current component, the origin phase and the amplitude of the signal according to the summation result of the three real parts of the three frequency domain signals and the instantaneous frequency obtained by calculation.

The implementation principle and technical effect of the system are the same as those of the method, and are not described herein again.

It must be noted that in any of the above embodiments, the methods are not necessarily executed in order of sequence number, and as long as it cannot be assumed from the execution logic that they are necessarily executed in a certain order, it means that they can be executed in any other possible order.

It will be understood by those skilled in the art that the foregoing is only a preferred embodiment of the present invention, and is not intended to limit the invention, and that any modification, equivalent replacement, or improvement made within the spirit and principle of the present invention should be included in the scope of the present invention.