Short-time adaptive time-varying filtering method based on Gaussian Newton algorithm

Abstract

The invention provides a short-time adaptive time-varying filtering method based on a Gauss-Newton algorithm. The method mainly solves the problems that an adaptive filter in the prior art can not process unknown interference due to the fact that reference signals need to be introduced in advance, and useful signals cannot be correctly recognized due to the fact that the interference and the useful signals are located in the same frequency band. According to the scheme adopted by the invention, the instantaneous angular frequency of interference can be automatically identified and filtered, so that the interference is removed, reference signals are not needed, meanwhile, the accuracy of estimating the instantaneous angular frequency of the interference is improved, the filtering effect is enhanced, and when the interference and the signals are in the same frequency band, the filtering efficiency is improved. According to the method, the instantaneous angular frequency of interference can be automatically identified for specified filtering, and useful signals cannot be lost.

Description

Technical Field

[0001] The invention belongs to the technical field of signal and information processing, and in particular relates to an adaptive time-varying filtering method. Background Art

[0002] Signal processing is a very important field in modern science and engineering. Whether in communications, radar, biomedical imaging or other fields, signal processing is indispensable. However, signals in actual environments are often subject to various interferences, which can seriously affect the quality and reliability of signals. Therefore, an efficient filter is needed to remove these interferences to ensure the accuracy and stability of the signal.

[0003] In the prior art, filters are classified according to their working principles and methods, and common ones include convolution filters and adaptive filters; convolution filters are simple and easy to understand, but they show high sensitivity when facing noise and sudden interference; adaptive filters can adaptively adjust filter parameters and are suitable for nonlinear systems, but they need to introduce reference signals in advance and cannot handle unknown interference; filters are classified according to their frequency response characteristics, and common ones include low-pass filters, high-pass filters, band-pass filters and band-stop filters, which can effectively suppress interference in a specified frequency band, but the range of the interference frequency needs to be known in advance. If the interference and the useful signal are in the same frequency band, while suppressing the interference, the useful signal will also be lost, and the useful signal cannot be correctly identified. In addition, traditional filters have certain limitations when dealing with time-varying interference and cannot effectively adaptively identify and remove interference. Therefore, a filter that can overcome the above shortcomings is now needed. Summary of the invention

[0004] In order to overcome the shortcomings of the prior art, the present invention provides a short-time adaptive time-varying filtering method based on the Gauss-Newton algorithm, which can automatically identify the instantaneous angular frequency of interference and perform filtering to remove the interference.

[0005] A short-time adaptive time-varying filtering method based on Gauss-Newton algorithm comprises the following steps:

[0006] S1: According to the input interference s(n) and sampling rate F s , set the instantaneous angular frequency parameter vector θ U , and initialize;

[0007] Set the instantaneous angular frequency parameter vector θ U The method is a fixed value method or a time-frequency analysis method; the time-frequency analysis method includes a short-time Fourier transform time-frequency analysis method and a Wigner-Wiley time-frequency analysis method;

[0008] S2: Divide the input interference s(n) into several short-term interferences, with a total number of L segments, ls is the current segment number;

[0009] S3: Traverse all short-term interferences and determine the current segment number l s Is it equal to the total number of segments L?

[0010] If l s If it is equal to the total number of segments L, the final instantaneous angular frequency estimation of all short-term interferences has been completed, the loop ends, and the process jumps to step S5;

[0011] If l s If it is not equal to the total number of segments L, execute l s =Count ls +1, Count ls = l s , and jump to step S4 to extract the first s The final instantaneous angular frequency of the short-term interference;

[0012] S4: Through the first s The short-term interference and the instantaneous angular frequency parameter vector θ U , using the basis function space expansion and Gauss-Newton algorithm to extract the lth s The final instantaneous angular frequency ω(n) of the short-term interference is obtained, and then the process goes to step S3;

[0013] S5: Calculate the filtering result y(n) according to the final instantaneous angular frequency ω(n) of the interference and the time-varying filter transfer function.

[0014] Furthermore, the steps of setting the instantaneous angular frequency parameter vector by the Wigner-Wiley time-frequency analysis method are as follows:

[0015] S1-1.1: According to the input interference s(n) and sampling rate F s , calculate the Wigner-Wiley distribution WVD(n1,ω s1 );

[0016] Wigner-Wiley distribution WVD(n1,ω s1 )The specific formula is:

[0017]

[0018] The sampling point number n1 ranges from [1,2,...,N w ],N w Indicates the length of short-term interference, angular frequency ω s1 The value is z(n1) represents the analytical signal, and the specific formula is:

[0019]

[0020] Wherein, j represents the imaginary number symbol;

[0021] S1-1.2: Calculate the Wigner-Wiley distribution WVD(n1,ω s1 ) time vector t w1 The angular frequency vector ω with the maximum energy corresponding to each moment w1 The specific formula is

[0022]

[0023] S1-1.3: Using the time vector t w1 The angular frequency vector ω with the maximum energy corresponding to each moment w1 Calculate the estimated instantaneous angular frequency parameter vector θ U , the specific formula is:

[0024]

[0025] in, Represents the time vector t w1 The i-th power of each element in, Q represents the maximum order of the basis function space, U w1 is the time vector t w1 The corresponding basis function space is represents the Moore–Penrose pseudoinverse.

[0026] Furthermore, the steps of the short-time Fourier transform time-frequency analysis method are:

[0027] S1-2.1: Input disturbance s(n) and sampling rate F s , calculate the short-time Fourier transform S(n2,ω s2 );

