Introduction

Digital audio systems are widely applied today in many areas of human activity. One of the biggest class of them are acoustic arrays that perform measurement and processing of sound field. One of the key branch in digital sound processing is Noise Source Identification (NSI) techniques. They play an important role in the acoustic camera, which is aimed to locate and characterize sound sources. It fuses images, received by the sensors in a complex image. This image consists of a camera image as a background and contour lines describing sound field as a foreground [Bil76].

Most of NSI techniques have a "narrowband" nature. This means that the sound map must be recalculated for each frequency of interest. Usually calculations are performed for the predefined set of bounds and corresponding center frequencies. The difference between actual frequency and defined one leads to inaccurate results. Therefore, knowledge of the original signal frequency is important for precise estimation of sound field parameters.

The estimate can be easily obtained in the steady state from long enough Fast Fourier Transform, but in non-stationary case it must be estimated instantly in real time. For fast, effective and precise estimation of any instantaneous parameter, the data sample must be as short as possible; therefore, appearance of impulsive noise has the most substantial influence on the estimation in such situation. Hence there is need to recover an original data from the degraded observations before main processing stage.

Filtering of an impulsive noise generally gives positive effect not only on frequency estimation, but also for NSI algorithms. The acoustical array is a multichannel system and almost simultaneous occurrence of impulses in many different channels can completely corrupt calculation results.

Problem Statement

In view of NSI application, attention in the paper is paid to a class of narrowband signals which parameters are changed slowly in time. The problem of the instantaneous frequency estimation can be interpreted as estimation on the limited observation interval (usually less than two periods of signal) during which the parameters are changed slowly and a narrowband signal is considered as a harmonic one. For description of a mixture of a digital narrowband signal s i with a white Gaussian noise η i of power σ 2 and impulsive noise ζ the following typical additive model of a data sample is used:

ζ υ + η + ϕ + − τ ω = ζ υ + η + = i i 0 i i i i i i ) ) 1 ( sin( ρ i s x , N i , 1 = ,(1)

where υ i is a sign function υ i =sgn(p i ), that can get values -1, 1, 0 with corresponding probabilities p ζ /2, p ζ /2, 1-p ζ /2; p ζ is an impulsive noise appearance probability; ζ is an impulsive noise amplitude (considered as constant for all impulses); ρ i is a signal instantaneous amplitude; ϕ 0 is an initial phase; ω i is an instantaneous angular frequency; τ is a sampling interval, N is a sample size. Further, the normalized frequency γ=ωτ is used to omit τ. Such the mixture model is close to a Bernoulli-Gaussian model of an impulsive noise process [Vas08]. The mixture sample (1) can be represented by the next probability density function

( ) ( ) ( ) ( ). , ζ σ, , , , | x σ , , , | x 1 , ζ , σ, , , , | x 0 ηζ ζ 0 η ζ ζ 0 n n n n n f p f p p f υ ϕ ρ γ ⋅ + ϕ ρ γ ⋅ − = υ ϕ ρ γ Ξ (2)

It corresponds to a Tukey-Huber model, which is the mostly used one for investigation of robust methods. It supposes that majority of sample counts have an expected distribution f η and some of counts belong to f ηζ . The latter ones are outliers that produce "tails" of the distribution. The aim of the paper is describing of the robust instantaneous frequency estimation procedure with removing of distorting outliers in (2).

Methodological Basis

Brief description of NSI techniques

The NSI techniques can be divided onto two categories: near-field acoustic holography (NAH) and beamforming [Bai13]. NAH is aimed to reconstruction of a sound field in the 3D space. Beamforming gives a map of a sound intensity by measurement of the signal response from a variety of directions. The main concept of the both NSI techniques is the next: the sound pressures measured by the microphones (more rarelyparticle velocities, captured by probes) are processed by an imaging algorithm of either type to calculate an acoustic map of the sound pressure or sound intensity with a snap to physical coordinates.

When NAH performs near-field imaging of noise sources, beamforming generally works in far field. Unlike beamforming, that carries out spatial filtering and maximizes the signal power from certain direction, NAH provides, based on measurements over a twodimensional aperture, a reconstruction of the threedimensional sound field from the source's boundary out to the far field. Precision of reconstruction depends mostly on microphone spacing distance, sound frequency and distance between source surface and measurement plane. NAH operates in a low frequency range, upper boundary of which is limited by a distance between microphones. On the other hand, beamforming works in the full frequency range but its use is reasonable only at high frequencies when it gives better resolution than holography.

Beamforming provide the best performance on the irregular arrays, which geometry (positions of microphones) is optimized in order to get the lowest possible level of false responses. As opposite to beamforming, NAH is usually performed on a regular array grid that is also the requirement for the classic Fourier NAH. Many modern NAH approaches perform calculation in the time domain and usually do not require location of microphones on some equidistant positions. Hence, it allows reconstructing of sound fields even with irregular array geometries.

Impulse noise filtering techniques

Conventional global filtering approach like a lowpass filtering assumes that both corrupted and uncorrupted samples must be processed. Median filters and other order statistics filters, that process a localized area, also typically modify uncorrupted samples as the transversal filtering is applied uniformly over the whole signal [Vas08]. In addition, median filter eliminates changes in the input signal with a duration less than a half size of the filter window and does not properly filter a set of consequent impulses longer than a half size. Some modifications of the classic median filter, that eliminate some its disadvantages, were developed [Geo11].

Some detection methods perform processing in time and frequency domains simultaneously, for example wavelet-based method in [Non08]. Yet another filtering approach uses fuzzy impulse detection, but mostly for images [Sch06].

Impulsive noise usually distorts a relatively small amount of total counts in the sample. Since usually a relatively large part of the signal counts remain unaffected by the impulsive noise. Hence, it is advantageous to replace only the noisy ones, leaving the uncorrupted counts unchanged. This ideology is implemented by the another approach that performs model-based two-stage filtering by using a linear prediction system [Esq02]. For audio signals, the most often used models are autoregressive (AR) or autoregressive moving average (ARMA) [Oud14]. In this case, the system consists of two main parts: detector and interpolator that perform individual processing of the each element of the data set.

Another independent class of methods is stochastic ones [Bas88]. They are close to model-based ones, but the statistical properties of data samples are used. A similar approach is proposed in the paper for prefiltering before frequency estimation.

Frequency Estimation

In this paper frequency estimation is carried out with using of the algorithm mentioned in the previous work [Pro12]. It has been synthesized with using an AR model for a single-tone harmonic signal mixed with a white Gaussian noise:

2 i 1 i 2 i 1 i i ) cos( 2 − − − − − γ = − α = s s s s s ,(3)

where α=2cos(γ) is a parameter of auto-regression.

The algorithm is based on the solution of the quadratic equation 0

2 2 2 = − α − α B ,(4)

with coefficient B calculated on the basis of the input signal sample ( )

x as: ( ) ( ) [ ] ( )       + − + = ∑ ∑ − = − + − = − + 1 2 1 1 1 2 2 2 1 1 2 / 2 N i i i i i N i i i i x x x x x x x x B .(5)

The equation (4) has two roots 2 2 ) , (

+ ± = α − + * B B

. Finally, frequency is estimated as

( ) 2 / arccos ) ( * + * α = γ

. Generally, this algorithm is robust in many cases with only the impulsive noise and without the Gaussian one [Pro09], but it becomes highly sensitive when Gaussian and impulsive noises occur simultaneously and are multiplied during calculations. Therefore, removing of impulses is a necessary condition before estimation procedure starts.

Impulsive noise detection

The NSI Robust parameters estimation usually consists of the next processing stages [Hub09]:

1) Data "errors" or corrupted counts detection; 2) Processing of detected counts by removing them from the sample or their restoration with the help of neighbor ones;

3) Pure robust estimation using the restored sample.

Here steps 1) -2) belong to impulse noise filtering.

The new stochastic approach based on the cluster analysis theory [Eve11] is proposed for performing impulses detection. The key idea is to apply clusterization methods for analysis of the input sample probability distribution and separate all counts in the sample onto (two) independent clusters: the subsample with normal counts and relatively small subsample with impulsive noise counts. The important assumption here is that the signal amplitude is relatively small in comparison to impulse amplitude, that allows to distinguish these clusters. The similar approach was proposed earlier by the author in application to electroencephalogram signals analysis [Pro13].

Generally, clusterization is carried out on the basis of some optimization criterion or an objective function. Usually it is an extent of objects density inside the cluster or an extent of distance between different clusters. Actually, most of cluster separation procedures can be considered as exact or approximate algorithm of some objective function optimization and finding a threshold.

