A time-frequency analysis method based on frequency-doubling symmetric wavelet

Abstract

The application provides a time-frequency analysis method based on frequency multiplication symmetric wavelet. The time-frequency analysis is an important analysis method for non-stationary signals, and its essence is to filter the signal by using a band-pass filter family. In the method provided by the application, the way of constructing the band-pass filter family is different from the previous method, specifically: the minimum frequency and the maximum frequency of the signal to be analyzed are obtained by Fourier transform, and the minimum octave frequency and the maximum octave frequency are obtained by taking the logarithm with base 2, and then a plurality of octave symmetric wavelets with half-interval overlap in the octave frequency domain are generated according to the preset octave number of a single filter, which are used as a band-pass filter family to filter and decompose the signal, and after obtaining the time domain decomposition signal, the Hilbert transform is performed and the instantaneous amplitude spectrum is obtained, and the time-frequency spectrum is generated based on the instantaneous amplitude spectrum of each time domain decomposition signal. The application can enrich the time-frequency analysis theory, and is expected to expand the application range of the time-frequency analysis method.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of signal processing, in particular to a time-frequency analysis method based on frequency-doubling symmetric wavelet. BACKGROUND

[0002] Fourier transform (FT) is widely used in various fields, and the basic idea is to decompose a signal into a series of superimposed continuous sinusoidal waves with different frequencies. The defect of FT is that the time information is lost, and it is impossible to determine when a specific signal occurs according to the transform result, that is, the FT transform is only suitable for analyzing stationary signals. In practice, most signals (such as seismic signals, electrocardiogram signals, etc.) are non-stationary signals, containing a large number of non-stationary characteristics such as mutations, singularities, event initiation and termination, etc. That is, the frequency domain characteristics of non-stationary signals change with time, and it is necessary to extract the frequency domain information of a certain time period or the time information corresponding to a certain frequency band. Time-frequency analysis is an important analysis method for non-stationary signals, which can analyze signals from time domain and frequency domain at the same time, and the specific methods include short-time Fourier transform (STFT), S transform, generalized S transform, wavelet transform, HHT transform (Hilbert-Huang transform), etc.

[0003] One of the problems to be solved by time-frequency analysis method is how to analyze only a small part of data. The basic idea of STFT is to add a small window to the signal and transform the signal in the small window, so as to reflect the local characteristics of the signal. However, its defects are: the size and shape of the window function are fixed and do not change with time and frequency, which is not conducive to the analysis of time-varying signals! High-frequency signals have short duration, and low-frequency signals have long duration. For high-frequency signals, a small time window is used for analysis, and for low-frequency signals, a large time window is used for analysis. Wavelet transform inherits and develops the localization idea of STFT, and overcomes the defect that the window size does not change with frequency. Wavelet is a special wave shape with limited length and average value of 0. "Small" refers to having a compact support set or approximate compact support set in time domain, and "wave" refers to having positive and negative alternating "undulation", that is, the direct current component is zero.

[0004] Decomposition of signal is the way to realize the analysis of only a small part of data in time-frequency analysis method. The essence of signal decomposition is to filter the signal by using a family of band-pass filters. As long as the frequencies of these band-pass filters can cover the entire frequency spectrum range of the original signal, the original signal can be obtained by combining these signals with different frequencies according to their component sizes when inverse transforming. There are various ways to select the family of band-pass filters, and different time-frequency analysis methods adopt different frequency domain segmentation methods when selecting the family of filters.

[0005] HHT(Hilbert-Huang) transform is to decompose a signal into a finite number of intrinsic mode functions(IMF) by empirical mode decomposition(EMD). The EMD decomposition process can be compared to a sifting process: the upper envelope of the original data sequence is obtained by fitting all the local maximum points, the lower envelope of the original data sequence is obtained by fitting all the local minimum points, the average of the upper envelope and the lower envelope is taken as the mean envelope, the original data sequence is subtracted from the mean envelope to obtain a new data sequence, if the new data still has negative local maximum and positive local minimum, it is not an intrinsic mode component, and needs to continue sifting. Each sifting is equivalent to a filtering process of the data, but not the traditional frequency filtering.

[0006] Wavelet transform is to decompose a signal into a series of superimposed wavelet functions. These wavelet functions are obtained by translating and scaling a basic wavelet. In wavelet transform, the wavelet functions obtained by scaling and translating the basic wavelet are equivalent to different filters. Wavelet transform is to filter the original signal by using the wavelet functions obtained by scaling and translating the basic wavelet, and the result of filtering is the wavelet coefficients.

[0007] Continuous wavelet transform(CWT) theoretically calculates wavelet coefficients under each possible scale factor and translation parameter, which has a very large amount of calculation and data, and many data are useless. Discrete wavelet transform(DWT) only selects part of the scale factors and translation parameters for calculation, which greatly reduces the amount of data for analysis. If the scale factor and the translation parameter are both selected as integer multiples of 2, such DWT is called dyadic wavelet transform. The usual discrete wavelet transform refers to dyadic wavelet transform. The Mallat algorithm for performing discrete wavelet transform uses filters to decompose the signal, which is commonly referred to as two-channel subband coding in digital signal processing. One filter is a low-pass filter, through which the approximate value A of the signal can be obtained, and the other is a high-pass filter, through which the detail value D of the signal can be obtained. Multi-level signal decomposition is called multi-resolution analysis(MRA). MRA divides the signal bandwidth into high-pass and low-pass two parts(corresponding to two filters, while the filter spectrum shape corresponding to different wavelets will also be different), and then continues to divide the low-pass part. Figure 1 Figure 1 As shown in the figure, the wavelet transform MRA spectrum decomposition diagram, the filter result of the high-pass part is reserved each time the signal is divided(the corresponding detail), because the energy of the signal is mostly distributed in the low frequency part, so the high frequency part is not divided, and the low frequency part is continuously divided into its corresponding detail A and approximation D, and more details are obtained by multi-level division.

[0008] ​The essence of HHT and wavelet transform used as time-frequency analysis is to decompose the signal based on the designed filter, only that the wavelet in the wavelet transform can be expressed in time domain, while the filtered signal in HHT is directly obtained by using empirical mode decomposition (EMD), and the time and frequency form of the filter has no specific rule and is not easy to be visualized.

[0009] The precision of the time-frequency analysis method is limited by the similarity between the filter family and the signal. For different signals, different characteristic filters should be selected to decompose the signal. For wavelet transform, the filter family is the wavelet after scaling and shifting, and for HHT, the filter family corresponds to the screening process of different orders of empirical mode decomposition (EMD). For a certain signal, the time-frequency analysis result using a certain wavelet as a filter has sufficient precision, but for another signal, the result using the wavelet as a filter may not be ideal. In order to adapt to different signal characteristics, it is necessary and important to design different filters to decompose the signal and do time-frequency analysis, but the form of the signal in each application field is different, and the requirement for the form of the wavelet as a filter and the precision of time-frequency analysis is endless. SUMMARY

[0010] The purpose of the present application is to propose a new time-frequency analysis method based on octave symmetric wavelet, which uses multiple octave symmetric wavelets with half interval overlap in octave frequency domain to segment the frequency domain to form a filter family, uses the constructed octave symmetric wavelet filter family to perform band-pass filtering and decomposition on the signal, calculates the instantaneous amplitude of the decomposed signal and extracts the corresponding main frequency, and then obtains the time-frequency analysis spectrum. The present application enriches the theory of time-frequency analysis and is expected to expand the application range of time-frequency analysis method.

[0011] The present application provides a time-frequency analysis method based on octave symmetric wavelet, which comprises:

[0012] Obtaining the minimum frequency and the maximum frequency of the signal to be analyzed in the normal frequency domain;

[0013] According to the minimum frequency and the maximum frequency, the minimum octave frequency and the maximum octave frequency of the signal to be analyzed in the octave frequency domain are calculated;

[0014] According to the minimum octave frequency, the maximum octave frequency and the preset number of single wavelet octave, the number of octave symmetric wavelets contained in the entire signal frequency range and the center frequency of each octave symmetric wavelet are determined;

[0015] According to the number of single wavelet octave and the center frequency of each octave symmetric wavelet, the amplitude spectrum form of each octave symmetric wavelet is designed in the octave frequency domain;

[0016] According to the amplitude spectrum form of each octave symmetric wavelet in the octave frequency domain, the frequency spectrum of each octave symmetric wavelet in the normal frequency domain is calculated;

[0017] The octave symmetric wavelets in the normal frequency domain are taken as a filter bank, and the frequency spectrum of the filter bank is used to filter the signal to be analyzed to obtain a decomposition signal spectrum and a time domain decomposition signal.

[0018] The time domain decomposition signal is subjected to Hilbert transform to obtain an instantaneous amplitude spectrum.

[0019] A time-frequency spectrum of the signal to be analyzed is generated based on the instantaneous amplitude spectrum of the time domain decomposition signal.

[0020] Optionally, the minimum frequency and the maximum frequency of the signal to be analyzed in the normal frequency domain are obtained by:

[0021] The frequency range of the signal to be analyzed is observed after Fourier transform, and the minimum frequency and the maximum frequency of the signal to be analyzed in the normal frequency domain are denoted as fmi and fma.

[0022] Optionally, the minimum octave frequency and the maximum octave frequency of the signal to be analyzed in the octave frequency domain are calculated based on the minimum frequency and the maximum frequency, and the calculation comprises:

[0023] The minimum frequency fmi and the maximum frequency fma of the signal to be analyzed are subjected to logarithm with base 2 to obtain an operation result, and the operation result is taken as the minimum octave frequency lfmi and the maximum octave frequency lfma of the signal to be analyzed in the octave frequency domain, and the calculation formula is specifically as follows:

[0024] lfmi = log2(fmi), and lfma = log2(fma).

[0025] Optionally, the number of octave symmetric wavelets contained in the entire signal frequency range and the center frequency of each octave symmetric wavelet are determined based on the minimum octave frequency, the maximum octave frequency and a preset number of single wavelet octaves, and the determination comprises:

[0026] If the preset number of single wavelet octaves is n, the number of non-overlapping filters with n octaves that can be accommodated in the signal frequency range is obtained by m = (lfma-lfmi) / n, and considering the form characteristics of the filter to be designed and the fact that the amplitude spectrum of all filters should be 1 after superposition, one filter is added between every two non-overlapping filters, that is, m-1 filters are added, so the total number of filters with n octaves in the octave frequency domain is 2m-1; according to the above design, the center frequency of the 2m-1 filters with n octaves in the octave domain is 2m-1 numbers with interval n / 2 from lfmi+n / 2 to lfma-n / 2.

[0027] Optionally, the step of designing the amplitude spectrum of each octave-symmetry wavelet in the octave frequency domain according to the single wavelet octave number and the center frequency of each octave-symmetry wavelet comprises:

[0028] In order to avoid Gibbs effect and to have many good characteristics, the cosine function in the range of [-π, π] is selected as the amplitude spectrum of the octave-symmetry wavelet in the octave frequency domain, and the formula is:

[0029]

[0030] In the formula, A is the amplitude spectrum of the octave-symmetry wavelet in the octave frequency domain, n is the octave number, In the formula, A is the amplitude spectrum of the octave-symmetry wavelet in the octave frequency domain, n is the octave number, In the formula, f0 is the center frequency of the spectrum of the octave-symmetry wavelet, In the formula, f0 is the center frequency of the spectrum of the octave-symmetry wavelet, In the formula, f0 is the center frequency of the spectrum of the octave-symmetry wavelet.

[0031] Optionally, the step of obtaining the spectrum of each octave-symmetry wavelet in the normal frequency domain according to the amplitude spectrum of each octave-symmetry wavelet in the octave frequency domain comprises:

[0032] For the octave-symmetry wavelet in the 2m-1 octave frequency domain, the corresponding octave frequency is transformed back to the frequency value, and the amplitude value corresponding to the uniformly sampled frequency value in the frequency domain is calculated by using the interpolation method to obtain the amplitude spectrum of the octave-symmetry wavelet in the normal frequency domain. The phase spectrum of the octave-symmetry wavelet is set to zero, and the spectrum of the octave-symmetry wavelet in the normal frequency domain can be obtained from the amplitude spectrum and the phase spectrum.

[0033] Optionally, the step of filtering the signal to be analyzed by using the filter bank spectrum of the filter bank comprising the octave-symmetry wavelets in the normal frequency domain to obtain the decomposition signal spectrum and the time-domain decomposition signal comprises:

