Pulse interference detection and suppression method for vector sensor based on trimmed median filter

Abstract

The present invention discloses a method for detecting and suppressing pulse interference of a vector sensor based on trimmed median filtering, designs a comparison template, calculates the correlation between the comparison template and the signal in the current sliding area by sliding the comparison template in a search area near the pulse interference, and selects the wavelet coefficients of the most similar part near the pulse interference as replacement wavelet coefficients of the pulse interference. The method of the present invention utilizes discrete wavelet transform to decompose the signal at multiple scales, ensures the calculation efficiency, and optimizes the threshold function of pulse interference detection based on trimmed median filtering, reduces the influence of strong pulses on threshold selection, and improves the detection rate of weak pulses.

Description

Technical Field

[0001] The invention belongs to the technical field of sensors, and in particular relates to a vector sensor pulse interference detection and suppression method based on trimmed median filtering. Background Art

[0002] As an advanced underwater sound signal detection device, the vector sensor can simultaneously and accurately measure the sound pressure and particle velocity of a specific point in the sound field. It has a stable dipole pointing characteristic and is not affected by frequency. The co-oscillation vector sensor indirectly obtains the particle velocity by measuring the fluid particle acceleration. Due to its good low-frequency characteristics, it has been widely used in the field of underwater acoustic measurement and detection. However, in actual applications, the co-oscillation vector sensor will produce a series of vibrations under the instigation of the water flow and the obstruction of the cable, resulting in instantaneous and strong pulse interference. Therefore, detecting and suppressing the pulse interference mixed in the output signal of the vector sensor is an important part of improving its receiving performance.

[0003] Short-time Fourier transform (STFT) in the time-frequency domain is a commonly used signal analysis method. This method obtains the frequency domain characteristics of the local signal by selecting a window function, and then slides the window function along the time axis, so that the signal spectrum at each local moment is presented to achieve joint analysis of the signal in the time domain and frequency domain. The length of the short-time Fourier transform window determines the time-frequency resolution. When processing short pulse signals, too high requirements for time resolution will lead to a decrease in frequency domain resolution. Wavelet analysis can change the size of the analysis "window" according to the decomposition scale, so that the wavelet analysis has a larger time resolution at a small scale and a larger frequency resolution at a large scale, which can better capture the characteristics of the signal. Continuous wavelet transform (CWT) can perform global analysis of the signal, provide a complete time-frequency diagram, and observe the changes of the signal at any scale. However, CWT needs to continuously calculate multiple scale transforms in the entire signal domain, which has high computational complexity and too much redundant information. Discrete wavelet transform (DWT) uses a filter array to construct multi-resolution time-frequency to perform the Mallat algorithm commonly used in wavelet decomposition. DWT reduces the amount of calculation by sampling and retains most of the important features of the signal. At the same time, due to the convenience of discrete wavelet transform in convolution, the Mallat algorithm can be easily implemented using digital integrated circuits, so DWT is widely used in practice.

[0004] Accurately detecting pulse interference is the prerequisite for suppressing interference. DWT decomposes the signal at different scales to obtain wavelet coefficients of each scale, which can be used to accurately reconstruct the original signal. The key to detecting pulse interference is to select a suitable threshold function to filter out the wavelet coefficients corresponding to the pulse interference for subsequent suppression processing. It is a commonly used pulse detection method to perform wavelet transform on singular signals such as pulses and perform median filtering on the wavelet coefficients to determine the pulse detection threshold and then detect the location of the pulse. However, when a weak pulse appears near a strong pulse, this method is prone to miss detection of the weak pulse. Therefore, choosing a suitable threshold function is crucial to simultaneously detect strong pulse interference and nearby weak pulses.

[0005] After completing the detection of pulse interference, the wavelet coefficients corresponding to the pulse interference are processed. The common processing methods are hard and soft thresholds and their improved threshold algorithms, that is, adjusting the wavelet coefficients corresponding to the pulse interference to the specified range of the set threshold function. This method is simple and efficient, but it destroys the statistical properties of the original signal and distorts it. Therefore, the focus of the research is to appropriately adjust the wavelet coefficients corresponding to the pulse interference to ensure that the original statistical properties of the signal are not destroyed while suppressing the pulse interference. Summary of the invention

[0006] In order to overcome the shortcomings of the prior art, the present invention provides a vector sensor pulse interference detection and suppression method based on trimmed median filtering, designs a comparison template, calculates the correlation between the comparison template and the signal in the current sliding area by sliding the comparison template in the search area near the pulse interference, and selects the wavelet coefficients of the most similar part near the pulse interference as the replacement wavelet coefficients of the pulse interference. The method of the present invention utilizes discrete wavelet transform to decompose the signal at multiple scales, ensures the calculation efficiency, and optimizes the threshold function of pulse interference detection based on trimmed median filtering, reduces the influence of strong pulses on threshold selection, and improves the detection rate of weak pulses.