[0028]

[0029] The sampling point number Angular frequency represents the window function; N w2 Indicates the length of the window function, the window function w(m w2 ) is expressed as:

[0030]

[0031] S1-2.2: Calculate the short-time Fourier transform S(n2,ω s2 ) time vector t w2 The angular frequency vector ω with the maximum energy corresponding to each moment w2 , the specific formula is:

[0032]

[0033] S1-2.3: Using the time vector t w2 The angular frequency vector ω with the maximum energy corresponding to each moment w2 Calculate the estimated instantaneous angular frequency parameter vector θ U , the specific formula is:

[0034]

[0035] in, Represents the time vector t w2 The i-th power of each element in, Q represents the maximum order of the basis function space, U w2 is the time vector t w2 The corresponding basis function space is represents the Moore–Penrose pseudoinverse.

[0036] Furthermore, in step S2, the specific steps of segmenting the input interference s(n) into a number of short-term interferences are:

[0037] The input interference s(n) is divided into several short-term interferences, the total number of segments is L, and the starting sequence number of each short-term interference is N. s (l s ), the end sequence number N of each short-term interference e (l s ), l s is the current segment number, Count ls for l s counter;

[0038] Among them, the total number of segments is L and the starting sequence number is N s (l s ), end sequence number N e (l s ), segment number l s Initial value, Count ls The initial values ​​are:

[0039]

[0040] N s (l s )=(l s -1)N m +1;

[0041]

[0042] l s =0;

[0043] Count ls =0;

[0044] Where N represents the total length of the input interference s(n), N m Indicates the length of movement when extracting each short-term interference; Indicates rounding down.

[0045] Furthermore, in step S4, by the first s The short-term interference and the instantaneous angular frequency parameter vector θ U , using the basis function space expansion and Gauss-Newton algorithm to extract the lth s The specific steps to calculate the final instantaneous angular frequency ω(n) of a short-term interference are:

[0046] S4.1: Extract the first s The short-term interference is used to establish the basis function space U, the pole contraction coefficient ρ(l) and the regularization matrix Γ (l) , and initialize the number of iterations of the Newton-Gauss algorithm, l, the initial value of l is l = 1; Count l is a counter; Count l =1

[0047] s s (n3) = s(n3 + l s N w -N w ),n3=1,2,...,N e (l s )-N s (l s ), s s (n3) is the first s Short-term interference;

[0048]

[0049] Γ (l) =10 -Q (10-1)I Q

[0050]

[0051] Among them I Q represents the Q-order identity matrix;

[0052] S4.2: Transform the instantaneous angular frequency parameter vector θ U Multiply the basis function space U to obtain the estimated instantaneous angular frequency ω of the short-term interference s (n4):

[0053] ω s (n4) = θ U U(n4),n4=1,2,...,N w

[0054] S4.3: Estimated instantaneous angular frequency ω s (n4) performing cosine and multiplication operations to obtain time-varying filter coefficients a(n4);

[0055] a(n4)=-2cos[ω s (n4)],n4=1,2,...,N w

[0056] S4.4: Calculate the prediction error ε(n4) and the instantaneous angular frequency parameter vector gradient ψ(n4) through the time-varying filter transfer function;

[0057] In step S4.4, the prediction error ε(n4) and the instantaneous angular frequency parameter vector gradient ψ(n4) are calculated by the time-varying filter transfer function. The specific formula is:

[0058] ε(n4)=s s (n4)+a(n4)s s (n4-1)+s s (n4-2)-ρ(n4)a(n4)ε(n4-1)-ρ 2 (n4)ε(n4-2),n4=1,2,...,N w

[0059]

[0060] in is the gradient of the time-varying filter coefficient a(n4) relative to the prediction error ε(n4), and the specific formula is:

[0061] s F (n4) = s s (n4)-ρ 2 (n4)s F (n4-2)-ρ(n4)s F (n4-1)a(n4),n4=1,2,...,N w

[0062] ε F (n4)=ε(n4)-ρ 2 (n4)ε F (n4-2)-ρ(n4)ε F (n4-1)a(n4),n4=1,2,...,N w

[0063]

[0064] Among them, s F (n4) and ε F (n4) is the intermediate variable;

[0065] S4.5: Determine the cost value β of the first cycle (l) and the cost value of the l-1th cycle β (l-1) Is the absolute error less than the threshold δ:

[0066] If the cost of the first cycle is β (l) and the cost value of the l-1th cycle β (l-1) If the absolute error is less than the threshold δ, then the current instantaneous angular frequency ω s (n4) is the final instantaneous angular frequency, go to step S3;

[0067] If the cost of the first cycle is β (l) and the cost value of the l-1th cycle β (l-1) If the absolute error of is greater than or equal to the threshold δ, execute step S4.6;

[0068] Cost value β (l) The specific formula is:

[0069]

[0070] The threshold value δ is 1×10 -6 ;

[0071] The formula for extracting the final instantaneous angular frequency ω(n) is:

[0072]

[0073]

[0074] ω(n)=ω s (nl s *N m +N m ),n=N ss ,N ss +1,...,N se

[0075] Where N ss Indicates extracting the first s The starting number of the final instantaneous angular frequency of the short-term interference, N se Indicates extracting the first s The ending sequence number of the final instantaneous angular frequency of the short-term interference segment;

[0076] S4.6: Update the instantaneous angular frequency parameter vector θ using the Gauss-Newton algorithm U As the instantaneous angular frequency parameter vector of the l+1th cycle, execute step S4.2;