[0034] The signal to be analyzed is transformed from the time domain to the frequency domain by using the Fourier transform to obtain the signal spectrum to be analyzed;

[0035] The signal spectrum to be analyzed is multiplied by the filter bank spectrum to obtain the decomposition signal spectrum in the normal frequency domain;

[0036] The decomposition signal spectrum is subjected to the Fourier inverse transform to obtain the time-domain expression of the decomposition signal.

[0037] Optionally, the Hilbert transform is performed on the time-domain decomposed signals and the instantaneous amplitude spectrum thereof is obtained, comprising:

[0038] An analytic signal corresponding to each decomposed signal is obtained by using the Hilbert transform formula, and the instantaneous amplitude corresponding to the analytic signal is obtained;

[0039] The specific formula of the Hilbert transform is as follows:

[0040]

[0041] wherein s i (t) represents the decomposed signal, represents the Hilbert transform of the decomposed signal, i=1, 2, …2m-1 represents the i th decomposed signal;

[0042] The specific formula for obtaining the instantaneous amplitude is as follows:

[0043]

[0044] wherein a i (t) represents the instantaneous amplitude data of the i th decomposed signal, and t is the time sample value after discrete sampling.

[0045] Optionally, the time-frequency spectrum of the signal to be analyzed is generated based on the instantaneous amplitude spectrum of the time-domain decomposed signal, comprising:

[0046] All the decomposed signals correspond to the same time axis but different main frequencies. Taking the same time axis corresponding to all the decomposed signals as the horizontal coordinate, taking the different main frequencies corresponding to different decomposed signals as the vertical coordinate, and taking the instantaneous amplitude spectrum corresponding to different decomposed signals as the numerical value of the two-dimensional time-frequency spectrum, the time-frequency analysis spectrum of the signal to be analyzed based on the frequency-symmetrical wavelet is obtained.

[0047] From the above technical solutions, it can be seen that the embodiments of the present application have the following advantages:

[0048] The application provides a time-frequency analysis method based on frequency multiplication symmetric wavelet. The time-frequency analysis is an important analysis method for non-stationary signals, and its essence is to filter a signal by using a band-pass filter family. The time-frequency analysis method based on frequency multiplication symmetric wavelet provided by the application is different from the previous method in the way of constructing a band-pass filter family, specifically: the minimum frequency and the maximum frequency of a signal to be analyzed are obtained by Fourier transform, the minimum frequency multiplication range frequency and the maximum frequency multiplication range frequency are obtained by taking the logarithm with 2 as the base, and then a plurality of frequency multiplication symmetric wavelets with half-interval overlap in the frequency multiplication range frequency domain are generated according to the preset frequency multiplication range number of a single filter, which are used as a band-pass filter family to filter and decompose the signal, the instantaneous amplitude of the decomposed signal is calculated, and the main frequency corresponding to each decomposed signal is extracted to draw a time-frequency spectrum. The application can enrich the time-frequency analysis theory and is expected to expand the application range of the time-frequency analysis method. BRIEF DESCRIPTION OF DRAWINGS

[0049] In order to more clearly illustrate the technical solutions in the embodiments of the application or the prior art, the following will briefly introduce the drawings needed to be used in the embodiments or the prior art description. Obviously, the drawings in the following description only constitute some embodiments of the application, and for those skilled in the art, other drawings can also be obtained without creative labor.

[0050] Figure 1 It is a wavelet transform MRA spectrum decomposition schematic diagram and filter shape of two different wavelets;

[0051] Figure 2 It is a flowchart of the signal time-frequency analysis method based on frequency multiplication symmetric wavelet provided by the embodiment of the application;

[0052] Figure 3 It is a comparison schematic diagram of the amplitude spectrum shape of each filter in the frequency multiplication range frequency domain and the frequency domain and the superposition of the amplitude spectrum shape when n=2 of the frequency multiplication range;

[0053] Figure 4 It is a time domain frequency multiplication symmetric wavelet schematic diagram when n=2 of the frequency multiplication range and different main frequencies provided by the embodiment of the application;

[0054] Figure 5 It is a decomposed signal spectrum schematic diagram after the signal is filtered by the frequency division filter with the frequency multiplication symmetric wavelet as the filter provided by the embodiment of the application;

[0055] Figure 6 It is a time domain decomposed signal schematic diagram after the signal is filtered by the frequency division filter with the frequency multiplication symmetric wavelet as the filter provided by the embodiment of the application;

[0056] Figure 7This is a schematic diagram of the time-frequency spectrum obtained by using a frequency-doubled symmetrical wavelet as a filter for the signal to be analyzed, as provided in an embodiment of this application. Detailed Implementation

[0057] The technical solutions of the embodiments of this application will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of this application, and not all embodiments. Based on the embodiments of this application, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of this application.

[0058] The accuracy of time-frequency analysis methods is limited by the similarity between the filter family and the signal. Different filters with different characteristics should be selected to decompose the signal for different signals. For wavelet transform, the filter family consists of wavelets after scaling and translation; for HHT, the filter family corresponds to the selection process of different orders of Empirical Mode Decomposition (EMD). For a certain signal, the time-frequency analysis results using a certain wavelet as a filter may have sufficient accuracy, but for another signal, the results may be unsatisfactory. To adapt to different signal characteristics, designing different filters to decompose the signal and perform time-frequency analysis is necessary and important. However, the forms of signals in various application fields vary greatly, and the requirements for wavelet morphology and time-frequency analysis accuracy as filters are endless.

[0059] Based on this, this application proposes the following technical solution, as detailed below:

[0060] In one embodiment, such as Figure 2 As shown, Figure 2 This is a flowchart illustrating a signal time-frequency analysis method provided in an embodiment of this application. The present application provides a signal time-frequency analysis method, specifically including the following:

[0061] S110: Obtain the minimum and maximum frequencies of the signal to be analyzed in the normal frequency domain.

[0062] Specifically, in this step, the signal to be analyzed can be acquired first, and then the minimum and maximum frequencies of the signal to be analyzed in the ordinary frequency domain can be obtained. The signal to be analyzed is a non-stationary signal.

[0063] S120: Based on the minimum and maximum frequencies, calculate the minimum and maximum octave frequencies of the signal to be analyzed in the octave frequency domain.