[0007] The technical solution adopted by the present invention to solve the technical problem is as follows:

[0008] Step 1: Pulse interference detection;

[0009] Step 1-1: Discrete wavelet transform;

[0010] Assume that the one-dimensional continuous time signal is f(t)∈L 2 (R),L 2 (R) is the function space of square integrable functions on the real number set R, that is:

[0011]

[0012] In the formula, t represents a continuous time variable, and the continuous wavelet transform formula corresponding to f(t) is:

[0013]

[0014] In the formula, ψ a,b (t) is obtained by scaling and translating the mother wavelet ψ(t), which is called a wavelet sequence; a,b∈R and a≠0, a is the scale factor, b is the translation factor; ψ * (.) is the complex conjugate of ψ(.);

[0015] The scale factor a and translation factor b of the continuous wavelet are transformed along the binary sequence {2 j}, j∈Z, and perform discretization sampling to obtain discrete wavelet transform:

[0016]

[0017] Among them, j represents the decomposition scale, n represents the discrete time point, and ψ j,n (t) represents the wavelet sequence after discretization sampling of scale factor and translation factor;

[0018] The Mallat algorithm is given by formula (4):

[0019]

[0020] Where g(m) and h(m) are the responses of the low-pass filter and high-pass filter, respectively, J represents the total number of scales of wavelet decomposition, m represents the index of the filter, M is the length of the filter, and cA is j (.) and cD j (.) represent the approximate coefficient and detail coefficient of discrete wavelet decomposition at the jth scale; cA j+1 (.) and cD j+1 (.) represent the approximate coefficient and detail coefficient of discrete wavelet decomposition at the j+1th scale respectively;

[0021] Step 1-2: trimmed median filtering;

[0022] A threshold function is set to detect pulse interference. The threshold function is proportional to the smooth envelope of the absolute value of the wavelet coefficients, but is not affected by the pulse interference. The threshold function is defined as:

[0023] Γ(j,n)=r f ·Env[DWT f (j,n)] (5)

[0024] Where Env[·] is the smooth envelope of the signal not affected by pulse interference, r f is an adjustable factor set empirically;

[0025] If the absolute value of the wavelet coefficient at a certain position is greater than the set threshold, then this position is determined to be a pulse interference, and the position index of all pulse interferences is obtained, recorded as η(j,n):

[0026]

[0027] Threshold function based on trimmed median filter for:

[0028]

[0029] Where α∈(0,1) is an adjustable factor; MED[.] represents median filtering of the absolute value of the wavelet coefficients, and l is the number of sample points of a single pulse interference;

[0030] The threshold function calculated by trimmed median filtering is used to detect pulse interference on the wavelet coefficients of each scale after DWT decomposition, and the detection results η(j,n) of each decomposition scale are obtained. η(j,n)=1 indicates that pulse interference is detected, and η(j,n)=0 indicates that pulse interference is not detected.

[0031] Step 1-3: Jointly determine pulse interference at multiple scales;

[0032] Perform differential operation on the detected pulse interference index:

[0033] q(j,n)=η(j,n)-η(j,n-1) (8)

[0034] Among them, q(j,n) is the result of the difference operation, and the index of the adjacent q(j,n)=1 and q(j,n)=-1 is recorded as the starting point n of the i-th pulse interference of the j-th scale. s (j,i) and end point n e (j,i); if the length of the j-th scale wavelet coefficient is recorded as L j c , then the length of the j+1th scale wavelet coefficient is recorded as

[0035] The wavelet coefficients of each scale are linearly interpolated as shown in formula (9) so that the length of the wavelet coefficients of all scales is the same as the length of the wavelet coefficients of the j=1 scale; the detection results after interpolation of each scale are recorded as: The starting and ending point indices are recorded as:

[0036]

[0037] Combine the detection results of all scales and perform voting statistics, recorded as V(n), then:

[0038]

[0039] Set the threshold T. If the voting result V(n) ≥ T, it is determined that pulse interference occurs at discrete time point n. Set T = J / 2 to indicate that at least half of the scales detect that there is a pulse here, which is effective pulse interference. The final detection result of the comprehensive judgment of each scale is:

[0040]

[0041] The start and end indexes of the pulse interference in the comprehensive judgment result are: The comprehensive judgment result is downsampled so that the length of the wavelet coefficients at each scale is the same as the length of the initial wavelet coefficients at each scale; the detection result of each scale after updating is recorded as η(j,n), and the relationship between it and the detection result η(n) of the comprehensive judgment of each scale is as follows:

[0042]

[0043] The starting and ending indexes of the pulse are recorded as: and

[0044] Step 2: Pulse interference suppression;

