CN108680786B - Adaptive filtering envelope extraction method for pulse signal frequency domain - Google Patents

Abstract

The invention discloses a method for extracting a pulse signal frequency domain self-adaptive filtering envelope, which comprises the following steps: the first step is as follows: acquiring a sampling data sequence of a pulse signal to be processed; the second step is that: calculating a power spectrum of the sampling data sequence; the third step: iterative smoothing and filter bandwidth parameter initialization; the fourth step: carrying out iterative smoothing processing on the power spectrum to obtain a smooth power spectrum; the fifth step: extracting power spectrum main lobe characteristic parameters; and a sixth step: generating adaptive filtering parameters; the seventh step: self-adaptive filtering in a frequency domain; eighth step: and extracting the envelope of the pulse signal. The extraction method generates the self-adaptive filtering parameters by extracting the main lobe characteristics of the pulse signals in real time, carries out self-adaptive filtering on the pulse signals, can realize an envelope extraction method matched with the type and the frequency of the pulse signals under the non-cooperative condition, has small operand and strong practicability, and is suitable for carrying out real-time processing on the signals.

Description

Adaptive filtering envelope extraction method for pulse signal frequency domain

Technical Field

The invention belongs to the technical field of signal processing, and particularly relates to a pulse signal frequency domain adaptive filtering envelope extraction method.

Background

The signal envelope is an important parameter for representing the time domain characteristics of the pulse signals, the accurate extraction of the pulse signal envelope can provide a basis for further accurate estimation of the arrival time and the pulse width of the pulse signals, the time delay difference of the two paths of pulse signal envelopes can be used for estimating the incoming wave direction of the signals, and the method has important theoretical and application values in the fields of sonar, radar, electronic warfare and the like, and is particularly prominent in the aspect of radar and sonar signal processing.

The conventional signal envelope extraction methods mainly comprise four types, namely (1) a Hilbert transform method, (2) a wavelet transform method, (3) an improved Hilbert transform method, such as an envelope extraction method based on Hilbert transform (or called Hilbert-yellow transform) of empirical mode decomposition and a Hilbert transform method based on blind source extraction, and (4) an envelope extraction method based on mathematical morphology processing, wherein the Hilbert transform method firstly converts an original signal into a complex analysis signal, and then takes a modulus value of the complex analysis signal as an envelope of the signal, and the method is based on FFT, simple in implementation, small in operand and good in envelope extraction effect under high signal-to-noise ratio, but is performed in the whole frequency interval and poor in noise suppression capability, and is very useful for envelope extraction of a narrow-band carrier signal, but has two inevitable defects, namely ① extraction of high-frequency components in the signal, so that time-varying adaptive analysis is not achieved ②.

The wavelet transformation method firstly carries out wavelet transformation on signals to obtain analytic signals, wherein the wavelet transformation is equivalent to multi-passband filtering on the signals by utilizing a series of bandpass filters, so that the obtained analytic signals are subjected to the bandpass filtering, and the real parts and the imaginary parts of the obtained analytic signals are orthogonal, and then the modulus values of the analytic signals are taken to obtain the envelopes of the signals. Compared with the Hilbert transform method, the envelope extracted by the method is more ideal, but the method has the difficulty that the selection of the wavelet transform scale is proper, the scale selection is proper, the frequency band covered by the band-pass filter can be matched with the signal frequency, and the envelope extraction effect is good; the scale selection is improper, so that the frequency band covered by the band-pass filter is mismatched with the signal frequency, the envelope extraction effect is rapidly deteriorated, the selection of the scale depends on the experience of a user to a great extent, and the universality is poor. The improved Hilbert transform method firstly utilizes a certain specific preprocessing method to carry out filtering and denoising processing on signals, and then further utilizes Hilbert transform to extract envelopes, such as the empirical mode decomposition and blind source extraction, and the like. The processing method based on mathematical morphology comprises the steps of preprocessing signals by utilizing morphological filtering, extracting signal envelopes by utilizing morphological closed operation, and eliminating noise envelopes by utilizing morphological open operation.

Disclosure of Invention

The purpose of the invention is as follows: aiming at the problems and the defects in the prior art, the invention provides a pulse signal frequency domain adaptive filtering envelope extraction method to meet the requirements of radar and sonar signal processing. In addition, the method has the same operation amount as the traditional Hilbert transform envelope extraction method, but has better noise suppression capability and self-adaption capability compared with the Hilbert transform method, clear principle, simple realization and strong engineering applicability.

The technical scheme is as follows: in order to achieve the purpose, the invention adopts the following technical scheme: a method for extracting a pulse signal frequency domain adaptive filtering envelope comprises the following steps:

(1) acquiring a pulse signal sampling data sequence to be processed: receiving real-time acquisition data of N sampling points from a sensor as a data sequence x (N) to be processed, wherein N is 0,1, …, N-1, or extracting N sampling point data containing the whole pulse signal from a memory as the data sequence x (N) to be processed, wherein N is 0,1, …, N-1, the number of the sampling points containing the whole pulse signal is N, and the value is an integer power of 2.

(2) Calculating a pulse signal power spectrum P (k) according to the data sequence x (n).

(3) Initialization of parameters related to the bandwidth of the iterative smoothing and filter: setting a maximum iteration number threshold of power spectrum iterative smoothing, wherein the precision control threshold and the iterative smoothing window length are respectively I₁，₁And M₁The maximum iteration number threshold, the precision control threshold and the iteration smoothing window length of the initial envelope iteration smoothing are respectively I₂，₂And M₂Bandwidth coefficient omega of narrow-band signal₁And bandwidth factor omega of broadband signal₂And judging characteristic parameter thresholds η of main lobes of the narrow-band pulse and the wide-band pulse.