[0064] Specifically, after obtaining the minimum and maximum frequencies of the signal to be analyzed in the ordinary frequency domain through step 110, the minimum and maximum octave frequencies of the signal to be analyzed in the octave frequency domain can be calculated based on the relationship between the ordinary frequency and the octave frequency.

[0065] wherein, if f1, f2 are the minimum frequency and the maximum frequency of the signal to be analyzed in the frequency domain, then is the minimum octave frequency and the maximum octave frequency of the signal to be analyzed in the octave frequency domain; is called octave.

[0066] S130: determining the number of octave symmetric wavelets contained in the entire signal frequency range and the center frequency of each octave symmetric wavelet according to the minimum octave frequency, the maximum octave frequency and the preset number of octaves contained in a single wavelet.

[0067] In this step, after obtaining the minimum octave frequency and the maximum octave frequency of the signal to be analyzed in the octave frequency domain through step S120, the present application can utilize the minimum octave frequency, the maximum octave frequency and the preset number of octaves contained in a single wavelet, and calculate the number of octave symmetric wavelets contained in the entire signal frequency range and the center frequency of each octave symmetric wavelet according to the morphological characteristics of the cos function.

[0068] Specifically, when calculating the number of octave symmetric wavelets in the octave frequency domain, the number of octaves contained in a single wavelet, the minimum octave frequency and the maximum octave frequency can be used to calculate the number of non-overlapping wavelets that can be accommodated in the octave frequency domain in the frequency band range of the signal to be analyzed. Due to the morphological characteristics of the cos function and the fact that the spectrum of all octave symmetric wavelets should be 1 after superposition, an octave symmetric wavelet can be supplemented between every two octave symmetric wavelets, thereby obtaining the total number of octave symmetric wavelets.

[0069] S140: designing the amplitude spectrum morphology of each octave symmetric wavelet in the octave frequency domain according to the number of octaves contained in a single wavelet and the center frequency of each octave symmetric wavelet.

[0070] Specifically, after obtaining the number of octave symmetric wavelets contained in the entire signal frequency range and the center frequency of each octave symmetric wavelet through step S130, the amplitude spectrum of each octave symmetric wavelet can be designed in the octave frequency domain. The reason for choosing to design the amplitude spectrum of the octave symmetric wavelet in the octave frequency domain rather than in the frequency domain is that the concept of "octave" and "symmetry" is simple and easy to implement in the octave frequency domain.

[0071] Illustratively, please refer to Figure 3 , Figure 3 the amplitude spectrum morphology of each filter in the octave frequency domain and the frequency domain and the superposition of the amplitude spectrum morphology provided by the embodiment of the present application when the octave n = 2, Figure 3 In the octave frequency domain, the filter amplitude spectrum morphology designed directly in the octave frequency domain highlights the characteristics of "octave" and "symmetry" more than the filter amplitude spectrum morphology converted into the frequency domain.

[0072] S150: Obtain the frequency spectrum of each frequency-symmetrical wavelet in the frequency domain according to the amplitude spectrum form of each frequency-symmetrical wavelet in the octave frequency domain.

[0073] Specifically, after the amplitude spectrum form of each frequency-symmetrical wavelet in the octave frequency domain is designed in step S140, the frequency spectrum form of each frequency-symmetrical wavelet in the frequency domain can be obtained. Since the filtering operation in the frequency domain can be simply expressed as the product of the filter spectrum and the signal spectrum, after the multiple frequency-symmetrical wavelets are designed in the octave frequency domain, the amplitude spectrum form of each frequency-symmetrical wavelet is converted from the octave frequency domain to the frequency domain, so as to realize the filtering process with the signal to be analyzed.

[0074] S160: Filter the signal to be analyzed by taking each frequency-symmetrical wavelet in the normal frequency domain as a filter family and taking the spectrum of each frequency-symmetrical wavelet as the spectrum of the filter family, to obtain the decomposition signal spectrum and the time-domain decomposition signal.

[0075] In this step, after the frequency spectrum form of each frequency-symmetrical wavelet in the frequency domain is obtained in step S150, the signal to be analyzed can be filtered by taking each frequency-symmetrical wavelet as a filter family and taking the spectrum of each frequency-symmetrical wavelet as the spectrum of the filter family, and then the spectrum of the filtered decomposition signal and the time-domain decomposition signal can be obtained. The time-domain band-pass filter family composed of the frequency-symmetrical wavelets changes the scale of the wavelets when the main frequency changes in the same octave, which is similar to the effect of the scale factor in the wavelet transform, and the translation is reflected in the filtering process expressed by the time-domain convolution operation (equivalent to the multiplication operation in the frequency domain).

[0076] It can be understood that the signal needs to be decomposed for time-frequency analysis, and the essence of signal decomposition is to filter the signal by using a band-pass filter family. As long as the frequency of the band-pass filter can cover the entire frequency spectrum range of the original signal, the original signal can be obtained by combining different frequency signals according to their component sizes in the inverse transform. In the present application, each frequency-symmetrical wavelet is used as a filter family to filter the non-stationary signal.

[0077] S170: Perform Hilbert transform on the time-domain decomposition signal and obtain the instantaneous amplitude spectrum thereof.

[0078] In this step, after the time-domain decomposition signal is obtained in step S160, the analytic signal corresponding to each decomposition signal can be obtained by using the Hilbert transform formula, and the corresponding instantaneous amplitude can be obtained.

[0079] The specific formula of the Hilbert transform is as follows:

[0080]

[0081] where s i(t) represents a decomposition signal, represents a Hilbert transform of a decomposition signal, i = 1, 2, … 2m-1 represents the i-th decomposition signal in the decomposition signal.

[0082] The formula for calculating the instantaneous amplitude is as follows:

[0083]

[0084] where a i (t) represents the instantaneous amplitude data of the i-th decomposition signal, t is the time sample value after discrete sampling.

[0085] S180: generating a time-frequency spectrum diagram of the signal to be analyzed based on the instantaneous amplitude spectrum of the time-domain decomposition signal.