[0077] Update the instantaneous angular frequency parameter vector θ UThe specific formula is:

[0078]

[0079] l=Count l +1

[0080] Count l = l

[0081] Where R (l) is the pseudo-Hessian matrix of the cost function of the lth cycle, is the gradient of the cost function at the lth cycle.

[0082] Furthermore, in step S5, the specific steps of calculating the filtering result y(n) according to the final instantaneous angular frequency ω(n) of the interference and the time-varying filter transfer function are:

[0083] y(n)=s(n)-2cos[ω(n)]s(n-1)+s(n-2)+2ρ e cos[ω(n)]y(n-1)-ρ e 2 y(n-2),n=1,2,...,N

[0084] where ρ e is the filter coefficient; filter coefficient ρ e The value is 0.995.

[0085] The beneficial effects of the present invention are:

[0086] 1. The present invention uses the Gauss-Newton algorithm to iteratively update the instantaneous angular frequency parameter vector, thereby obtaining accurate instantaneous angular frequency and filtering, without the need for a reference signal, and is completely driven by data;

[0087] 2. Compared with the prior art, the short-time adaptive time-varying filtering method based on the fixed initial value Gauss-Newton algorithm improves the accuracy of estimating the instantaneous angular frequency of the interference, thereby enhancing the filtering effect, and when the interference and the signal are in the same frequency band, the present invention can automatically identify the instantaneous angular frequency of the interference for designated filtering without causing loss to the useful signal; when the instantaneous angular frequency parameter vector is set to estimate the instantaneous angular frequency of the first short-time interference using the short-time Fourier transform time-frequency analysis method or the Wigner-Wiley time-frequency analysis method, the number of iterations is further reduced on the basis of the fixed value method, thereby reducing the amount of calculation;

[0088] 3. The present invention adopts a time-varying filter to reduce the influence of noise on the instantaneous angular frequency estimation, and has a strong anti-noise performance;

[0089] 4. For multi-component signals, the present invention obtains the instantaneous angular frequency of the highest angular frequency component thereof, and thus can further process multi-component interference by filtering the filtered signal again. BRIEF DESCRIPTION OF THE DRAWINGS

[0090] Figure 1 It is the flow chart of short-time adaptive time-varying filter based on Gauss-Newton algorithm;

[0091] Figure 2 Flowchart for setting the instantaneous angular frequency parameter vector for the short-time Fourier transform time-frequency analysis method;

[0092] Figure 3 Flowchart for setting the instantaneous angular frequency parameter vector for the Wigner-Wiley time-frequency analysis method;

[0093] Figure 4 The short-time Fourier transform time-frequency distribution diagram of the sine frequency modulation interference simulated in Example 1 of the present invention;

[0094] Figure 5 A comparison diagram of a real instantaneous frequency curve and a filtered instantaneous frequency curve of a simulated sinusoidal frequency modulation interference according to Embodiment 1 of the present invention;

[0095] Figure 6 This is a comparison diagram of the simulated sinusoidal frequency modulation interference and filtering results of Example 1 of the present invention;

[0096] Figure 7 A Wigner-Wiley time-frequency distribution diagram of the simulated secondary frequency modulation interference of Embodiment 2 of the present invention;

[0097] Figure 8 A comparison diagram of a real instantaneous frequency curve and a filtered instantaneous frequency curve of a simulated secondary frequency modulation interference in Example 2 of the present invention;

[0098] Fig. 9 This is a comparison diagram of simulated secondary frequency modulation interference and filtering results in Example 2 of the present invention. DETAILED DESCRIPTION

[0099] The present invention is further described below in conjunction with the accompanying drawings and embodiments.

[0100] A short-time adaptive time-varying filtering method based on Gauss-Newton algorithm comprises the following steps:

[0101] S1: Input interference s(n) and sampling rate F s , set the instantaneous angular frequency parameter vector θ U , and initialize;

[0102] Instantaneous angular frequency parameter vector θ U for:

[0103] θU =[ω1,ω2,.ω i ..,ω Q ]

[0104] where ω i is the i-th fitting parameter for estimating the instantaneous angular frequency, and Q represents the maximum order of the basis function space;

[0105] θ U The initial value of is a fixed value or obtained by time-frequency analysis method;

[0106] θ U When the initial value of is fixed, Where 0 Q-1 is a 1×(Q-1) vector;

[0107] In order to take into account the computational complexity of the present invention and the fitting effect of the final instantaneous angular frequency, the preferred Q value is 4, and the instantaneous angular frequency parameter vector is initialized

[0108] Methods of time-frequency analysis include short-time Fourier transform and Wigner-Wiley distribution method;

[0109] Using the Wigner-Wiley time-frequency analysis method, the instantaneous angular frequency parameter vector θ is set U The initial value of is as follows:

[0110] S1-1.1: According to the input interference s(n) and sampling rate F s , calculate the Wigner-Wiley distribution WVD(n1,ω s1 );

[0111] Wigner-Wiley distribution WVD(n1,ω s1 )The specific formula is:

[0112]

[0113] The sampling point number n1 ranges from [1,2,...,N w ],N w Indicates the length of short-term interference, the value is 256, the angular frequency ω s1 The value is z(n1) represents the analytical signal, and the specific formula is:

[0114]

[0115] Wherein, j represents the imaginary number symbol;