[0045] Step 2-1: pulse interference wavelet coefficient processing;

[0046] Wavelet coefficients DWT of the original signal f Subtract the wavelet coefficients at the pulse interference from (j,n) to obtain the wavelet coefficients DWT without pulse interference components x (j,n), that is:

[0047]

[0048] Calculate the signal envelope after removing the pulse interference, denoted as ξ x (j,n):

[0049]

[0050] in:

[0051]

[0052] Similarly, calculate the signal envelope of the pulse interference, denoted as ξ η (j,n);

[0053] The following processing is performed on the wavelet coefficients of the pulse interference and the pulse interference vicinity to obtain the processing template Ξ(j,n) of the wavelet coefficients in the pulse interference region:

[0054]

[0055] Multiply the processing template Ξ(j,n) by the original signal wavelet coefficient to achieve zero wavelet coefficients at the pulse interference and reset the wavelet coefficients near the pulse interference in proportion:

[0056]

[0057] Step 2-2: Replacement of pulse interference coefficient;

[0058] The starting point and end point of the pulse interference below no longer distinguish the specific decomposition scale j, and j is deleted from the symbol. The starting point of the i-th pulse interference at the j-th scale is recorded as n s (j,i) is End point e (j,i) is

[0059] (a) Design a comparison template for the current pulse interference. The starting range of the comparison template is the same as the length of the processing template. L2 represents the distance between the starting point and the end point of the pulse interference.

[0060] Used to define the scope of comparison templates and processing templates;

[0061] (b) Define a search area. The starting range of the search area is L1>L2 and the processing template slides in the search area, and the energy in the area is calculated after multiplying it with the current wavelet coefficient; L1 represents the distance between the starting point and the end point of the pulse interference and extends to both sides, which is used to define the range of the search area;

[0062] Through (a) and (b), the energy value in the current pulse interference area calculated by the comparison template is obtained, which is recorded as is a constant; use the processing template to slide in the search area and calculate the energy value in the current template, recorded as Compare and Similarity degree, the similarity calculation formula is:

[0063]

[0064] in:

[0065]

[0066] In the formula is the energy of the current search window; To compare the energy of the template, for the i-th pulse interference is a constant; n1 is the kth point of the sliding window, n2 is the kth point of the comparison template where the pulse is located, K i is the number of sample points of the comparison template;

[0067] After the above similarity calculation, the index of the start of the search window with the maximum similarity is obtained Then the index of the start position to be replaced for impulse interference and is:

[0068]

[0069] wherein represents the start point of the optimal search window;

[0070] If L3 sample points are extended on both sides of the replacement part and L3 < L2, the new area of the replacement part is

[0071]

[0072] For the processing template in Equation (16), find the minimum value, and record the index of the first point as The index of the last point is Then the replacement smoothing template at this time is:

[0073]

[0074] In the formula

[0075]

[0076] The wavelet coefficients of each scale after final impulse interference suppression are:

[0077]

[0078] Preferably, the Env[·] is implemented using a median filter.

[0079] Preferably, the calculation definition of the median filter is:

[0080] Γ(j,n) = r f MED[DWT f (j,n - l), …, DWT f (j,n), …, DWT f (j,n + l)];

[0081] The length L of the median filter med needs to satisfy:

[0082] L med > 2l + 1.

[0083] A computer program that causes a computer to execute the above sensor pulse self-noise detection and suppression method.

[0084] An electronic device comprises: a processor and a memory; the memory is used to store a computer program, and the processor is used to execute the computer program stored in the memory, so that the electronic device executes the above-mentioned sensor pulse self-noise detection and suppression method.

[0085] A computer-readable storage medium stores a computer program, which, when executed by a processor, implements the above-mentioned sensor pulse self-noise detection and suppression method.

[0086] A chip includes: a processor for calling and running a computer program from a memory, so that a device equipped with the chip executes the above-mentioned sensor pulse self-noise detection and suppression method.

[0087] A computer program product, the computer program product comprising a computer storage medium, the computer storage medium storing a computer program, the computer program comprising instructions executable by at least one processor, and the above-mentioned sensor pulse self-noise detection and suppression method is implemented when the instructions are executed by the at least one processor.

[0088] The beneficial effects of the present invention are as follows:

[0089] The method of the present invention utilizes discrete wavelet transform to decompose the signal at multiple scales, thereby ensuring computational efficiency, and optimizes the threshold function of pulse interference detection based on trimmed median filtering, thereby reducing the influence of strong pulses on threshold selection and improving the detection rate of weak pulses. BRIEF DESCRIPTION OF THE DRAWINGS

[0090] Figure 1 Schematic diagram of Mallat algorithm.

[0091] Figure 2 A diagram showing the relationship between processing templates, replacement templates, and search areas.