(4) And performing iterative smoothing processing on the P (k) to obtain a smoothed power spectrum S (k).

(5) Extracting the original power spectrum and the main lobe characteristic parameter gamma of the smoothed power spectrum_pAnd gamma_S。

(6) Adaptive filtering parameters h (k) are generated.

(7) And carrying out adaptive filtering on the pulse signal according to the filtering parameters H (k).

(8) The pulse signal envelope y (n) is extracted.

Preferably, in the step (2), fast fourier transform is performed on the data sequence x (n), and discrete fourier transform x (l) and power spectrum p (k) of the data sequence are obtained through calculation, specifically including the following steps:

step 2-1: calculating the discrete Fourier transform of x (n):

where l is the discrete frequency index of X (l), and j represents an imaginary unit, i.e.

This equation can be implemented by fast fourier transform.

Step 2-2: calculating the power spectrum of x (n) according to x (l):

where k is the discrete frequency index of the power spectrum p (k), | represents the modulo operation, and | x (k) | is the signal amplitude spectrum.

Preferably, in the step (2), the discrete fourier transform x (l) of the step 2-1 calculating x (n) can be realized by fast fourier transform;

preferably, in step (3), a maximum iteration number threshold of power spectrum iterative smoothing is set, and the precision control threshold and the iterative smoothing window length are respectively I₁，₁And M₁The maximum iteration number threshold, the precision control threshold and the iteration smoothing window length of the initial envelope iteration smoothing are respectively I₂，₂And M₂Bandwidth coefficient omega of narrow-band signal₁And bandwidth factor omega of broadband signal₂Judging the characteristic parameter threshold η of the main lobe of the narrow-band pulse and the wide-band pulse, wherein I₁And I₂A value is a positive integer greater than 2,₁and₂the value is positive number less than 1, the length of the iterative smoothing window of the power spectrum satisfies 3 and is less than or equal to M₁An odd number less than or equal to N/2-2, and an envelope iteration smoothing window length of 5 less than or equal to M₂Odd number, omega, not more than 49₁To satisfy 0 ≤ ω₁Any number, omega, not more than 0.05₂To satisfy 0.05 < omega₂Any number not more than 0.1, η is any number satisfying 1.5 not more than η not more than 2.5₁The value is 10, and the number is,₁value of 0.01, M₁A value of 5, I₂The value is 5, and the number of the grooves,₂a value of 0.02, M₂Value of 11, omega₁Has a value of 0.05, omega₂Is 0.1 and η is 2.0.

Preferably, in the step (4), the iterative smoothing processing is performed on the power spectrum p (k) to obtain a smoothed power spectrum s (k), and the method specifically includes the following steps:

step 4-1: the iterative smoothing times i of the initialized power spectrum is 0, and the iterative smoothing result S of the power spectrum₀(k) Comprises the following steps:

S₀(k)＝P(k),k＝0,1,2…,N/2-1

step 4-2: let the iteration number i be i +1, and smooth the last iteration of the power spectrum to obtain a result S_i-1(k) Performing smoothing to obtain current smoothing result S_i(k):

Step 4-3: judging whether the maximum iterative smoothing times of the power spectrum is reached, namely judging that I is less than or equal to I₁If yes, entering the step 4-4, otherwise jumping to the step 4-6;

step 4-4: respectively calculating last smoothing results S of power spectrum_i-1(k) And current smoothing result S_i(k) Sum of squares J of residuals from original power spectrum p (k)_i-1And J_iAre respectively as

And 4-5: judgment of | J_i-1-J_i|≤₁J_iIf yes, entering the step 4-6, otherwise returning to the step 4-2; wherein₁A threshold is controlled for power spectrum iterative smoothing precision;

and 4-6: let S (k) become S_i(k) And k is 0,1,2 …, N/2-1, and a smoothed power spectrum s (k) is obtained.

Preferably, in the step (5), the original power spectrum and the smoothed power spectrum main lobe characteristic parameter γ are extracted_pAnd gamma_SThe method specifically comprises the following steps:

step 5-1: searching discrete frequency index k corresponding to maximum value of smoothed power spectrum S (k)_M

Wherein

Representing that the discrete frequency index corresponding to the maximum value of S (k) is searched within the range of 1 ≦ k ≦ N/2-1;

step 5-2: subtracting half of the maximum value from s (k) to obtain a half-power difference power spectrum z (k):

Z(k)＝S(k)-S(k_M)/2,k＝0,1,2…,N/2-1

step 5-3: searching for the discrete frequency index k of the first zero crossing on the left of the Z (k) peak_LThe searching process comprises the following steps:

step 5-3-1: initializing a frequency start search index k of a first zero crossing on the left side of a Z (k) peak value_l＝k_M-1；

Step 5-3-2: judgment of k _l0 or Z (k)_l) If the result is less than 0, jumping to the step 5-3-4 if the result is true, and otherwise, entering the step 5-3-3;

step 5-3-3: let k_l＝k _l1 and returning to step 5-3-2;

step 5-3-4: let k_L＝k_lObtaining the discrete frequency index k of the first zero-crossing point on the left side of the Z (k) peak value_L；

Step 5-4: searching for the first zero-crossing discrete frequency index k to the right of the Z (k) peak_RThe searching process comprises the following steps:

step 5-4-1: initializing a frequency start search index k of a first zero crossing on the right side of a Z (k) peak value_r＝k_M+1；

Step 5-4-2: judgment of k_rN/2-1 or Z (k)_r) If the value is less than 0, jumping to the step 5-4-4 if the value is less than 0, otherwise, entering the step 5-4-3;