[0116] S1-1.2: Calculate the Wigner-Wiley distribution WVD(n1,ω s1 ) time vector t w1The angular frequency vector ω with the maximum energy corresponding to each moment w1 The specific formula is

[0117]

[0118] S1-1.3: Using the time vector t w1 The angular frequency vector ω with the maximum energy corresponding to each moment w1 Calculate the estimated instantaneous angular frequency parameter vector θ U , the specific formula is:

[0119]

[0120] in, Represents the time vector t w1 The i-th power of each element in, Q represents the maximum order of the basis function space, U w1 is the time vector t w1 The corresponding basis function space is represents Moore-Penrose pseudoinverse;

[0121] The instantaneous angular frequency parameter vector θ is set using the short-time Fourier transform time-frequency analysis method. U The initial value of is as follows:

[0122] S1-2.1: Input disturbance s(n) and sampling rate F s , calculate the short-time Fourier transform S(n2,ω s2 );

[0123]

[0124] The sampling point number Angular frequency represents the window function; N w2 Indicates the length of the window function, the value is 32, the window function w(m w2 ) is expressed as:

[0125]

[0126] S1-2.2: Calculate the short-time Fourier transform S(n2,ω s2 ) time vector t w2 The angular frequency vector ω with the maximum energy corresponding to each moment w2 , the specific formula is:

[0127]

[0128] S1-2.3: Using the time vector t w2 The angular frequency vector ω with the maximum energy corresponding to each momentw2 Calculate the estimated instantaneous angular frequency parameter vector θ U , the specific formula is:

[0129]

[0130] in, Represents the time vector t w2 The i-th power of each element in, Q represents the maximum order of the basis function space, U w2 is the time vector t w2 The corresponding basis function space is represents Moore-Penrose pseudoinverse;

[0131] S2: Divide the input interference s(n) into several short-term interferences, the total number of segments is L, and the starting sequence number of each short-term interference is N s (l s ), the end sequence number N of each short-term interference e (l s ), l s is the current segment number, Count ls for l s counter;

[0132] Among them, the total number of segments is L and the starting sequence number is N s (l s ), end sequence number N e (l s ), segment number l s Initial value, Count ls The initial values ​​are:

[0133]

[0134] N s (l s )=(l s -1)N m +1;

[0135]

[0136] l s =0;

[0137] Count ls =0;

[0138] Where N represents the total length of the input interference s(n), m Indicates the length of movement when extracting each short-term interference; Indicates rounding down;

[0139] The present invention extracts the moving length value of 32 when each short-term interference is extracted;

[0140] S3: Determine segment number l s Is it equal to the total number of segments L?

[0141] If l s If it is equal to the total number of segments L, the final instantaneous angular frequency estimation of all short-term interferences has been completed, the loop ends, and the process jumps to step S5;

[0142] If l s If it is not equal to the total number of segments L, execute l s =Count ls +1, Count ls = l s , and jump to step S4 to extract the first s The final instantaneous angular frequency of the short-term interference;

[0143] S4: Through the first s The short-term interference and θ U , using the basis function space expansion and Gauss-Newton algorithm to extract the lth s The final instantaneous angular frequency of a short-term interference:

[0144] S4.1: Extract the first s The short-term interference is used to establish the basis function space U, the pole contraction coefficient ρ(l) and the regularization matrix Γ (l) , and initialize the number of iterations of the Newton-Gauss algorithm, l, the initial value of l is l = 1; Count l is a counter; Count l =1

[0145] s s (n3) = s(n3 + l s N w -N w ),n3=1,2,...,N e (l s )-N s (l s ), s s (n3) is the first s Short-term interference;

[0146]

[0147] Γ (l) =10 -Q (10-1)I Q

[0148]

[0149] Among them I Q represents the Q-order identity matrix;

[0150] S4.2: Transform the instantaneous angular frequency parameter vector θ U Multiply the basis function space U to obtain the estimated instantaneous angular frequency ω of the short-term interference s (n4):

[0151] ω s (n4) = θ U U(n4),n4=1,2,...,N w

[0152] S4.3: Estimated instantaneous angular frequency ω s (n4) performing cosine and multiplication operations to obtain time-varying filter coefficients a(n4);

[0153] a(n4)=-2cos[ω s (n4)],n4=1,2,...,N w

[0154] S4.4: Calculate the prediction error ε(n4) and the instantaneous angular frequency parameter vector gradient ψ(n4) through the time-varying filter transfer function;

[0155] In step S4.4, the prediction error ε(n4) and the instantaneous angular frequency parameter vector gradient ψ(n4) are calculated by the time-varying filter transfer function. The specific formula is:

[0156] ε(n4)=s s (n4)+a(n4)s s (n4-1)+s s (n4-2)-ρ(n4)a(n4)ε(n4-1)-ρ 2 (n4)ε(n4-2),n4=1,2,...,N w

[0157]

[0158] in is the gradient of the time-varying filter coefficient a(n4) relative to the prediction error ε(n4), and the specific formula is:

[0159] s F (n4) = s s (n4)-ρ 2 (n4)s F (n4-2)-ρ(n4)s F (n4-1)a(n4),n4=1,2,...,N w

[0160] ε F (n4)=ε(n4)-ρ2 (n4)ε F (n4-2)-ρ(n4)ε F (n4-1)a(n4),n4=1,2,...,N w

[0161]

[0162] Among them, s F (n4) and ε F (n4) is the intermediate variable;

[0163] S4.5: Determine the cost value β of the first cycle (l) and the cost value of the l-1th cycle β (l-1) Is the absolute error less than the threshold δ:

[0164] If the cost of the first cycle is β (l) and the cost value of the l-1th cycle β (l-1) If the absolute error is less than the threshold δ, then the current instantaneous angular frequency ω s (n4) is the final instantaneous angular frequency, go to step S3;

[0165] If the cost of the first cycle is β (l) and the cost value of the l-1th cycle β (l-1) If the absolute error of is greater than or equal to the threshold δ, execute step S4.6;

[0166] Cost value β (l) The specific formula is:

[0167]

[0168] The threshold value δ is 1×10 -6 , cost value β (0) =0;

[0169] The formula for extracting the final instantaneous angular frequency ω(n) is:

[0170]

[0171] ω(n)=ω s (nl s *N m +N m ),n=N ss ,N ss +1,...,N se

[0172] Where N ss Indicates extracting the first s The starting number of the final instantaneous angular frequency of the short-term interference, Nse Indicates extracting the first s The ending sequence number of the final instantaneous angular frequency of the short-term interference segment;

[0173] S4.6: Update the instantaneous angular frequency parameter vector θ using the Gauss-Newton algorithm U As the instantaneous angular frequency parameter vector of the l+1th cycle, execute step S4.2;

[0174] Update the instantaneous angular frequency parameter vector θ U The specific formula is:

[0175]

[0176] l=Count l +1

[0177] Count l = l

[0178] Where R (l) is the pseudo-Hessian matrix of the cost function of the lth cycle, is the gradient of the cost function of the lth cycle;

[0179] S5: Calculate the filtering result y(n) according to the final instantaneous angular frequency ω(n) of the signal and the time-varying filter transfer function;

[0180] y(n)=s(n)-2cos[ω(n)]s(n-1)+s(n-2)+2ρ e cos[ω(n)]y(n-1)-ρ e 2 y(n-2),n=1,2,...,N

[0181] where ρ e is the filter coefficient; filter coefficient ρ e The value is 0.995.

[0182] In the embodiment of the present invention, the model for simulating the sinusoidal frequency modulation interference s1(n) and the quadratic frequency modulation interference s2(n) is:

[0183]

[0184] Where A is the interference amplitude, f c is the carrier frequency, A f is the amplitude of the instantaneous frequency change of the sinusoidal FM interference, f f is the frequency of the instantaneous frequency change of the sinusoidal FM interference, is the initial phase of the instantaneous frequency change of the sinusoidal FM interference, F s is the sampling rate of interference, f1 is the starting frequency of secondary FM interference, f2 is the ending frequency of secondary FM interference, is the initial phase of the interference, N is the total length of the interference. w(n) has a mean of 0 and a variance of σ 2 The additive white Gaussian noise is determined by the interference-to-noise ratio JNR. The specific formula is:

[0185]

[0186] In order to accurately evaluate the filter effect, the instantaneous frequency mean square error is used as the evaluation index. The specific formula is:

[0187]

[0188] Among them, f r (n) is the real instantaneous frequency of the interference, f(n) is the instantaneous frequency to be filtered out, and the specific formula is:

[0189]

[0190] Embodiment 1:

[0191] In this embodiment, sinusoidal frequency modulation interference is used, and the parameters are set as follows: interference amplitude A=1, sampling rate F s =1000Hz, carrier frequency f c =320, the amplitude of the instantaneous frequency change of the interference is A f =40, the frequency of the instantaneous frequency change of the interference f f =2, the initial phase of the interference instantaneous frequency change Interference initial phase The total length of interference N = 4000, the interference-to-noise ratio JNR = 30dB, Figure 4 The short-time Fourier transform time-frequency distribution diagram of sinusoidal FM interference is given.

[0192] According to step 1, input interference s(n) and sampling rate F s , set the instantaneous angular frequency parameter vector θ U , and initialize the instantaneous angular frequency parameter vector

[0193] According to step 2, the input interference s(n) is segmented into several short-term interferences, and the total number of segments L=117 and the starting sequence number N of each short-term interference are calculated. s (l s ) and the end sequence number N of each short-term interference e (l s ), and initialize the paragraph number l s =0, counter Count ls =0.

[0194] According to step 3, determine the current segment number l sIs it equal to the total number of segments L?

[0195] If l s If it is equal to the total number of segments L, the final instantaneous angular frequency estimation of all short-term interferences has been completed, the loop ends, and jumps to step 5.

[0196] If l s If it is not equal to the total number of segments L, execute l s =Count ls +1, Count ls = l s , and jump to step S4 to extract the first s The final instantaneous angular frequency of a short-term interference.

[0197] According to step 4, through the first s The short-term interference and the instantaneous angular frequency parameter vector θ U , using the basis function space expansion and Gauss-Newton algorithm to extract the lth s The final instantaneous angular frequency of the short-term interference, go to step 3.

[0198] According to step 5, the filtering result y(n) is calculated according to the final instantaneous angular frequency ω(n) and the time-varying filter transfer function.

[0199] The comparison of the real instantaneous frequency curve of the simulated sinusoidal FM interference and the filtered instantaneous frequency curve is shown in the figure below. Figure 5 As shown in the figure, the mean square error of the instantaneous frequency is 0.328, and the filtered instantaneous frequency curve basically fits the real instantaneous frequency curve with a small error. This result also proves that the algorithm of the present invention can effectively use the Gauss-Newton algorithm to iteratively update the instantaneous angular frequency parameter vector, thereby obtaining an accurate instantaneous angular frequency and accurately filtering out interference. The comparison between the simulated sinusoidal frequency modulation interference and the filtering result is shown in the figure. Figure 6 As shown, the interference has been basically filtered out.

[0200] Embodiment 2:

[0201] In this embodiment, secondary frequency modulation interference is used, and the parameters are set as follows: interference amplitude A=1, sampling rate F s =1000Hz, the starting frequency of interference f1 = 200, the ending frequency of interference f2 = 410, the initial phase The total length of interference N = 4000, the interference-to-noise ratio JNR = 30dB, Figure 7 The Wigner-Wiley time-frequency distribution diagram of secondary FM interference is given.

[0202] According to step 1, input interference s(n) and sampling rate F s , set the instantaneous angular frequency parameter vector θ U, and use the short-time Fourier transform time-frequency analysis method to initialize the instantaneous angular frequency parameter vector θ U =[1.26,0.16,-1.89,5.91].

[0203] According to step 2, the input interference s(n) is segmented into several short-term interferences, and the total number of segments L=117 and the starting sequence number N of each short-term interference are calculated. s (l s ) and the end sequence number N of each short-term interference e (l s ), and initialize the paragraph number l s =0, counter Count ls =0.

[0204] According to step 3, determine the current segment number l s Is it equal to the total number of segments L?

[0205] If l s If it is equal to the total number of segments L, the final instantaneous angular frequency estimation of all short-term interferences has been completed, the loop ends, and jumps to step 5.

[0206] If l s If it is not equal to the total number of segments L, execute l s =Count ls +1, Count ls = l s , and jump to step S4 to extract the first s The final instantaneous angular frequency of a short-term interference.

[0207] According to step 4, through the first s The short-term interference and the instantaneous angular frequency parameter vector θ U , using the basis function space expansion and Gauss-Newton algorithm to extract the lth s The final instantaneous angular frequency of the short-term interference, go to step 3.

[0208] According to step 5, the filtering result y(n) is calculated according to the final instantaneous angular frequency ω(n) and the time-varying filter transfer function.

[0209] The comparison diagram of the real instantaneous frequency curve of the simulated secondary FM interference and the filtered instantaneous frequency curve is as follows: Figure 8As shown in the figure, the instantaneous frequency at the start of filtering is 199.63 Hz, the real start frequency is 200 Hz, the instantaneous frequency at the end of filtering is 409.80 Hz, the real end frequency is 410 Hz, the mean square error of the instantaneous frequency is 0.295, and the filtered instantaneous frequency curve basically fits the real instantaneous frequency curve with a small error. Such results also confirm that the algorithm of the present invention can effectively use the Gauss-Newton algorithm to iteratively update the instantaneous angular frequency parameter vector, thereby obtaining an accurate instantaneous angular frequency and accurately filtering out interference. The comparison diagram of the simulated secondary frequency modulation interference and the filtering result is shown in the figure. Fig. 9 As shown, the interference has been basically filtered out.

[0210] Embodiment 3:

[0211] In this embodiment, 100 independent experiments were performed under the conditions of interference-to-noise ratio of 30 dB, 40 dB and 50 dB respectively.

[0212] Among them, the secondary frequency modulation interference parameters are set as: interference amplitude A = 1, sampling rate F s =1000Hz, the starting frequency f1 of the interference is a random number in [200,300], the ending frequency f2 of the interference is a random number in [350,450], the initial phase is a random number in [0,π], and the total length of interference is N=4000.

[0213] The parameters of the sinusoidal FM interference are set as follows: interference amplitude A = 1, sampling rate F s =1000Hz, carrier frequency f c is a random number in [250,350], interfering with the amplitude of the instantaneous frequency change A f =40, the frequency of the instantaneous frequency change of the interference f f is a random number in [0,2], interfering with the initial phase of the instantaneous frequency change is a random number in [0,π], interfering with the initial phase is a random number in [0,2π], and the total length of interference is N=4000.

[0214] Under different signal-to-noise ratios, short-time Fourier transform and Wigner-Wiley time-frequency analysis methods are used to estimate the instantaneous angular frequency of interference, and a time-varying filter is designed according to the instantaneous angular frequency of interference for filtering, and a short-time adaptive time-varying filter based on the Gauss-Newton algorithm is used to calculate the interference suppression amount. At the same time, when the statistical interference-to-noise ratio JNR=30dB, the fixed value method, short-time Fourier transform time-frequency analysis method, and Wigner-Wiley time-frequency analysis method are used to set the instantaneous angular frequency parameter vector θ U The impact on the number of iterations when estimating the instantaneous angular frequency of the first short-term interference.

[0215] The interference suppression amount is calculated as follows:

[0216] △J=JNR1-JNR2

[0217] Wherein, JNR1 represents the interference-to-noise ratio before filtering, and JNR2 represents the interference-to-noise ratio after filtering.

[0218] The results are summarized in the following table:

[0219] Table 1 Average interference suppression amount of different methods

[0220]

[0221]

[0222] Table 2 Average number of iterations for different initial methods to estimate the instantaneous angular frequency of the first short-term interference

[0223]

[0224] As can be seen from Table 1, the short-time adaptive time-varying filter based on the Gauss-Newton algorithm has the best effect. When the interference-to-noise ratio is 20dB and 30dB, this method can almost completely filter out the interference; when the interference-to-noise ratio is 40dB, it can also effectively filter out most of the interference. In contrast, the two methods of estimating the instantaneous angular frequency of the interference using short-time Fourier transform and Wigner-Wiley time-frequency analysis method, and designing a time-varying filter based on the frequency for filtering, have large errors in estimating the instantaneous angular frequency of the interference, and there are more residual interferences after filtering. In the case of sinusoidal frequency modulation interference, the Wigner-Wiley time-frequency analysis method is easily affected by cross terms, thus misestimating the instantaneous angular frequency of the interference, resulting in the failure of this method. However, in the case of secondary frequency modulation interference, the Wigner-Wiley time-frequency analysis method is more accurate in estimating the instantaneous angular frequency of the interference, and its filtering effect is better than the time-varying filter based on short-time Fourier transform.

[0225] As shown in Table 2, the instantaneous angular frequency parameter vector θ is set by using the short-time Fourier transform time-frequency analysis method or the Wigner-Wiley time-frequency analysis method. U When estimating the instantaneous angular frequency of the first short-term interference, compared with the fixed value method, the number of iterations can be reduced, the accuracy can be improved, and the incorrect estimation of the instantaneous angular frequency of the interference can be prevented.

[0226] The above description is only a specific implementation of the present application, so that those skilled in the art can understand or implement 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 range consistent with the principles and novel features applied for herein.

Claims

1. A short-time adaptive time-varying filtering method based on Gauss-Newton algorithm, characterized in that: The steps include: S1: According to the input interference s(n) and sampling rate F s , set the instantaneous angular frequency parameter vector θ U , and initialize; Set the instantaneous angular frequency parameter vector θ U The method is a fixed value method or a time-frequency analysis method; the time-frequency analysis method includes a short-time Fourier transform time-frequency analysis method and a Wigner-Wiley time-frequency analysis method; S2: Divide the input interference s(n) into several short-term interferences, with a total number of segments of L, l s is the current segment number; S3: Traverse all short-term interferences and determine the current segment number l s Is it equal to the total number of segments L? If l s If it is equal to the total number of segments L, the final instantaneous angular frequency estimation of all short-term interferences has been completed, the loop ends, and the process jumps to step S5; If l s If it is not equal to the total number of segments L, execute l s =Count ls +1, Count ls = l s , and jump to step S4 to extract the first s The final instantaneous angular frequency of the short-term interference; S4: Through the first s The short-term interference and the instantaneous angular frequency parameter vector θ U , using the basis function space expansion and Gauss-Newton algorithm to extract the lth s The final instantaneous angular frequency ω(n) of the short-term interference is obtained, and then the process goes to step S3; S5: Calculate the filtering result y(n) according to the final instantaneous angular frequency ω(n) of the interference and the time-varying filter transfer function.

2. The short-time adaptive time-varying filtering method based on Gauss-Newton algorithm according to claim 1, characterized in that: The steps of setting the instantaneous angular frequency parameter vector in the Wigner-Wiley time-frequency analysis method are as follows: S1-1.1: According to the input interference s(n) and sampling rate F s , calculate the Wigner-Wiley distribution WVD(n1,ω s1 ); S1-1.2: Calculate the Wigner-Wiley distribution WVD(n1,ω s1 ) time vector t w1 The angular frequency vector ω with the maximum energy corresponding to each moment w1 ; S1-1.3: Using the time vector t w1 The angular frequency vector ω with the maximum energy corresponding to each moment w1 Calculate the estimated instantaneous angular frequency parameter vector.

3. The short-time adaptive time-varying filtering method based on Gauss-Newton algorithm according to claim 1 is characterized in that: The steps of the short-time Fourier transform time-frequency analysis method are: S1-2.1: Input disturbance s(n) and sampling rate F s , calculate the short-time Fourier transform S(n2,ω s2 ); The sampling point number Angular frequency w(m w2 ) represents the window function; N w2 Indicates the length of the window function, the window function w(m w2 ) is expressed as: S1-2.2: Calculate the short-time Fourier transform S(n2,ω s2 ) time vector t w2 The angular frequency vector ω with the maximum energy corresponding to each moment w2 , the specific formula is: S1-2.3: Using the time vector t w2 The angular frequency vector ω with the maximum energy corresponding to each moment w2 Calculate the estimated instantaneous angular frequency parameter vector θ U , the specific formula is: in, Represents the time vector t w2 The i-th power of each element in, Q represents the maximum order of the basis function space, U w2 is the time vector t w2 The corresponding basis function space is represents the Moore–Penrose pseudoinverse.

4. The short-time adaptive time-varying filtering method based on Gauss-Newton algorithm according to claim 1, characterized in that: In step S2, the specific steps of segmenting the input interference s(n) into a number of short-term interferences are: The input interference s(n) is divided into several short-term interferences, the total number of segments is L, and the starting sequence number of each short-term interference is N. s (l s ), the end sequence number N of each short-term interference e (l s ), l s is the current segment number, Count ls for l s counter; Among them, the total number of segments is L and the starting sequence number is N s (l s ), end sequence number N e (l s ), segment number l s Initial value, Count ls The initial values ​​are: N s (L s )=(l s -1)N m +1; l s ＝0； Count ls =0; Where N represents the total length of the input interference s(n), N m Indicates the length of movement when extracting each short-term interference; Indicates rounding down.

