On a nonlocal spectrogram for denoising one-dimensional signals

Applied Mathematics and Computation 244 (2014) 859–869

Contents lists available at ScienceDirect

Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc

On a nonlocal spectrogram for denoising one-dimensional signals q G. Galiano ⇑, J. Velasco Dpt. of Mathematics, Universidad de Oviedo, c/ Calvo Sotelo, 33007-Oviedo, Spain

a r t i c l e

i n f o

Keywords: Spectrogram Image restoration Local and nonlocal filtering Reassigned spectrogram Instantaneous frequency Biomedical signals

a b s t r a c t In previous works, we investigated the use of local filters based on partial differential equations (PDE) to denoise one-dimensional signals through the image processing of time– frequency representations, such as the spectrogram. In these image denoising algorithms, the particularity of the image was hardly taken into account. We turn, in this paper, to study the performance of nonlocal filters, like Neighborhood or Yaroslavsky filters, in the same problem. We show that, for certain iterative schemes involving the Neighborhood filter, the computational time is drastically reduced with respect to Yaroslavsky or nonlinear PDE based filters, while the outputs of the filtering processes are similar. This is heuristically justified by the connection between the fast Neighborhood filter applied to a spectrogram and the corresponding accurate reassigned spectrogram of the image. This correspondence holds only for time–frequency representations of one-dimensional signals, not for usual images, and in this sense the particularity of the image is exploited. We compare through a series of experiments on synthetic and biomedical signals the performance of local and nonlocal filters. Ó 2014 Elsevier Inc. All rights reserved.

1. Introduction Denoising one-dimensional signals is an important topic which is usually addressed from filter theory in time or frequency domains. In some applications in which processing speed is not a fundamental issue, filters defined in the joint time–frequency domain may be considered, usually improving the filtering process. Examples of this situation are found in Electrocardiogram (ECG) and other biomedical signals [23], human voice analysis [25] or animal sound analysis [13]. With respect to the latter, in previous works [7–9] we investigated the use of time–frequency distributions to estimate the number of wolves howling in a giving recording to provide an estimation of the number of individuals in a pack. This estimation is the basis for counting regional wolf populations which is of interest for both ecological and economic purposes, since authorities must reimburse the cost of cattle killed by this protected specie [21]. Of course, and despite the quality of recording devices, field recordings are usually affected for a variety of undesirable signals which range from low amplitude broad spectrum long duration signals, like wind, to signals localized in time, like cattle bells, or localized in spectrum, like car engines. Clearly, the addition of all these signals generates an unstructured noise in the background of the wolves chorus which must be treated for a proper signal analysis.

q

Supported by Spanish DGI Project MTM2010-18427.

⇑ Corresponding author.

E-mail addresses: [email protected] (G. Galiano), [email protected] (J. Velasco). http://dx.doi.org/10.1016/j.amc.2014.07.003 0096-3003/Ó 2014 Elsevier Inc. All rights reserved.

860

G. Galiano, J. Velasco / Applied Mathematics and Computation 244 (2014) 859–869

Medical signals are another good example of this situation. Due to the electromagnetic fields created by measuring devices, the usual low frequency signals to be acquired are contaminated by a background noise which is usually in the same frequency band that the signal of interest. Therefore, fine denoising techniques must be applied to segregate the signal of interest from the noise. In general, the denoising procedure is not aimed to recovering a clean signal but to produce a clean time–frequency representation of the signal which allows further analysis techniques. For instance, for the first example given above, the instantaneous frequency (IF) estimation allows counting the number of different individuals howling in the recording (each individual being identified with an IF line). In this article, we also provide an example in which the filtering of the spectrogram energy content of an ECG signal is used to identify an arrhythmia episode. In [7,8,10], we considered nonlinear diffusion image denoising techniques applied to the spectrogram of a sound signal, a wolf chorus. Although, as above mentioned, execution time is not the main issue for this type of problems, we found that nonlinear diffusion algorithms require a high computational time, making their use not operative in many situations. In addition, these filters do not take advantage of the special characteristics of the image produced from the spectrogram, i.e. they operate on the spectrogram as in any other image. In this article we show that nonlocal filters, specially the Neighborhood filter [15], are computationally more efficient to deal with these images and give similar results. Moreover, we point out a relationship between nonlocal image filtering and an effective technique for one-dimensional signals known as reassigned spectrogram, see [14,3]. In few words, the reassigned spectrogram operator takes the energy contributions represented by the spectrogram, initially placed in the geometric centers of the smoothing window used to define the short time Fourier transform of the signal, and reallocates these energy contributions to their gravity centers, producing in this way an energy concentration on the IF lines. As we indicate in Section 2.1, a similar result is obtained when applying nonlocal filters to the spectrogram of the signal, thus demonstrating the singular effect of nonlocal filtering on this kind of images. The outline of the article is as follows. We present in Section 2 the mathematical framework of the problem and the filtering techniques proposed in this article for one-dimensional signal denoising. In particular, we justify our choice of the Neighborhood filter as an inexpensive approximation to other well known denoising algorithms, while producing similar results on spectrogram images. We also give an argument linking nonlocal filtering to the reassigned spectrogram, providing some insight on the special effect of these filters when applied to spectrogram images. In Section 3, we introduce the discrete problem and deduce the corresponding formulas for algorithm implementation. Apart from the Neighborhood filter, we consider the Yaroslavsky-SUSAN [26,22] filter and a nonlinear diffusion filter based in the Total Variation norm [19,10], for comparison purposes. Then, we demonstrate the performance of these filters by applying them to three noisy signals (synthetic, wolf chorus and ECG) and give quantitative comparisons based on the Mean Square Error (MSE), and the visual inspection of the processed spectrograms and other related magnitudes. Finally, in Section 4 we give our conclusions. 2. Mathematical framework Let f 2 L2 ðRÞ denote a one-dimensional signal and WVðf ; ; Þ be its Wigner–Ville distribution, defined as

Z

 s  s f t þ f t  eiws ds; 2 2 R

WVðf ; t; xÞ ¼

where f denotes the complex conjugate of f. The Wigner–Ville distribution has received much attention for instantaneous frequency estimation due to its excellent concentration for mono-signals and many other desirable mathematical properties, see [17]. However, it is well known that it presents high amplitude sign-varying cross-terms for multi-component signals which makes its interpretation difficult. For attenuating these interference terms several approaches have been followed, mainly based on the smoothing of WV by convolution with a suitable regularizing kernel. Special mention is due to the spectrogram, which may be defined either as the energy density function of the short time Fourier transform

G u ðf ; t; xÞ ¼

Z

f ðsÞuðs  tÞeixs ds

ð1Þ

R

for some real, symmetric and normalized window u 2 L2 ðRÞ, i.e.

 u ðf ; t; xÞ; Su ðf ; t; xÞ ¼ G u ðf ; t; xÞG

ð2Þ

or as the convolution product of the Wigner–Ville distributions of the signal and the window

Su ðf ; t; xÞ ¼

Z R2

~ ÞWVðf ; t  ~t; x  x ~ Þd~tdx ~: WVðu; ~t; x

ð3Þ

In practice, for finite time duration and band-limited one-dimensional signals the time–frequency domain is bounded. Let us denote this domain by X, and write a generic time–frequency point ðt; xÞ 2 X by x. Dropping the dependence on the signal f and the window u in the notation for the spectrogram, the conservation of energy implies kSkL2 ðXÞ ¼ kf kL2 ð0;TÞ , and the additional assumption f ; u 2 L1 ð0; TÞ implies S 2 L1 ðXÞ as may be easily deduced from (1) and (2). Therefore, after a suitable normalization, we may consider the spectrogram of a signal as an image and apply to it well-known techniques of image

G. Galiano, J. Velasco / Applied Mathematics and Computation 244 (2014) 859–869

861

denoising and enhancing for the estimation of the underlaying signal. In addition, the definition of S as a convolution product involving the window u ensures further regularity when u is smooth. For instance, for Gaussian windows, S 2 C 1 ðXÞ. First attempts to denoising and enhancing spectrogram images were made via local differential filters [7,8,10] based on corresponding well established methods, see [1,5,19], in which the filtering process of a intensity image S : X ! ½0; 1 at x 2 X is based only on the intensities in a neighborhood of x. Although the resulting denoised spectrogram greatly improves the IF estimation of the signal, both computational time and low energy harmonic removing were drawbacks which motivated other approaches, see [9]. In this paper we show that nonlocal filters, specially the Neighborhood filter, are computationally more efficient to deal with these images and give similar results. The main idea in Neighborhood filters, see [4] for a detailed description, is to take into account in the filtering process the global gray scale values to define neighboring pixels. Thus, the denoised value at pixel x is defined as an average of values at pixels having a gray scale value close to SðxÞ. The resulting nonlocal algorithm is written as

NF h SðxÞ ¼ where CðxÞ ¼

R X

1 CðxÞ

Z

exp 

X

jSðxÞ  SðyÞj2 h

2

! SðyÞdy;

ð4Þ

  2 exp jSðxÞ  SðyÞj2 h dy is a normalization factor.

Since one application of the Neighborhood filter is usually not enough for image denoising, an iterative scheme based on formula (4) is carried out,

1 Siþ1 ðxÞ ¼ C m ðxÞ with C m ðxÞ ¼

R X

Z X

exp 

jSm ðxÞ  Sm ðyÞj2 h

2

! Si ðyÞdy;

ð5Þ

  2 dy, where m ¼ 0, for the linear fixed-kernel iterative scheme, or m ¼ i for the exp jSm ðxÞ  Sm ðyÞj2 h

