CN111681667B - Underwater sound signal denoising method based on adaptive window filtering and wavelet threshold optimization - Google Patents

Underwater sound signal denoising method based on adaptive window filtering and wavelet threshold optimization Download PDF

Info

Publication number
CN111681667B
CN111681667B CN202010509512.6A CN202010509512A CN111681667B CN 111681667 B CN111681667 B CN 111681667B CN 202010509512 A CN202010509512 A CN 202010509512A CN 111681667 B CN111681667 B CN 111681667B
Authority
CN
China
Prior art keywords
noise
signal
gaussian
new
threshold
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202010509512.6A
Other languages
Chinese (zh)
Other versions
CN111681667A (en
Inventor
王景景
李嘉恒
董新利
杨星海
施威
徐凌伟
郭瑛
李海涛
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Qingdao University of Science and Technology
Original Assignee
Qingdao University of Science and Technology
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Qingdao University of Science and Technology filed Critical Qingdao University of Science and Technology
Priority to CN202010509512.6A priority Critical patent/CN111681667B/en
Publication of CN111681667A publication Critical patent/CN111681667A/en
Priority to PCT/CN2021/088635 priority patent/WO2021258832A1/en
Application granted granted Critical
Publication of CN111681667B publication Critical patent/CN111681667B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G10—MUSICAL INSTRUMENTS; ACOUSTICS
    • G10L—SPEECH ANALYSIS TECHNIQUES OR SPEECH SYNTHESIS; SPEECH RECOGNITION; SPEECH OR VOICE PROCESSING TECHNIQUES; SPEECH OR AUDIO CODING OR DECODING
    • G10L21/00—Speech or voice signal processing techniques to produce another audible or non-audible signal, e.g. visual or tactile, in order to modify its quality or its intelligibility
    • G10L21/02—Speech enhancement, e.g. noise reduction or echo cancellation
    • G10L21/0208—Noise filtering
    • G10L21/0216—Noise filtering characterised by the method used for estimating noise
    • H—ELECTRICITY
    • H04—ELECTRIC COMMUNICATION TECHNIQUE
    • H04B—TRANSMISSION
    • H04B11/00—Transmission systems employing sonic, ultrasonic or infrasonic waves
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2218/00—Aspects of pattern recognition specially adapted for signal processing
    • G06F2218/02—Preprocessing
    • G06F2218/04—Denoising
    • G06F2218/06—Denoising by applying a scale-space analysis, e.g. using wavelet analysis
    • G—PHYSICS
    • G10—MUSICAL INSTRUMENTS; ACOUSTICS
    • G10L—SPEECH ANALYSIS TECHNIQUES OR SPEECH SYNTHESIS; SPEECH RECOGNITION; SPEECH OR VOICE PROCESSING TECHNIQUES; SPEECH OR AUDIO CODING OR DECODING
    • G10L25/00—Speech or voice analysis techniques not restricted to a single one of groups G10L15/00 - G10L21/00
    • G10L25/27—Speech or voice analysis techniques not restricted to a single one of groups G10L15/00 - G10L21/00 characterised by the analysis technique
    • G—PHYSICS
    • G10—MUSICAL INSTRUMENTS; ACOUSTICS
    • G10L—SPEECH ANALYSIS TECHNIQUES OR SPEECH SYNTHESIS; SPEECH RECOGNITION; SPEECH OR VOICE PROCESSING TECHNIQUES; SPEECH OR AUDIO CODING OR DECODING
    • G10L25/00—Speech or voice analysis techniques not restricted to a single one of groups G10L15/00 - G10L21/00
    • G10L25/45—Speech or voice analysis techniques not restricted to a single one of groups G10L15/00 - G10L21/00 characterised by the type of analysis window

Landscapes

  • Engineering & Computer Science (AREA)
  • Signal Processing (AREA)
  • Computer Networks & Wireless Communication (AREA)
  • Computational Linguistics (AREA)
  • Quality & Reliability (AREA)
  • Health & Medical Sciences (AREA)
  • Audiology, Speech & Language Pathology (AREA)
  • Human Computer Interaction (AREA)
  • Physics & Mathematics (AREA)
  • Acoustics & Sound (AREA)
  • Multimedia (AREA)
  • Measurement Of Velocity Or Position Using Acoustic Or Ultrasonic Waves (AREA)

Abstract

The invention discloses an underwater sound signal denoising method based on adaptive window filtering and wavelet threshold optimization. Firstly, describing Gaussian/non-Gaussian pulse noise in an underwater acoustic channel by combining an S alpha S distribution model and a normal distribution model; a median filtering method based on a self-adaptive window is designed, the size of a filtering window is corrected according to the number of noise points in the window, and non-Gaussian pulse noise is suppressed; then, based on an improved artificial bee colony method GDES-ABC, threshold parameters of a wavelet threshold denoising method are optimized, and the suppression capability of Gaussian noise is improved. The invention can effectively inhibit Gaussian/non-Gaussian pulse noise in a complex underwater acoustic environment, improve the receiving capability of 2FSK, QPSK, 16QAM and other underwater acoustic communication signals and obtain higher output signal-to-noise ratio and noise suppression ratio.

Description

Underwater sound signal denoising method based on adaptive window filtering and wavelet threshold optimization
Technical Field
The invention belongs to the technical field of underwater acoustic signal denoising, and particularly relates to an underwater acoustic signal denoising method based on Adaptive Window Filtering (AWFM) and wavelet threshold optimization (GDES) in a Gaussian/non-Gaussian impulse noise environment.
Background
The sound wave is widely applied in the field of underwater communication, and in the process of underwater transmission and processing, sound wave signals are influenced by underwater complex Gaussian/non-Gaussian pulse noise, so that the sound wave signals are degraded and distorted, and the communication quality is reduced. The signal denoising technology is a signal processing method for improving signal quality and reducing noise influence, and is widely applied to the fields of underwater acoustic communication and the like.
The sudden non-gaussian impulse noise emitted by underwater submarine exploration, marine organisms, sea surface waves, submarine earthquakes and the like can be suppressed by a Standard Media Filter (SMF) method. However, the SMF simultaneously processes the useful signal portion of the received signal, resulting in distortion of the useful signal.
Gaussian noise can be effectively suppressed by a filtering method, a wavelet transformation method, an empirical mode decomposition method and the like. The denoising method based on the wavelet threshold can obtain the progressive optimal estimation of the original signal, and is widely applied. The main factors that determine the performance of the method are the accurate estimation of the threshold and the reasonable construction of the threshold function. The unified threshold value based on the multidimensional independent normal variable decision theory is commonly used under the assumption that a noise model is a Gaussian noise model; however, the uniform threshold depends on an accurate estimate of the noise variance and is difficult to apply in the case of an actually unknown noise variance. Common threshold functions comprise a hard threshold, a soft threshold, a semi-soft threshold and the like, the method processes wavelet coefficients according to a fixed structure, the self-adaptability is lacked, and the flexibility of signal processing is reduced. To overcome the above limitations, a method based on group intelligence optimization is introduced to improve the wavelet threshold denoising performance.
However, in an actual marine background noise environment, both gaussian and non-gaussian impulse noise are included, and the wavelet threshold denoising method based on intelligent optimization is mainly used for gaussian noise processing, is difficult to be directly applied to overall processing of marine underwater acoustic noise, still has a lot of disadvantages, and is specifically embodied in that: firstly, a general principle of uniformly establishing a threshold function is lacked, so that the construction of the threshold function is difficult; secondly, the determination of the threshold parameter is an iterative process, and usually a suboptimal value rather than an optimal value is achieved; and as the number of method iterations increases, population diversity decreases, the above optimization method may fall into local minima.
Disclosure of Invention
The invention provides an underwater acoustic signal denoising method based on Adaptive Window Filtering (AWFM) and wavelet threshold optimization (GDES) in a Gaussian/non-Gaussian pulse noise environment, which is used for overcoming the defects of the prior art.
The invention firstly describes Gaussian/non-Gaussian impulse noise in an underwater acoustic channel by combining an S alpha S distribution model and a normal distribution model, and the model is a limit distribution which can keep a generation mechanism and a propagation condition of complex marine background noise and can better describe the marine environment background noise. A median filtering method based on a self-adaptive window is designed, the size of the filtering window is corrected according to the number of noise points in the window, non-Gaussian pulse noise is suppressed, the processing of the non-noise points is avoided, the distortion of useful signals is effectively reduced, meanwhile, the size of the filtering window is self-adaptively adjusted based on the noise content, and the filtering performance and the calculation complexity of the method are effectively balanced; then, based on an improved artificial bee colony method GDES-ABC, threshold parameters of a wavelet threshold denoising method are optimized, and the suppression capability of Gaussian noise is improved.
In order to realize the purpose of the invention, the invention is realized by adopting the following specific technical scheme:
an underwater sound signal denoising method based on adaptive window filtering and wavelet threshold optimization comprises the following steps:
s1: acquiring underwater acoustic signal data, and describing Gaussian/non-Gaussian pulse noise in an underwater acoustic channel by using the conventional data receiving model;
s2: suppressing non-Gaussian pulse noise based on a median filtering method of a self-adaptive window to obtain underwater sound signal data without the non-Gaussian pulse noise;
s3: and then, the underwater acoustic signal data without the non-Gaussian pulse noise obtained in the step S2 is processed based on an improved artificial bee colony wavelet threshold optimization method (recorded as GDES-ABC), Gaussian noise is suppressed, and finally the water acoustic signal data after denoising is obtained.
The step S1 is specifically as follows:
s1-1: a signal receiving model, wherein the signal noise model adopts an S alpha S distribution model and a normal distribution model:
for a single-transmitting single-receiving underwater acoustic communication system, a time domain signal y (t) received by a receiving end is represented in a digital form, and a group of discrete samples are represented as follows:
y(i)＝s(i)+e(i),i＝1,2,...,N
where s (i) is a noise-free desired signal having random amplitude and phase; e (i) additive ocean background noise; n is the number of samples;
s1-2: gaussian/non-gaussian pulse noise model:
the probability density function selected for the gaussian noise model is as follows: wherein x is an instantaneous value of the noise pressure;
Figure GDA0002551939170000021
the signal-to-noise ratio SNR is defined as follows:
Figure GDA0002551939170000022
wherein
Figure GDA0002551939170000023
Variance of the desired signal and gaussian noise, respectively;
underwater non-gaussian noise sources including sound waves emitted by submarine exploration, marine organisms, sea surface waves, submarine earthquakes, and the like; the probability density function of the sound wave signals emitted by the noise source is similar to normal distribution, but the probability of tailing and strong amplitude is higher, the duration is shorter, the spike pulse characteristic is realized, and the method belongs to a burst non-Gaussian pulse signal;
the alpha stable distribution is adopted to describe the underwater spike pulse noise and has greater advantage than the Gaussian distribution if the characteristic function of the random variable X
Figure GDA0002551939170000024
Can be expressed as:
Figure GDA0002551939170000025
wherein
Figure GDA0002551939170000026
- ∞ < a <infinityis a real number, represents a position parameter, and
Figure GDA0002551939170000027
Figure GDA0002551939170000028
the random variable X obeys alpha stable distribution; wherein alpha is more than 0 and less than or equal to 2, the characteristic index is represented, the pulse characteristic degree is determined, the smaller alpha is, the stronger the pulse is, the larger alpha is, the closer alpha is to the Gaussian process, and when alpha is 2, the Gaussian distribution is obtained; -1 ≦ β ≦ 1 is a symmetry parameter for determining the slope of the distribution; gamma > 0 is a dispersion coefficient, the meaning of which is similar to the variance of Gaussian distribution;
when β is 0, α stable distribution characteristic function is expressed as
Figure GDA0002551939170000031
At this time, the distribution is called symmetrical alpha stable distribution and is recorded as X-S alpha S; assuming that the position parameter a is 0, the probability density function of the S α S distribution is:
Figure GDA0002551939170000032
non-gaussian impulse noise based on the S α S distribution cannot calculate the variance, so the noise magnitude is described by using a mixed signal-to-noise ratio MSNR, which is defined as follows:
Figure GDA0002551939170000033
wherein
Figure GDA0002551939170000034
And γ represents a variance of a desired signal and a dispersion coefficient of non-gaussian impulse noise, respectively;
to better describe the noise in the actual environment, the underwater acoustic noise model is assumed to be obtained by superposing Gaussian noise and non-Gaussian pulse noise models, so that the underwater acoustic noise e (i) is defined as
e(i)＝eGauss(i)+eSαS(i)
Wherein eGauss(i)、eSαS(i) Are respectively composed of
Figure GDA0002551939170000035
And
Figure GDA0002551939170000036
and (4) generating.
The step S2 is as follows:
s2-1: and (3) noise point detection:
assume that the signal received at the receiving end is y ═ y (1), y (2),.., y (n)]The initial sliding window W has a length LWTaking out a sample W (i) corresponding to the time-th reception signal y excluding the central point y (i) by using the initial sliding window W, where 2n + 1:
w(i)＝[w1(i),w2(i),...,w2n(i)]
＝[y(i-n),...,y(i-1),y(i+1),...,y(i+n)]
ordering the signal points in w (i) from small to large
r(i)＝sort(w(i))
＝[r1(i),r2(i),...r2n(i)]
Wherein sort (·) is a ranking function; let Med mean (r (i)), mean () means taking the median; defining a differential noise identifier as
Figure GDA0002551939170000037
For a given predetermined pulse threshold TnoiseIf d (i) > TnoiseIf yes, y (i) is determined to be an impulse noise point, and n (i) is made equal to 1, otherwise y (i) is an expected signal, and n (i) is made equal to 0, wherein n (i) represents an impulse mark; in the underwater sound receiving signal, the sound velocity is c, the amplitude is A, and the sampling frequency is fsCarrier frequency of fcThe rate of change of any adjacent sampling point will not exceed
Figure GDA0002551939170000041
And the sampling length of the underwater sound receiving signal is
Figure GDA0002551939170000042
Thereby setting the pulse threshold value to
Figure GDA0002551939170000043
S2-2: adaptive window size determination:
w length L for initial sliding windowWWhen the center point y (i) is an impulse noise point, calculating the number of noise points in the window:
Figure GDA0002551939170000048
the new window is marked as WnewThe length is as follows:
Figure GDA0002551939170000044
s2-3: noise filtering:
according to the new window size
Figure GDA0002551939170000045
And pulse mark n (i) filtering the received signal; wherein the non-noise points in the new window are unchanged and the noise points are replaced by the median value of the signal in the new window; suppose yip(i) Is WnewInner noise point, taking out new window and removing yip(i) All signal samples w ofnew(i)：
Figure GDA0002551939170000046
To wnew(i) The medium signal points are ordered from small to large to obtain
Figure GDA0002551939170000047
The noise point yip(i) Substituted by the formula:
y′ip(i)＝median(rnew(i))
wherein y'ip(i) Is the filtered signal.
The details of S3 above are as follows:
to y 'first'ip(i) Performing wavelet transformation to obtain wavelet coefficient wj,kWherein j, k represents the kth coefficient of the jth layer;
and constructing a new threshold function:
s3-1: a new threshold function construction:
the invention proposes a new adaptive threshold function:
Figure GDA0002551939170000051
wherein
Figure GDA0002551939170000052
Is an exponential factor, taking a non-negative number, lambdajA j-th layer threshold, wherein j is 1, 2.., and L is the number of decomposition layers;
s3-2: the new threshold function property proves that:
from the definition of continuity, it is easy to prove that the new adaptive threshold function is at (- ∞, - λj)，(-λj,+λj) And (+ λ)j, + ∞) is continuous; when w isj,k＞λjThen, the new adaptive threshold function can be written as:
Figure GDA0002551939170000053
then
Figure GDA0002551939170000054
When in use
Figure GDA0002551939170000055
When | wj,k|＜λjNew adaptive threshold function
Figure GDA0002551939170000056
Namely, it is
Figure GDA0002551939170000057
Thus, it is possible to provide
Figure GDA0002551939170000058
A new adaptive threshold function is available at point wj,k＝λjContinuously; when w isj,k＜-λjThen, the new adaptive threshold function can be written as:
Figure GDA0002551939170000059
then
Figure GDA00025519391700000510
When | wj,k|≤λjWhen the temperature of the water is higher than the set temperature,
Figure GDA00025519391700000511
namely, it is
Figure GDA00025519391700000512
Therefore, the temperature of the molten metal is controlled,
Figure GDA0002551939170000061
a new adaptive threshold function is available at point wj,k＝-λjContinuously;
when w isj,kThe time → + ∞ times,
Figure GDA0002551939170000062
when w isj,kThe time is → -infinity,
Figure GDA0002551939170000063
namely, it is
Figure GDA0002551939170000064
Thus, it is possible to provide
Figure GDA0002551939170000065
Is an asymptote of the new adaptive threshold function;
s3-3: determining a threshold parameter to be optimized:
the invention relates to a threshold lambda in a new adaptive threshold functionjAnd an exponential factor
Figure GDA0002551939170000067
The estimation method is regarded as an unknown threshold parameter, then the threshold parameter is optimized and estimated based on a GDES-ABC method, the estimation precision and speed are improved, and the denoising performance of the method is guaranteed;
s3-4: population initialization based on a good point set:
population initialization based on the good point set can effectively improve population diversity and avoid the method from falling into local optimum too early; the construction method of the good point comprises the following steps:
Figure GDA0002551939170000066
where p is the smallest prime number satisfying (p-3)/2 ≧ D, D is the dimension of the solution, deci {. cndot. } denotes the fractional part, rkIs a good point; thus, the set of good points [ PSN(1),PSN(2),...,PSN(SN)]TThe construction method of (2) is as follows:
PSN(i)＝{deci{r1*i},...,deci{rD*i}},i＝1,2,...,SN
wherein [. ]]TIndicating transposition, SN is the population size; the initial population is
X＝Lb+(Ub-Lb)*PSN
Wherein Ub and Lb are the upper and lower bounds of the solution, respectively, D is the current dimensional solution space, D is 1, 2.
S3-5: field search strategy guided by dynamic elite population:
the dynamic elite population contains better solutions in the population, and the scale of the solutions is different along with the difference of the iteration times; the dynamic smart English group-based neighborhood search can effectively accelerate the convergence speed of the method and improve the search efficiency; first calculate each X i1,2, the fitness value of SN, and then taking the better Telite ceil (p' SN) bees to form a dynamic elite population DXE i′1, 2., Telite, where ceil (·) denotes rounding up and p' is the proportion of elite population in all populations, determined according to the following formula:
Figure GDA0002551939170000071
wherein p ismaxAnd pminRespectively representing the maximum value and the minimum value of p'; t is the current number of iterations, tmaxIs the maximum number of iterations; from the above formula, it can be seen that in the initial stage of the method, t and p' are very small, and at this time, the dynamic elite population contains several optimal solutions, so that based on the elite population, the neighborhood search is more targeted, and the convergence rate can be greatly accelerated; at the later stage of the method, t and p' are larger, at the moment, the dynamic elite population contains more solutions, wherein some relatively poor solutions may be contained, so that the population has more diversity, and the capability of the method for jumping out of local optimum and finding out global optimum is enhanced;
the improved neighborhood search method based on the dynamic elite population of the GDES-ABC method is as follows:
vid＝DXECd+φid(Gbestd-xkd),d＝1,2,...,D
wherein phi isidIs [ -1,1 [ ]]Random real numbers in between; gbest is the global optimal solution, xkdIs xidRandom neighborhood of (D) DXECdIs a dynamic elite population center, the expression is as follows:
Figure GDA0002551939170000072
the neighborhood search strategy based on the dynamic elite population is described as follows: for each neighborhood search, the hiring bee randomly searches neighborhoods with the same probability and generates new solutions from the improved neighborhood search method; observing the bee to randomly search neighborhoods from the elite population, and generating a new solution by an improved neighborhood searching method; for the observation bees, if the new solution is better than the current solution, selecting the new solution to perform next neighborhood search; otherwise, in the next neighborhood searching, randomly searching the neighborhood from the elite population again to perform the next neighborhood searching; the neighborhood search strategy is carried out randomly, so that the population diversity is ensured, and invalid search is avoided;
s3-6: simulated annealing selection mechanism:
assume that at the T-th iteration, the current temperature is TtThe annealing parameter is K, and a new solution is obtained by an improved neighborhood search method and is VtWith a fitness value of fitv(ii) a The simulated annealing selection mechanism is as follows: if fitv＞fitiThen accept a new solution, where fitiIs the current solution fitness value; otherwise, accepting the new solution by using the probability P, wherein the acceptance probability is defined as
Figure GDA0002551939170000081
Wherein
Figure GDA0002551939170000082
Changing with the number of iterations; wherein beta is not more than 1 and is a constant with the value of 0.7 and sigmafitThe standard deviation of the fitness of all solutions;
from the acceptance probability, the method is initially characterized by small T and TtThe size is larger, so that P is larger, the method accepts a certain difference value, and the bee colony has larger development capability; at the end of the process, T becomes larger, TtAnd P is reduced, and the method rejects the difference value with higher probability, thereby ensuring the capability of searching the optimal solution and avoiding invalid search.
S3-7: taking the mean square error between the signal to be denoised and the denoised signal as a fitness function of the improved artificial bee colony method in S3, and obtaining an optimal threshold parameter under the condition of obtaining the minimum mean square error; the fitness function is expressed as:
Figure GDA0002551939170000083
where s (i) is a training signal,
Figure GDA0002551939170000089
to reconstruct the training signal, N is the training signal length. Known from the new adaptive threshold function, the threshold function is the parameter λjAnd
Figure GDA0002551939170000085
a function of (a); once λjAnd
Figure GDA0002551939170000086
determining a new self-adaptive threshold function of the threshold function, and obtaining a wavelet coefficient after the threshold is shrunk, thereby reconstructing a denoised signal
Figure GDA0002551939170000087
Therefore, in the GDES-ABC method, the parameter λ can be determinedjAnd
Figure GDA0002551939170000088
the composed vector is regarded as the position of the honey source, and the optimal threshold parameter is obtained by minimizing the fitness function.
And finally, performing contraction processing on the wavelet coefficient by constructing a new threshold function and an optimal threshold parameter to obtain a new wavelet coefficient, and performing inverse wavelet transformation to obtain a de-noising signal.
Compared with the prior art, the invention has the following advantages and technical effects:
the method inhibits non-Gaussian pulse noise based on a median filtering method of a self-adaptive window, optimizes threshold parameters of a wavelet threshold denoising method based on an improved artificial bee colony method GDES-ABC, and improves the inhibition capability of Gaussian noise. The invention can obtain higher output signal-to-noise ratio (SNR) and Noise Suppression Ratio (NSR), and effectively improves the receiving capability of the underwater acoustic communication machine to 2FSK, QPSK, 16QAM and other underwater acoustic communication signals in Gaussian/non-Gaussian pulse noise environment.
Drawings
FIG. 1 is a schematic flow diagram of the present invention;
FIG. 2 is a schematic diagram of an underwater artificial and natural non-Gaussian noise source;
FIG. 3LWWhen 2n +1An initial sliding window schematic;
FIG. 4 is a comparison of different threshold functions of the present invention;
FIG. 5 is a comparison graph of output SNR after denoising different underwater acoustic communication signals based on different denoising methods, which changes with input SNR (5-1, 5-2, 5-3 correspond to 2FSK, QPSK, 16QAM, respectively);
FIG. 6 is a comparison graph of the output NSR after denoising different underwater acoustic communication signals based on different denoising methods and changing with the input SNR (6-1, 6-2, 6-3 correspond to 2FSK, QPSK, 16QAM, respectively);
FIG. 7 is a comparison graph of output SNR after denoising different underwater acoustic communication signals based on different denoising methods, which changes with input MSNR (7-1, 7-2, 7-3 correspond to 2FSK, QPSK, 16QAM, respectively);
FIG. 8 is a comparison graph of the output NSR after denoising different underwater acoustic communication signals based on different denoising methods and changing with the input MSNR (8-1, 8-2, 8-3 correspond to 2FSK, QPSK, 16QAM, respectively);
FIG. 9 is a time domain waveform contrast diagram (9-1, 9-2, 9-3 correspond to 2FSK, QPSK, 16QAM, respectively) after denoising different underwater acoustic communication signals based on the invention.
Detailed Description
In order to make the objects, technical solutions and advantages of the present invention more apparent, the present invention will be described in further detail with reference to the accompanying drawings and examples.
Example 1:
referring to fig. 1, the underwater acoustic signal denoising method based on AWFM + GDES in the gaussian/non-gaussian impulse noise environment according to the embodiment includes the following steps:
s1: describing Gaussian/non-Gaussian pulse noise in an underwater acoustic channel by combining an S alpha S distribution model and a normal distribution model; the method comprises the following specific steps:
s1-1: the signal receiving model comprises the following steps:
for a single-transmitting single-receiving underwater acoustic communication system, a time domain signal y (t) received by a receiving end is represented in a digital form, and a group of discrete samples are represented as follows:
y(i)＝s(i)+e(i),i＝1,2,...,N
where s (i) is a noise-free desired signal having random amplitude and phase; e (i) additive ocean background noise; n is the number of samples;
s1-2: gaussian/non-gaussian pulse noise model:
the probability density function of the instantaneous value x of the sound pressure of the Gaussian noise source is
Figure GDA0002551939170000091
The SNR is defined as follows:
Figure GDA0002551939170000092
wherein
Figure GDA0002551939170000093
Variance of the desired signal and gaussian noise, respectively;
referring to fig. 2, there are artificial and natural non-gaussian noise sources under water, including sound waves emitted by undersea exploration, marine life, sea surface waves, and undersea earthquakes. The probability density function of the sound wave signals emitted by the noise source is similar to normal distribution, but the probability of tailing and strong amplitude is higher, the duration is shorter, the spike pulse characteristic is realized, and the method belongs to a burst non-Gaussian pulse signal;
it is more advantageous to describe the underwater spike impulse noise with the stable α distribution than the gaussian distribution if the characteristic function of the random variable X can be expressed as:
Figure GDA0002551939170000094
wherein
Figure GDA0002551939170000095
- ∞ < a <infinityis a real number, represents a position parameter, and
Figure GDA0002551939170000101
Figure GDA0002551939170000102
the random variable X obeys alpha stable distribution; wherein, alpha is more than 0 and less than or equal to 2, which represents a characteristic index and determines the characteristic degree of the pulse, the smaller alpha is, the stronger the pulse is, the larger alpha is, the closer alpha is to the Gaussian process, and when alpha is 2, the Gaussian distribution is obtained; -1 ≦ β ≦ 1 is a symmetry parameter for determining the slope of the distribution; gamma > 0 is a dispersion coefficient, the meaning of which is similar to the variance of Gaussian distribution;
when β is 0, α stable distribution characteristic function is expressed as
Figure GDA0002551939170000103
At this time, the distribution is called symmetrical alpha stable distribution and is recorded as X-S alpha S; assuming that a is 0, the probability density function of the S α S distribution is:
Figure GDA0002551939170000104
non-gaussian impulse noise based on the S α S distribution cannot calculate the variance, so the noise magnitude is described using MSNR, which is defined as follows:
Figure GDA0002551939170000105
wherein
Figure GDA0002551939170000106
And γ represents a variance of a desired signal and a dispersion coefficient of non-gaussian impulse noise, respectively; to better describe the noise in the actual environment, the underwater acoustic noise model is assumed to be obtained by superposing Gaussian noise and non-Gaussian pulse noise models, so that the underwater acoustic noise e (i) is defined as
e(i)＝eGauss(i)+eSαS(i)
Wherein eGauss(i)、eSαS(i) Are respectively composed of
Figure GDA0002551939170000107
And
Figure GDA0002551939170000108
and (4) generating.
S2: based on a median filtering method of a self-adaptive window, recording the median filtering method as AWFM, and inhibiting non-Gaussian pulse noise; the method comprises the following specific steps:
s2-1: and (3) noise point detection:
assume that the signal received at the receiving end is y ═ y (1), y (2),.., y (n)]The initial sliding window W has a length LWReferring to fig. 3, a sample W (i) corresponding to the time-th reception signal y excluding the center point y (i) is extracted using the initial sliding window W as 2n + 1:
w(i)＝[w1(i),w2(i),...,w2n(i)]
＝[y(i-n),...,y(i-1),y(i+1),...,y(i+n)]
ordering the signal points in w (i) from small to large
r(i)＝sort(w(i))
＝[r1(i),r2(i),...r2n(i)]
Where sort (·) is a ranking function. Let Med mean (r (i)), mean () denotes taking the median. Defining a differential noise identifier as
Figure GDA0002551939170000111
For a given predetermined pulse threshold TnoiseIf d (i) > TnoiseThen, y (i) is determined as the impulse noise point, and n (i) is made equal to 1, otherwise y (i) is the desired signal, and n (i) is equal to 0, where n (i) represents the impulse mark. In the underwater sound receiving signal, the sound velocity is c, the amplitude is A, and the sampling frequency is fsCarrier frequency of fcThe rate of change of any adjacent sampling point will not exceed
Figure GDA0002551939170000112
And the sampling length of the underwater sound receiving signal is
Figure GDA0002551939170000113
Thereby setting the impact threshold value to
Figure GDA0002551939170000114
S2-2: adaptive window size determination:
w length L for initial sliding windowWWhen the center point y (i) is an impulse noise point, calculating the number of noise points in the window:
Figure GDA0002551939170000118
the new window is marked as WnewThe length is as follows:
Figure GDA0002551939170000115
s2-3: noise filtering:
according to the new window size LWnew(y (i)) and pulse labels n (i) filtering the received signal. Wherein the non-noise points within the new window are unchanged and the noise points are replaced by the median value of the signal within the new window. Suppose yip(i) Is WnewInner noise point, taking out new window and removing yip(i) All signal samples w ofnew(i)：
Figure GDA0002551939170000116
To wnew(i) The medium signal points are ordered from small to large to obtain
Figure GDA0002551939170000117
The noise point yip(i) Substituted by the formula:
y′ip(i)＝median(rnew(i))
wherein y'ip(i) Is the filtered signal.
S3: optimizing a threshold parameter of a wavelet threshold denoising method based on a GDES-ABC method, and inhibiting Gaussian noise; the method comprises the following specific steps:
s3-1: a new threshold function construction:
the invention proposes a new adaptive threshold function:
Figure GDA0002551939170000121
wherein
Figure GDA0002551939170000122
Is an exponential factor, taking a non-negative number, lambdajA j-th layer threshold, j 1,2, and L is the number of decomposition layers.
S3-2: the new threshold function property proves that:
from the definition of continuity, it is easy to prove that the new adaptive threshold function is at (- ∞, - λj)，(-λj,+λj) And (+ λ)j, + ∞) is continuous. When w isj,k＞λjThen, the new adaptive threshold function can be written as:
Figure GDA0002551939170000123
then
Figure GDA0002551939170000124
When in use
Figure GDA0002551939170000125
When | wj,k|＜λjNew adaptive threshold function
Figure GDA0002551939170000126
Namely, it is
Figure GDA0002551939170000127
Thus, it is possible to provide
Figure GDA0002551939170000128
A new adaptive threshold function is available at point wj,k＝λjAnd (4) continuous. When w isj,k＜-λjThen, the new adaptive threshold function can be written as:
Figure GDA0002551939170000129
then
Figure GDA00025519391700001210
When | wj,k|≤λjWhen the temperature of the water is higher than the set temperature,
Figure GDA0002551939170000131
namely, it is
Figure GDA0002551939170000132
Therefore, the temperature of the molten metal is controlled,
Figure GDA0002551939170000133
a new adaptive threshold function is available at point wj,k＝-λjAnd (4) continuous.
When w isj,kThe time → + ∞ times,
Figure GDA0002551939170000134
when w isj,kThe time is → -infinity,
Figure GDA0002551939170000135
namely, it is
Figure GDA0002551939170000136
Thus, it is possible to provide
Figure GDA0002551939170000137
Is an asymptote of the new adaptive threshold function, as shown in fig. 4, where threshold λ is 5 and the horizontal axis represents wavelet coefficients, ranging from-10, 10]And the vertical axis is a wavelet coefficient obtained after threshold shrinkage. As can be seen from the figure, the threshold function shown in the new adaptive threshold function is a compromise strategy between a soft threshold and a hard threshold, has better continuity and smoothness, retains a large wavelet coefficient and has stronger fidelity to a target signal.
S3-3: determining a threshold parameter to be optimized:
the invention relates to a threshold lambda in a new adaptive threshold functionjAnd an exponential factor
Figure GDA0002551939170000138
The method is regarded as unknown threshold parameters, then optimized estimation is carried out on the threshold parameters based on a GDES-ABC method, and estimation precision and speed are improved, so that denoising performance of the method is guaranteed.
S3-4: population initialization based on a good point set:
population initialization based on the good point set can effectively improve population diversity and avoid the method from falling into local optimum too early. The construction method of the good point comprises the following steps:
Figure GDA0002551939170000141
where p is the smallest prime number satisfying (p-3)/2 ≧ D, D is the dimension of the solution, deci {. cndot. } denotes the fractional part, rkIs the best point. Thus, the set of good points [ PSN(1),PSN(2),...,PSN(SN)]TThe construction method of (2) is as follows:
PSN(i)＝{deci{r1*i},...,deci{rD*i}},i＝1,2,...,SN
wherein [. ]]TIndicating transposition, SN is the population size; the initial population is
X＝Lb+(Ub-Lb)*PSN
Wherein Ub and Lb are the upper and lower bounds of the solution, respectively, D is the current dimensional solution space, D is 1, 2.
S3-5: field search strategy guided by dynamic elite population:
the dynamic elite population contains better solutions in the population, the scale of which varies with the number of iterations. The dynamic smart English group-based neighborhood search can effectively accelerate the convergence speed of the method and improve the search efficiency; first calculate each X i1,2, the fitness value of SN, and then taking the better Telite ceil (p' SN) bees to form a dynamic elite population DXE i′1, 2., Telite, where ceil (·) denotes rounding up and p' is the proportion of elite population in all populations, determined according to the following formula:
Figure GDA0002551939170000142
wherein p ismaxAnd pminRespectively representing the maximum value and the minimum value of p'; t is the current number of iterations, tmaxIs the maximum number of iterations. As can be seen from the above formula, t and p' are small in the early stages of the method, where the dynamic elite population contains several optimal solutions, and thusBased on the elite population, the neighborhood search is more targeted, and the convergence rate can be greatly accelerated. At the later stage of the method, t and p' are very small and large, and at the moment, the dynamic elite population contains more solutions, wherein relatively poor solutions may be contained, so that the population has more diversity, and the capability of the method for jumping out of local optimum and finding global optimum is enhanced.
The improved neighborhood search method based on the dynamic elite population of the GDES-ABC method is as follows:
vid＝DXECd+φid(Gbestd-xkd),d＝1,2,...,D
wherein phi isidIs [ -1,1 [ ]]Random real numbers in between; gbest is the global optimal solution, xkdIs xidRandom neighborhood of (D) DXECdIs a dynamic elite population center, the expression is as follows:
Figure GDA0002551939170000143
the neighborhood search strategy based on the dynamic elite population is described as follows: for each neighborhood search, the hiring bee randomly searches neighborhoods with the same probability, and new solutions are generated by the improved neighborhood search method. And observing the bee to randomly search neighborhoods from the elite population, and generating a new solution by an improved neighborhood searching method. For the observer bees, if the new solution is better than the current solution, the new solution is selected for the next neighborhood search. Otherwise, in the next neighborhood searching, randomly searching the neighborhood from the elite population again to perform the next neighborhood searching. The neighborhood search strategy is carried out randomly, so that the population diversity is ensured, and invalid search is avoided.
S3-6: simulated annealing selection mechanism:
assume that at the T-th iteration, the current temperature is TtThe annealing parameter is K, and a new solution is obtained by an improved neighborhood search method and is VtWith a fitness value of fitv. The simulated annealing selection mechanism is as follows: if fitv＞fitiThen accept a new solution, where fitiIs the current solution fitness value; whether or notThe new solution is accepted with a probability P, the acceptance probability being defined as
Figure GDA0002551939170000151
Wherein
Figure GDA0002551939170000152
Varying with the number of iterations. Wherein β ≦ 1 is a constant, usually 0.7, σfitIs the standard deviation of fitness of all solutions.
From the acceptance probability, the method is initially characterized by small T and TtLarger, and therefore larger P, the method accepts a certain difference, and the bee colony has larger exploitation capability. At the end of the process, T becomes larger, TtAnd P is reduced, and the method rejects the difference value with higher probability, thereby ensuring the capability of searching the optimal solution and avoiding invalid search.
S3-7: and taking the mean square error between the signal to be denoised and the denoised signal as a fitness function of the improved artificial bee colony method in S3, and obtaining the optimal threshold parameter under the condition of obtaining the minimum mean square error. The fitness function is expressed as:
Figure GDA0002551939170000153
where s (i) is a training signal,
Figure GDA0002551939170000154
to reconstruct the training signal, N is the training signal length. Known from the new adaptive threshold function, the threshold function is the parameter λjAnd
Figure GDA0002551939170000155
as a function of (c). Once λjAnd
Figure GDA0002551939170000156
determining a new self-adaptive threshold function of the threshold function, and obtaining a wavelet coefficient after the threshold is shrunk, thereby reconstructing a denoised signal
Figure GDA0002551939170000157
Therefore, in the GDES-ABC method, the parameter λ can be determinedjAnd
Figure GDA0002551939170000158
the composed vector can be regarded as the position of the honey source, and the optimal threshold parameter is obtained by minimizing the fitness function.
Example 2: comparative analysis of simulation test results
In this embodiment, common underwater acoustic communication signals such as 2FSK, QPSK, and 16QAM signals are regarded as SOI, and additive white gaussian noise and non-gaussian impulse noise are combined into underwater acoustic noise, so as to verify the performance of the present invention. Wherein the invention is denoted as AWFM + GDES; the computer configuration used for the simulation was: intel i5-4570 processor, Windows 7 operating system, 4G memory, MATLAB R2015 b.
The output SNR is defined as follows:
Figure GDA0002551939170000161
the Noise Suppression Ratio (NSR) is defined as follows:
Figure GDA0002551939170000162
wherein s (i) and
Figure GDA0002551939170000163
respectively a desired signal and an estimated signal;
Figure GDA0002551939170000164
and
Figure GDA0002551939170000165
respectively, the desired signal and the estimated signal mean, N being the signal length.
Fig. 5 shows the variation curve of output SNR with input SNR after AWFM + GDES and other methods denoise 2FSK, QPSK, 16QAM signals, respectively. Where the input MSNR is 20dB and the other parameter settings are as in table 1. As can be seen from the figure, for each signal, as the input SNR increases, the output SNR obtained by the 6 methods increases first and then stabilizes. Furthermore, the proposed AWMF method achieves a higher output SNR than the SMF method, and the proposed AWMF + GDES method achieves the highest output SNR than the other 5 methods. For 2FSK and QPSK signals, the AWFM + GDES method achieves a higher output SNR than the other three methods when the input SNR > -5 dB. For 16QAM signal, when the input SNR > -10dB, the AWFM + GDES method can obtain higher output SNR than the other three methods. This means that the proposed AWMF + GDES method is more suitable for de-noising of 16QAM signals. When the input SNR is 20dB, for 2FSK, QPSK, 16QAM signals, the AWFM + GDES method obtains output SNR at least 2.3dB, 1.3dB, 1.2dB higher than the other methods, respectively.
Method parameter set-up as presented in Table 1
Figure GDA0002551939170000166
Figure GDA0002551939170000171
Fig. 6 shows the curves of the changes of the output NSR with the input SNR after the 2FSK, QPSK, and 16QAM signals are denoised by different denoising methods. The parameters are set as in table 1, and the input MSNR is 20 dB. It can be seen from the graph that, as the input SNR increases, the output NSR obtained by the 6 methods increases first and then becomes stable, and the trend thereof is basically consistent with the result shown in fig. 5, which illustrates the effectiveness of the proposed AWMF + GDES method from another aspect.
Fig. 7 and 8 show output SNR and NSR results after denoising 2FSK, QPSK, and 16QAM signals when different denoising methods change with input MSNR, respectively. Where the input SNR is 5dB and the other parameter settings are as in table 1. As can be seen from the figure, the 6 methods obtain a gradual increase in output SNR and NSR as the input MSNR increases. Compared with the other five methods, the proposed AWMF + GDES method obtains the highest output SNR and NSR, which shows that the proposed AWMF + GDES method can obtain better underwater sound signal denoising performance in terms of both the output SNR and the NSR.
Fig. 9 shows the results of denoising 2FSK, QPSK, and 16QAM signals by the underwater acoustic signal denoising method based on AWFM + GDES, respectively. The parameter setting section is shown in table 1, where input SNR is 5dB and MSNR is 20 dB. As can be seen from the figure, all desired signals are contaminated by underwater acoustic noise. However, most of the noise is eliminated after the AWMF + GDES-based processing. And the denoised signal retains more detail information of the desired signal. Therefore, the underwater sound signal denoising method based on the AWFM + GDES is more suitable for the actual underwater environment.
The above-mentioned embodiments are merely preferred embodiments of the present invention, and the scope of the present invention is not limited thereto, so that variations based on the shape and principle of the present invention should be covered within the scope of the present invention.

Claims (3)

1. A denoising method of underwater sound signals based on adaptive window filtering and wavelet threshold optimization is characterized by comprising the following steps:
s1: acquiring underwater acoustic signal data, and describing Gaussian/non-Gaussian pulse noise in an underwater acoustic channel by using the conventional data receiving model;
s2: suppressing non-Gaussian pulse noise based on a median filtering method of a self-adaptive window to obtain underwater sound signal data without the non-Gaussian pulse noise;
s3: then, the underwater acoustic signal data without the non-Gaussian pulse noise obtained in the step S2 is processed based on an improved artificial bee colony wavelet threshold optimization method to suppress Gaussian noise, and finally the underwater acoustic signal data after denoising is obtained;
the step S1 is specifically as follows:
s1-1: a signal receiving model, wherein the signal noise model adopts an S alpha S distribution model and a normal distribution model:
for a single-transmitting single-receiving underwater acoustic communication system, a time domain signal y (t) received by a receiving end is represented in a digital form, and a group of discrete samples are represented as follows:
y(i)＝s(i)+e(i),i＝1,2,...,N
where s (i) is a noise-free desired signal having random amplitude and phase; e (i) additive ocean background noise; n is the number of samples;
s1-2: gaussian/non-gaussian pulse noise model:
the probability density function selected for the gaussian noise model is as follows: wherein x is an instantaneous value of the noise pressure;
Figure FDA0002987289830000011
the signal-to-noise ratio SNR is defined as follows:
Figure FDA0002987289830000012
wherein
Figure FDA0002987289830000013
Variance of the desired signal and gaussian noise, respectively;
the alpha stable distribution is adopted to describe the underwater spike pulse noise and has greater advantage than the Gaussian distribution if the characteristic function of the random variable X
Figure FDA0002987289830000014
Can be expressed as:
Figure FDA0002987289830000015
wherein
Figure FDA0002987289830000016
- ∞ < a <infinityis a real number, represents a position parameter, and
Figure FDA0002987289830000017
Figure FDA0002987289830000018
the random variable X obeys alpha stable distribution; wherein alpha is more than 0 and less than or equal to 2, the characteristic index is represented, the pulse characteristic degree is determined, the smaller alpha is, the stronger the pulse is, the larger alpha is, the closer alpha is to the Gaussian process, and when alpha is 2, the Gaussian distribution is obtained; -1 ≦ β ≦ 1 is a symmetry parameter for determining the slope of the distribution; gamma > 0 is a dispersion coefficient, the meaning of which is similar to the variance of Gaussian distribution;
when β is 0, α stable distribution characteristic function is expressed as
Figure FDA0002987289830000021
At this time, the distribution is called symmetrical alpha stable distribution and is recorded as X-S alpha S; assuming that the position parameter a is 0, the probability density function of the S α S distribution is:
Figure FDA0002987289830000022
non-gaussian impulse noise based on the S α S distribution cannot calculate the variance, so the noise magnitude is described by using a mixed signal-to-noise ratio MSNR, which is defined as follows:
Figure FDA0002987289830000023
wherein
Figure FDA0002987289830000024
And gamma represents the variance and non-gaussian pulse of the desired signal, respectivelyDispersion coefficient of impulse noise;
assuming that the underwater acoustic noise model is obtained by superposing Gaussian noise and non-Gaussian pulse noise models, defining the underwater acoustic noise e (i) as
e(i)＝eGauss(i)+eSαS(i)
Wherein eGauss(i)、eSαS(i) Are respectively composed of
Figure FDA0002987289830000025
And
Figure FDA0002987289830000026
and (4) generating.
2. The denoising method according to claim 1, wherein the step S2 is specifically as follows:
s2-1: and (3) noise point detection:
assume that the signal received at the receiving end is y ═ y (1), y (2),.., y (n)]The initial sliding window W has a length LWTaking out a sample W (i) corresponding to the time-th reception signal y excluding the central point y (i) by using the initial sliding window W, where 2n + 1:
w(i)＝[w1(i),w2(i),...,w2n(i)]
＝[y(i-n),...,y(i-1),y(i+1),...,y(i+n)]
ordering the signal points in w (i) from small to large
r(i)＝sort(w(i))
＝[r1(i),r2(i),...r2n(i)]
Wherein sort (·) is a ranking function; let Med mean (r (i)), mean () means taking the median; defining a differential noise identifier as
Figure FDA0002987289830000027
For a given predetermined pulse threshold TnoiseIf d (i) > TnoiseThen, y (i) is determined as the impulse noise point, andlet n (i) be 1, otherwise y (i) be the desired signal, and n (i) be 0, where n (i) represents the pulse marker; in the underwater sound receiving signal, the sound velocity is c, the amplitude is A, and the sampling frequency is fsCarrier frequency of fcThe rate of change of any adjacent sampling point will not exceed
Figure FDA0002987289830000028
And the sampling length of the underwater sound receiving signal is
Figure FDA0002987289830000031
Thereby setting the pulse threshold value to
Figure FDA0002987289830000032
S2-2: adaptive window size determination:
w length L for initial sliding windowWWhen the center point y (i) is an impulse noise point, calculating the number of noise points in the window:
Figure FDA0002987289830000033
the new window is marked as WnewThe length is as follows:
Figure FDA0002987289830000034
s2-3: noise filtering:
according to the new window size
Figure FDA0002987289830000035
And pulse mark n (i) filtering the received signal; wherein the non-noise points in the new window are unchanged and the noise points are replaced by the median value of the signal in the new window; suppose yip(i) Is WnewInner noise point, taking out new window and removing yip(i) All signal samples w ofnew(i)：
Figure FDA0002987289830000036
To wnew(i) The medium signal points are ordered from small to large to obtain
Figure FDA0002987289830000037
The noise point yip(i) Substituted by the formula:
y′ip(i)＝median(rnew(i))
wherein y'ip(i) Is the filtered signal.
3. The denoising method of claim 1, wherein the step S3 is specifically as follows:
to y 'first'ip(i) Performing wavelet transformation to obtain wavelet coefficient wj,kWherein j, k represents the kth coefficient of the jth layer; and constructing a new threshold function:
s3-1: a new threshold function construction:
Figure FDA0002987289830000038
wherein
Figure FDA0002987289830000039
Is an exponential factor, taking a non-negative number, lambdajA j-th layer threshold, wherein j is 1, 2.., and L is the number of decomposition layers;
s3-2: determining a threshold parameter to be optimized:
selecting a threshold lambda in the new adaptive threshold functionjAnd an exponential factor
Figure FDA0002987289830000044
The method is regarded as an unknown threshold parameter, and then optimized estimation is carried out on the threshold parameter based on a GDES-ABC method, so that the estimation precision and speed are improved;
s3-3: population initialization based on a good point set:
population initialization based on the good point set can effectively improve population diversity and avoid the method from falling into local optimum too early; the construction method of the good point comprises the following steps:
Figure FDA0002987289830000041
where p is the smallest prime number satisfying (p-3)/2 ≧ D, D is the dimension of the solution, deci {. cndot. } denotes the fractional part, rkIs a good point; thus, the set of good points [ PSN(1),PSN(2),...,PSN(SN)]TThe construction method of (2) is as follows:
PSN(i)＝{deci{r1*i},...,deci{rD*i}},i＝1,2,...,SN
wherein [. ]]TIndicating transposition, SN is the population size; the initial population is
X＝Lb+(Ub-Lb)*PSN
Wherein Ub and Lb are the upper and lower bounds of the solution, respectively, D is the current dimensional solution space, D is 1, 2.
S3-4: field search strategy guided by dynamic elite population:
the dynamic elite population contains better solutions in the population, and the scale of the solutions is different along with the difference of the iteration times; the dynamic smart English group-based neighborhood search can effectively accelerate the convergence speed of the method and improve the search efficiency; first calculate each Xi1,2, the fitness value of SN, and then taking the better Telite ceil (p' SN) bees to form a dynamic elite population DXEi′1, 2., Telite, where ceil (·) denotes rounding up and p' is the proportion of elite population in all populations, determined according to the following formula:
Figure FDA0002987289830000042
wherein p ismaxAnd pminRespectively representing the maximum value and the minimum value of p'; t is the current number of iterations, tmaxIs the maximum number of iterations; the improved neighborhood search method based on the dynamic elite population of the GDES-ABC method is as follows:
vid＝DXECd+φid(Gbestd-xkd),d＝1,2,...,D
wherein phi isidIs [ -1,1 [ ]]Random real numbers in between; gbest is the global optimal solution, xkdIs xidRandom neighborhood of (D) DXECdIs a dynamic elite population center, the expression is as follows:
Figure FDA0002987289830000043
the neighborhood search strategy based on the dynamic elite population is described as follows: for each neighborhood search, the hiring bee randomly searches neighborhoods with the same probability and generates new solutions from the improved neighborhood search method; observing the bee to randomly search neighborhoods from the elite population, and generating a new solution by an improved neighborhood searching method; for the observation bees, if the new solution is better than the current solution, selecting the new solution to perform next neighborhood search; otherwise, in the next neighborhood searching, randomly searching the neighborhood from the elite population again to perform the next neighborhood searching; the neighborhood search strategy is performed randomly;
s3-5: simulated annealing selection mechanism:
assume that at the T-th iteration, the current temperature is TtThe annealing parameter is K, and a new solution is obtained by an improved neighborhood search method and is VtWith a fitness value of fitv(ii) a The simulated annealing selection mechanism is as follows: if fitv＞fitiThen accept a new solution, where fitiIs the current solution fitness value; otherwise, accepting the new solution by using the probability P, wherein the acceptance probability is defined as
Figure FDA0002987289830000051
Wherein
Figure FDA0002987289830000052
Changing with the number of iterations; wherein beta is not more than 1 and is a constant with the value of 0.7 and sigmafitThe standard deviation of the fitness of all solutions;
from the acceptance probability, the method is initially characterized by small T and TtThe size is larger, so that P is larger, the method accepts a certain difference value, and the bee colony has larger development capability; at the end of the process, T becomes larger, TtDecreasing, P decreases, the method rejects the worse value with a greater probability;
s3-6: taking the mean square error between the signal to be denoised and the denoised signal as a fitness function of the improved artificial bee colony method in S3, and obtaining an optimal threshold parameter under the condition of obtaining the minimum mean square error; the fitness function is expressed as:
Figure FDA0002987289830000053
where s (i) is a training signal,
Figure FDA0002987289830000054
for reconstructing the training signal, N is the training signal length; known from the new adaptive threshold function, the threshold function is the parameter λjAnd
Figure FDA0002987289830000055
a function of (a); once λjAnd
Figure FDA0002987289830000056
determining, determining new adaptive threshold function of threshold function, and shrinking threshold to obtain wavelet coefficientReconstructing the denoised signal
Figure FDA0002987289830000057
And finally, performing contraction processing on the wavelet coefficient by constructing a new threshold function and an optimal threshold parameter to obtain a new wavelet coefficient, and performing inverse wavelet transformation to obtain a de-noising signal.
CN202010509512.6A 2020-06-23 2020-06-23 Underwater sound signal denoising method based on adaptive window filtering and wavelet threshold optimization Active CN111681667B (en)

Priority Applications (2)

Application Number Priority Date Filing Date Title
CN202010509512.6A CN111681667B (en) 2020-06-23 2020-06-23 Underwater sound signal denoising method based on adaptive window filtering and wavelet threshold optimization
PCT/CN2021/088635 WO2021258832A1 (en) 2020-06-23 2021-04-21 Method for denoising underwater acoustic signal on the basis of adaptive window filtering and wavelet threshold optimization

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010509512.6A CN111681667B (en) 2020-06-23 2020-06-23 Underwater sound signal denoising method based on adaptive window filtering and wavelet threshold optimization

Publications (2)

Publication Number Publication Date
CN111681667A CN111681667A (en) 2020-09-18
CN111681667B true CN111681667B (en) 2021-05-04

Family

ID=72454897

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010509512.6A Active CN111681667B (en) 2020-06-23 2020-06-23 Underwater sound signal denoising method based on adaptive window filtering and wavelet threshold optimization

Country Status (2)

Country Link
CN (1) CN111681667B (en)
WO (1) WO2021258832A1 (en)

Families Citing this family (25)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111681667B (en) * 2020-06-23 2021-05-04 青岛科技大学 Underwater sound signal denoising method based on adaptive window filtering and wavelet threshold optimization
CN112133321B (en) * 2020-09-23 2021-05-14 青岛科技大学 Underwater acoustic signal Gaussian/non-Gaussian noise suppression method based on blind source separation
CN112530449B (en) * 2020-10-20 2022-09-23 国网黑龙江省电力有限公司伊春供电公司 Speech enhancement method based on bionic wavelet transform
CN112633225B (en) * 2020-12-31 2023-07-18 矿冶科技集团有限公司 Mining microseism signal filtering method
CN113822329B (en) * 2021-08-10 2023-11-03 国网新源控股有限公司 Method and device for processing main shaft swing degree signal of hydroelectric generating set
CN114462458B (en) * 2022-04-11 2022-07-08 自然资源部第一海洋研究所 Ship underwater signal noise reduction and target enhancement method
CN114734440B (en) * 2022-04-15 2023-09-05 同济大学 Precise calibration method for kinematic parameters of hybrid double-arm transfer robot
CN115100864B (en) * 2022-06-24 2023-06-06 北京联合大学 Traffic signal control optimization method based on improved sparrow search algorithm
CN115953790B (en) * 2022-09-29 2024-04-02 江苏智联天地科技有限公司 Label detection and identification method and system
CN115816514B (en) * 2023-02-16 2023-06-27 极限人工智能有限公司 Robot joint brake fault detection method and system based on measurement of electric variables
CN116559421B (en) * 2023-04-03 2024-05-31 杭州臻稀生物科技有限公司 Automatic fluorescence immunoassay analyzer and analysis method
CN116304581B (en) * 2023-05-10 2023-07-21 佛山市钒音科技有限公司 Intelligent electric control system for air conditioner
CN116859711B (en) * 2023-06-27 2024-04-12 中科沃业江苏生物有限公司 IABC (advanced integrated circuit controller) -based fuzzy PID (proportion integration differentiation) optimized constant flow pump control system
CN116485884B (en) * 2023-06-28 2023-09-12 四川君安天源精酿啤酒有限公司 Real-time positioning method and system for finish brewing beer bottle mouth based on computer vision
CN116597856B (en) * 2023-07-18 2023-09-22 山东贝宁电子科技开发有限公司 Voice quality enhancement method based on frogman intercom
CN116776749B (en) * 2023-08-22 2023-11-17 吉林大学 Uncertainty-considered optimization method for wireless charging device of electric automobile
CN117370737B (en) * 2023-12-08 2024-02-06 成都信息工程大学 Unsteady state non-Gaussian noise removing method based on self-adaptive Gaussian filter
CN117395626B (en) * 2023-12-11 2024-02-09 厦门大学 Underwater acoustic network water quality monitoring data collection method based on meta learning and NOMA
CN117633554B (en) * 2023-12-14 2024-05-14 艾信智慧医疗科技发展（苏州）有限公司 Medical box type logistics transmission monitoring and early warning system
CN117668468B (en) * 2024-01-31 2024-04-26 山东省舜天化工集团有限公司 Intelligent analysis management system for chemical preparation data
CN117665788B (en) * 2024-02-01 2024-04-05 湖南科技大学 Noise processing method based on microwave measurement data
CN117786325B (en) * 2024-02-27 2024-04-30 山东晨晖电子科技有限公司 Ambient temperature wisdom early warning system of heavy-calibre thing networking water gauge
CN117828282B (en) * 2024-03-06 2024-06-04 山东泰霖信息工程有限公司 Data efficient processing method based on adaptive filtering
CN117973429B (en) * 2024-04-01 2024-06-07 南京信息工程大学 Model parameter ratio estimation method applied to non-Gaussian noise filtering
CN118114031A (en) * 2024-04-28 2024-05-31 长鹰恒容电磁科技(成都)有限公司 Radio waveform prediction method and system based on machine learning

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101322323A (en) * 2005-12-05 2008-12-10 艾利森电话股份有限公司 Echo detection
CN107800659A (en) * 2017-10-12 2018-03-13 西安电子科技大学 LFM signal modulation method for parameter estimation under Alpha Stable distritation noises
CN109063676A (en) * 2018-08-24 2018-12-21 广东石油化工学院 A kind of adaptive time-frequency method method and system for power signal
CN111144318A (en) * 2019-12-27 2020-05-12 苏州联视泰电子信息技术有限公司 Point cloud data noise reduction method for underwater sonar system

Family Cites Families (11)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CA2507901A1 (en) * 2004-05-21 2005-11-21 Imaging Dynamics Company Ltd. De-noising digital radiological images
JP5641186B2 (en) * 2010-01-13 2014-12-17 ヤマハ株式会社 Noise suppression device and program
US9053367B2 (en) * 2012-11-09 2015-06-09 Seiko Epson Corporation Detector evolution with multi-order contextual co-occurrence
EP3134853B1 (en) * 2014-04-25 2020-07-01 King Abdullah University Of Science And Technology System and method for image reconstruction, analysis, and/or de-noising
CN105528768A (en) * 2015-12-10 2016-04-27 国网四川省电力公司天府新区供电公司 Image denoising method
CN106530254A (en) * 2016-11-14 2017-03-22 山东理工大学 Algorithm for inhibiting mixed noise of images based on wavelet threshold function and improved median filtering fusion
CN107705260B (en) * 2017-10-03 2019-02-12 深圳市贝斯达医疗股份有限公司 The denoising system and calculation method of medical X-ray image
CN109816596B (en) * 2017-11-21 2020-12-22 中移（杭州）信息技术有限公司 Image denoising method and device
CN110515063A (en) * 2019-08-13 2019-11-29 青岛海洋科学与技术国家实验室发展中心 Underwater acoustic signal processing method and apparatus based on the steady wavelet transform of iteration
CN110765834B (en) * 2019-08-25 2020-07-17 青岛科技大学 Parameter wavelet threshold signal denoising method based on improved artificial bee colony algorithm
CN111681667B (en) * 2020-06-23 2021-05-04 青岛科技大学 Underwater sound signal denoising method based on adaptive window filtering and wavelet threshold optimization

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101322323A (en) * 2005-12-05 2008-12-10 艾利森电话股份有限公司 Echo detection
CN107800659A (en) * 2017-10-12 2018-03-13 西安电子科技大学 LFM signal modulation method for parameter estimation under Alpha Stable distritation noises
CN109063676A (en) * 2018-08-24 2018-12-21 广东石油化工学院 A kind of adaptive time-frequency method method and system for power signal
CN111144318A (en) * 2019-12-27 2020-05-12 苏州联视泰电子信息技术有限公司 Point cloud data noise reduction method for underwater sonar system

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
Characteristic Function Based Parameter Estimation for Ocean Ambient Noise;Xuebo Zhang;《2018 IEEE 3rd International Conference on Image, Vision and Computing (ICIVC)》;20181018;817-820 *
盲源分离技术在水声信号中的应用研究;王成刚;《舰船科学技术》;20150630;196-198 *

Also Published As

Publication number Publication date
CN111681667A (en) 2020-09-18
WO2021258832A1 (en) 2021-12-30

Similar Documents

Publication Publication Date Title
CN111681667B (en) Underwater sound signal denoising method based on adaptive window filtering and wavelet threshold optimization
Wang et al. A novel underwater acoustic signal denoising algorithm for Gaussian/non-Gaussian impulsive noise
CN110765834B (en) Parameter wavelet threshold signal denoising method based on improved artificial bee colony algorithm
CN112735456B (en) Speech enhancement method based on DNN-CLSTM network
CN110490816B (en) Underwater heterogeneous information data noise reduction method
CN115442191B (en) Communication signal noise reduction method and system based on relative average generation countermeasure network
JP5115952B2 (en) Noise suppression device and noise suppression method
CN113222854B (en) Wavelet transform CT image denoising method
CN113238190A (en) Ground penetrating radar echo signal denoising method based on EMD combined wavelet threshold
CN114446314A (en) Voice enhancement method for deeply generating confrontation network
CN110211598A (en) Intelligent sound noise reduction communication means and device
CN114842867A (en) DFT-based audio sinusoidal signal frequency estimation method and system
CN116203634A (en) Ghost wave removing method based on low-rank constraint
CN116520333A (en) Multi-base buoy fusion detection method and system based on complex interference
Sun et al. Wavelet denoising method based on improved threshold function
CN115065578B (en) DFT channel estimation method based on improved self-adaptive threshold
CN114299938B (en) Intelligent voice recognition method and system based on deep learning
CN115828120B (en) Self-adaptive identification method, system and computer equipment for ship traffic behavior mode
CN117493776B (en) Geophysical exploration data denoising method and device and electronic equipment
CN116757969B (en) Image blind denoising method and system based on self-adaptive curvature feature fusion
CN116070094B (en) Underwater sound signal processing method based on adaptive wavelet threshold function
CN117527570B (en) Sensor cluster position optimization method based on edge reinforcement learning
CN114936622B (en) Underwater sound target identification method and device based on cyclic generation countermeasure network
CN118016088A (en) Whale signal enhancement method based on learning difference
CN118248158A (en) Conference scene non-stationary noise elimination method based on wavelet threshold optimization

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant