A signal frequency identification method based on active aliasing and dealiasing

Abstract

The present disclosure discloses a signal frequency identification method based on active aliasing and dealiasing, comprising: 1. sampling a signal to be identified, obtaining dual-channel sampling sample data, and intercepting two groups of data of the same length therefrom; 2. recording the relative delay time of the two groups of data and performing Fourier transform to obtain two groups of Fourier transform coefficient vectors, obtaining a spectrum based on the two groups of vectors, and obtaining an index position through peak search; 3. extracting corresponding complex numbers from the two groups of vectors according to the index position, obtaining intermediate variables based on the extracted complex numbers and calculating the phase angle value of the intermediate variables, and calculating an equivalent aliasing frequency and an equivalent sampling frequency according to the phase angle value and the relative delay time of the two groups of data; 4. changing the interception position of the two groups of data, repeating steps 2 and 3 to obtain a series of aliasing frequencies and equivalent sampling frequencies and calculating the remainder frequency; 5. solving a group of congruence equations constructed based on the aliasing frequency and the remainder frequency to obtain the frequency of the signal to be identified.

Description

Technical Field

[0001] The present disclosure belongs to the technical field of signal parameter identification and non-destructive testing, and particularly relates to a signal frequency identification method based on active aliasing and dealiasing. Background Art

[0002] In the field of signal processing, frequency estimation of undersampled signals is a common problem. Delay estimation uses the phase change during a certain time delay to estimate the signal's frequency. To accurately identify the signal's frequency, the delay interval must be less than the Nyquist interval. When the target frequency is high, an extremely short delay interval is required to accurately obtain the signal frequency. This increases the requirements for sampling time accuracy and significantly reduces the system's robustness to noise and sampling time fluctuations, rendering traditional delay estimation methods ineffective.

[0003] The remainder theorem uses the idea of ​​coprime to recover the signal frequency from undersampled signals. However, the remainder theorem requires at least two samplers with coprime sampling rates to perform uniform sampling. Since multiple sampling devices with different processing rates are required, this places higher requirements on the hardware. In addition, for passive sampling, such as direction of arrival estimation and blade tip timing vibration measurement, the implementation of different sampling rates depends on the number and layout of sensors. For example, to achieve two coprime sampling frequencies in blade tip timing, at least four sensors are required, which greatly increases the measurement cost. Therefore, the application of the remainder theorem in this type of passive sampling is hindered. Summary of the Invention

[0004] In response to the deficiencies in the prior art, the purpose of the present disclosure is to provide a signal frequency recognition method based on active aliasing and dealiasing. This method uses active aliasing delay estimation to obtain a series of equivalent sampling rates and residue frequencies from a severely undersampled signal, and then uses a robust closed remainder theorem algorithm to dealias the residue frequencies to obtain a frequency estimate of the signal to be identified.

[0005] To achieve the above objectives, the present disclosure provides the following technical solutions:

[0006] A signal frequency identification method based on active aliasing and dealiasing comprises the following steps:

[0007] S100: performing delayed coprime sampling on the signal to be identified to obtain dual-channel sampling sample data, and intercepting two groups of data of the same length from the dual-channel sampling sample data;

[0008] S200: Record the relative delay time of the two sets of data, perform Fourier transform on the two sets of data respectively, obtain two sets of Fourier transform coefficient vectors, take the absolute value of the Fourier transform coefficient vector to obtain the spectrum, and obtain the index position corresponding to the frequency peak through peak search;

[0009] S300: extracting corresponding complex numbers from dual-path Fourier transform coefficients according to index positions, multiplying the first path complex number by the conjugate of the second path complex number to obtain an intermediate variable, calculating a phase angle value of the intermediate variable, and calculating an equivalent aliasing frequency and an equivalent sampling frequency according to the phase angle value and a relative delay time of two sets of data;

[0010] S400: changing the interception positions of the two sets of data so that the relative delay of the two sets of data changes, repeating steps S200 and S300 to obtain a series of aliasing frequencies and equivalent sampling frequencies, and calculating the remainder frequency based on the series of aliasing frequencies;

[0011] S500: Solve a series of congruence equations with the equivalent sampling frequency as the modulus and the remainder frequency as the remainder by using the remainder theorem to obtain the frequency of the signal to be identified.

[0012] Preferably, in step S100, two samplers with the same sampling rate are used to perform delayed coprime sampling on the signal to be identified, and the sampling rate is recorded as f s , the sampling period is recorded as T0, f s =1 / T0, the delay time of the two sampling channels is recorded as τ.

[0013] Preferably, in step S200, the two groups of Fourier transform coefficient vectors are respectively expressed as:

[0014]

[0015] Wherein, d1(n) represents the signal obtained by sampling the first channel; j represents the imaginary number symbol, N represents the length of the signal selected for analysis, which is determined manually. In this example, N=700; n represents the number of iterations, from 1 to N; k represents an integer from 1 to N; D1(k) represents the pair [d1(1), d1(2), ..., d1(N)] T The kth data after discrete Fourier transform;

[0016]

[0017] Where d2(n) represents the signal obtained by sampling the second channel; D2(k) represents the pair [d2(1), d2(2), ..., d2(N)] T The kth data after discrete Fourier transform.

[0018] Preferably, in step S300, the equivalent aliasing frequency is calculated by the following formula:

[0019]

[0020] in,(·) *represents the conjugate operation, angle[·] represents the operation of finding the complex phase angle, φ represents the phase angle value, D1 and D2 represent the Fourier coefficients obtained by Fourier transforming the two intercepted sets of data, D1(c) and D2(c) represent the frequencies of the c-th values ​​in the two Fourier coefficient vectors, π represents pi, and τ1 represents the relative delay time between the two intercepted sets of data.

[0021] Preferably, the complex phase angle operation is performed by the following formula:

[0022]

[0023] Here, α+bi represents a complex number with a real part and b imaginary part.

[0024] Preferably, in step S300, the equivalent sampling frequency is expressed as:

[0025]

[0026] in, Indicates the delay is τ v The equivalent sampling frequency in this case is v=1, 2, 3,….

[0027] Preferably, in step S400, when the signal to be identified is a real-valued signal, each equivalent aliasing frequency can obtain two possible remainder frequencies, namely or Wherein, v represents an integer, that is, v = 1, 2, 3, ..., represents the vth equivalent sampling frequency, represents the vth equivalent aliasing frequency, Represents the vth remainder frequency; when the signal to be identified is a complex signal, the remainder frequency is equal to the equivalent aliasing frequency.

[0028] Preferably, in step S500, the congruence equations are expressed as follows:

[0029]

[0030]

[0031]

[0032] Among them, n1, n2, n3 represent the folding coefficient, f n represents the frequency of the signal to be solved.

[0033] Compared with the prior art, the present invention has the following beneficial effects:

[0034] This method significantly reduces the average sampling rate of samples, enabling frequency estimation of the signal to be identified using only two low-speed sampling channels. Unlike traditional multi-rate sampling, which uses circuit hardware to uniformly sample multiple rates, this method uses two sampling channels with the same sampling rate to sample the signal. Based on the acquired data, the delay between the two channels is used to generate multiple equivalent sampling frequencies. This algorithm effectively achieves the effect of multi-rate sampling, significantly simplifying the hardware and reducing the cost of frequency identification. BRIEF DESCRIPTION OF THE DRAWINGS

[0035] Figure 1 A flowchart of a signal frequency identification method based on active aliasing and dealiasing provided in one embodiment of the present disclosure;

[0036] Figure 2 This is the principle diagram of dual-path delayed coprime sampling;

[0037] Figure 3 This is the time domain diagram of the simulated signal acquired through dual-path delayed mutual prime sampling;

[0038] Figure 4 is the result of the first equivalent aliasing frequency calculation;

[0039] Figure 5 is the result of the second equivalent aliasing frequency calculation;

[0040] Figure 6 Calculate the result for the third equivalent aliasing frequency:

[0041] FIG7( a ) is a schematic diagram of a first interception of dual-channel sampling sample data;

[0042] FIG7( b ) is a second interception diagram for dual-channel sampling sample data;

[0043] FIG7( c ) is a schematic diagram of a third interception of dual-channel sampling sample data. DETAILED DESCRIPTION

[0044] The following will refer to the attached Figures 1 to 7(c) Specific embodiments of the present disclosure are described in detail. Although specific embodiments of the present disclosure are shown in the accompanying drawings, it should be understood that the present disclosure can be implemented in various forms and should not be limited by the embodiments described herein. Instead, these embodiments are provided to enable a more thorough understanding of the present disclosure and to fully convey the scope of the present disclosure to those skilled in the art.

[0045] It should be noted that certain words are used in the specification and claims to refer to specific components. Those skilled in the art should understand that technicians may use different nouns to refer to the same component. This specification and claims do not use the difference in nouns as a way to distinguish components, but use the difference in the functions of the components as the criterion for distinction. As mentioned throughout the specification and claims, "including" or "comprising" is an open term, so it should be interpreted as "including but not limited to". The subsequent description of the specification is a preferred embodiment of the present disclosure, but the description is based on the general principles of the specification and is not used to limit the scope of the present disclosure. The scope of protection of the present disclosure shall be as defined by the attached claims.

[0046] To facilitate understanding of the embodiments of the present disclosure, further explanation will be given below using specific embodiments as examples in conjunction with the accompanying drawings, and the accompanying drawings do not constitute a limitation on the embodiments of the present disclosure.

[0047] like Figure 1 As shown, the present disclosure proposes a method for frequency identification of undersampled signals based on active aliasing and dealiasing, comprising the following steps:

[0048] S100: performing delayed coprime sampling on the signal to be identified to obtain dual-channel sampling sample data, and intercepting two groups of data of the same length from the dual-channel sampling sample data;

[0049] In this step, first generate the simulation signal using the following formula:

[0050]

[0051] Where d(t) represents the generated simulation signal, f, A and Respectively represent the frequency, amplitude and phase of the simulation signal d(t), N(t) represents the noise term, for example, f = 875 Hz, A = 1 mm, The signal-to-noise ratio in the simulation is 20 dB and the simulation duration is 4 seconds.

[0052] Then, using two sampling rates f s = 180Hz sampler performs delayed coprime sampling on the simulation signal d(t), with a sampling period of T0 = 1 / 180s. The sampling time of the second sampler is delayed by τ1 = 1 / 450s compared to the sampling time of the first sampler. At this time, the ratio of the time delay of the sampling time of the two samplers to the sampling period is: τ1 / T0 = 2 / 5, 2 and 5 are two coprime numbers. Each sampler samples 720 sample points, for a total of 1440 samples. The collected dual-channel sampling sample data is as follows Figure 3 shown.

[0053] The two discrete signal sample vectors collected by the two samplers are:

[0054] d1=[d1(1), d1(2), d1(3),…, d1(720)] T =[0.5371, -0.2738, -0.8966, ..., 1.1168] T

[0055] d2=[d2(1), d2(2), d2(3),…, d2(720)] T =[0.5202, -0.1472, -0.4603, ..., 1.4022] T

[0056] S200: Record the relative delay time of the two sets of data, perform Fourier transform on the two sets of data respectively, obtain two Fourier transform coefficient vectors, take the absolute value of the Fourier transform coefficient vector to obtain the spectrum, and obtain the index position corresponding to the frequency peak through peak search;

[0057] In this step, as shown in Figure 7(a), two sets of data of the same length are intercepted from the dual-channel sample data. In Figure 7(a), the two data interception windows used are aligned, and the two sets of sampled data are represented by gray squares, with each square representing one data. First, the first N data of the first channel and the first N data of the second channel are selected for analysis. It can be seen that the relative delay between the two sets of data is τ1 = 1 / 450s. The two intercepted sets of data are Fourier transformed according to the following formula to obtain two sets of Fourier transform coefficient vectors D1 and D2:

[0058]

[0059] Wherein, d1(n) represents the signal obtained by sampling the first channel; j represents the imaginary number symbol, N represents the length of the signal selected for analysis, which is determined manually. In this example, N=700; n represents the number of iterations, from 1 to N; k represents an integer from 1 to N; D1(k) represents the pair [d1(1), d1(2), ..., d1(N)] T The kth data after discrete Fourier transform.

[0060] D1 is the Fourier transform coefficient vector of the first channel of data, D1=[D1(1), D1(2), …, D2(N)] T .

[0061]

[0062] Where d2(n) represents the signal obtained by sampling the second channel; D2(k) represents the pair [d2(1), d2(2), ..., d2(N)] TThe kth data after discrete Fourier transform.

[0063] D2 is the Fourier transform coefficient vector of the second data, D2 = [D2(1), D2(2), ..., D2(N)] T .

[0064] By finding the absolute value of D1 or D2, an amplitude spectrum can be obtained. The position index of the aliasing frequency can be located by peak search, which is recorded as c. For example, the absolute value of D1 is found and the index of the peak position is found. The index position c = 101, such as Figure 4 Shown in the upper part.

[0065] S300: extracting corresponding complex numbers from dual-path Fourier transform coefficients according to index positions, multiplying the first path complex number by the conjugate of the second path complex number to obtain an intermediate variable, calculating a phase angle value of the intermediate variable, and calculating an equivalent aliasing frequency and an equivalent sampling frequency according to the phase angle value and a relative delay time of two sets of data;

[0066] In this step, two corresponding Fourier coefficients D1(c) and D2(c) are extracted from the Fourier coefficient vectors of the first and second paths of data according to the index position c. Figure 4 In the example, when the index position c=101, D1(101)=0.1481+0.4729jj is an imaginary unit. D2(101)=-0.0224+0.4962j

[0067] Multiply D1(c) by the conjugate of D2(c) to get the intermediate variable, that is,

[0068]

[0069]

[0070] angel[0.2312-0.0841j]=-0.349rad, which is equal to Figure 4 The phase value extracted from .

[0071] Find the phase angle φ of the intermediate variable, that is And divide by 2πτ1 to get the first equivalent aliasing frequency The calculation formula is as follows:

[0072]

[0073] in,(·) * represents the conjugate operation, and angle[·] represents the operation of finding the complex phase angle, which is defined as follows:

[0074]

[0075] Here, a+bi represents a complex number with a real part and b imaginary part.

[0076] For example, when c = 101, τ1 = 1 / 450s, the first equivalent aliasing frequency is

[0077]

[0078] S400: changing the interception positions of the two sets of data so that the relative delay of the two sets of data changes, repeating steps S200 and S300 to obtain a series of aliasing frequencies and equivalent sampling frequencies, and calculating the remainder frequency based on the series of aliasing frequencies;

[0079] In this step, as shown in Figure 7(b), the window interception position of the dual-channel sample data is changed, that is, the interception window of the first channel data is shifted so that the first interception window position lags behind the second interception window position by one data bit, and the 2nd to N+1th data of the first channel data are selected for Fourier transform to obtain the new Fourier transform coefficient vector D1 of the first channel data:

[0080]

[0081] The window cut-off position of the second data channel remains unchanged, and Fourier transform is still performed on the 1st to Nth data channels. The time delay of the first sampling point of the two data channels is T0-τ1=1 / 300s, recorded as τ2. According to the index position 101, the corresponding complex numbers are extracted from the 1st and 2nd Fourier transform coefficient vectors, and the complex number of the 1st channel is multiplied by the conjugate of the complex number of the 2nd channel to obtain the extracted phase angle φ=-0.5234rad, as shown in Figure 5 shown.

[0082] The second equivalent aliasing frequency is calculated according to the following formula:

[0083]

[0084] As shown in Figure 7(c), the window interception position of the dual-channel sample data is changed again, that is, the first interception window position is one data bit ahead of the second interception window position, that is, the 1st to Nth data of the first channel data and the 2nd to N+1th data of the same channel data are selected for Fourier transform respectively, and two sets of Fourier transform coefficient vectors D are obtained. 1a and D 1b , according to index position c = 101 from D 1a and D 1b Extract the corresponding complex number from the Fourier transform coefficient vector and convert D 1a (101) multiplied by D 1b(101) and the extracted phase angle is φ = -0.8728 rad, as Figure 6 shown.

[0085]

[0086]

[0087] In this case, the time delay between the first sampling points of the two channels of data is T0, recorded as τ3. The third equivalent aliasing frequency is calculated according to the following formula:

[0088]

[0089] By continuously changing the window interception position of the dual-channel sampling sample data, different degrees of delay time can be achieved, namely:

[0090] τ=[τ1, τ2, τ3, τ4, τ5, τ6,…] T =[τ1,T0-τ1,T0,T0+τ1,2T0-τ1,2T0,…] T .

[0091] By taking the reciprocal of the delay time, we can get the equivalent sampling frequency under the corresponding delay condition, that is,

[0092]

[0093] in, Indicates the delay is τ v The equivalent sampling frequency in this case is v=1, 2, 3,….

[0094] Furthermore, if the simulated signal is a real-valued signal, there are two possibilities for each remainder frequency.

[0095] or

[0096] or

[0097] or

[0098] Wherein, v represents an integer, v = 1, 2, 3, represents the vth remainder frequency.

[0099] For example, The corresponding remainder frequencies obtained have the following 8 possibilities.

[0100]

[0101] S500: Solve a series of congruence equations with the equivalent sampling frequency as the modulus and the remainder frequency as the remainder by using the remainder theorem to obtain the frequency of the signal to be identified.

[0102] In this step, first take the first three equivalent sampling frequencies and the corresponding remainder frequencies and combine them in pairs. Use the equivalent sampling frequency as the modulus and the remainder frequency as the remainder to establish a congruence equation system:

[0103]

[0104]

[0105]

[0106] Among them, n1, n2, n3 represent the folding coefficient, f n represents the frequency of the signal to be solved.

[0107] Substitute all possible cases of the remainder frequency into the above equations in turn. Each possible case has a corresponding frequency solution. According to the consistency of the approximate frequency range a priori and the obtained frequency value, the correct case is determined, and the corresponding frequency solutions are averaged to obtain the final frequency estimate.

[0108] The final frequency estimation solution process is as follows:

[0109] Substituting the value of the remainder frequency in the first case and the equivalent sampling frequency into the equation, we have:

[0110]

[0111]

[0112]

[0113] Use the closed robust remainder theorem algorithm to solve each system of equations. Take the first system of equations as an example:

[0114] M=gcd(450,300)=150, where gcd represents the greatest common divisor.

[0115]

[0116] in is the modular inverse of Γ1 with respect to Γ2,

[0117] Where γ1=Γ1·Γ2 / Γ1=2, b 2，1 It is the modular inverse of γ1 / Γ2 with respect to r2.

[0118]

[0119]

[0120] Solving all equations, we can get the following results:

[0121]

[0122] Theoretically, when the residual frequencies are correctly chosen, the frequencies recovered by different systems of equations are very similar. Based on the approximate range of frequencies and the consistency of the obtained results, we know that Case 7 corresponds to the correct case.

[0123] So there is

[0124] In addition, it should be noted that the above situation corresponds to the analysis process of a real-valued signal. If the simulation signal is a complex-valued signal, for example, the simulation signal in the embodiment is modified to:

[0125]

[0126] Then the correspondence between the remainder frequency and the equivalent aliasing frequency is as follows:

[0127]

[0128] At this time, the calculation result of the equivalent aliasing frequency is: So the corresponding remainder frequency Solving this case yields the same result as above.

[0129] Finally, the obtained signal frequency estimates are averaged to obtain the final frequency estimate:

[0130]