nonlinear varying-kernel iterative scheme. In Singer et al. [20], a probabilistic approach is used for the understanding of the fixed-kernel iterative Neighborhood filter (m ¼ 0 in (5)). Its rate of convergence towards the constant steady state and their deblurring properties were analyzed interpreting the successive iterations as a random walk. The iterated composition of the corresponding transition probability matrix was then shown to produce the noise reduction at a certain asymptotic rate, acting specially on the high frequency components. More recently [11], both the fixed and the varying-kernels schemes have been reformulated in terms of the so-called decreasing rearrangement of the image S, which is closely related to the image histogram. Apart from a large complexity reduction of the algorithm explaining its low computational effort, some qualitative properties have been proved, among which:  The asymptotic behavior of the NF as a shock filter of the type introduced by Álvarez et al. [2], inducing pixel mass accumulation in local maxima of the histogram. In this way, the instantaneous frequency lines which are initially broaden by the spectrogram, become thinner and more concentrated in the final image. In addition, since the noise intensity distribution is expected to be highly varying, it will be attenuated in the final image.  The contrast change character of the NF, i.e. the existence of a continuous and increasing function g : R ! R such that NF h ðSðxÞÞ ¼ gðSðxÞÞ, leaving invariant the structure of the level lines of S. The latter property is the fundament for the fast algorithm implemented in this article, since the algorithm complexity only depends on the number of level lines. Indeed, if N denotes the image size in pixels, the complexity of local PDE based methods, see for instance (18), or of other nonlocal filters, such as the Bilateral (11) or the Nonlocal Means [4] filter, are of the order C  N, where C depends on different parameters used in those models: number of time iterations in PDE models, or kernel window sizes in nonlocal methods. A typical value of C may be as high as of the order 105 . However, the complexity of the rearranged version of the NF depends on the number of the initial image intensity levels, which we assume quantized in Q levels, and on some small constant related to the window size, h, resulting in a complexity of the order c  Q , being a typical value of c of the order 102 . Thus, the complexity of the NF is independent of the image size, once the level lines of the initial image have been identified. However, denoising quality of the NF is, in general, poorer. The staircasing effect, always present in algorithms reducing the Total Variation of the initial image, is especially strong for the NF. In fact, after several iterations of the varying-kernel algorithm, and depending on the window size h, the NF output image concentrates most of its pixels mass on few level sets, producing a segmentation-like effect on the initial image. Although this may represent a severe inconvenient for natural images, in which the image smoothness is primary for visual perception, it is not that important for spectrogram analysis, since we are mainly interested in the extraction of time–frequency lines representing the instantaneous frequencies of the underlying 1D signal.

862

G. Galiano, J. Velasco / Applied Mathematics and Computation 244 (2014) 859–869

2.1. Relationship between nonlocal image filtering and time–frequency signal processing We finish this section by describing an important property which is related to the special character of the spectrogram interpretation as an image, and which is in close relationship with the so-called reassigned spectrogram. As it is well known, an important drawback of the spectrogram with respect to the WV distribution is the broadening of the IF lines produced as a direct consequence of the convolution with a smooth window, see (3). To override this inconvenient, it is suggested in [14] that instead of assigning the averaged energy to the geometric center of the smoothing kernel, ^ ¼ ð^t; x ^ Þ, x ¼ ðt; xÞ, as it is done for the spectrogram, one assigns it to the center of gravity of these energy contributions, x which is certainly more representative of the local energy distribution of the signal, see Fig. 1. As deduced in [3], the gravity center may be computed by the following formulas (component-wise)

^ðxÞ ¼ x

1 SðxÞ

Z

yWVðf ; yÞWVðu; y  xÞ

R2

dy : 2p

^ ðxÞ, leading to define the The spectrogram value is then moved from the point x where it has been computed to this centroid x reassigned spectrogram as

1 2p

RSðxÞ ¼

Z

^ ðyÞÞdy: SðyÞdðx  x

R2

^ðyÞ  y, may be expressed through a For our discussion, it is of special interest the fact that the reallocation vector, rðyÞ ¼ x potential related to the spectrogram [6],

rðyÞ ¼

1 r log SðyÞ; 2

when u is a Gaussian window of unit variance, thus obtaining

1 2p

RSðxÞ ¼

  1 SðyÞd x  y  r log SðyÞ dy: 2 R2

Z

ð6Þ

Therefore, the spectrogram energy content at y will contribute to the reassigned spectrogram energy content at x if

xy ¼ Let

1 rSðyÞ: 2SðyÞ

ð7Þ

s 2 R be a small parameter and consider the set of y 2 R2 such that xy ¼

r 1 rSðyÞ; jrSðyÞj 2SðyÞ

with r ¼ jrSðyÞj þ s. Observe that for hand

  x  y  

s ¼ 0; y satisfies (7) and therefore contributes to RSðxÞ. Also from (7) we get, on one

2  2  1 r 1 rSðyÞ ¼ jrSðyÞj2 1 2SðyÞ jrSðyÞj 4S2 ðyÞ

ð8Þ

and, on the other hand, using Taylor’s expansion of order one,

jSðxÞ  SðyÞj  jðx  yÞ  rSðyÞj ¼

r jrSðyÞj; 2SðyÞ

Spectrogram

Reassigned Spectrogram

1300

1300

1200

1200

1200

1100 1000 900 800

Frequency

1300

Frequency

Frequency

Wigner−Ville Distribution

ð9Þ

1100 1000 900 800

0.4

0.5

0.6

Time

0.7

0.8

1100 1000 900 800

0.4

0.5

0.6

Time

0.7

0.8

0.4

0.5

0.6

0.7

0.8

Time

Fig. 1. Detail of the Wigner–Ville distribution, spectrogram and reassigned spectrogram corresponding to the clean signal of Experiment 1.

G. Galiano, J. Velasco / Applied Mathematics and Computation 244 (2014) 859–869

863

if jx  yj is small. Using the definition of r and gathering (8) and (9), we obtain

  x  y  

2  1 s2 1 rSðyÞ  2 jSðxÞ  SðyÞj2 : 2SðyÞ r jrSðyÞj2

Let us consider a Gaussian approximation to the Dirac delta with variance s. Then the integrand in (6) may be approximated as

 d xy

!  1 jSðxÞ  SðyÞj2 rSðyÞ  exp  : 2SðyÞ r 2 jrSðyÞj2

Therefore, setting h ¼ rjrSðyÞj  jrSðyÞj2 and taking into account that the computations were local, i.e. for jx  yj small, we get that the reassigned spectrogram of S, given by (6), may be approximated by

RSðxÞ  CðxÞYF h;q SðxÞ:

ð10Þ

Here, YF h;q stands for the Yaroslavsky filter [26], an interesting variation of the Neighborhood filter in which the domain of integration in (4) is restricted to a spatial neighborhood of the base pixel x, i.e. to a ball Bq ðxÞ ¼ fy 2 X : jx  yj < qg. The term CðxÞ is the usual normalization factor. Similarly, the reassigned spectrogram behavior can be linked to the SUSAN [22] and the Bilateral [24] filters, which incorporate the domain restriction as a weighted distance to the reference pixel

SNF h;q SðxÞ ¼

1 CðxÞ

Z

exp 

X

jy  xj2

q2



jSðxÞ  SðyÞj2 h

!

2

SðyÞdy:

ð11Þ

Observe that we recover the Neighborhood filter from the Yaroslavsky, the SUSAN or the Bilateral filters by performing the limit q ! 1, i.e. allowing jy  xj to be arbitrarily large. Let us finally remark that the approximation (10) just gives an idea of the behavior and relationship between the nonlocal filters and the reassigned spectrogram, in the sense that both tend to concentrate the spectrogram energy on the lines of maxima (ridges) of the spectrogram. However, the assumption h  jrSðyÞj2 should be taken with care, since h is taken constant in the nonlocal filtering algorithms. 3. Discretization and numerical experiments In this section we present numerical demonstrations on spectrograms of synthetic and biomedical signals of the nonlocal Neighborhood filter (5), the Yaroslavsky filter (11) iterated as in (5), and a Total Variation based PDE local filter, see (15)–(17) below. In order to compare the results, we use an objective measure, the Mean Square Error (MSE) between the denoised and the clean signals, and a subjective measure based on the visual inspection of the spectrogram itself and on some related magnitudes: the IF lines and the time-energy content. 3.1. Discretization Computation of the spectrogram is a standard operation performed by applying the discrete fast Fourier transform (dfft) to the convolution of the signal with the window, which we choose as a zero mean unit variance Gaussian window. Once the spectrogram is produced, it is normalized in some range ½0; Q , with Q > 0, obtaining in this way a discrete function S : X ! ½0; Q . The choice Q ¼ 255 corresponds to the usual digital image range. Let us consider the regular mesh of X ¼ ð0; TÞ  ð0; FÞ, as produced by the spectrogram computation, given by fðtm ; xn Þg, with t m ¼ mT=ðM þ 1Þ, for m ¼ 0; 1; . . . ; M þ 1 and xn ¼ nF=ðN þ 1Þ, for n ¼ 0; 1; . . . ; N þ 1, and consider the corresponding barycenters xm;n for m ¼ 0; 1; . . . ; M and n ¼ 0; 1; . . . ; N. To each xm;n we associate a surface element dA ¼ TF=MN. For Q 2 N, we consider the Q quantized spectrogram S0 ðxÞ ¼ k if SðxÞ 2 ½k  12 ; k þ 12Þ \ ½0; Q , for k ¼ 0; 1; . . . ; Q , inducing the partition of X given by

X¼

Q [

Bk ;

with Bk ¼ fx 2 X : S0 ðxÞ ¼ kg:

ð12Þ

k¼0

According to Eq. (5) for the iteration of the Neighborhood filter, and using (12), we define for j ¼ 0; 1; . . . ; Q , and for the iteration index i ¼ 0; 1; . . .

  SNF iþ1 ð x 2 Bj Þ ¼

1 C 0 ðxÞ

Z X

exp 

jS0 ðxÞ  S0 ðyÞj2 2

h

!

PQ SNF i ðyÞdy ¼

with SNF 0 ¼ S0 . We consider the zero-order approximations

k¼0 e

jk 2 Þ h

ð

PQ

k¼0 e

R Bk

SNF i ðyÞdy

jk 2 Þ h

ð

jBk j

;

ð13Þ

864

G. Galiano, J. Velasco / Applied Mathematics and Computation 244 (2014) 859–869

Z Bk

X

SNF i ðyÞdy 

SNF i ðy m;n ÞdA;

and jBk j  CðkÞdA;

ð14Þ

ym;n 2Bk

where CðkÞ denotes the number of nodes contained in Bk . Relabeling SNF iþ1 in terms of its level sets, the approximation to formula (13) (for which we keep the same notation) reads

PQ SNF iþ1 ðjÞ

k¼0 e

¼

jk 2 Þ h

ð

P

NF ym;n 2Bk Si ðy m;n Þ

PQ

ðjk Þ h k¼0 e

2

;

CðkÞ

or, in a more compact form

SNF iþ1 ðjÞ ¼

haðj; Þ; S i i ; haðj; Þ; Ci

with jk 2

X

S i ðkÞ ¼

aðj; kÞ ¼ eð h Þ ;

SNF i ðy m;n Þ

ym;n 2Bk

and h; i denoting the scalar product in RQ þ1 . In a similar way, we may find the following expression for the iteration of the Yaroslavsky type filter given by (11). However, for this filter we can not operate on terms of level sets, and a particular expression for each xm;n 2 Bj must be given, explaining the large time execution differences between Neighborhood and Yaroslavsky filtering:

PQ SYiþ1 ðxm;n Þ

¼

jk 2 Þ h

R

yxm;n 2

eð q Þ SYi ðyÞdy : PQ ðjkÞ2 R ðyxm;n Þ2 h e q dy k¼0 e B

ð k¼0 e

Bk

k

Approximating the integral terms as in (14), we obtain

PQ SYiþ1 ðxm;n Þ

¼

ð k¼0 e

jk 2 Þ h

PQ

k¼0 e

P

ð yp;q 2Bk e

jk 2 Þ h

ð

P

yp;q xm;n 2

Þ

q

yp;q 2Bk e

SYi ðyp;q Þ

yp;q xm;n 2

ð

q

:

Þ

Finally, for comparison purposes, we used the following denoising algorithm based on a convection–diffusion nonlinear PDE related to that introduced by Rudin et al. [19]. The problem is written as: find S : 0; 1e  X ! ½0; Q  such that

   @S e rS 1 ¼ 0 in 0;  X; þ r logðS0 Þ  rS  div @s 2 jrSj e 

rS  n ¼ 0 on 0; Sð0; Þ ¼ S0

1

e

ð15Þ

 @ X;

ð16Þ

in X;

ð17Þ

where S0 ¼ Gr  S0 , with Gr a Gaussian regularizing kernel of variance r. Here, e > 0 controls the strength of transport towards the maxima of S0 , in relation to the mean curvature diffusion term. The addition of the transport term in (15), as well as the bounded time domain ½0; 1e , which are not usual in image denoising algorithms, are motivated by the notion of differential reassignment introduced by Auger et al. [5], which is proper of one-dimensional signals, see [8] for further details. The discretization of (15)–(17) follows the standard finite element methodology for the total variation denoising problem, see [10]. We consider a time semi-implicit Euler scheme and a P1 continuous finite element approximation in space. In particular, for the approximation, Skþ1 , at time t ¼ t kþ1 , we use the diffusive coefficients obtained in the previous time step tk , that is,

0

1

Skþ1  Sk  1 B C þ r logðS0 Þ  rSk  div@qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi rSkþ1 A ¼ 0 dt 2 2 2 ~ jrSk j þ 

ð18Þ

for some ~  > 0. The convective term is computed by an upwind scheme, evaluating for each node the derivatives in the element corresponding to the upwind direction. Therefore, for each time step, we solve an elliptic problem with a diffusive term coefficient bounded from above. As mentioned in the introduction to this section, we use two measures to compare the denoising results:

865

G. Galiano, J. Velasco / Applied Mathematics and Computation 244 (2014) 859–869

(1) The Mean Square Error between the clean spectrogram, Sc ðxÞ, and the filtered spectrogram after I iterations, SI ðxÞ

MSE ¼

kSc  SI k ; kSc k

where k  k denotes the L2 norm of a matrix. (2) The visual inspection of the IF lines. We use a simple algorithm to produce candidates to IF lines of the corresponding spectrograms. We consider the truncation v ðt; xÞ ¼ Sðt; xÞ if Sðt; xÞ P b, and v ðt; xÞ ¼ 0 elsewhere, with b ¼ MeanX ðSÞ in the experiments. For each t 2 ½0; T we find the N connected components of the set fx 2 ð0; FÞ : v ðt; xÞ > 0g, say C n ðtÞ, for n ¼ 1; . . . ; NðtÞ, and define the function

R

xv ðt; xÞdx ; v ðt; xÞdx C n ðtÞ

C ðtÞ IFðt; nÞ ¼ Rn

ð19Þ

which is the frequency gravity center of the component C n ðtÞ. In this way, we shrink each connected component to one point to which we assign the average image intensity through the function INTðt; nÞ ¼ MeanC n ðtÞ ðv ðt; ÞÞ. Finally, we plot function IF only for components with averaged intensity, INT, greater than a certain threshold, i 2 ½0; Q . This final image does not seem to be very sensible under small perturbations of the parameters b and i. In all the experiments we have chosen the usual image quantization into 256 levels (Q ¼ 255). For the Neighborhood and Yaroslavsky filters, the variance of the Gaussian neighborhood for the level lines is taken as h ¼ 10. For the latter, we fix the variance of the local neighborhood as q ¼ 10 so, according to [4], the nonlinear diffusion regime applies. For the Neighborhood algorithm, the iteration stopping criterium is established in terms of the relative difference between two consecutive iterations, i.e. we stop the algorithm if

kSiþ1  Si kL2 < tol; kSi kL2 with tol empirically chosen to perform few iterations. For the PDE based filter algorithm we stop when the artificial time

send ¼ 1=e is reached. The parameter e, controlling the ratio transport-diffusion is empirically chosen. For the Yaroslavsky filter we use the same number of iterations as for the PDE algorithm. In Table 1, the execution time expressed in seconds refers to the tic-toc Matlab command for a script run under a Core i7 (quad) processor, built on a standard laptop. Experiment 1. We use a one second 4 KHz synthetic signal composed by the addition of two signals. The first is the addition of pure tones and chirps:

  x1 ðtÞ ¼ c1 sinð2p500tÞ þ sinð2p700tÞ þ sinð2p1000t 2 Þ þ sinð2p600t3 Þ ; while the second, x2 , is a uniformly distributed real random variable. We normalize them to have kxi kL2 ¼ 1 (so the constant c1 ) and define the test signal as x ¼ x1 þ x2 , i.e., with SNR ¼ 0. We then compute the spectrogram, S, for a Gaussian window of unit variance. The size of the corresponding image is 2048  512 pixels. Other parameters are chosen as follows: tol ¼ 0:04; e ¼ 0:02 and ds ¼ 2:5, implying send ¼ 50. In Fig. 2 we show the graphic results. The left panel contains the spectrogram of the noisy signal and the corresponding images filtered with the three described algorithms. In comparison with the original spectrogram, all filtered images are very closed from each other. We observe that the image obtained by Neighborhood filtering keeps more noise than the others, but also experiments less diffusive effects. This is specially noticeable in the broken IF segments appearing at high frequencies, where the Yaroslavsky and the TV filtering diffuse locally while the Neighborhood filter just remove the information of this area. The similarity of these denoised images is also checked by comparing their corresponding MSE, see Table 1. However, the most impressive result is that concerning the execution times which gives a ratio larger than 200 in favor of the Neighborhood algorithm. Experiment 2. We use a recording done in wilderness [16], from where we extracted a signal of 2:25 s which is affected by a strong background noise. Since field data recorders are set to 44.1 kHz meanwhile wolves signals are rarely out of the range 200–3000 Hz, we start by filtering and down-sampling the signal to speed up computations. The resulting spectrogram image size is 172  690 pixels. Other parameters are chosen as follows: tol ¼ 0:01; e ¼ 0:2 and ds ¼ 0:25, implying send ¼ 5. Table 1 Mean Square Error (MSE) and execution time (s), of the different algorithms and experiments. A smaller MSE indicates that the estimate is closer to the original image. The numbers have to be compared on each row.

Experiment 1 Experiment 2 Experiment 3

MSE Ex. time Ex. time Ex. time

Neighborhood

Yaroslavsky

Nonlinear diff.

0.28 2.33 0.24 1.03

0.16 874 95 125

0.17 255 22 32

866

G. Galiano, J. Velasco / Applied Mathematics and Computation 244 (2014) 859–869

Instantaneous frequency 2000

1400

1400

Frequency

Frequency

Spectrogram 2000

800

800

200

200 0

0

0.5

Time

1400

1400

Frequency

Frequency

2000

800

200

200 0

0.5

0

Time 2000

1400

1400

Frequency

Frequency

1

Instantaneous frequency

2000

800

0.5

Time

Yaroslavsky filter

800

200

200 0

0.5

0

Time

1

Instantaneous frequency 2000

1400

1400

Frequency

2000

800

0.5

Time

Total Variation

Frequency

1

Instantaneous frequency

Neighborhood filter 2000

800

0.5

Time

800

200

200 0

0.5

Time

0

0.5

1

Time

Fig. 2. Experiment 1. Left: spectrogram of the noisy signal and filtered images obtained with the different algorithms. Right: corresponding IF lines according to formula (19).

In Fig. 3 we show the graphic results. Again, the left panel contains the spectrogram of the signal and the results of filtering. In this case, the noise is not Gaussian and the filtering algorithms are less adapted to the problem. The Neighborhood filter keeps most of the relevant information (the fairly straight lines) but leaves part of the low frequency

867

G. Galiano, J. Velasco / Applied Mathematics and Computation 244 (2014) 859–869

Instantaneous frequency

2600

2600

2000

2000

Frequency

Frequency

Spectrogram

1400

800

200

1400

800

200 0

0.5

1

1.5

2

0

0.5

Time

2600

2600

2000

2000

Frequency

Frequency

2

Instantaneous frequency

1400

800

200

1400

800

200 0

0.5

1

1.5

2

0

0.5

Time

1

1.5

2

Time

Yaroslavsky filter

Instantaneous frequency

2600

2600

2000

2000

Frequency

Frequency

1.5

Time

Neighborhood filter

1400

800

200

1400

800

200 0

0.5

1

1.5

2

0

0.5

Time

1

1.5

2

Time

Total Variation

Instantaneous frequency

2600

2600

2000

2000

Frequency

Frequency

1

1400

800

200

1400

800

200 0

0.5

1

Time

1.5

2

0

0.5

1

1.5

2

Time

Fig. 3. Experiment 2. Left: spectrogram of the noisy signal and filtered images obtained with the different algorithms. Right: corresponding IF lines according to formula (19).

868

G. Galiano, J. Velasco / Applied Mathematics and Computation 244 (2014) 859–869

Clean signal

Noisy signal

1.2

0.9

0.6

1.5

Frequency

1.5

Frequency

Frequency

1.5

1.2

0.9

0.6 250

300

350

400

450

−3

4

300

350

400

450

250

x 10

Noisy signal

1

350

400

450

400

450

x 10

400

450

3

2

1

Time

350

Extraction

−3

4

Energy

2

300

300

Time

3

Energy

Energy

3

250

0.9

Time

Clean signal

−3

x 10

1.2

0.6 250

Time 4

Extraction

2

1

250

300

350

400

450

250

300

Time

350

Time

Fig. 4. Experiment 3. Recovery of an arrhythmia episode detection from a noisy signal.

noise unfiltered. The Yaroslavsky filter removes the noise, but also the IF lines with less energy content. In a middle way, the TV algorithm removes most of the noise and keeps a large part of the IF lines, being for this example the most adequate algorithm. However, execution times show again large differences in favor of the Neighborhood filter. Experiment 3. We used ECG signals obtained from the MIT-BIH noise stress test database [18,12], from which we extracted a segment of 240 s (minutes 4 to 8), containing an arrhythmia episode in its central part. The recording was made using physically active volunteers and standard ECG recorders. Noise was added beginning after the first 5 min of each record, during two-minute segments alternating with two-minute clean segments. The signal-to-noise ratio was 12 dB. The added noise was of the type electrode motion artifact, which is generally considered the most troublesome, since it can mimic the appearance of ectopic beats and cannot be removed easily by simple filters, as can noise of other types. For the clean image, see Fig. 4 (first row), the spectrogram (size 236  676) shows the usual energy spread of arrhythmia episodes around t ¼ 360. Integration of the energy in the frequency domain shows a pick at the same time instant, which allows to detect the episode, see Fig. 4 (second row). In the central column of the same figure we repeat these calculations for the noisy signal. As it may be observed, the energy pick has been hidden due to the time coincidence with the high energy noise. In the last column, we show the result of a filtering-extraction procedure. Let S0 be the spectrogram of the noisy signal and Sn the image resulting from the application of the Neighborhood filter. Then we subtract Sn > a from S0 , for some small a 2 ½0; Q  (a ¼ 1 in the experiment). The corresponding image and its integration in the frequency domain are shown in the last column of Fig. 4. As we see, after the subtraction of the most energetic components of the noisy signal, the noise and the normal beat frequency line, the arrhythmia episode may be again detected. Although not shown, the results with the other algorithms are similar, having again in this experiment much larger execution times than that of the Neighborhood filter algorithm.

4. Conclusions One dimensional signal denoising and enhancing is a topic of interest in may fields of applications, among which the Biosignals is of particular importance. In previous works we addressed the problem by filtering the image produced by some time–frequency representation of the signal, usually the spectrogram, through PDE based algorithms. Although the filtering results on the spectrogram image

G. Galiano, J. Velasco / Applied Mathematics and Computation 244 (2014) 859–869

869

greatly improved those obtained from the direct filtering of the one-dimensional signal, the execution times were prohibitive for long duration signals. In this article we replaced PDE based methods by nonlocal filtering image techniques, observing a similar filtering behavior for both methodologies. However, nonlocal filtering by the Neighborhood filter offers execution times several orders of magnitude lower than the other techniques, while the results are similar when applied to the specific images produced by time–frequency representations. The reason of this suited behavior of nonlocal filtering on spectrogram images has been motivated by the link between these filters and the reassigned spectrogram, an improved version of the spectrogram which concentrates the energy in the instantaneous frequency lines. We believe that this approach can be further exploited by introducing weighting terms in the nonlocal filters which bring even closer the results of both operations. This topic will be the subject of future research. References [1] L. Álvarez, P.L. Lions, J.M. Morel, Image selective smoothing and edge detection by nonlinear diffusion. II, SIAM J. Numer. Anal. 29 (3) (1992) 845–866. [2] L. Álvarez, L. Mazorra, Signal and image restoration using shock filters and anisotropic diffusion, SIAM J. Numer. Anal. 31 (2) (1994) 590–605. [3] F. Auger, P. Flandrin, Improving the readability of time–frequency and time-scale representations by the method of reassignment, IEEE Trans. Signal Process. 43 (5) (1995) 1068–1089. [4] A. Buades, B. Coll, J.M. Morel, A review of image denoising algorithms, with a new one, Multiscale Model. Simul. 4 (2005) 490–530. [5] E. Chassandre-Mottin, I. Daubechies, F. Auger, P. Flandrin, Differential reassignment, IEEE Signal Process. Lett. 4 (10) (1997) 293–294. [6] E. Chassandre-Mottin, I. Daubechies, F. Auger, P. Flandrin, Differential reassignment, IEEE Signal Process. Lett. 4 (10) (1997) 293–294. [7] B. Dugnol, C. Fernández, G. Galiano, Wolves counting by spectrogram image processing, Appl. Math. Comput. 186 (1) (2007) 820–830. [8] B. Dugnol, C. Fernández, G. Galiano, J. Velasco, Implementation of a diffusive differential reassignment method for signal enhancement: an application to wolf population counting, Appl. Math. Comput. 193 (2007) 374–384. [9] B. Dugnol, C. Fernández, G. Galiano, J. Velasco, On a chirplet transform based method applied to separating and counting wolf howls, Signal Process. 88 (2008) 1817–1826. [10] B. Dugnol, C. Fernández, G. Galiano, J. Velasco, Evolution nonlinear diffusion–convection PDE models for spectrogram enhancement, Numer. Anal. Appl. Math., AIP Conf. Proc. 1048 (2008) 166–169. [11] G. Galiano, J. Velasco, Neighborhood filters and the decreasing rearrangement, J. Math. Imaging Vision, in press. http://dx.doi.org/10.1007/s10851-0140522-3. [12] A.L. Goldberger, L.A.N. Amaral, L. Glass, J.M. Hausdorff, P.Ch. Ivanov, R.G. Mark, J.E. Mietus, G.B. Moody, C.K. Peng, H.E. Stanley, Physiobank, physiotoolkit, and physionet: components of a new research resource for complex physiologic signals, Circulation 101 (2000) e215–e220. (June 13). [13] S.L. Hopp, M.J. Owren, C.S. Evans (Eds.), Animal Acoustic Communication. Sound Analysis and Research Methods, Springer, Berlin, 1998. [14] K. Kodera, R. Gendrin, C. de Villedary, Analysis of time-varying signals with small BT values, IEEE Trans. Acoust. Speech Signal Process., ASSP-26 1 (1978) 64–76. [15] J.S. Lee, Digital image smoothing and the sigma filter, Comput. Vision Graph. 24 (1983) 255–269. [16] L. LLaneza, V. Palacios, Asesores en Recursos Naturales, S.L., Field recordings obtained in wilderness in Asturias (Spain) in the 2003 campaign. [17] S. Mallat, A Wavelet Tour of Signal Processing, Academic Press, London, 1998. [18] G.B. Moody, W.E. Muldrow, R.G. Mark, A noise stress test for arrhythmia detectors, Comput. Cardiol. 11 (1984) 381–384. [19] L. Rudin, S. Osher, E. Fatemi, Nonlinear total variation based noise removal algorithms, Physica D 60 (1992) 259–268. [20] A. Singer, Y. Shkolnisky, B. Nadler, Diffusion interpretation of nonlocal neighborhood filters for signal denoising, SIAM J. Imaging Sci. 2 (2009) 118–139. [21] A. Skonhoft, The costs and benefits of animal predation: an analysis of Scandinavian wolf re-colonization, Ecol. Econ. 58 (4) (2006) 830–841. [22] S.M. Smith, J.M. Brady, SUSAN–a new approach to low level image processing, Int. J. Comput. Vision 23 (1997) 45–78. [23] L. Sörnmo, P. Laguna, Bioelectrical Signal Processing in Cardiac and Neurological Applications, Academic Press, 2005. [24] C. Tomasi, R. Manduchi, Bilateral filtering for gray and color images, in: Proceedings of the 6th International Conference on Computer Vision, 1998, pp. 839–846. [25] H. Xu, Z.-H. Tan, P. Dalsgaard, B. Lindberg, Robust speech recognition by nonlocal means denoising processing, IEEE Signal Process. Lett. 15 (2008) 701– 704. [26] L.P. Yaroslavsky, Digital Picture Processing. An Introduction, Springer-Verlag, Berlin, 1985.