5. The short-time adaptive time-varying filtering method based on Gauss-Newton algorithm according to claim 1, characterized in that: In the step S4, the first s The short-term interference and the instantaneous angular frequency parameter vector θ U , using the basis function space expansion and Gauss-Newton algorithm to extract the lth s The specific steps to calculate the final instantaneous angular frequency ω(n) of a short-term interference are: S4.1: Extract the first s The short-term interference is used to establish the basis function space U, the pole contraction coefficient ρ(l) and the regularization matrix Γ (l) , and initialize the number of iterations of the Newton-Gauss algorithm, l, the initial value of l is l = 1; Count l is a counter; Count l =1s s (n3) = s(n3 + l s N w -N w ),n3=1,2,...,N e (l s )-N s (l s ), s s (n3) is the first s Short-term interference; C (l) =10 -Q (10-l)I Q Among them I Q represents the Q-order identity matrix; S4.2: Transform the instantaneous angular frequency parameter vector θ U Multiply the basis function space U to obtain the estimated instantaneous angular frequency ω of the short-term interference s (n4): ω s (n4)=θ U U(n4),n4=1,2,...,N w S4.3: Estimated instantaneous angular frequency ω s (n4) performing cosine and multiplication operations to obtain time-varying filter coefficients a(n4); a(n4)=-2cos[ω s (n4)],n4=1,2,...,N w S4.4: Calculate the prediction error ε(n4) and the instantaneous angular frequency parameter vector gradient ψ(n4) through the time-varying filter transfer function; The prediction error ε(n4) and the instantaneous angular frequency parameter vector gradient ψ(n4) are calculated by the time-varying filter transfer function in step S4.

4. The specific formula is: ε(n4)=s s (n4)+a(n4)s s (n4-1)+s s (n4-2)-ρ(n4)a(n4)ε(n4-1)-ρ 2 (n4)ε(n4-2),n4=1,2,...,N w in is the gradient of the time-varying filter coefficient a(n4) relative to the prediction error ε(n4), and the specific formula is: s F (n4)=s s (n4)-ρ 2 (n4)s F (n4-2)-ρ(n4)s F (n4-1)a(n4),n4=1,2,...,N w e F (n4)=ε(n4)-ρ 2 (n4)e F (n4-2)-p(n4)e F (n4-1)a(n4),n4=1,2,...,N w Among them, s F (n4) and ε F (n4) is the intermediate variable; S4.5: Determine the cost value β of the first cycle (l) and the cost value of the l-1th cycle β (l-1) Is the absolute error less than the threshold δ: If the cost of the first cycle is β (l) and the cost value of the l-1th cycle β (l-1) If the absolute error is less than the threshold δ, then the current instantaneous angular frequency ω s (n4) is the final instantaneous angular frequency, go to step S3; If the cost of the first cycle is β (l) and the cost value of the l-1th cycle β (l-1) If the absolute error of is greater than or equal to the threshold δ, execute step S4.6; Cost value β (l) The specific formula is: The threshold value δ is 1×10 -6 ; The formula for extracting the final instantaneous angular frequency ω(n) is: ω(n)=ω s (nl) s *N m +N m ),n=N ss ,N ss +1,...,N se Where N ss Indicates extracting the first s The starting number of the final instantaneous angular frequency of the short-term interference, N se Indicates extracting the first s The ending sequence number of the final instantaneous angular frequency of the short-term interference segment; S4.6: Update the instantaneous angular frequency parameter vector θ using the Gauss-Newton algorithm U As the instantaneous angular frequency parameter vector of the l+1th cycle, execute step S4.2; Update the instantaneous angular frequency parameter vector θ U The specific formula is: l=Count l +1 Count l =l Where R (l) is the pseudo-Hessian matrix of the cost function of the lth cycle, ▽β (l) is the gradient of the cost function at the lth cycle.

6. The short-time adaptive time-varying filtering method based on Gauss-Newton algorithm according to claim 1, characterized in that: In step S5, the specific steps of calculating the filtering result y(n) according to the final instantaneous angular frequency ω(n) of the interference and the time-varying filter transfer function are: y(n)=s(n)-2cos[ω(n)]s(n-1)+s(n-2)+2ρ e cos[ω(n)]y(n-1)-ρ e 2 y(n-2),n=1,2,...,N where ρ e is the filter coefficient; filter coefficient ρ e The value is 0.

995.

7. The short-time adaptive time-varying filtering method based on Gauss-Newton algorithm according to claim 2, characterized in that: The calculation of the Wigner-Wiley distribution WVD(n1,ω s1 ) steps are: Wigner-Wiley distribution WVD(n1,ω s1 )The specific formula is: The sampling point number n1 ranges from [1,2,...,N w ],N w Indicates the length of short-term interference, angular frequency ω s1 The value is z(n1) represents the analytical signal, and the specific formula is: Here, j represents the imaginary number symbol.

8. The short-time adaptive time-varying filtering method based on Gauss-Newton algorithm according to claim 2, characterized in that: Calculate the Wigner-Wiley distribution WVD(n1,ω s1 ) time vector t w1 The angular frequency vector ω with the maximum energy corresponding to each moment w1 The specific steps are: Among them, ω w1 is the angular frequency vector with the maximum energy corresponding to each moment.

9. The short-time adaptive time-varying filtering method based on Gauss-Newton algorithm according to claim 2, characterized in that: Using the time vector t w1 The angular frequency vector ω with the maximum energy corresponding to each moment w1 Calculate the estimated instantaneous angular frequency parameter vector θ U , the specific steps are: in, Represents the time vector t w1 The i-th power of each element in, Q represents the maximum order of the basis function space, U w1 is the time vector t w1 The corresponding basis function space is represents the Moore–Penrose pseudoinverse.

Images