step 5-4-3: let k_r＝k_r+1 and returning to step 5-4-2;

step 5-4-4: let k_R＝k_rObtaining the discrete frequency index k of the first zero-crossing point on the right side of the Z (k) peak value_R。

Step (ii) of5-5: calculating the primary power spectrum P (k) main lobe characteristic parameter gamma_p

And 5-6: calculating a main lobe characteristic parameter gamma of the smoothed power spectrum S (k)_S

Preferably, in the step (6), generating the adaptive filtering parameter h (k) specifically includes the following steps:

step 6-1: initializing adaptive filter parameters

H(k)＝0,k＝0,1,2…,N-1

Step 6-2: respectively calculating the discrete frequency indexes T of the start and the end of the transition band of the narrow-band pulse signal filter_N1And T_N2And a passband start and end discrete frequency index P_N1And P_N2：

T_N1＝max(round[(1-0.5ω₁)k_M],0)，T_N2＝min(round[(1+0.5ω₁)k_M],N/2-1)

P_N1＝max(k_M-1,T_N1)，P_N2＝min(k_M+1,T_N2)

Where round [ ] represents a rounding operation, and max (,) and min (,) represent a large and small value operation, respectively.

Step 6-3: respectively calculating the start and end discrete frequency indexes T of the transition band of the broadband pulse signal filter_W1And T_W2And a passband start and end discrete frequency index P_W1And P_W2：

T_W1＝max(round[(1-0.5ω₂)k_M],0)，T_W2＝min(round[(1+0.5ω₂)k_M],N/2-1)

P_W1＝max(k_L-1,T_W1)，P_W2＝max(k_R+1,T_W2)

Step 6-4: obtaining discrete frequency index T from adaptive filter transition band start and end according to signal type₁And T₂And a passband start and end discrete frequency index P₁And P₂：

Namely, max (gamma)_p,γ_s)/min(γ_p,γ_s) If > η is true, if true, the signal is judged to be a narrow-band pulse signal and the following steps are carried out:

T₁＝T_N1,T₂＝T_N2,P₁＝P_N1,P₂＝P_N2

otherwise, judging the signal to be a broadband pulse signal and ordering:

T₁＝T_W1,T₂＝T_W2,P₁＝P_W1,P₂＝P_W2

step 6-5: the transition band part of the adaptive filter coefficients H (k) is modified so that k satisfies T₁≤k≤T₂The part of (A) is:

step 6-6: the pass band part of the adaptive filter coefficients H (k) is modified such that k satisfies P₁≤k≤P₂The part of (A) is: h (k) ═ 1.

Preferably, in step (7), the result of discrete fourier transform x (k) of the pulse signal is adaptively filtered to obtain a filtered discrete fourier transform y (k):

Y(k)＝X(k)H(k),k＝0,1,2…,N-1

preferably, in the step (8), the extracting of the pulse signal envelope y (N), where N is 0,1 …, N-1 specifically includes the following steps:

step 8-1: performing inverse discrete Fourier transform on the filtered discrete Fourier transform Y (k) to obtain a pulse signal filtered time-domain complex signal y_c(n)：

Step 8-2: taking a complex number y_c(n) obtaining the initial envelope data sequence y of the pulse signal by the modulus value₀(n)：

y₀(n)＝|y_c(n)|,n＝0,1…,N-1

Step 8-3: initializing the iterative smoothing times q of the envelope data to be 0;

step 8-4: let the iteration number q be q +1 and smooth the result y for the last iteration of the envelope data sequence_q-1(n) smoothing to obtain current smoothing result y_q(n)：

And 8-5: judging whether the maximum iterative smoothing times of the envelope is reached, namely judging that q is less than or equal to I₂If yes, entering the step 8-6, otherwise jumping to the step 8-8;

and 8-6: respectively calculating last smoothing result y of envelope data sequence_q-1(n) and current smoothing result y_q(n) and the initial envelope data sequence y₀Residual sum of squares Ψ of (n)_q-1And Ψ_qRespectively is as follows:

and 8-7: determine | Ψ_q-1-Ψ_q|≤₂Ψ_qIf yes, entering the step 8-8, otherwise returning to the step 8-4;

and 8-8: let y (n) be y_q(N), N is 0,1,2 …, N-1, resulting in a sequence of signal envelope data:

y(n),n＝0,1,2…,N-1。

preferably, in the step (8), the step 8-1 performs inverse discrete fourier transform on the filtered discrete fourier transform y (k) to obtain the pulse signal filtered time-domain complex signal y_c(n) may be implemented by a fast fourier transform.

Has the advantages that: compared with the prior art, the method has the following beneficial effects:

(1) the envelope extraction method of the invention automatically judges the type of the pulse signal under the non-cooperative condition by extracting the main lobe characteristic of the pulse signal in real time, estimates the bandwidth parameter of the signal in real time, generates the adaptive filter parameter matched with the type and the frequency of the pulse signal, and carries out adaptive filtering on the pulse signal, thereby obtaining the signal processing gain of approximate matched filtering.

(2) According to the envelope extraction method, after the discrete Fourier transform result of the filtered signal is subjected to the inverse discrete Fourier transform and a modulus value is taken to obtain a signal initial envelope data sequence, the initial envelope is further subjected to iterative smoothing processing, burrs existing in the initial envelope are reduced, and the problem that the burrs in the envelope are more due to the existence of high-frequency components in the conventional Hilbert transform method is effectively solved.

Drawings

FIG. 1 is a schematic flow diagram of the process of the present invention;

FIG. 2 is the original power spectrum of the simulated pulse signal of example 1;