[0131] The frequency estimation value is very close to the actual frequency value of 875 Hz set by the simulation signal, with an error of only 0.0033 Hz, which illustrates the effectiveness of the method described in the present disclosure.

[0132] The present disclosure realizes a series of equivalent sampling rates by changing the phase difference, allowing the identification of frequencies from signals with extremely low sampling rates. The method is simple and feasible, has low computational complexity, and can be used for extracting the frequencies of undersampled signals.

[0133] The above description is only a preferred example of the present invention and is not intended to limit the present invention. Any modifications, equivalent substitutions and improvements made within the spirit and principles of the present invention should be included in the scope of protection of the present invention.

Claims

1. A signal frequency identification method based on active aliasing and dealiasing, comprising the following steps: S100: performing delayed coprime sampling on the signal to be identified to obtain dual-channel sampling sample data, and intercepting two groups of data of the same length from the dual-channel sampling sample data; S200: Record the relative delay time of the two sets of data, perform Fourier transform on the two sets of data respectively, obtain two sets of Fourier transform coefficient vectors, take the absolute value of the Fourier transform coefficient vector to obtain the spectrum, and obtain the index position corresponding to the frequency peak through peak search; S300: extracting corresponding complex numbers from dual-path Fourier transform coefficients according to index positions, multiplying the first path complex number by the conjugate of the second path complex number to obtain an intermediate variable, calculating a phase angle value of the intermediate variable, and calculating an equivalent aliasing frequency and an equivalent sampling frequency according to the phase angle value and a relative delay time of two sets of data; S400: changing the interception positions of the two sets of data so that the relative delay of the two sets of data changes, repeating steps S200 and S300 to obtain a series of aliasing frequencies and equivalent sampling frequencies, and calculating the remainder frequency based on the series of aliasing frequencies; S500: Solving a series of congruence equations with the equivalent sampling frequency as the modulus and the remainder frequency as the remainder using the remainder theorem to obtain the frequency of the signal to be identified; in, In step S400, when the signal to be identified is a real-valued signal, two remainder frequencies can be obtained for each equivalent aliasing frequency. or ,in, represents an integer, , represents the vth equivalent sampling frequency, represents the vth equivalent aliasing frequency, Represents the vth remainder frequency; when the signal to be identified is a complex signal, the remainder frequency is equal to the equivalent aliasing frequency.

2. The method according to claim 1, wherein In step S100, two samplers with the same sampling rate are used to perform delayed coprime sampling on the signal to be identified. The sampling rate is recorded as , the sampling period is recorded as , , the delay time of the two sampling channels is recorded as .

3. The method according to claim 1, wherein In step S200, the two groups of Fourier transform coefficient vectors are respectively expressed as: ， in, Indicates the signal obtained by sampling the first channel; Represents the imaginary number symbol, ; Indicates the length of the signal selected for analysis, which is determined manually. In this example, N=700; Indicates the number of iterations, from 1 to N; Indicates 1 to integer; Express After discrete Fourier transform individual data; ， Wherein, d2(n) represents the signal obtained by sampling the second channel; Express After discrete Fourier transform data.

4. The method according to claim 1, wherein In step S300, the equivalent aliasing frequency is calculated by the following formula: ， in, represents the conjugate operation, represents the operation of finding the complex phase angle, represents the phase angle value, and Respectively represent the Fourier coefficients obtained after Fourier transform of the two sets of intercepted data, and Respectively represent the frequency as the c-th value in the two Fourier coefficient vectors, represents pi, Indicates the relative delay time between two sets of intercepted data.

5. The method according to claim 4, wherein The complex phase angle operation is performed by the following formula: ， Here, a+bi represents a complex number with a real part and b imaginary part.

6. The method according to claim 1, wherein In step S300, the equivalent sampling frequency is expressed as: in, Indicates the delay is The equivalent sampling frequency in this case is .

7. The method according to claim 1, wherein In step S500, the congruence equations are expressed as follows: ， ， ， in, represents the folding factor, represents the frequency of the signal to be solved.

Images