A Narrowband Sound Signal Frequency Estimation with Impulsive Noise Filtering Iurii Chyrka Mathematical Methods for Sensor Information Processing Department of Institute of Information and Communication Technologies of Bulgarian Academy of Sciences 25 A, Acad. G. Bonchev str., 1113 Sofia, Bulgaria Yurasyk88@gmail.com Abstract estimation in such situation. Hence there is need to recover an original data from the degraded observations The new stochastic approach of impulsive before main processing stage. noise filtering, based on the input sample Filtering of an impulsive noise generally gives separation by modified clusterization criterion, positive effect not only on frequency estimation, but is proposed. It is shown that preliminary also for NSI algorithms. The acoustical array is a filtration by the proposed procedure provides multichannel system and almost simultaneous robust narrowband sound frequency estimation occurrence of impulses in many different channels can and eliminates failures of the estimation completely corrupt calculation results. algorithm caused by the impulsive noise. 2. Problem Statement 1. Introduction In view of NSI application, attention in the paper is Digital audio systems are widely applied today in paid to a class of narrowband signals which parameters many areas of human activity. One of the biggest class are changed slowly in time. The problem of the of them are acoustic arrays that perform measurement instantaneous frequency estimation can be interpreted and processing of sound field. One of the key branch in as estimation on the limited observation interval digital sound processing is Noise Source Identification (usually less than two periods of signal) during which (NSI) techniques. They play an important role in the the parameters are changed slowly and a narrowband acoustic camera, which is aimed to locate and signal is considered as a harmonic one. characterize sound sources. It fuses images, received by For description of a mixture of a digital narrowband the sensors in a complex image. This image consists of signal si with a white Gaussian noise ηi of power σ2 and a camera image as a background and contour lines impulsive noise ζ the following typical additive model describing sound field as a foreground [Bil76]. of a data sample is used: Most of NSI techniques have a “narrowband” nature. xi = si + ηi + υ i ζ = ρ i sin(ωi τ(i − 1) + ϕ 0 ) + ηi + υi ζ , This means that the sound map must be recalculated for each frequency of interest. Usually calculations are i = 1, N , (1) performed for the predefined set of bounds and where υi is a sign function υi =sgn(pi), that can get corresponding center frequencies. The difference values –1, 1, 0 with corresponding probabilities pζ/2, between actual frequency and defined one leads to pζ/2, 1–pζ/2; pζ is an impulsive noise appearance inaccurate results. Therefore, knowledge of the original probability; ζ is an impulsive noise amplitude signal frequency is important for precise estimation of (considered as constant for all impulses); ρi is a signal sound field parameters. instantaneous amplitude; ϕ0 is an initial phase; ωi is an The estimate can be easily obtained in the steady instantaneous angular frequency; τ is a sampling state from long enough Fast Fourier Transform, but in interval, N is a sample size. Further, the normalized non-stationary case it must be estimated instantly in frequency γ=ωτ is used to omit τ. Such the mixture real time. For fast, effective and precise estimation of model is close to a Bernoulli–Gaussian model of an any instantaneous parameter, the data sample must be impulsive noise process [Vas08]. as short as possible; therefore, appearance of impulsive The mixture sample (1) can be represented by the next noise has the most substantial influence on the probability density function f Ξ (x n | γ , ρ, ϕ 0 , σ, p ζ , ζ , υ n ) = (1 − p ζ ) ⋅ f η (x n | γ , ρ, ϕ 0 , σ ) Copyright © 2015 for the individual papers by the papers' authors. (2) Copying permitted only for private and academic purposes. + p ζ ⋅ f ηζ (x n | γ , ρ, ϕ 0 , σ, ζ , υ n ). In: A. Bădică, M. Colhon (eds.): Proceedings of the 2015 Balkan It corresponds to a Tukey-Huber model, which is the Conference on Informatics: Advances in ICT mostly used one for investigation of robust methods. It 40 supposes that majority of sample counts have an 3.2 Impulse noise filtering techniques expected distribution fη and some of counts belong to fηζ. The latter ones are outliers that produce “tails” of Conventional global filtering approach like a low- the distribution. The aim of the paper is describing of pass filtering assumes that both corrupted and the robust instantaneous frequency estimation uncorrupted samples must be processed. Median filters procedure with removing of distorting outliers in (2). 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 3. Methodological Basis signal [Vas08]. In addition, median filter eliminates changes in the input signal with a duration less than a 3.1 Brief description of NSI techniques half size of the filter window and does not properly filter a set of consequent impulses longer than a half The NSI techniques can be divided onto two size. Some modifications of the classic median filter, categories: near-field acoustic holography (NAH) and that eliminate some its disadvantages, were developed beamforming [Bai13]. NAH is aimed to reconstruction [Geo11]. of a sound field in the 3D space. Beamforming gives a Some detection methods perform processing in time map of a sound intensity by measurement of the signal and frequency domains simultaneously, for example response from a variety of directions. The main concept wavelet-based method in [Non08]. Yet another filtering of the both NSI techniques is the next: the sound approach uses fuzzy impulse detection, but mostly for pressures measured by the microphones (more rarely – images [Sch06]. particle velocities, captured by probes) are processed Impulsive noise usually distorts a relatively small by an imaging algorithm of either type to calculate an amount of total counts in the sample. Since usually a acoustic map of the sound pressure or sound intensity relatively large part of the signal counts remain with a snap to physical coordinates. unaffected by the impulsive noise. Hence, it is When NAH performs near-field imaging of noise advantageous to replace only the noisy ones, leaving sources, beamforming generally works in far field. the uncorrupted counts unchanged. This ideology is Unlike beamforming, that carries out spatial filtering implemented by the another approach that performs and maximizes the signal power from certain direction, model-based two-stage filtering by using a linear NAH provides, based on measurements over a two- prediction system [Esq02]. For audio signals, the most dimensional aperture, a reconstruction of the three- often used models are autoregressive (AR) or dimensional sound field from the source’s boundary out autoregressive moving average (ARMA) [Oud14]. In to the far field. Precision of reconstruction depends this case, the system consists of two main parts: mostly on microphone spacing distance, sound detector and interpolator that perform individual frequency and distance between source surface and processing of the each element of the data set. measurement plane. NAH operates in a low frequency Another independent class of methods is stochastic range, upper boundary of which is limited by a distance ones [Bas88]. They are close to model-based ones, but between microphones. On the other hand, beamforming the statistical properties of data samples are used. A works in the full frequency range but its use is similar approach is proposed in the paper for pre- reasonable only at high frequencies when it gives better filtering before frequency estimation. resolution than holography. Beamforming provide the best performance on the irregular arrays, which geometry (positions of 3.3 Frequency Estimation microphones) is optimized in order to get the lowest In this paper frequency estimation is carried out with possible level of false responses. As opposite to using of the algorithm mentioned in the previous work beamforming, NAH is usually performed on a regular [Pro12]. It has been synthesized with using an AR array grid that is also the requirement for the classic model for a single-tone harmonic signal mixed with a Fourier NAH. Many modern NAH approaches perform white Gaussian noise: calculation in the time domain and usually do not si = αsi −1 − si− 2 = 2 cos( γ ) si −1 − si− 2 , (3) require location of microphones on some equidistant positions. Hence, it allows reconstructing of sound where α=2cos(γ) is a parameter of auto-regression. fields even with irregular array geometries. The algorithm is based on the solution of the quadratic equation 41 α 2 − 2Bα − 2 = 0 , (4) minimal sum of cluster variances, which for a two- with coefficient B calculated on the basis of the input component sample can be written as: signal sample (x ) as: { ( ) ( )} y ( opt ) = arg min σ l2 y ( v −1), N + σ r2 y ( v ), N . v (6) (xi+1 + xi −1 )2 − 2xi2 / 2 (xi+1 xi + xi xi−1 ) .(5) N −1 N −1 B(x ) = ∑[ i =2 ] ∑  i =2  In order to increase sensitivity to presence of the second “impulsive” subsample the modified criterion is The equation (4) has two roots considered: α ∗ ( + , − ) = B ± B 2 + 2 . Finally, frequency is estimated as γ * = arccos (α ∗( + ) / 2) . { ( ) ( )} y ( opt ) = arg min σ l4 y ( v −1), N + σ r4 y ( v ), N . v (7) Generally, this algorithm is robust in many cases It allows correct treatment of mixtures that contain with only the impulsive noise and without the Gaussian small amount of impulses. one [Pro09], but it becomes highly sensitive when Finally, the proposed impulse filtering procedure can Gaussian and impulsive noises occur simultaneously be represented in the form of four consequent steps of and are multiplied during calculations. Therefore, calculations. On the first step of processing in order to removing of impulses is a necessary condition before facilitate detection of the sample with bipolar pulse it is estimation procedure starts. necessary to take the absolute value z n = xn , n = 0, N − 1 . Taking into account that time order of 4. Impulsive noise detection discrete values (1) is insignificant, it can be () transformed via y = ℜ z into an ordered statistics The NSI Robust parameters estimation usually y = ( y(1), N , y ( 2 ), N ,... y( N ), N ) , (8) consists of the next processing stages [Hub09]: 1) Data “errors” or corrupted counts detection; which must be separated onto two parts with 2) Processing of detected counts by removing them corresponding standard deviation values: σl, σr. from the sample or their restoration with the help of The next step includes the calculation of values of neighbor ones; the above mentioned criterion for each discrete rank 3) Pure robust estimation using the restored sample. (v),N. (the function (7) is used here) The separation Here steps 1) – 2) belong to impulse noise filtering. threshold is simply found as argmin among all these The new stochastic approach based on the cluster values. When the threshold is obtained, all counts that analysis theory [Eve11] is proposed for performing exceed its value are marked as defective or i. e. impulses detection. The key idea is to apply containing a noise impulse. clusterization methods for analysis of the input sample On the last stage of processing the detected impulses probability distribution and separate all counts in the must be replaced by interpolation or extrapolation sample onto (two) independent clusters: the subsample using “good” neighbor counts. In view of assumption with normal counts and relatively small subsample with on small duration of the impulses, only the simplest impulsive noise counts. The important assumption here linear interpolation is used in this paper. is that the signal amplitude is relatively small in comparison to impulse amplitude, that allows to 5. Simulation results distinguish these clusters. The similar approach was The effectiveness of impulsive noise filtering was proposed earlier by the author in application to analyzed in connection to its influence on the frequency electroencephalogram signals analysis [Pro13]. estimation process by the aforementioned algorithm. Generally, clusterization is carried out on the basis The flowchart of the frequency estimation process of some optimization criterion or an objective function. including the impulse filtering one can see in Fig. 1. Usually it is an extent of objects density inside the Statistical simulations by the Monte-Carlo approach cluster or an extent of distance between different were done under the next conditions: a signal sample clusters. Actually, most of cluster separation size N=50, the sample contains one period of the signal, procedures can be considered as exact or approximate hence the normalized frequency γ = 2π / 50 ≈ 0.126 , algorithm of some objective function optimization and SNR=20 dB, signal to impulsive noise ratio is –20 dB, finding a threshold. number of numerical simulations for each plot is The most widespread methods for optimal threshold 10000. v calculation [Kor89] are based on the criterion of the 42 Figure 2: Comparison of algorithm failure probabilities for different filtration methods 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 Figure 3: Comparison of precision indicators of high enough. Longer filter window allows getting additional suppression of many impulses but there must frequency estimations for different filtration methods be tradeoff between noise filtering and the signal Similarly to situation with failures, the estimation deformation due influence of the filter. The proposed algorithm does not give reliable estimates, when the separation procedure for impulse noise filtration makes probability of noise appearance is high. The median failures caused by impulses action almost impossible. filter reduces the deviation of a mean in about two In addition, the comparison of the precision of times. The separation procedure provides a mean and a frequency estimation was carried out in three standard deviation close to constant values. aforementioned cases without any filtering, with The carried out statistical studies have shown that filtration by the classic median filter and by using the the proposed method works well with relatively big separation. The plots of mean and standard deviation impulses, when signal to impulsive noise ratio is –10 are shown in Fig. 3. 43 dB or lower. When the amplitude of impulses is low, USA, March 31 - April 4, 2008). ICASSP '08. two clusters of good and corrupted counts are located IEEE, Las Vegas, NV, 1593–1596, 2008. very close to each other and the polymodality of [Sch06] S. Schulte, M. Nachtegael, V. De Witte, D. impulse distribution vanishes. It is hard to find the threshold in this case. Van der Weken, E. Kerre. A Fuzzy Impulse Noise Detection and Reduction Method. IEEE Transactions on Image Processing. 15(5): 6. Conclusions 1153–1162, May 2006. The proposed impulse noise filtering procedure [Esq02] P. Esquef, M. Karjalainen and V. Valimaki. provides better frequency estimation quality in Detection of clicks in audio signals using comparison to conventional median filter when the warped linear prediction. In Proceedings of impulse to signal amplitude ratio is 10 dB or more. the 14th International Conference on Digital This is caused by better sensitivity of the procedure to appearance of big impulses in the sample. At the same Signal Processing (New York, USA, March time, it provides constant estimation precision, when 31 - April 4, 2008). ICDSP '02. IEEE, New impulse appearance probability is not bigger than 0.15. York, NY, 2: 1085–1088, 2002. This is a consequence of bigger number of points that [Oud14] L. Oudre. Automatic detection and removal of have to be interpolated. Therefore, linear interpolation impulsive noise in audio signals. Image distorts the signal structure and worsen frequency Processing On Line. Preprint. February 2014. estimation at bigger appearance probabilities. [Bas88] M. Basseville. Detecting Changes in Signals and Systems-A Survey. Automatica. 24(3): Acknowledgments 309–326, May 1988. The research work reported in the paper was partly [Pro12] I. G. Prokopenko, I. P. Omelchuk, and Y. D. supported by the Project AComIn "Advanced Chyrka. Radar signal parameters estimation in Computing for Innovation"; grant 316087, funded by the MTD tasks. International Journal of the FP7 Capacity Programme. Electronics and Telecommunications (JET). 58(2): 159–164, June 2012. References [Pro09] I. G. Prokopenko, I. P. Omelchuk and G. Y. Sokolov. Robustness of quasioptimal [Bil76] J. Billingsley, R. Kinns. 1976. The acoustic frequency estimator to impulsive noise. telescope, Journal of Sound and Vibration, 48, Electronics and Control Systems. 20(2): 69– 485–510, 1976. 74, 2009. (In Ukrainian). [Vas08] Saeed V. Vaseghi. Advanced Digital Signal [Hub09] P. Huber, E. Ronchett. Robust Statistics – 2nd Processing and Noise Reduction – 4th ed. A ed. A John Wiley and Sons Ltd, 18–20, 2009. John Wiley and Sons Ltd., Singapore, 2008. [Eve11] B. S. Everitt, S. Landau, M. Leese, D. Stahl. [Bai13] R. M. Bai, J.G. Ih, J. Benesty. Acoustic Array Cluster Analysis – 5th ed. A John Wiley and Systems. John Wiley & Sons Singapore Pte. Sons Ltd., London. 2011. Ltd., Singapore, 2013. [Pro13] I. G. Prokopenko, I. P. Omelchuk, and Y. D. [Geo11] Geoffrine Judith M. C., N. Kumarasabapathy. Chyrka. Pat. № 84833, UA, IPC G06F7/06. Study and analysis of impulse noise reduction Method for electroencephalogram spike filters. Signal & Image Processing: An detection. Publ. date: 11.11.2013 (In International Journal (SIPIJ). 2(1): 82–92, Ukrainian). March 2011. [Kor89] E. A. Korniliev, I. G. Prokopenko, V. M. [Non08] R. C. Nongpiur. Impulse noise removal in Chuprin. Robust algorithms in automated speech using wavelets. In Proceedings of the systems of information processing. Technics, IEEE International Conference on Acoustics, Kyiv, 1989 (In Russian). Speech and Signal Processing (Las Vegas, 44