FIG. 3 is the power spectrum of the simulation pulse signal after smoothing in example 1;

FIG. 4 is the adaptive filter coefficients of example 1;

FIG. 5 is a simulated signal envelope extracted in example 1;

FIG. 6 is the original power spectrum of the simulated pulse signal of example 2;

FIG. 7 is the power spectrum of the simulation pulse signal after smoothing in example 2;

FIG. 8 is the adaptive filter coefficients of example 2;

fig. 9 is a simulated signal envelope extracted in example 2.

Detailed Description

The invention is further described with reference to the following figures and examples:

as shown in fig. 1, a method for extracting a frequency domain adaptive filtering envelope of an impulse signal includes the following steps:

(1) acquiring a pulse signal sampling data sequence to be processed: receiving real-time acquisition data of N sampling points from a sensor as a data sequence x (N) to be processed, wherein N is 0,1, …, N-1, or extracting N sampling point data containing the whole pulse signal from a memory as the data sequence x (N) to be processed, wherein N is 0,1, …, N-1, the number of the sampling points containing the whole pulse signal is N, and the value is an integer power of 2.

(2) Calculating a pulse signal power spectrum P (k) from the data sequence x (n): performing fast Fourier transform on the data sequence x (n), and calculating to obtain discrete Fourier transform X (l) and a power spectrum P (k) of the data sequence, wherein the method specifically comprises the following steps:

step 2-1: calculating the discrete Fourier transform of x (n):

where l is the discrete frequency index of X (l), and j represents an imaginary unit, i.e.

This equation can be implemented by fast fourier transform.

Step 2-2: calculating the power spectrum of x (n) according to x (l):

where k is the discrete frequency index of the power spectrum p (k), | represents the modulo operation, and | x (k) | is the signal amplitude spectrum.

In the step (2), the discrete fourier transform x (l) of the step 2-1 for calculating x (n) can be realized by fast fourier transform, so that the calculation efficiency of the algorithm is improved.

(3) Initialization of parameters related to the bandwidth of the iterative smoothing and filter: setting a maximum iteration number threshold of power spectrum iterative smoothing, wherein the precision control threshold and the iterative smoothing window length are respectively I₁，₁And M₁Threshold of maximum iteration number of initial envelope iteration smoothing, precision controlThe system threshold and the iterative smoothing window are respectively I₂，₂And M₂Bandwidth coefficient omega of narrow-band signal₁And bandwidth factor omega of broadband signal₂Judging characteristic parameter threshold η of main lobe of narrow-band pulse and wide-band pulse, where I₁And I₂A value is a positive integer greater than 2,₁and₂the value is positive number less than 1, the length of the iterative smoothing window of the power spectrum satisfies 3 and is less than or equal to M₁An odd number less than or equal to N/2-2, and an envelope iteration smoothing window length of 5 less than or equal to M₂An odd number of not more than 49, η is an arbitrary number satisfying 1.5 not more than η not more than 2.5, omega₁To satisfy 0 ≤ ω₁Any number, omega, not more than 0.05₂To satisfy 0.05 < omega₂Any number less than or equal to 0.1. Preferred is I₁The value is 10, and the number is,₁value of 0.01, M₁A value of 5, I₂The value is 5, and the number of the grooves,₂a value of 0.02, M₂Value of 11, omega₁Has a value of 0.05, omega₂Is 0.1 and η is 2.0.

(4) Performing iterative smoothing processing on P (k) to obtain a smoothed power spectrum S (k): performing iterative smoothing processing on the power spectrum P (k) to obtain a smoothed power spectrum S (k), and specifically comprising the following steps:

step 4-1: the iterative smoothing times i of the initialized power spectrum is 0, and the iterative smoothing result S of the power spectrum₀(k) Comprises the following steps:

S₀(k)＝P(k),k＝0,1,2…,N2-1

step 4-2: let the iteration number i be i +1, and smooth the last iteration of the power spectrum to obtain a result S_i-1(k) Performing smoothing to obtain current smoothing result S_i(k)：

Step 4-3: judging whether the maximum iterative smoothing times of the power spectrum is reached, namely judging that I is less than or equal to I₁If yes, entering the step 4-4, otherwise jumping to the step 4-6;

step 4-4: respectively calculating last smoothing results S of power spectrum_i-1(k) And current smoothing result S_i(k) Residue of the original power spectrum P (k)Sum of squared differences J_i-1And J_iRespectively is as follows:

and 4-5: judgment of | J_i-1-J_i|≤₁J_iIf yes, entering the step 4-6, otherwise returning to the step 4-2; wherein₁A threshold is controlled for power spectrum iterative smoothing precision;

and 4-6: let S (k) become S_i(k) And k is 0,1,2 … and N2-1, and a smoothed power spectrum s (k) is obtained.

(5) Extracting the original power spectrum and the main lobe characteristic parameter gamma of the smoothed power spectrum_pAnd gamma_SThe method specifically comprises the following steps:

step 5-1: searching discrete frequency index k corresponding to maximum value of smoothed power spectrum S (k)_M

Wherein

Representing that the discrete frequency index corresponding to the maximum value of S (k) is searched within the range of 1 ≦ k ≦ N/2-1;

step 5-2: subtracting half of the maximum value from s (k) to obtain a half-power difference power spectrum z (k):

Z(k)＝S(k)-S(k_M)/2,k＝0,1,2…,N/2-1

step 5-3: searching for the discrete frequency index k of the first zero crossing on the left of the Z (k) peak_LThe searching process comprises the following steps:

step 5-3-1: initializing a frequency start search index k of a first zero crossing on the left side of a Z (k) peak value_l＝k_M-1；