[0092] FIG3( a ) is a time domain waveform diagram of a Y channel signal according to an embodiment of the present invention.

[0093] FIG3( b ) is a time-frequency diagram of the Y channel signal according to an embodiment of the present invention.

[0094] FIG. 4( a ) is a graph showing the median filter cA4 threshold curve and detection results in an embodiment of the present invention.

[0095] FIG. 4( b ) is a diagram showing a trimmed median filter cA4 threshold curve and a detection result diagram according to an embodiment of the present invention.

[0096] FIG. 5( a ) is a diagram showing a trimmed median filter cA4 threshold curve and a detection result diagram according to an embodiment of the present invention.

[0097] FIG. 5( b ) is a diagram showing a trimmed median filter cD4 threshold curve and a detection result diagram according to an embodiment of the present invention.

[0098] FIG. 5( c ) is a diagram showing a trimmed median filter cD3 threshold curve and a detection result diagram according to an embodiment of the present invention.

[0099] FIG. 5( d ) is a diagram showing a trimmed median filter cD2 threshold curve and a detection result diagram according to an embodiment of the present invention.

[0100] Figure 6 Detection results at each scale and multi-scale superposition results of the embodiment of the present invention.

[0101] Figure 7 A diagram of multi-scale comprehensive decision detection results according to an embodiment of the present invention.

[0102] Figure 8 The signal envelope and pulse interference envelope of the embodiment of the present invention.

[0103] Fig. 9 Schematic diagram of a processing template according to an embodiment of the present invention.

[0104] Fig.10 Schematic diagram of replacing a template according to an embodiment of the present invention.

[0105] Fig.11 The time domain diagram of the original signal in the embodiment of the present invention and the time domain diagram after median filtering and processing by the method.

[0106] Fig.12 The time-frequency diagram of the original signal in the embodiment of the present invention and the time-frequency diagram after being processed by the method. DETAILED DESCRIPTION

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

[0108] In order to solve the problem of instantaneous and strong pulse interference caused by the sudden vibration of the vector sensor under the instigation of water flow and the obstruction of the cable, an improved median filtering method - trimmed median filtering is used to detect strong and weak pulses, improve the detection rate of weak pulses, and solve the problem of missed detection of weak pulses near strong pulses.

[0109] In order to improve the signal quality after pulse interference suppression, the present invention designs a comparison template, which calculates the correlation between the comparison template and the signal in the current sliding area by sliding the comparison template in the search area near the pulse interference, and selects the wavelet coefficients of the most similar part near the pulse interference as the replacement wavelet coefficients of the pulse interference.

[0110] The technical solution of the present invention is composed of two processing modules, namely a pulse interference detection module and a pulse interference suppression module.

[0111] 1) Pulse interference detection;

[0112] Pulse interference detection can be divided into the following two steps:

[0113] 1.1 Discrete Wavelet Transform;

[0114] Assume a one-dimensional continuous time signal f(t)∈L 2 (R),L 2 (R) The function space of square integrable functions on the real number set R, that is:

[0115]

[0116] The corresponding continuous wavelet transform formula is:

[0117]

[0118] The scale factor a and translation factor b of the continuous wavelet are transformed along the binary sequence {2 j}, j∈Z, and perform discretization sampling to obtain discrete wavelet transform:

[0119]

[0120] It can be seen that the discrete wavelet transform performs a secondary extraction of the translation factor b. The Mallat algorithm is given by formula (4):

[0121]

[0122] Mallat algorithm Figure 1 shown.

[0123] 1.2 Pruning median filtering;

[0124] A threshold function is set to detect pulse interference. The threshold function is proportional to the smooth envelope of the absolute value of the wavelet coefficients, but is not affected by the pulse interference. The threshold function is defined as:

[0125] Γ(j,n)=r f ·Env[DWT f (j,n)]

[0126] If the absolute value of the wavelet coefficient at a certain position is greater than the set threshold, this position is considered to be pulse interference, and the position index of all pulses is obtained, recorded as η(n):

[0127]

[0128] Env[·] used to calculate the threshold function in equation (5) is usually implemented using a median filter, because the median filter retains the step function type signal and is robust to pulse interference. The traditional median filter calculation is defined as:

[0129] Γ(j,n)=r f MED[DWT f (j,nl),…,DWT f (j,n),…,DWT f (j,n+l)]

[0130] The length of the median filter is L med Need to meet:

[0131] L med >2l+1

[0132] Traditional median filtering can remove strong pulse interference with large amplitude. However, when there is a weak pulse near the strong pulse, this method is prone to miss the weak pulse. The present invention proposes a threshold selection method based on trimmed median filtering. First, a certain proportion of extreme value data in the wavelet coefficients is removed to reduce the influence of strong pulses on threshold selection. Then, median filtering is performed on the wavelet coefficients to determine the detection threshold, thereby improving the detection rate of weak pulses. Threshold function based on trimmed median filtering for:

[0133]

[0134] Where α∈(0,1) is an adjustable factor. It is used to remove the wavelet coefficients corresponding to the strong pulse interference of the tail α proportion after the median filter is sorted in ascending order, reducing the influence of strong pulses on the threshold selection, so that weak pulses can be detected and the detection rate of weak pulses can be improved. The threshold function calculated by the trimmed median filter is used to calculate the wavelet coefficients of each scale after DWT decomposition to realize pulse interference detection, and the detection results η(j,n) of each decomposition scale are obtained. η(j,n)=1 means that the pulse interference is detected, and η(j,n)=0 means that the pulse interference is not detected.

[0135] 1.3 Multi-scale joint determination of pulse interference;

[0136] In order to reduce the missed detection rate and effectively suppress weak pulses, the following joint processing is performed:

[0137] Perform differential operation on the detected pulse interference index:

[0138] q(j,n)=η(j,n)-η(j,n-1)

[0139] Due to the characteristics of the Mallat algorithm, the length of the wavelet coefficients is halved each time the decomposition is performed. If the length of the j-th scale wavelet coefficient is recorded as Then the length of the j+1th scale wavelet coefficient is recorded as

[0140] Since the lengths of wavelet coefficients at different scales are different, the wavelet coefficients at each scale are linearly interpolated as shown in formula (9) so that the lengths of wavelet coefficients at all scales are the same as the length of wavelet coefficients at scale j = 1. The detection results after interpolation at each scale are recorded as: The starting and ending point indices are recorded as:

[0141]

[0142] Combine the detection results of all scales and perform voting statistics, recorded as V(n), then:

[0143]

[0144] Set the threshold T. If the voting result V(n) ≥ T, it is determined that pulse interference occurs at discrete time point n. Set T = J / 2 to indicate that at least half of the scales detect the presence of pulses here to be considered effective pulse interference. Therefore, the final detection result of the comprehensive judgment of each scale is:

[0145]

[0146] The start and end indexes of the pulse interference in the comprehensive judgment result are: The comprehensive judgment result is downsampled so that the length of the wavelet coefficients at each scale is the same as the length of the initial wavelet coefficients at each scale. The updated detection result of each scale is recorded as η(j,n), and the relationship between it and the detection result η(n) of the comprehensive judgment of each scale is:

[0147]

[0148] The starting and ending indexes of the pulse are recorded as and

[0149] 2) Pulse interference suppression;

[0150] Pulse interference suppression can be divided into the following two steps:

[0151] 2.1 Pulse interference wavelet coefficient processing;

[0152] After completing the detection of pulse interference in 1), the detection results of pulse interference at various scales η(j,n) are obtained. In order to avoid simply setting the wavelet coefficients at the location of the pulse interference to zero, which will cause artifacts near the pulse interference (mutation signals or noise introduced due to improper processing), affecting the effect of pulse interference suppression, the wavelet coefficients near the pulse interference also need to be replaced. fSubtract the wavelet coefficients at the pulse interference from (j,n) to obtain the wavelet coefficients DWT without pulse interference components x (j,n), that is:

[0153]

[0154] Calculate the signal envelope after removing the pulse interference, denoted as ξ x (j,n):

[0155] ξ x (j,n)=min{φ1(j,n),φ2(j,n)}

[0156] in

[0157] φ1(j,n)=max{DWT x (j,nl),…,DWT x (j,n)}

[0158] φ2(j,n)=max{DWT x (j,n),…,DWT x (j,n+l)}

[0159] Similarly, calculate the signal envelope of the pulse interference, denoted as ξ η The following processing is performed on the wavelet coefficients of the pulse interference and the area near the pulse interference to obtain the processing template Ξ(j,n) of the wavelet coefficients in the pulse interference area:

[0160]

[0161] Multiply the processing template with the original signal wavelet coefficient to achieve zero wavelet coefficients at the pulse interference and reset the wavelet coefficients near the pulse interference in proportion:

[0162]

[0163] 2.2 Replacement of pulse interference coefficient;

[0164] In order to find the most suitable wavelet coefficients to replace the wavelet coefficients of the pulse interference suppressed in 2.1, the following processing steps are designed. For the sake of simplicity in notation, the starting point and end point of the pulse interference below no longer distinguish the specific decomposition scale j, and j is deleted in the notation, and the starting point of the i-th pulse interference at the j-th scale is recorded as n s (j,i) is End point e (j,i) is

[0165] (a) Design a comparison template for the current pulse interference. The starting range of the comparison template has the same length as the processing template above, both of which are This step aims to calculate the energy within the area after removing impulse interference, which serves as a reference for calculating the energy of the replacement area in the subsequent calculation.