[0086] In this step, after obtaining the instantaneous amplitude spectrum of the time-domain decomposition signal in step S170, the time-frequency spectrum diagram of the signal to be analyzed can be further obtained. All decomposition signals correspond to the same time axis, but correspond to different main frequencies. Taking the same time axis corresponding to all decomposition signals as the horizontal coordinate, taking the different main frequencies corresponding to different decomposition signals as the vertical coordinate, and taking the instantaneous amplitude spectrum corresponding to different decomposition signals as the numerical value of the two-dimensional time-frequency spectrum diagram, the time-frequency analysis spectrum diagram of the signal to be analyzed based on the frequency doubling symmetric wavelet is obtained.

[0087] In the above embodiment, time-frequency analysis is an important analysis method for non-stationary signals, and its essence is to filter the signal by using a band-pass filter family. The time-frequency analysis method based on the frequency doubling symmetric wavelet proposed in the present application is different from the previous method in the way of constructing the filter family, specifically: by Fourier transform, the minimum frequency and the maximum frequency of the signal to be analyzed are obtained, the logarithm with base 2 is taken to obtain the minimum octave frequency and the maximum octave frequency, and then according to the preset octave number of a single filter, a plurality of frequency doubling symmetric wavelets with overlapping half-interval in the octave frequency domain are generated, which are used as a band-pass filter family to filter and decompose the signal, the instantaneous amplitude of the decomposition signal is calculated, and the main frequency corresponding to each decomposition signal is extracted to draw a time-frequency spectrum diagram. The present application can enrich the theory of time-frequency analysis, and is expected to expand the application range of time-frequency analysis method.

[0088] In one embodiment, the minimum frequency and the maximum frequency of the signal to be analyzed in the frequency domain obtained in step S110 can include:

[0089] S111: observing the frequency range of the signal to be analyzed after Fourier transform, and recording fmi and fma as the minimum frequency and the maximum frequency of the signal to be analyzed in the normal frequency domain.

[0090] In the embodiment, when performing time-frequency analysis on the signal to be analyzed, Fourier transform can be performed on the signal to be analyzed first to obtain an amplitude spectrum of the signal to be analyzed, and the frequency range is observed to obtain the minimum frequency and the maximum frequency of the signal to be analyzed in the frequency domain, and the minimum frequency is recorded as fmi and the maximum frequency is recorded as fma.

[0091] In one embodiment, the calculation of the minimum octave frequency and the maximum octave frequency of the signal to be analyzed in the octave frequency domain according to the minimum frequency and the maximum frequency in step S120 can include:

[0092] S121: taking the minimum frequency fmi and the maximum frequency fma of the signal to be analyzed as logarithms with base 2 to obtain an operation result, and taking the operation result as the minimum octave frequency lfmi and the maximum octave frequency lfma of the signal to be analyzed in the octave frequency domain, and the calculation formula is specifically as follows:

[0093] lfmi=log2(fmi), lfma=log2(fma)

[0094] In the embodiment, after the minimum frequency fmi and the maximum frequency fma of the signal to be analyzed are obtained, the minimum octave frequency lfmi and the maximum octave frequency lfma of the signal to be analyzed in the octave frequency domain can be calculated according to the relationship between the octave frequency and the signal frequency. Taking the frequency as a logarithm with base 2 has the physical meaning of octave. After the logarithmic operation of the frequency, the high frequency is compressed to a smaller range, which is convenient for investigating the relationship between the frequency and the energy from the perspective of the octave.

[0095] In one embodiment, the determination of the number of octave symmetric wavelets contained in the entire signal frequency range and the center frequencies of the octave symmetric wavelets according to the minimum octave frequency, the maximum octave frequency and the preset number of wavelets in step S130 can include:

[0096] S131: if the preset number of wavelets is n, the number of non-overlapping filters with an octave number of n that can be accommodated in the signal frequency range is obtained by m=(lfma-lfmi) / n. Considering the morphological characteristics of the filter to be designed and the fact that the amplitude spectrum of all filters should be 1 after superposition, one filter is added between every two non-overlapping filters, that is, m-1 filters are added, so the total number of filters with an octave number of n in the octave frequency domain is 2m-1. According to the above design, the center frequencies of the 2m-1 filters with an octave number of n in the octave domain are 2m-1 values starting from lfmi+n / 2 to lfma-n / 2 with an interval of n / 2.

[0097] For example, a non-stationary signal to be analyzed varies in frequency as follows at different time periods:

[0098] x(t) = x1(t) + x2(t) + x3(t), where,

[0099] where the frequencies of the signal to be analyzed at different time periods are 40 Hz, 130 Hz and 220 Hz respectively, the minimum frequency of the frequency distribution range is 2 Hz, the maximum frequency is 512 Hz, i.e. fmi = 2, fma = 512, then the minimum octave frequency lfmi = 1, the maximum octave frequency lfma = 9, if each octave symmetric wavelet contains n = 2 octaves, then the number of octave symmetric wavelets without overlap in the octave frequency domain can be calculated as m = (lfma-lfmi) / n = 4. The amplitude spectrum of all filters should be a horizontal line with the same amplitude after summation, and the form of the amplitude spectrum of the octave symmetric wavelet in the octave frequency domain should adopt the cos function form in the range of [-π, π], considering the above two factors, m-1 octave symmetric wavelets are inserted between every two octave symmetric wavelets, so that the octave symmetric wavelets in the octave frequency domain present a form of half-interval overlap. For this embodiment, there are 2m-1 = 7 octave symmetric wavelets.

[0100] Specifically, in this embodiment, when the filters do not overlap, the minimum octave frequency of the first filter is 1, the maximum octave frequency is 3, and the center octave frequency is 2. If the case of infinite expansion of the minimum frequency to the low frequency end is considered, the first filter can be set as a low-pass filter without being limited to the center frequency and the minimum frequency. For details, please refer to Figure 3 (a) and Figure 3 (c); the minimum octave frequency of the second filter is 3, the maximum octave frequency is 5, and the center octave frequency is 4; the minimum octave frequency of the third filter is 5, the maximum octave frequency is 7, and the center octave frequency is 6; the minimum octave frequency of the fourth filter is 7, the maximum octave frequency is 9, and the center octave frequency is 8. If the case of infinite expansion of the maximum frequency to the high frequency end is considered, the last filter can be set as a high-pass filter without being limited to the center frequency and the maximum frequency. For details, please refer to Figure 3 (a) and Figure 3 (c). Because the spectrum of each octave symmetric wavelet in the octave frequency domain is in the form of a cos function, the spectrum summation of the four octave symmetric wavelets thus constructed cannot cover the entire frequency axis with equal amplitude, so a filter can be supplemented between every two filters in a half-interval overlap manner, and a total of 3 octave symmetric wavelets with center octave frequencies of 3, 5 and 7 are supplemented to form a family of octave symmetric wavelet filters with half-interval overlap in the octave frequency domain.