Step 5-3-2: judgment of k _l0 or Z (k)_l) If the result is less than 0, jumping to the step 5-3-4 if the result is true, and otherwise, entering the step 5-3-3;

step 5-3-3: let k_l＝k _l1 and returning to step 5-3-2;

step 5-3-4: let k_L＝k_lObtaining the discrete frequency index k of the first zero-crossing point on the left side of the Z (k) peak value_L。

Step 5-4: searching for the first zero-crossing discrete frequency index k to the right of the Z (k) peak_RThe searching process comprises the following steps:

step 5-4-1: initializing a frequency start search index k of a first zero crossing on the right side of a Z (k) peak value_r＝k_M+1；

Step 5-4-2: judgment of k_rN/2-1 or Z (k)_r) If the value is less than 0, jumping to the step 5-4-4 if the value is less than 0, otherwise, entering the step 5-4-3;

step 5-4-3: let k_r＝k_r+1 and returning to step 5-4-2;

step 5-4-4: let k_R＝k_rObtaining the discrete frequency index k of the first zero-crossing point on the right side of the Z (k) peak value_R。

Step 5-5: calculating the primary power spectrum P (k) main lobe characteristic parameter gamma_p：

And 5-6: calculating a main lobe characteristic parameter gamma of the smoothed power spectrum S (k)_S：

(6) Generating adaptive filtering parameters h (k), specifically including the following steps:

step 6-1: initializing adaptive filter parameters

H(k)＝0,k＝0,1,2…,N-1

Step 6-2: respectively calculating the discrete frequency indexes T of the start and the end of the transition band of the narrow-band pulse signal filter_N1And T_N2And a passband start and end discrete frequency index P_N1And P_N2：

T_N1＝max(round[(1-0.5ω₁)k_M],0)，T_N2＝min(round[(1+0.5ω₁)k_M],N/2-1)

P_N1＝max(k_M-1,T_N1)，P_N2＝min(k_M+1,T_N2)

Where round [ ] represents a rounding operation, and max (,) and min (,) represent a large and small value operation, respectively.

Step 6-3: respectively calculating the start and end discrete frequency indexes T of the transition band of the broadband pulse signal filter_W1And T_W2And a passband start and end discrete frequency index P_W1And P_W2：

T_W1＝max(round[(1-0.5ω₂)k_M],0)，T_W2＝min(round[(1+0.5ω₂)k_M],N/2-1)

P_W1＝max(k_L-1,T_W1)，P_W2＝max(k_R+1,T_W2)

Step 6-4: obtaining discrete frequency index T from adaptive filter transition band start and end according to signal type₁And T₂And a passband start and end discrete frequency index P₁And P₂：

Namely, max (gamma)_p,γ_s)/min(γ_p,γ_s) If > η is true, if true, the signal is judged to be a narrow-band pulse signal and the following steps are carried out:

T₁＝T_N1,T₂＝T_N2,P₁＝P_N1,P₂＝P_N2

otherwise, judging the signal to be a broadband pulse signal and ordering:

T₁＝T_W1,T₂＝T_W2,P₁＝P_W1,P₂＝P_W2

step 6-5: the transition band part of the adaptive filter coefficients H (k) is modified so that k satisfies T₁≤k≤T₂The part of (A) is:

step 6-6: the pass band part of the adaptive filter coefficients H (k) is modified such that k satisfies P₁≤k≤P₂The part of (A) is: h (k) ═ 1;

(7) adaptively filtering the pulse signal according to the filter parameter h (k), that is, adaptively filtering the result x (k) of the discrete fourier transform of the pulse signal to obtain a filtered discrete fourier transform y (k):

Y(k)＝X(k)H(k),k＝0,1,2…,N-1

(8) the method for extracting the pulse signal envelope y (n) specifically comprises the following steps:

step 8-1: performing inverse discrete Fourier transform on the filtered discrete Fourier transform Y (k) to obtain a pulse signal filtered time-domain complex signal y_c(n)：

Step 8-2: taking a complex number y_c(n) obtaining the initial envelope data sequence y of the pulse signal by the modulus value₀(n)：

y₀(n)＝|y_c(n)|,n＝0,1…,N-1

Step 8-3: initializing the iterative smoothing times q of the envelope data to be 0;

step 8-4: let the iteration number q be q +1 and smooth the result y for the last iteration of the envelope data sequence_q-1(n) smoothing to obtain current smoothing result y_q(n)

And 8-5: judging whether the maximum number of iterative smoothing times of the envelope is reached, namely judging that q is less than or equal to I₂If yes, entering the step 8-6, otherwise jumping to the step 8-8;

and 8-6: respectively calculating last smoothing result y of envelope data sequence_q-1(n) and current smoothing result y_q(n) and the initial envelope data sequence y₀Residual sum of squares Ψ of (n)_q-1And Ψ_qRespectively is as follows:

and 8-7: determine | Ψ_q-1-Ψ_q|≤₂Ψ_qIf yes, entering the step 8-8, otherwise returning to the step 8-4;

and 8-8: let y (n) be y_q(N), N is 0,1,2 …, N-1, resulting in a sequence of signal envelope data:

y(n),n＝0,1,2…,N-1。

wherein, the step 8-1 can be realized by fast Fourier transform, and the calculation efficiency of the algorithm is improved.

In the embodiment of the invention, the simulation rectangular envelope pulse signal model is as follows:

where a is the amplitude of the signal and,

to an initial phase, τ₀For the pulse signal arrival time, τ is the pulse width, T is the received signal duration, f₁Is the signal starting frequency, f₂For signal termination frequency, when f₂＝f₁Time is a narrow band signal, otherwise, is a wide band signal, w (t) is a mean value of 0, and variance is σ²White gaussian noise, variance σ²Is determined by the signal-to-noise ratio SNR: SNR of 10log₁₀[A²/(2σ²)]。

At a sampling frequency f_sThe pulse signal is subjected to discrete sampling to obtain a pulse signal sampling data sequence:

wherein n is₀＝int(τ₀f_s)，M＝int(τf_s),N＝int(Tf_s)。

Example 1:

the simulation signal parameters are respectively set as: signal amplitude a 1, initial phase

Pulse width tau is 0.512s, received signal duration T is 1.024s, and pulse signal arrival time tau₀0.256s, signal start frequency f₁400Hz, signal termination frequency f₂400Hz, i.e. the artificial signal is a narrow-band pulse signal, the sampling frequency f_s2000Hz, 2048 observation data sequence points N, and 3dB SNR.

According to the step (2), calculating the power spectrum P (k) of the data sequence x (n), as shown in FIG. 2.

In the step (3), setting a power spectrum iteration smoothing maximum iteration number threshold I₁10, precision control threshold₁0.01, iterative smoothing window length M₁Envelope iteration smoothing maximum iteration number threshold I of 5₂5, the precision control threshold₂0.02, iterative smoothing window length M₂Narrow band signal bandwidth factor ω 11₁0.05, bandwidth factor ω of the wideband signal₂And (4) judging that the main lobe characteristic parameter threshold η of the narrow-band pulse and the wide-band pulse is 2 when the main lobe characteristic parameter threshold is 0.1.

According to the step (4), performing iterative smoothing processing on p (k) to obtain a smoothed power spectrum s (k), as shown in fig. 3, it can be seen from comparing fig. 2 and fig. 3 that: after iterative smoothing processing, the main lobe of the power spectrum of the narrow-band pulse signal is obviously widened, and the whole frequency band spectrum is smoother.

According to the step (5), extracting the original power spectrum and the main lobe characteristic parameter gamma of the smoothed power spectrum_pAnd gamma_SAre respectively as

γ_p＝45.6772，γ_S＝1.4108

The adaptive filter parameters h (k) generated according to step (6) are shown in fig. 4.

According to the steps (7) and (8), the extracted pulse signal envelope is shown in fig. 5, and as can be seen from fig. 5, the front and back edge features of the extracted narrow-band pulse signal envelope are obvious, so that a basis can be provided for further extracting the arrival time and the end time of the narrow-band pulse signal.

Example 2:

the simulation signal parameters are respectively set as: signal amplitude a 2, initial phase

Pulse width τ is 0.5s, received signal duration T is 1s, and pulse signal arrival time τ is₀0.25s, signal start frequency f₁350Hz, signal termination frequency f₂450Hz, i.e. the simulated signal is a broadband pulse signal, the sampling frequency f_s2048Hz, 2048 points N of the observation data sequence, and 6dB SNR of the signal-to-noise ratio.

According to the step (2), calculating the power spectrum P (k) of the data sequence x (n), as shown in FIG. 6.

In the step (3), setting a power spectrum iteration smoothing maximum iteration number threshold I₁10, precision control threshold₁0.01, iterative smoothing window length M₁7, envelope iteration smoothing maximum iteration number threshold I₂10, precision control threshold₂0.02, iterative smoothing window length M₂13, narrow-band signal bandwidth factor ω₁0.05, bandwidth factor ω of the wideband signal₂And (4) judging that the main lobe characteristic parameter threshold η of the narrow-band pulse and the wide-band pulse is 2 when the main lobe characteristic parameter threshold is 0.1.

According to the step (4), performing iterative smoothing processing on p (k) to obtain a smoothed power spectrum s (k), as shown in fig. 7, comparing fig. 6 and fig. 7, it can be seen that: after iterative smoothing processing, the main lobe change of the power spectrum of the broadband pulse signal is not obvious, and the whole frequency band spectrum is smoother.

According to the step (5), extracting the original power spectrum and the main lobe characteristic parameter gamma of the smoothed power spectrum_pAnd gamma_SRespectively as follows:

γ_p＝1.1081，γ_S＝1.1974

the adaptive filter parameters h (k) generated according to step (6) are shown in fig. 8.

According to the steps (7) and (8), the extracted pulse signal envelope is shown in fig. 9, and as can be seen from fig. 9, the front and rear edge features of the extracted wideband pulse signal envelope are obvious, so that a basis can be provided for further extracting the arrival time and the end time of the wideband pulse signal.

Claims (8)

1. A method for extracting a self-adaptive filtering envelope of a pulse signal frequency domain is characterized by comprising the following steps:

(1) acquiring a pulse signal sampling data sequence x (N) to be processed, wherein N is 0,1, … and N-1, and N is the number of sampling points corresponding to the pulse width length of a detected pulse signal and is an integer power of 2;

(2) calculating an original power spectrum P (k) of the pulse signal according to the data sequence x (n);

(3) carrying out parameter initialization on iterative smoothing and filter bandwidth;

(4) carrying out iterative smoothing processing on the original power spectrum P (k) to obtain a smoothed power spectrum S (k);

(5) extracting main lobe characteristic parameters of the original power spectrum P (k) and the smooth power spectrum S (k) respectively;

(6) generating adaptive filtering parameters h (k);

(7) performing adaptive filtering on the discrete Fourier transform result X (k) of the pulse signal according to the adaptive filtering parameter H (k);

(8) extracting a pulse signal envelope y (n);

in the step (5), the main lobe characteristic parameters of the original power spectrum p (k) and the smoothed power spectrum s (k) are respectively extracted, which specifically includes the following steps:

step 5-1: searching discrete frequency index k corresponding to maximum value of smoothed power spectrum S (k)_M：

Wherein the content of the first and second substances,

representing that the discrete frequency index corresponding to the maximum value of S (k) is searched within the range of 1 ≦ k ≦ N/2-1;

step 5-2: subtracting half of the maximum value from s (k) to obtain a half-power difference power spectrum z (k):

Z(k)＝S(k)-S(k_M)/2,k＝0,1,2…,N/2-1

step 5-3: searching for the discrete frequency index k of the first zero crossing on the left of the Z (k) peak_LThe searching process comprises the following steps:

step 5-3-1: initializing a frequency start search index k of a first zero crossing on the left side of a Z (k) peak value_l＝k_M-1；

Step 5-3-2: judgment of k_l0 or Z (k)_l) If the result is less than 0, jumping to the step 5-3-4 if the result is true, and otherwise, entering the step 5-3-3;

step 5-3-3: let k_l＝k_l1 and returning to step 5-3-2;

step 5-3-4: let k_L＝k_lObtaining the discrete frequency index k of the first zero-crossing point on the left side of the Z (k) peak value_L；

Step 5-4: searching for the first zero-crossing discrete frequency index k to the right of the Z (k) peak_RThe searching process comprises the following steps:

step 5-4-1: initializing a frequency start search index k of a first zero crossing on the right side of a Z (k) peak value_r＝k_M+1；

Step 5-4-2: judgment of k_rN/2-1 or Z (k)_r) If the value is less than 0, jumping to the step 5-4-4 if the value is less than 0, otherwise, entering the step 5-4-3;

step 5-4-3: let k_r＝k_r+1 and returning to step 5-4-2;

step 5-4-4: let k_R＝k_rObtaining the discrete frequency index k of the first zero-crossing point on the right side of the Z (k) peak value_R；

Step 5-5: calculating the primary power spectrum P (k) main lobe characteristic parameter gamma_p

And 5-6: calculating a smoothed power spectrum S (k) main lobe characteristic parameter gamma_S：

2. The method according to claim 1, wherein in step (1), the processed pulse signal sample data sequence is obtained by the following method: receiving real-time acquisition data of N sampling points from a sensor as a data sequence x (N) to be processed; or extracting the N sampling point data containing the whole pulse signal from the memory as the data sequence x (N) to be processed.

3. The method for extracting adaptive filtering envelopes in frequency domain of impulse signal as claimed in claim 1, wherein in the step (2), the original power spectrum p (k) of the impulse signal is calculated for the data sequence x (n) by the following method: performing fast Fourier transform on the data sequence x (n), and calculating to obtain discrete Fourier transform X (k) and an original power spectrum P (k) of the data sequence, wherein the method specifically comprises the following steps:

step 2-1: calculating the discrete Fourier transform of x (n):

where k is the discrete frequency index of X (k) and j represents the unit of an imaginary number, i.e.

This formula is implemented by fast fourier transform;

step 2-2: calculating a power spectrum for x (n) according to x (k):

wherein k is the discrete frequency index of the original power spectrum p (k), | represents the modulo value operation, and | x (k) | is the signal amplitude spectrum.

4. The method for extracting the adaptive filtering envelope in the frequency domain of the impulse signal as claimed in claim 1, wherein in the step (3), the parameters of the iterative smoothing and the filter bandwidth are initialized by the following method: setting a maximum iteration number threshold of power spectrum iterative smoothing, wherein the precision control threshold and the iterative smoothing window length are respectively I₁，₁And M₁The maximum iteration number threshold, the precision control threshold and the iteration smoothing window length of the initial envelope iteration smoothing are respectively I₂，₂And M₂Bandwidth coefficient omega of narrow-band signal₁And bandwidth factor omega of broadband signal₂Judging the characteristic parameter threshold η of the main lobe of the narrow-band pulse and the wide-band pulse, wherein I₁And I₂A value is a positive integer greater than 2,₁and₂the value is positive number less than 1, the length of the iterative smoothing window of the power spectrum satisfies 3 and is less than or equal to M₁An odd number less than or equal to N/2-2, and an envelope iteration smoothing window length of 5 less than or equal to M₂Odd number, omega, not more than 49₁To satisfy 0 ≤ ω₁Any number, omega, not more than 0.05₂To satisfy 0.05 < omega₂Any number not less than 0, η is any number satisfying 1.5 not less than η not more than 2.5.

5. The method for extracting adaptive filtering envelopes of a pulse signal frequency domain according to claim 1, wherein in the step (4), an iterative smoothing process is performed on an original power spectrum p (k) by using a method as follows to obtain a smoothed power spectrum s (k), and the method specifically includes the following steps:

step 4-1: the iterative smoothing times i of the initialized power spectrum is 0, and the iterative smoothing result S of the power spectrum₀(k) Comprises the following steps:

S₀(k)＝P(k),k＝0,1,2…,N/2-1

step 4-2: let the iteration number i be i +1, and smooth the last iteration of the power spectrum to obtain a result S_i-1(k) Performing smoothing to obtain current smoothing result S_i(k)：

Wherein the content of the first and second substances,M₁iteratively smoothing the window length for the power spectrum;

step 4-3: judging whether the maximum iterative smoothing times of the power spectrum is reached, namely judging that I is less than or equal to I₁If yes, entering the step 4-4, otherwise jumping to the step 4-6; wherein I₁A maximum iteration number threshold for power spectrum iterative smoothing;

step 4-4: respectively calculating last smoothing results S of power spectrum_i-1(k) And current smoothing result S_i(k) Sum of squares J of residuals from original power spectrum p (k)_i-1And J_iRespectively is as follows:

and 4-5: judgment of | J_i-1-J_i|≤₁J_iIf yes, entering the step 4-6, otherwise returning to the step 4-2; wherein₁A threshold is controlled for power spectrum iterative smoothing precision;

and 4-6: let S (k) become S_i(k) And k is 0,1,2 …, N/2-1, and the smoothed power spectrum s (k) is obtained.

6. The method according to claim 1, wherein in the step (6), the adaptive filtering parameters h (k) are generated by using the following method, k being 0,1,2 …, N-1, and specifically comprising the following steps:

step 6-1: initializing adaptive filtering parameters:

H(k)＝0,k＝0,1,2…,N-1

step 6-2: respectively calculating the discrete frequency indexes T of the start and the end of the transition band of the narrow-band pulse signal filter_N1And T_N2And a passband start and end discrete frequency index P_N1And P_N2：

T_N1＝max(round[(1-0.5ω₁)k_M],0)，T_N2＝min(round[(1+0.5ω₁)k_M],N/2-1)

P_N1＝max(k_M-1,T_N1)，P_N2＝min(k_M+1,T_N2)

Wherein, ω is₁Is a narrow band signal bandwidth coefficient, round [ ]]Representing a rounding operation, max (,) and min (,) representing a large and small value operation, respectively, k_MRepresenting the discrete frequency index corresponding to the maximum value of the smoothed power spectrum S (k);

step 6-3: respectively calculating the start and end discrete frequency indexes T of the transition band of the broadband pulse signal filter_W1And T_W2And a passband start and end discrete frequency index P_W1And P_W2：

T_W1＝max(round[(1-0.5ω₂)k_M],0)，T_W2＝min(round[(1+0.5ω₂)k_M],N/2-1)

P_W1＝max(k_L-1,T_W1)，P_W2＝max(k_R+1,T_W2)

Wherein, ω is₂For the bandwidth factor, k, of the broadband signal_LDiscrete frequency index for the first zero crossing point on the left side of peak value of half-power difference power spectrum Z (k)_RIndexing the first zero-crossing discrete frequency on the right side of a half-power difference power spectrum Z (k) peak value;

step 6-4: obtaining discrete frequency index T from adaptive filter transition band start and end according to signal type₁And T₂And a passband start and end discrete frequency index P₁And P₂：

Namely, max (gamma)_p,γ_s)/min(γ_p,γ_s) If > η is true, if true, the signal is judged to be a narrow-band pulse signal and the following steps are carried out:

T₁＝T_N1,T₂＝T_N2,P₁＝P_N1,P₂＝P_N2

otherwise, judging the signal to be a broadband pulse signal and ordering:

T₁＝T_W1,T₂＝T_W2,P₁＝P_W1,P₂＝P_W2

wherein, γ_pIs the main lobe characteristic parameter, gamma, of the original power spectrum P (k)_STo smooth the power spectrum S: (k) η is a main lobe characteristic parameter threshold for judging the narrow-band pulse and the wide-band pulse;

step 6-5: modifying the transition band part of the adaptive filter parameters H (k) so that k satisfies T₁≤k≤T₂The part of (A) is:

step 6-6: modifying the passband portion of the adaptive filter parameters H (k) such that k satisfies P₁≤k≤P₂The part of (A) is: h (k) ═ 1.

7. The method for extracting adaptive filtering envelopes in frequency domain of impulse signals as claimed in claim 1, wherein in the step (7), the discrete fourier transform result x (k) of the impulse signals is adaptively filtered by the following method to obtain the filtered discrete fourier transform y (k):

Y(k)＝X(k)H(k),k＝0,1,2…,N-1

where h (k) is an adaptive filtering parameter.

8. The method for extracting the envelope of the pulse signal by frequency domain adaptive filtering according to claim 7, wherein in the step (8), the envelope y (N) of the pulse signal is extracted by the following method, where N is 0,1 …, and N-1, and the method specifically includes the following steps:

step 8-1: performing inverse discrete Fourier transform on the filtered discrete Fourier transform Y (k) to obtain a pulse signal filtered time-domain complex signal y_c(n)：

Step 8-2: taking a complex number y_c(n) obtaining the initial envelope data sequence y of the pulse signal by the modulus value₀(n)：

y₀(n)＝|y_c(n)|,n＝0,1…,N-1

Step 8-3: initializing the iterative smoothing times q of the envelope data to be 0;

step 8-4: let the iteration number q be q +1 and smooth the result y for the last iteration of the envelope data sequence_q-1(n) smoothing to obtain current smoothing result y_q(n)：

Wherein M is₂Iterating the smoothing window length for the initial envelope;

and 8-5: judging whether the maximum iterative smoothing times of the envelope is reached, namely judging that q is less than or equal to I₂If yes, entering the step 8-6, otherwise jumping to the step 8-8; wherein, I₂A maximum iteration number threshold for initial envelope iteration smoothing;

and 8-6: respectively calculating last smoothing result y of envelope data sequence_q-1(n) and current smoothing result y_q(n) and the initial envelope data sequence y₀Residual sum of squares Ψ of (n)_q-1And Ψ_qRespectively is as follows:

and 8-7: determine | Ψ_q-1-Ψ_q|≤₂Ψ_qIf yes, entering the step 8-8, otherwise returning to the step 8-4;

wherein the content of the first and second substances,₂a threshold is controlled for the initial envelope iteration smoothing precision;

and 8-8: let y (n) be y_q(N), N is 0,1,2 …, N-1, resulting in a sequence of signal envelope data:

y(n),n＝0,1,2…,N-1。

Images