[0166] (b) Define a search area. The starting range of the search area is L1 > L2, and the processing template slides within the search area. After multiplying with the current wavelet coefficients, calculate the energy within this area.

[0167] Through the above two steps, obtain the energy value within the current impulse interference area calculated by the comparison template, denoted as is a constant; use the processing template to slide within the search area and calculate the energy value within the current template, denoted as Compare and for the similarity degree. The similarity calculation formula is:

[0168]

[0169] where

[0170]

[0171] n1 = n + k

[0172]

[0173] After completing the above similarity calculation, obtain the index of the starting position of the search window with the maximum similarity Then the starting position index of the impulse interference to be replaced and are:

[0174]

[0175] However, in order to minimize artifacts as much as possible, the replaced part needs to be smoothly blended at the boundaries. Therefore, it is necessary to extend and overlap the replaced part on both sides to achieve smooth blending. It is stipulated that L3 sample points are extended on both sides of this part and L3 < L2, then the new area of the replaced part is

[0176]

[0177] To smoothly replace the wavelet coefficients in this impulse interference area, it is necessary to design a replacement smoothing template to achieve smooth replacement from to For the processing template in Equation (16), find the minimum value. The points where the minimum value is taken are usually not unique. Denote the index of the first point as The index of the last point is Then the replacement smoothing template at this time is:

[0178]

[0179] In the formula

[0180]

[0181] The relationship diagram of processing template, replacement template and search area is as follows Figure 2 As shown. Then the wavelet coefficients of each scale after the final pulse interference suppression are:

[0182]

[0183] Example:

[0184] In order to verify the effective performance of the present invention, the embodiments of the present invention verify the actual effect of the above-mentioned pulse interference detection and suppression method on actual data.

[0185] The data comes from the Y channel of a domestic isochronous vector sensor. The target signal in the data is a 275Hz CW signal, the original sampling rate of the signal f0 = 5000Hz, the signal duration is 102.4s, and the signal is downsampled to a sampling rate of 2500Hz. The duration of the pulse interference is 0.3s, and the number of sampling points corresponding to the time domain is l0 = 750.

[0186] The signal time domain waveform is shown in Figure 3(a), and the time-frequency diagram is shown in Figure 3(b). From the time-frequency diagram in Figure 3(b), it can be seen that the pulse interference continues from the low frequency band to the high frequency band. The signal is transformed by J=4 scale discrete wavelet transform. According to the discrete wavelet decomposition characteristics, the length of each scale pulse interference after decomposition is l j , l j =l0 / 2 j The signal frequency range corresponding to the wavelet coefficients decomposed at each scale is [0,f s / 2 j ], where the approximation coefficient cA j The corresponding frequency range is [0, f s / 2 j+1 ], detail coefficient cD j The corresponding frequency range is [f s / 2 j+1 ,f s / 2 j ]. The wavelet approximation coefficient cA is scaled with j = 4 jAs an example, the threshold functions of median filtering and pruned median filtering with a pruning coefficient of α=0.15 are calculated respectively, and the detection effects are shown in Figure 4(a) and Figure 4(b). The upper figure of Figure 4(a) is the threshold function calculated by median filtering, and the lower figure of Figure 4(a) is the result of its detection. The upper figure of Figure 4(b) is the threshold function calculated by pruned median filtering, and the lower figure of Figure 4(b) is the result of its detection. By comparison, it can be seen that the threshold function calculated by pruned median filtering can detect more weak pulse interference and improve the detection rate of weak pulse interference.

[0187] Figure 5(a), Figure 5(b), Figure 5(c), and Figure 5(d) respectively show the threshold function curves calculated by the trimmed median filtering method for cA4, cD4, cD3, and cD2 and the pulse interference detection results. Figure 6 The results of each scale detection and the superposition of all scale detections are shown. In the step of multi-scale joint determination of pulse interference, according to the scale of wavelet decomposition J=4, it can be known that the threshold of multi-scale joint is T=J / 2=2, that is, it is considered that voting statistics are performed on each position at four scales, and the statistical value reaches 2 times or more, which is considered to be pulse interference, and less than 2 times is considered a false alarm and is discarded. Figure 7 The figure above shows the result after superposition of all scale detections, and the figure below shows the pulse interference detection result after multi-scale joint determination.