[0101] Further, all the frequency-doubling symmetry wavelet amplitude spectrums in the octave frequency domain are superimposed to verify that the superimposed frequency-doubling symmetry wavelet amplitude spectrums cover the entire frequency spectrum of the signal to be analyzed with equal amplitude, which can be referred to in detail in Figure 3 , Figure 3 (b) an amplitude spectrum superposition diagram of each filter in the octave frequency domain when the octave n = 2, provided by the embodiment of the present application, Figure 4 (d) an amplitude spectrum superposition diagram of each filter in the ordinary frequency domain when the octave n = 2, provided by the embodiment of the present application.

[0102] In one embodiment, the step S140 of designing the amplitude spectrum of each frequency-doubling symmetry wavelet in the octave frequency domain according to the number of octaves of the single wavelet and the center frequency of each frequency-doubling symmetry wavelet can include:

[0103] S141: In order to avoid the Gibbs effect and because the symmetrical amplitude spectrum has many good characteristics, a cos function form in the range of [-π, π] is selected as the amplitude spectrum form of the frequency-doubling symmetry wavelet in the octave frequency domain, and the formula is:

[0104]

[0105] In the formula, A is the amplitude spectrum value of the frequency-doubling symmetry wavelet in the octave frequency domain, n is the number of octaves, is the center octave frequency of the frequency-doubling symmetry wavelet, is the minimum octave frequency of the frequency-doubling symmetry wavelet in the octave frequency domain, is the maximum octave frequency of the frequency-doubling symmetry wavelet in the octave frequency domain. In the embodiment, the frequency-doubling symmetry wavelet is used as the filter of the signal to be analyzed. The smoother the spectrum form of the frequency-doubling symmetry wavelet is, the more coordinated the relationship between the frequency and the energy is, thereby avoiding the Gibbs effect. Since the single-period cosine function curve has the properties of symmetry and smoothness, the single-period cosine function curve is selected as the amplitude spectrum form of the frequency-doubling symmetry wavelet in the octave frequency domain. According to the preset number of octaves, the horizontal axis of the octave coordinate system is divided. Therefore, the form of each frequency-doubling symmetry wavelet in the octave frequency domain is the same, but the corresponding main frequency is different. The cos function form in the range of [-π, π] is used as the amplitude spectrum form of each frequency-doubling symmetry wavelet in the octave frequency domain.

[0106] In one embodiment, the step 150 of obtaining the spectrum of each frequency-doubling symmetry wavelet in the frequency domain according to the amplitude spectrum form of each frequency-doubling symmetry wavelet in the octave frequency domain can include:

[0107]

[0108] ​S151: for the octave frequency domain of the octave symmetric wavelet, the corresponding octave frequency is inverse transformed back to the frequency value, and the amplitude value corresponding to the uniformly sampled frequency value in the frequency domain is calculated by using the interpolation method to obtain the amplitude spectrum form of the octave symmetric wavelet in the ordinary frequency domain, and the phase spectrum of the octave symmetric wavelet is set to zero, so that the frequency spectrum of the octave symmetric wavelet in the ordinary frequency domain can be obtained from the amplitude spectrum and the phase spectrum.

[0109] In the embodiment, before the filter family is used for filtering the signal to be analyzed, the amplitude spectrum form in the ordinary frequency domain can be obtained according to the amplitude spectrum form in the octave frequency domain, and then the phase spectrum can be set, and finally the frequency spectrum of each octave symmetric wavelet in the ordinary frequency domain is obtained. In the present application, the phase spectrum of the octave symmetric wavelet is set to 0.

[0110] It can be understood that it is easier to realize the characteristics of "octave" and "symmetry" by directly designing the amplitude spectrum form of the filter in the octave frequency domain. However, in the subsequent filtering operation of the filter and the signal, the frequency domain filter spectrum and the signal spectrum are usually multiplied, so after the multiple octave symmetric wavelets are designed in the octave frequency domain, the amplitude spectrum form of each octave symmetric wavelet is converted from the octave frequency domain to the frequency domain, so as to realize the filtering process of the filter on the signal to be analyzed.

[0111] As shown in Figure 4 , the time domain octave symmetric wavelet of the present application is shown in the following figure: Figure 4 , which is the time domain octave symmetric wavelet of the present application with n=2 and different main frequencies. Figure 5 In the figure, the higher the main frequency, the larger the frequency range, and the more compact the wavelet in the time domain, and the better the time domain resolution when used for time-frequency analysis.

[0112] In one embodiment, the step 160 uses each octave symmetric wavelet as a filter family, and the spectrum of the filter family is used for filtering the signal to be analyzed to obtain the filtering decomposition spectrum and the time domain decomposition signal, which can include:

[0113] S161: using Fourier transform, the signal to be analyzed is transformed from the time domain to the frequency domain to obtain the frequency spectrum of the signal to be analyzed.

[0114] S162: multiplying the frequency spectrum of the signal to be analyzed with the frequency spectrum of the filter family to obtain the decomposition signal spectrum in the frequency domain.

[0115] S163: performing Fourier inverse transform on the decomposition signal spectrum to obtain the time domain expression of the decomposition signal.

[0116] In the embodiment, when filtering the signal to be analyzed, Fourier transform is used to transform the time-domain signal to be analyzed into frequency domain, so as to obtain the frequency spectrum of the signal to be analyzed in the frequency domain. Then, the frequency spectrum of the signal to be analyzed is multiplied by the frequency spectrum of the filter family, so as to obtain the frequency spectrum of the decomposition signal in the frequency domain. Since the instantaneous amplitude spectrum of the decomposition signal needs to be calculated in the time domain, Fourier inverse transform is used to convert the frequency spectrum of the decomposition signal into the time-domain decomposition signal.

[0117] Specifically, the filter family is used to filter the signal to be analyzed, which can be performed in the time domain or in the frequency domain. In the time domain, the filtering operation is represented as the convolution of the time-domain form of the signal to be analyzed and the time-domain form of the filter. In the frequency domain, the filtering operation is represented as the product of the frequency spectrum of the signal to be analyzed and the frequency spectrum of the filter. The frequency-domain filtering operation is used in the application, and the specific calculation formula is as follows:

[0118] S i (f) = S0(f) · H i (f), i = 1, 2, … 2m-1

[0119] In the formula, S0(f) is the frequency spectrum of the signal to be analyzed, H i (f), i = 1, 2, … 2m-1 is the frequency spectrum of each frequency-doubling symmetric wavelet in the filter family, S i (f), i = 1, 2, … 2m-1 is the frequency spectrum of the decomposition signal obtained by filtering the signal to be analyzed through the frequency-doubling symmetric wavelet filter family.

[0120] In a specific implementation, the Fourier inverse transform calculation formula used when converting the decomposition signal spectrum into the time-domain decomposition signal in the application can be as follows:

[0121]

[0122] In the formula, s i (t) is the time-domain decomposition signal, S i (f), i = 1, 2, … 2m-1 is the frequency spectrum of the decomposition signal.

[0123] As shown in Figure 6 and Figure 5 , Figure 6 is a decomposition signal spectrum diagram of the embodiment of the application after filtering the signal through the frequency-doubling symmetric wavelet as the filter, Figure 7 is a time-domain decomposition signal diagram of the embodiment of the application after filtering the signal through the frequency-doubling symmetric wavelet as the filter.

[0124] In one embodiment, the step S170 of performing Hilbert transform on the time domain decomposed signals and obtaining the instantaneous amplitude spectrum thereof can include:

[0125] S171: performing Hilbert transform on each decomposed signal by using the Hilbert transform formula, and obtaining the corresponding instantaneous amplitude thereof.

[0126] In a specific implementation, the present application can perform Hilbert transform on each decomposed signal by using the Hilbert transform formula, and the calculation formula of the Hilbert transform is as follows:

[0127]

[0128] wherein s i (t) represents the decomposed signal, represents the Hilbert transform of the decomposed signal, i = 1, 2, … 2m-1 represents the i-th decomposed signal.

[0129] Further, the instantaneous amplitude is obtained, and the calculation formula of the instantaneous amplitude of the decomposed signal after the Hilbert transform in the present application is specifically as follows:

[0130]

[0131] wherein a i (t) represents the instantaneous amplitude data of the i-th decomposed signal, and t is the time sample value after the discrete sampling.

[0132] For example, a non-stationary signal to be analyzed is filtered by a filter family composed of 2m-1 frequency doubling symmetrical wavelets, and 2m-1 decomposed signals are obtained. The instantaneous amplitude is obtained by performing Hilbert transform on the 2m-1 decomposed signals, and thus 2m-1 one-dimensional arrays a i (t) are obtained. Each one-dimensional array a i (t) contains the instantaneous amplitude value at the same time sample value as the signal to be analyzed.

[0133] In one embodiment, the step S180 of generating the time-frequency spectrum of the signal to be analyzed based on the decomposed signal instantaneous amplitude spectrum can include:

[0134] S181: all the decomposed signals correspond to the same time axis but different main frequencies. Taking the same time axis corresponding to all the decomposed signals as the horizontal coordinate, taking the different main frequencies corresponding to different decomposed signals as the vertical coordinate, and taking the instantaneous amplitude spectrum corresponding to different decomposed signals as the numerical value of the two-dimensional time-frequency spectrum, the time-frequency analysis spectrum of the signal to be analyzed based on the frequency doubling symmetrical wavelet is obtained.

[0135] Specifically, the 2m-1 one-dimensional arrays ai (t) constitute a two-dimensional array, along the row direction is the discrete time sample, along the column direction is a i (t) the main frequency of each frequency-symmetrical wavelet in the corresponding filter family, after being plotted in the form of contour lines, is a time-frequency spectrum with the horizontal axis representing time, the vertical axis representing the main frequency of each frequency-symmetrical wavelet, and the contour line amplitude representing the instantaneous amplitude.

[0136] As shown schematically in Figure 7 , the time-frequency spectrum of the analyzed signal is obtained by using frequency-symmetrical wavelets as filters. Figure 7 When a single filter contains a large number of frequency octaves, the number of frequency-symmetrical wavelet filters generated is small, and the discrete values of the frequency vertical axis in the time-frequency spectrum of the analyzed signal generated accordingly are also small, and the accuracy of the time-frequency spectrum is low. In order to improve the accuracy, the number of frequency octaves contained by each filter can be reduced, and the number of filters is correspondingly increased, and the accuracy of the time-frequency spectrum of the analyzed signal is also correspondingly improved. ​ In the embodiment, the number of frequency octaves of each filter is 0.5.

[0137] Finally, it should be noted that the relational terms herein such as first and second and the like are used solely to distinguish one entity or action from another entity or action without necessarily requiring or implying any such actual relationship or order between such entities or actions. Moreover, the terms "comprises", "comprising", or any other variations thereof, are intended to cover a non-exclusive inclusion, such that a process, method, article, or apparatus that comprises a list of elements does not include only those elements but can include other elements not expressly listed or inherent to such process, method, article, or apparatus. Without further limitation, an element preceded by "comprises a" does not, without more constraints, foreclose the existence of additional identical elements in the process, method, article, or apparatus that comprises the recited element.

[0138] The embodiments in the specification are described in a progressive manner, and each embodiment focuses on the differences from other embodiments. The embodiments can be combined as needed, and the same and similar parts refer to each other.

[0139] The above description of the disclosed embodiments enables a person skilled in the art to implement or use the present application. Various modifications to these embodiments will be apparent to those skilled in the art, and the general principles defined herein can be implemented in other embodiments without departing from the spirit or scope of the present application. Therefore, the present application will not be limited to the embodiments shown herein, but will conform to the widest scope consistent with the principles and novel features disclosed herein.

Claims

1. A time-frequency analysis method based on frequency-doubled symmetric wavelets, characterized in that, The method comprises: acquiring minimum frequency and maximum frequency of a signal to be analyzed in a normal frequency domain; calculating minimum octave frequency and maximum octave frequency of the signal to be analyzed in an octave frequency domain according to the minimum frequency and the maximum frequency; determining the number of octave symmetric wavelets contained in the whole signal frequency range and the center frequency of each octave symmetric wavelet according to the minimum octave frequency, the maximum octave frequency and a preset single wavelet octave number; designing the amplitude spectrum form of each octave symmetric wavelet in the octave frequency domain according to the single wavelet octave number and the center frequency of each octave symmetric wavelet; calculating the frequency spectrum of each octave symmetric wavelet in the normal frequency domain according to the amplitude spectrum form of each octave symmetric wavelet in the octave frequency domain; filtering the signal to be analyzed by taking each octave symmetric wavelet in the normal frequency domain as a filter family and taking the frequency spectrum thereof as filter family frequency spectrum, to obtain a decomposition signal frequency spectrum and a time domain decomposition signal; performing Hilbert transform on the time domain decomposition signal and calculating the instantaneous amplitude spectrum thereof; generating a time-frequency spectrum of the signal to be analyzed based on the instantaneous amplitude spectrum of the time domain decomposition signal.

2. The time-frequency analysis method of claim 1, wherein, The acquisition of the minimum frequency and the maximum frequency of the signal to be analyzed in the normal frequency domain comprises: The frequency range of the signal to be analyzed is observed after Fourier transformation and the minimum frequency and maximum frequency in the normal frequency domain are recorded and are the minimum frequency and maximum frequency in the normal frequency domain of the signal to be analyzed.

3. The time-frequency analysis method of claim 1, wherein, The calculation of the minimum octave frequency and the maximum octave frequency of the signal to be analyzed in the octave frequency domain according to the minimum frequency and the maximum frequency comprises: a minimum frequency of the signal to be analyzed and a maximum frequency taking the logarithm to the base 2, obtaining an operation result, and taking the operation result as a minimum octave frequency of the signal to be analyzed in the octave frequency domain and a maximum octave frequency The calculation formula is as follows: 。 4. The time-frequency analysis method of claim 1, wherein, The determination of the number of octave symmetric wavelets contained in the whole signal frequency range and the center frequency of each octave symmetric wavelet according to the minimum octave frequency, the maximum octave frequency and a preset single wavelet octave number comprises: If the preset octave number of a single wavelet is ,use The octave band that can be accommodated within the signal frequency range is obtained as follows: The number of non-overlapping filters should be determined by considering the shape characteristics of the proposed filters and the fact that the sum of the amplitude spectra of all filters should equal 1. Therefore, one filter should be added between every two non-overlapping filters, i.e., increasing the number of filters by 1. There are one filter, therefore the octave band in the frequency domain is... The total number of filters is One; based on the above design, The number of octaves is The center frequency of the filter in the octave domain From Beginning, to Termination, interval is of A number.

5. The time-frequency analysis method of claim 1, wherein, The design of the amplitude spectrum form of each octave symmetric wavelet in the octave frequency domain according to the single wavelet octave number and the center frequency of each octave symmetric wavelet comprises: To avoid the Gibbs effect, and because the symmetric amplitude spectrum has many good properties, the cos function shape in the range The cos function shape in the range [0, 1] is chosen as the amplitude spectrum shape of the octave frequency domain multiplier symmetric wavelet, and the formula is: ， ； wherein is the amplitude spectrum shape of the octave-symmetric wavelet in the octave frequency domain, n is the number of octaves, is the center frequency of the octave-symmetric wavelet in the octave frequency domain, is the minimum octave frequency of the octave-symmetric wavelet in the octave frequency domain, is the maximum octave frequency of the octave-symmetric wavelet in the octave frequency domain.

6. The time-frequency analysis method of claim 1, wherein, The calculation of the frequency spectrum of each octave symmetric wavelet in the normal frequency domain according to the amplitude spectrum form of each octave symmetric wavelet in the octave frequency domain comprises: To The frequency of the frequency domain of the octave frequency symmetry wavelet is transformed back to the frequency value, and the amplitude value corresponding to the uniformly sampled frequency value in the frequency domain is calculated by using the interpolation method, so as to obtain the amplitude spectrum form of the octave frequency symmetry wavelet in the ordinary frequency domain. The phase spectrum of the octave frequency symmetry wavelet is set to zero, and the frequency spectrum of the octave frequency symmetry wavelet in the ordinary frequency domain can be obtained from the amplitude spectrum and the phase spectrum.

7. The time-frequency analysis method of claim 1, wherein, The filtering of the signal to be analyzed by taking each octave symmetric wavelet in the normal frequency domain as a filter family and taking the frequency spectrum thereof as filter family frequency spectrum, to obtain a decomposition signal frequency spectrum and a time domain decomposition signal, comprises: transforming the signal to be analyzed from a time domain to a frequency domain by Fourier transform to obtain a signal to be analyzed spectrum; multiplying the signal to be analyzed spectrum and the filter family frequency spectrum to obtain a decomposition signal spectrum in the normal frequency domain; performing Fourier inverse transform on the decomposition signal spectrum to obtain a time domain expression of the decomposition signal.

8. The time-frequency analysis method of claim 1, wherein, The Hilbert transform on the time domain decomposition signal and the calculation of the instantaneous amplitude spectrum thereof comprise: solving an analytic signal corresponding to each decomposition signal by a Hilbert transform formula and obtaining the instantaneous amplitude corresponding thereto; the specific Hilbert transform formula is as follows: ， ； wherein denotes the decomposition signal, denotes the Hilbert transform of the decomposition signal, denotes the i-th decomposition signal; the formula for solving the instantaneous amplitude is as follows: ， ； wherein denotes the instantaneous amplitude data of the i-th decomposition signal, t is the time sample value after discrete sampling.

9. The time-frequency analysis method of claim 1, wherein, The generation of a time-frequency spectrum of the signal to be analyzed based on the instantaneous amplitude spectrum of the time domain decomposition signal comprises: All the decomposition signals correspond to the same time axis but different main frequencies. With the same time axis corresponding to all the decomposition signals as the horizontal coordinate, with the different main frequencies corresponding to different decomposition signals as the vertical coordinate, and with the instantaneous amplitude spectrum corresponding to different decomposition signals as the numerical value of the two-dimensional time-frequency spectrum, a time-frequency analysis spectrum of the signal to be analyzed based on the frequency-doubling symmetric wavelet is obtained.

Images