The most widespread methods for optimal threshold v calculation [Kor89] are based on the criterion of the minimal sum of cluster variances, which for a twocomponent sample can be written as:

( ) ( ) { } N v r N v l v opt y y y ),( 2 ), 1 ( 2 )( min arg σ + σ = − .(6)

In order to increase sensitivity to presence of the second "impulsive" subsample the modified criterion is considered:

( ) ( ) { } N v r N v l v opt y y y ),( 4 ), 1 ( 4 )( min arg σ + σ = − .(7)

It allows correct treatment of mixtures that contain small amount of impulses.

Finally, the proposed impulse filtering procedure can be represented in the form of four consequent steps of calculations. On the first step of processing in order to facilitate detection of the sample with bipolar pulse it is necessary to take the absolute value

n n x z = , 1 , 0 − = N n

. Taking into account that time order of discrete values (1) is insignificant, it can be transformed via ( )

z y ℜ = into an ordered statistics ( ) ( ) ( ) ( ) N N N N y y y y , , 2 , 1 ,... , = ,(8)

which must be separated onto two parts with corresponding standard deviation values: σ l , σ r .

The next step includes the calculation of values of the above mentioned criterion for each discrete rank (v),N. (the function ( 7) is used here) The separation threshold is simply found as argmin among all these values. When the threshold is obtained, all counts that exceed its value are marked as defective or i. e. containing a noise impulse.

On the last stage of processing the detected impulses must be replaced by interpolation or extrapolation using "good" neighbor counts. In view of assumption on small duration of the impulses, only the simplest linear interpolation is used in this paper.

Simulation results

The effectiveness of impulsive noise filtering was analyzed in connection to its influence on the frequency estimation process by the aforementioned algorithm. The flowchart of the frequency estimation process including the impulse filtering one can see in Fig. 1.

Statistical simulations by the Monte-Carlo approach were done under the next conditions: a signal sample size N=50, the sample contains one period of the signal, hence the normalized frequency 126 . 0 50 / 2 ≈ π = γ , SNR=20 dB, signal to impulsive noise ratio is -20 dB, number of numerical simulations for each plot is 10000.

Figure 1: Flowchart of the frequency estimation process Fig. 2 shows plots of dependencies of estimation algorithm failure probability on impulsive noise appearance probability. Three cases are regarded: without filtering, with filtration by a classic median filter and by using the proposed separation. From the presented figure one can see that the appearance of the impulsive noise with big power significantly worsen performance of the frequency estimation algorithm. Small decreasing of failure rate at higher noise probability can be explained by the fact that some impulses start to appear one after another and such a sequence has less influence on the algorithm.

The median filter of size 5 makes a situation better and substantially decreases probability of failure, at low appearance probability in particular, but its rate is still high enough. Longer filter window allows getting additional suppression of many impulses but there must be tradeoff between noise filtering and the signal deformation due influence of the filter. The proposed separation procedure for impulse noise filtration makes failures caused by impulses action almost impossible.

In addition, the comparison of the precision of frequency estimation was carried out in three aforementioned cases without any filtering, with filtration by the classic median filter and by using the separation. The plots of mean and standard deviation are shown in Fig. 3. Similarly to situation with failures, the estimation algorithm does not give reliable estimates, when the probability of noise appearance is high. The median filter reduces the deviation of a mean in about two times. The separation procedure provides a mean and a standard deviation close to constant values.

The carried out statistical studies have shown that the proposed method works well with relatively big impulses, when signal to impulsive noise ratio is -10 dB or lower. When the amplitude of impulses is low, two clusters of good and corrupted counts are located very close to each other and the polymodality of impulse distribution vanishes. It is hard to find the threshold in this case.

Conclusions

The proposed impulse noise filtering procedure provides better frequency estimation quality in comparison to conventional median filter when the impulse to signal amplitude ratio is 10 dB or more. This is caused by better sensitivity of the procedure to appearance of big impulses in the sample. At the same time, it provides constant estimation precision, when impulse appearance probability is not bigger than 0.15. This is a consequence of bigger number of points that have to be interpolated. Therefore, linear interpolation distorts the signal structure and worsen frequency estimation at bigger appearance probabilities.