[0188] Suppress and replace the detected pulse interference. Draw the signal envelope and the pulse interference envelope together Figure 8 In the above example, the processing template of (16) is designed according to the envelope of the signal and the pulse interference. Taking the first pulse as an example, the processing template Ξ(j,n) of the wavelet coefficients in the pulse interference region is plotted. Fig. 9 , the extension length of the processing template L2 = l j The length of the pulse search area on both sides is L1 = 3l j . Pulse interference replacement template such as Fig.10 As shown, the extension length of the replacement template is L3 = l j / 10. The signal time domain diagram and the original signal waveform diagram after the pulse interference is detected and suppressed by two threshold functions are shown in the figure below. Fig.11 As shown, Fig.11 The above figure is the time domain diagram of the original signal. Fig.11 The middle figure is the time domain figure after the threshold function detection and suppression calculated by the median filter. Fig.11 The figure below is a time domain diagram of pulse interference detection and suppression by this method. Compared with the median filter method, this method suppresses 10 additional weak pulses when the pruning coefficient α=0.15, which are marked with dotted arrows. This proves that this method reduces the missed detection rate and effectively suppresses weak pulses. The time-frequency diagram of the signal before and after processing by this method is compared. Fig.12 As shown, Fig.12 The left picture is the time-frequency diagram of the original signal. Fig.12 The right figure is the time-frequency diagram after pulse suppression using this method.

[0189] According to the embodiment, it can be considered that the pulse interference detection and suppression method proposed in the present invention is feasible. The discrete wavelet transform is used to reduce the amount of calculation, which can meet the requirements for pulse interference detection and optimize the threshold function of pulse interference detection based on the trimmed median filter, reduce the influence of strong pulses on threshold selection, and improve the detection rate of weak pulses. Actual data processing shows that under the condition that the trimming ratio is between 10% and 20%, the method of the present invention significantly improves the effect of pulse interference detection.

Claims

1. A vector sensor pulse interference detection and suppression method based on trimmed median filtering, characterized in that: The steps include: Step 1: Pulse interference detection; Step 1-1: Discrete wavelet transform; Assume that the one-dimensional continuous time signal is f(t)∈L 2 (R),L 2 (R) is the function space of square integrable functions on the real number set R, that is: In the formula, t represents a continuous time variable, and the continuous wavelet transform formula corresponding to f(t) is: In the formula, ψ a,b (t) is obtained by scaling and translating the mother wavelet ψ(t), which is called a wavelet sequence; a,b∈R and a≠0, a is the scale factor, b is the translation factor; ψ * (.) is the complex conjugate of ψ(.); The scale factor a and translation factor b of the continuous wavelet are transformed along the binary sequence {2 j }, j∈Z, and perform discretization sampling to obtain discrete wavelet transform: Among them, j represents the decomposition scale, n represents the discrete time point, and ψ j,n (t) represents the wavelet sequence after discretization sampling of scale factor and translation factor; The Mallat algorithm is given by formula (4): Where g(m) and h(m) are the responses of the low-pass filter and high-pass filter, respectively, J represents the total number of scales of wavelet decomposition, m represents the index of the filter, M is the length of the filter, and cA is j (.) and cD j (.) represent the approximate coefficient and detail coefficient of discrete wavelet decomposition at the jth scale; cA j+1 (.) and cD j+1 (.) represent the approximate coefficient and detail coefficient of discrete wavelet decomposition at the j+1th scale respectively; Step 1-2: trimmed median filtering; A threshold function is set to detect pulse interference. The threshold function is proportional to the smooth envelope of the absolute value of the wavelet coefficients, but is not affected by the pulse interference. The threshold function is defined as: Γ(j,n)=r f ·Env[DWT f (j,n)] (5); Where Env[·] is the smooth envelope of the signal not affected by pulse interference, r f is an adjustable factor set empirically; If the absolute value of the wavelet coefficient at a certain position is greater than the set threshold, then this position is determined to be a pulse interference, and the position index of all pulse interferences is obtained, recorded as η(j,n): Threshold function based on trimmed median filter for: Where α∈(0,1) is an adjustable factor; MED[.] represents median filtering of the absolute value of the wavelet coefficients, and l is the number of sample points of a single pulse interference; The threshold function calculated by trimmed median filtering is used to detect pulse interference on the wavelet coefficients of each scale after DWT decomposition, and the detection results η(j,n) of each decomposition scale are obtained. η(j,n)=1 indicates that pulse interference is detected, and η(j,n)=0 indicates that pulse interference is not detected. Step 1-3: Jointly determine pulse interference at multiple scales; Perform differential operation on the detected pulse interference index: q(j,n)=η(j,n)-η(j,n-1) (8); Among them, q(j,n) is the result of the difference operation, and the index of the adjacent q(j,n)=1 and q(j,n)=-1 is recorded as the starting point n of the i-th pulse interference of the j-th scale. s (j,i) and end point n e (j,i); if the length of the j-th scale wavelet coefficient is recorded as Then the length of the j+1th scale wavelet coefficient is recorded as The wavelet coefficients of each scale are linearly interpolated as shown in formula (9) so that the length of the wavelet coefficients of all scales is the same as the length of the wavelet coefficients of the j=1 scale; the detection results after interpolation of each scale are recorded as: The starting and ending point indices are recorded as: Combine the detection results of all scales and perform voting statistics, recorded as V(n), then: Set the threshold T. If the voting result V(n) ≥ T, it is determined that pulse interference occurs at discrete time point n. Set T = J / 2 to indicate that at least half of the scales detect that there is a pulse here, which is effective pulse interference. The final detection result of the comprehensive judgment of each scale is: The start and end indexes of the pulse interference in the comprehensive judgment result are: The comprehensive judgment result is downsampled so that the length of the wavelet coefficients at each scale is the same as the length of the initial wavelet coefficients at each scale; the detection result of each scale after updating is recorded as η(j,n), and the relationship between it and the detection result η(n) of the comprehensive judgment of each scale is as follows: The starting and ending indexes of the pulse are recorded as: and Step 2: Pulse interference suppression; Step 2-1: pulse interference wavelet coefficient processing; Wavelet coefficients DWT of the original signal f Subtract the wavelet coefficients at the pulse interference from (j,n) to obtain the wavelet coefficients DWT without pulse interference components x (j,n), that is: Calculate the signal envelope after removing the pulse interference, denoted as ξ x (j,n): x x (j,n)=min{φ1(j,n),φ2(j,n)} (14); in: Similarly, calculate the signal envelope of the pulse interference, denoted as ξ η (j,n); The following processing is performed on the wavelet coefficients of the pulse interference and the pulse interference vicinity to obtain the processing template Ξ(j,n) of the wavelet coefficients in the pulse interference region: Multiply the processing template Ξ(j,n) by the original signal wavelet coefficient to achieve zero wavelet coefficients at the pulse interference and reset the wavelet coefficients near the pulse interference in proportion: Step 2-2: Replacement of pulse interference coefficients; The starting point and end point of the pulse interference below no longer distinguish the specific decomposition scale j, and j is deleted from the symbol. The starting point of the i-th pulse interference at the j-th scale is recorded as n s (j,i) is End point e (j,i) is (a) Design a comparison template for the current pulse interference. The starting range of the comparison template is the same as the length of the processing template. L2 represents the distance between the starting point and the end point of the pulse interference. Used to define the scope of comparison templates and processing templates; (b) Define a search area. The starting range of the search area is L1>L2 and the processing template slides in the search area, and the energy in the area is calculated after multiplying it with the current wavelet coefficient; L1 represents the distance between the starting point and the end point of the pulse interference and extends to both sides, which is used to define the range of the search area; Through (a) and (b), the energy value in the current pulse interference area calculated by the comparison template is obtained, which is recorded as is a constant; use the processing template to slide in the search area and calculate the energy value in the current template, recorded as Compare and Similarity degree, the similarity calculation formula is: in: In the formula is the energy of the current search window; To compare the energy of the template, for the i-th pulse interference is a constant; n1 is the kth point of the sliding window, n2 is the kth point of the comparison template where the pulse is located, K i is the number of sample points of the comparison template; After completing the above similarity calculation, the index of the start of the search window with the maximum similarity is obtained Then the starting position index of the pulse interference to be replaced and for: in, Indicates the starting point of the optimal search window; If L3 sample points are extended on both sides of the replacement part and L3 < L2, the new area of the replacement part is For the processing template in equation (16), find the minimum value and record the index of the first point as The index of the last point is Then the replacement smoothing template at this time is: The wavelet coefficients of each scale after the final pulse interference suppression are:

2. The method for detecting and suppressing pulse interference of a vector sensor based on trimmed median filtering according to claim 1 is characterized in that: The Env[·] is implemented using a median filter.

3. The method for detecting and suppressing pulse interference of a vector sensor based on trimmed median filtering according to claim 2 is characterized in that: The median filter calculation is defined as: Γ(j,n)=r f MED[DWT f (j,nl),…,DWT f (j,n),…,DWT f (j,n+l)]; The length of the median filter is L med Need to meet: L med >2l+1。 4. A computer program, characterized in that The computer program enables a computer to execute the method according to any one of claims 1 to 3.

5. An electronic device, characterized in that: include: Processor and memory; The memory is used to store a computer program, and the processor is used to execute the computer program stored in the memory, so that the electronic device executes the method as claimed in any one of claims 1 to 3.

6. A computer-readable storage medium having a computer program stored thereon, characterized in that: When the computer program is executed by a processor, the method according to any one of claims 1 to 3 is implemented.

7. A chip, characterized in that: include: A processor, configured to call and run a computer program from a memory, so that a device equipped with the chip executes a method as claimed in any one of claims 1 to 3.

8. A computer program product, characterized in that The computer program product comprises a computer storage medium storing a computer program, wherein the computer program comprises instructions executable by at least one processor, and when the instructions are executed by the at least one processor, the method according to any one of claims 1 to 3 is implemented.

Images