Robust Time-Frequency Distributions Based on ... - Semantic Scholar
1
Robust Time-Frequency Distributions Based on the Robust Short Time Fourier Transform Igor Djurovi´c1 , LJubiˇsa Stankovi´c1 and Braham Barkat2 1
Elektrotehniˇcki fakultet Podgorica, University of Montenegro, Podgorica, Montenegro. Tel./Fax. +381 81 244 921. E-mail: [email protected], [email protected].
2
Nanyang Technology University, School of EEE, Block S2, Nanyang Avenue, Singapore 639798. E-mail: [email protected]. Abstract Design of time-frequency distributions (TFDs) that are robust to the impulse noise influence is considered. The robust TFDs based on the robust short-time Fourier transform (STFT) are proposed. An efficient procedure to evaluate the robust STFT is given. Robust TFDs based on the robust STFT have better energy concentration around the signal instantaneous frequency (IF) than the robust STFT itself. Also, these TFDs are more resistant to higher impulse noise than the robust TFDs obtained using the local autocorrelation function (LAF) based minimization problem.
I. Introduction An overview of time-frequency distributions (TFDs) that have attracted significant attention in the signal processing community can be found in [1], [2]. Influence of Gaussian noise on TFDs has been extensively studied [3]. The instantaneous frequency (IF), a concept closely related to time-frequency (T-F) analysis, has attracted wide interest. The IF estimation of signals embedded in Gaussian noise has been addressed in several research works, among them [4], [5]. However, only few research works have dealt with other types of noise [6]. In numerous practical situations, environment cannot be modeled with Gaussian noise. Namely, due to natural or man-made phenomenon, the resulting noise could have rare but very strong values. This kind of noise is referred to as impulse noise. More details on the impulse noise and impulse noise modeling can be
2
found in [7]. TFDs are sensitive to this kind of noise. An intensive research effort has been put into design of TFDs that are robust to such a noise [8]-[10]. These TFDs are based on the robust statistics concept, introduced by Huber [11]. This concept has been successfully used in various signal processing applications. Katkovnik introduced this concept into spectral analysis by developing the robust M-Fourier transform (MFT) [8]. In the sequel, transforms that are defined on the basis of the robust statistics concept will be referred to as ‘robust’. The robust M-FT has been defined in an implicit form. An iterative procedure is required for its evaluation. In order to avoid iterative procedures, the marginal-median filter form of the robust FT was proposed in [9]. Both forms of the robust FTs yield accurate spectra estimates for stationary signals embedded in impulse noises. They are almost as accurate as the standard FT for Gaussian noise environment. Generalization to the T-F analysis of non-stationary signals corrupted by impulse noise has been proposed in [9], [10]. The robust Wigner distribution (WD) can be evaluated as a median of the modulated local autocorrelation function (LAF) samples. Note that all minimization problems based on the Hermitian LAF yield real-valued robust TFDs [9]. The L-filter form of robust TFDs is proposed for handling signals corrupted by a mixture of Gaussian and impulse noise in [10]. Higher-order TFDs yield a better energy concentration around their IF than the STFT in the case of noiseless signals. In this paper, we will show that the same fact holds for the robust TFDs. In addition, IF estimates based on peaks of the higher-order robust TFDs have smaller bias as compared to those based on the robust STFT peaks. However, higher-order TFDs are more sensitive to the impulsive noise than the robust STFT. To mitigate this, we propose a new form of the robust TFDs based on the robust STFT. In this way, a significant improvement, comparing to the results obtained by using the robust TFDs based on the LAF, is achieved. In addition, we propose a recursive technique to evaluate the robust STFT. As a consequence, the robust TFDs based on the robust STFT can be efficiently calculated in real-time. After reviewing the existing robust TFDs in Section II, we introduce robust TFDs based on the robust STFT, as well as the recursive procedure to evaluate them, in Section III. A comparison study of various robust TFDs, using IF estimation, is presented in Section IV.
3
II. A Review of the Robust FT and TFDs A. Robust FT Consider discretized version of the useful signal f(n) that is embedded in a white noise ν(n), x(n) = f(n) + ν(n). The standard FT of noisy signal x(n),
X(ω) =
1 N
XN−1 n=0
x(n)e−jωn = mean{x(n)e−jωn | n ∈ [0, N − 1]},
(1)
can be defined as a solution of the optimization problem1 :
X(ω) = arg min J(m; ω) m
(2)
where the penalty function can be written as
J(m; ω) =
XN−1 n=0
F (x(n)e−jωn − m),
(3)
for the loss-function F (e) = |e|2 . The standard DFT, X(ω), is the ML estimate of the useful signal spectra F (ω) for a Gaussian noise environment. It is sensitive to the impulse noise. For a known probability density function (pdf) of the noise the ML estimate of the FT can be determined. The ML estimates are sensitive to variations of the noise pdf. To avoid this difficulty, Huber introduced the M-estimation concept whose main idea is as follows. A class of noises is considered. The ‘worst’ noise is defined as the one whose pdf has the heaviest tail. The ML estimate for this noise is the M-estimate for the used class. The M-estimate produces reasonable accuracy for the entire class. The Laplacian noise was found as the worst one for many impulse noise classes considered in practice. The loss-function F (e) = |e|, that produces the ML estimate for the Laplace noise, results in robust behavior for almost all other impulse disturbances. The robust M-FT in implicit form, obtained as the solution of the 1
Minimizing the penalty function with respect to the complex-valued argument is done in the following manner. Penalty function
J (m) is considered as the 2D function of real-valued arguments J(mr , mi ; ω) where mr and mi are real and imaginary part of m. The searched solution is formed as X(ω) = R(ω) + jI(ω), where (R(ω), I(ω)) are solution of equation (R(ω), I(ω)) = arg J (mr , mi ; ω).
min
(mr ,mi )
4
minimization problem (2)-(3) for F (e) = |e|, is given by [8] XR (ω) =
PN−1
x(n)e−jωn n=0 |x(n)e−jωn −XR (ω)| PN−1 1 n=0 |x(n)e−jωn −XR (ω)|
.
(4)
To avoid computational problems associated with iterative procedures, required for solving (4), an alternative form of the robust M-FT, based on the loss function F (e) = | Re(e)| + | Im(e)|, was proposed in [9]. This alternative form, called the marginal median form of the M-FT, is given by [9] n
o
n
o
XM (ω) = median Re{x(n)e−jωn }| n ∈ [0, N − 1] + j median Im{x(n)e−jωn }| n ∈ [0, N − 1] .
(5)
Transforms (4) and (5) behave similarly. Their accuracy is high for impulse noise. For Gaussian noise they are slightly worse than the standard FT. The median-filter form does not suffer from convergence problem that may exist in (4). B. Robust TFDs Robust Short Time Fourier Transform: The robust STFT forms, which can be seen as an extension of the robust FT, are defined by using the optimization problem ST F T (n, ω) = arg min J(m; n, ω) = arg min m
m
XN−1 k=0
F (x(n + k)e−jωk − m).
(6)
For the loss functions F (e) = |e|2 , F (e) = |e| and F (e) = | Re(e)| + | Im(e)|, yield the standard STFT (ST F TS (n, ω)), the robust STFT in the implicit form (ST F TR (n, ω)), and the robust STFT in the marginalmedian form (ST F TM (n, ω)), respectively. These transforms are given by ST F TS (n, ω) = ST F TR (n, ω) =
XN−1
x(n + k)e−jωk k=0 PN−1 x(n+k)e−jωk k=0 |x(n+k)e−jωk −ST F TR (n,ω)| PN−1 1 k=0 |x(n+k)e−jωk −ST F TR (n,ω)| 1 N
(7) (8)
ST F TM (n, ω) = median{Re{x(n + k)e−jωk }|k ∈ [0, N − 1]} +jmedian{Im{x(n + k)e−jωk }|k ∈ [0, N − 1]}.
(9)
Robust TFDs Based on the LAF: The robust WD was introduced in [9] as a solution of the minimization problem W D(n, ω) = arg min J(m; n, ω) = arg min m
m
XN−1 k=0
F (rxx (n, k)e−2jωk − m)
(10)
5
where the LAF is given by rxx (n, k) = x(n + k)x∗ (n − k). For F (e) = |e|2 , we obtain the standard WD (W DS (n, ω)) and for F (e) = |e|, we obtain the robust WD (W DR (n, ω)), i.e., W DS (n, ω) =
1 N
XN−1 k=0
rxx (n, k)e−2jωk = mean{rxx (n, k)e−2jωk | k ∈ [0, N − 1]},
W DR (n, ω) = median{Re{rxx (n, k)e−2jωk }| k ∈ [0, N − 1]}.
(11) (12)
Here, the loss function F (e) = | Re(e)| + | Im(e)| yields the same solution as the one obtained for F (e) = |e|. The LAFs, rxx (n, k), of some common TFDs are given in Table I [12]-[15], where the acronyms CD, LWD and PWVD stand for the Cohen class of distributions, the L-Wigner distribution and the polynomial Wigner-Ville ∗ (n, −k), the distribution, respectively. It can be shown that for any Hermitian LAF, i.e., rxx (n, k) = rxx
optimization problem (10) results in a real solution for any even loss function. Using the appropriate LAF expression, we define the standard TFD (T F (n, ω)) as the solution of the above optimization problem when F (e) = |e|2 . Similarly, we define the robust TFD (T FR (n, ω)) as the solution obtained for F (e) = |e|. These two quantities are equal to n
T F (n, ω) = mean rxx (n, k)e−2jωk | k ∈ [0, N − 1] n
o
(13) o
T FR (n, ω) = median Re{rxx (n, k)e−2jωk }| k ∈ [0, N − 1] .
Representation
LAF rxx (n, k)
WD
x(n + k)x∗ (n − k)
CD1
cT (n, k) ∗n x(n + k)x∗ (n − k)
LWD
xL (n + k/L)xL∗ (n − k/L)
PWVD2 q = 4
x2 (n + 0.675k)x2∗ (n − 0.675k)x∗ (n + 0.85k)x(n − 0.85k)
(14)
TABLE I 1
cT (n, k) real and symmetric kernel function in the time-lag domain.
2
the PWVD form used in the
paper.
L-filter Forms: The L-filters are effective in filtering stationary signals embedded in a mixture of Gaussian
6
and impulsive noise. They are also successfully applied to non-stationary signals [10] and for design of the robust filters of high-pass data [16]. The L-form of the STFT is given as
ST F TL (n, ω) =
XN−1 l=0
al [r(n, ω; l) + j · i(n, ω; l)],
(15)
where r(n, ω; l) and i(n, ω; l) are elements of the sets R(n, ω) = {Re{x(n + k)e−jωk }| k ∈ [0, N − 1]} and I(n, ω) = {Im{x(n + k)e−jωk }| k ∈ [0, N − 1]}, i.e., r(n, ω; l) ∈ R(n, ω), and i(n, ω; l) ∈ I(n, ω), sorted into non-decreasing order as r(n, ω; l) ≤ r(n, ω; l + 1), i(n, ω; l) ≤ i(n, ω; l + 1). The L-filter forms of the TFDs given in Table I, can be written as
T FL (n, ω) =
XN−1 l=0
al r(n, ω; l),
(16)
where r(n, ω; l) are sorted elements of the set R(n, ω) = {Re{rxx (n, k)e−2jωk }|k ∈ [0, N − 1]}. The filter coefficients al , l = 0, 1, ..., N − 1, are usually determined to satisfy the constraints
PN−1 l=0
al = 1 and al =
aN−l−1 . We will use the α-trimmed form of the L-filter whereby the coefficients are given as: al = 1/[N − α(N − 2)], l ∈ [(N − 2)α, (−N + 2)α + N − 1], α ∈ [0, 0.5], and al = 0 elsewhere. Note that, for α = 0 and α = 0.5, the standard and the median forms follow, respectively. III. Robust TFDs Based on the Robust STFT Realization of standard TFDs, based on the standard STFT, is a well known and described in numerous research papers [13]-[17]. This realization technique is used for many reasons, such as: computational efficiency, to avoid aliasing effects in TFDs, or to produce highly concentrated TFDs with reduced interferences. Since the STFT can be calculated recursively, it means that different TFDs can be realized in real-time. In addition to the just stated benefits, our main objective in this section is to introduce robust TFDs based on the robust STFT, that have stronger resistance to the impulse noise than the previously reviewed robust TFDs. To this aim, a recursive realization of the robust STFT is proposed in this section. As we have previously seen, the robust STFT can be realized by using expressions (8), (9) or (15). In the sequel, we can use any of these forms to realize the robust TFDs. These robust TFDs will be more resilient to impulsive noise than those realized by using the minimization problems discussed earlier. The robust STFT used to realize the robust TFDs will
7
be denoted by ST F Tρ (t, ω), where ρ can be any of the letters M, R or L, used to refer to the robust STFT forms in (8), (9) and (15). The notation of the various robust TFDs, obtained by using the robust STFT, will be preceded by r. A. Proposed Transforms In the light of what we discussed above, and using the same technique to realize standard TFDs based on the standard STFT, we itemize below the different forms of the robust TFDs based on the robust STFT. •
The windowed robust STFT can be calculated by rST F T (n, ω) = ST F Tρ (n, ω) ∗ω W (ω), where W (ω) is the FT of the desired window function and ∗ω represents the convolution in frequency. For instance, for a Hanning window (discrete domain), this convolution reduces to rST F T (n, k∆ω) =
P1
i=−1 ci ST F Tρ (n, (k+
i)∆ω), where c0 = 1/2, c±1 = 1/4, and ∆ω is the frequency step. Similar expressions can be derived for other window types. •
The robust forms of the WD based on the robust STFT are determined as [12] 1 rW D(n, ω) = π
Z
θ
ST F Tρ (n, ω + θ)ST F Tρ∗ (n, ω − θ)dθ.
(17)
The pseudo form of the rWD can be determined in the same manner as the windowed form of ST F Tρ (n, ω). •
The S-method (SM) [12] combines favorable properties of the STFT (non-aliasing for the Nyquist rate sampled signal and cross-terms reduction) and the WD (high concentration). We define its robust counterpart as rSM (n, k) =
1 π
Z
θ
P (θ) ST F Tρ (n, ω + θ) ST F Tρ∗ (n, ω − θ)dθ,
(18)
where P (θ) is the window function in the frequency domain. Wider window means better T-F concentration, while narrower window means better cross-terms interference suppression. Full details on the window length determination can be found in [12]. •
Realization of the quadratic Cohen’s class of TFDs by using the standard STFT was introduced in [17]. We define their corresponding robust Cohen distributions as rCD0 (n, k∆ω) =
X
i
λi |ST F Tρ (n, k∆ω) ∗ω Qi (k∆ω)|2 .
(19)
8
Here, λi represents the eigenvalues of the rotated kernel function in the time-lag domain of the particular member of the CD class, while Qi (k∆ω) represents the FT of the corresponding eigenvectors. Details about this decomposition and its associated quantities can be found in [17]. •
The LWD, proposed in [13], is an appropriate tool for nonlinear FM signals. This distribution can be evaluated recursively by using the WD or the SM. In a similar way, we propose to evaluate its robust counterpart based on the rWD or the rSM defined in (17) and (18). We define the robust LWD by using the recursive expression rLW DL (n, ω) =
Z
1 π
θ
P (θ)rLW DL/2 (n, ω + θ)rLW DL/2 (n, ω − θ)dθ,
(20)
where rLW D1 (n, ω) = rW D(n, ω) or rLW D1 (n, ω) = rSM(n, ω). •
The PWVD is a T-F distribution proposed to deal with polynomial FM signals. For multicomponent signals, this distribution suffers from the presence of cross-terms. To mitigate this problem, an implementation procedure for the fourth order PWVD was proposed in [14]. Its robust counterpart can be defined as rP W V D(n, k∆ω) =
X
l
P (l∆ω)rLW D2 (n, (k + l)∆ω)rW D(n, (k + [l/A])∆ω),
(21)
where A = 0.85/1.35 and [·] denotes rounding to the closest integer value. B. Example In this subsection, we illustrate the reason for introducing robust TFDs based on the robust STFT. Let us consider an FM signal given by n
o
f(t) = exp{jφ(t)} = exp j112π[t arctan(2πt) − ln(1 + 4π2 t2 )/4π] , for −1/2 ≤ t ≤ 1/2, sampled at ∆t = 1/512. Performance of the TFD used in the analysis will be measured based on its accuracy in estimating the signal’s IF given by2 ω(t) = φ0 (t) = 112π arctan(2πt). 2
(22)
Relation (22) represents semi-intuitive definition of the instantaneous frequency for complex-valued signals. We assume that
the complex-valued signal is produced as an analytical extension of the corresponding real-valued signal. The analytical extension is valid mapping to complex domain signals that satisfy conditions defined by Vakman [18]. In order to meet these requirements
9
Peak of the TFD is used as an IF estimator3 ,
ω ˆ (t) = arg max T F (t, ω). ω
(23)
The impulsive noise is added to the noiseless signal, x(t) = f(t) + ε3 (ν13 (t) + jν23 (t)), where νi (t), i = 1, 2, are zero-mean Gaussian noises whose variance is equal to 1. We assume ν1 (t) and ν2 (t) to be uncorrelated. This impulsive noise model was proposed in [8]. We compared the robust TFDs realized by using the robust STFT with those realized by using the LAF optimization problem. Some of these results are depicted in Figs.1,2, for ε = 0.65 and ε = 1, respectively. The following TFDs are shown: standard and median-filter based STFTs (the squared magnitude of the robust STFT is depicted, i.e., robust spectrogram), rWDs and rPWVDs calculated based on the LAF and the rSTFT. It can be seen that in a moderate impulse noise environment (Fig.1) both rWD forms are of the same accuracy. However, the rPWVD produced by using the LAF is inaccurate. In a heavier noise (Fig.2) both robust TFDs calculated by using the LAF are inaccurate, while the rTFDs calculated by using the rSTFT in initial stage are accurate, indicating that these robust TFDs are more resistant to the high impulsive noise. To validate and quantify these results, a statistical performance analysis will be presented in the next section. An explanation of the above behavior can be given as follows. The performance of all robust techniques based on the robust statistics concept degrade as the noise impulse rate increases [11]. For instance, the number of noise impulses in the LAF expression, used for the robust WD, is approximately double of that in the signal used for the robust STFT. For the robust PWVD, the number is much higher. Thus, this representation yields the worst results when the noise we assume that the amplitude in our signal is slow varying comparing to the signal phase |A0 (t)| ¿ |φ0 (t)|. Other forms of complex-domain signal representations are studied in [19] and references therein. 3
TFDs are commonly used as the IF estimators. The simplest used approach is based on the maxima of the T-F representation.
Namely, the standard T-F representations are usually highly concentrated around the IF of the FM signals. The robust TFDs used as the IF estimators should meet the same condition. Katkovnik studied the robust periodogram as an IF estimator [24]. He developed asymptotic expressions for bias and variance for this transform. He has also shown that the robust periodogram has the same bias as the standard one. It can be proven that the same conclusion holds in the case of robust TFDs. This topic is briefly addressed in the Appendix. It can be concluded that the robust TFDs are also highly concentrated along the IF, and that they can be good IF estimators. Detailed numerical study of the robust WD as an IF estimator has been done in [20].
10
Fig. 1. TFDs for moderate impulse noise environment ε = 0.65: (a) Standard STFT; (b) Median-filter based rSTFT; (e) rWD calculated by using the median filter and LAF; (f) rWD calculated by using median filter form of STFT; (i) rPWVD calculated by using median filter and LAF; (j) rPWVD calculated by using median filter form of STFT. Third and fourth columns - IF estimates obtained by using the TFD depicted in the first and second column.
increases. Therefore, we propose elimination of the impulse noise influence by using the rSTFT in initial stage, and after that applying well known procedures for realization of the higher order TFDs based on the STFT. C. Recursive Procedure for the rSTFT Realization Evaluation of various forms of the robust STFT requires computationally demanding iterative or sorting procedures for each point in the T-F domain. To avoid this situation, a recursive procedure to evaluate the robust STFT is proposed [21]. Consequently, the different robust TFDs discussed previously can be implemented in real-time. Let us first recall the recursive procedure to evaluate the standard STFT. A (discrete domain) recursive expression of the standard STFT is given by: ST F TS (n + 1, k∆ω) = [ N1 x(n + N) − N1 x(n) + ST F TS (n, k∆ω)]ej2πk/N . After a number of recursions, and due to accumulation of the roundoff error in the registers, the calculation of the STFT by using the FFT algorithm should be repeated. This
11
Fig. 2. TFDs for high impulse noise environment ε = 1: (a) Standard STFT; (b) Median-filter based rSTFT; (e) rWD calculated by using median filter and LAF; (f) rWD calculated by using median filter form of STFT; (i) rPWVD calculated by using median filter and LAF; (j) rPWVD calculated by using median filter form of STFT. Third and fourth columns - IF estimates obtained by using the TFD depicted in the first and second column.
recursive procedure can be used for a large number of consecutive samples with an acceptable small error. For the robust STFT, which can be evaluated by using any of the forms discussed earlier, the recursive expression cannot be directly used because, in this case, the sample x(n + N) may be heavily corrupted by the impulse noise. To remedy this situation, we need to estimate the noiseless signal samples from the noise observation. Common tool for denoising in the impulse noise environment are median and other related filters [22]. These filters are lowpass and they are not suitable for signals with high spectral content. To this aim robust filters based on the median and myriad filters, admitting negative weights, are proposed in [23]. Here, de-noising is performed on the robust filter designed in the spectral domain [16]. The de-noised (enhanced) signal can be calculated as [16] fρ (n + k) =
N 2π
Z
π
−π
ST F Tρ (n, ω)ejωk dω, k = 0, 1, ..., N − 1
(24)
12
1. Evaluate, for the first sample, the initial robust STFT ST F Tρ (n, ω) by using either (8), (9) or (15). 2. Determine the inverse function fρ (n + k) using (24). Set r = 0 and p = 0. 3. If |x(n + N)| ≤ (1 + η) max |fρ (n + k)| k ∈ [0, N) is satisfied, then, set fρ (n + N) = x(n + N) and p = 0. Otherwise, set fρ (n + N) = fρ (n + N − 1) and p = p + 1. 4. Evaluate the robust STFT at the new time instant n + 1 as ST F Tρ (n + 1, k∆ω) = [ N1 fρ (n + N) −
1 j2πk/N N fρ (n) + ST F Tρ (n, k∆ω)]e
and set r = r + 1 and n = n + 1.
5. If r ≤ R or p ≤ P then and go to Step 3, else go to Step 1. TABLE II Recursive procedure for the evaluation of the robust short-time Fourier transform.
where it is assumed that this inversion significantly reduces the impulsive noise influence. Now, if the magnitude of x(n + N) is smaller than, or close to, the maximum of |fρ (n + k)| (i.e., |x(n + N)| ≤ max |fρ (n + k)| or |x(n + N)| ≈ max |fρ (n + k)|), then we can assume that x(n + N) is not heavily corrupted by the impulsive noise. Otherwise, we assume that it is corrupted by the noise, and we should take the previous known sample fρ (n + N − 1) as an estimate of fρ (n + N). The recursive procedure is given in Table II. Note that R, the maximal number of recurrences used to compute ST F Tρ (n, ω), is now smaller than in the case of the standard STFT due to the constraint imposed by Step 3. In the previous examples, we set R = N/4. In some situations, the value of |x(n + N)| could be larger than magnitude of the samples in the interval [n, n + N − 1]. To discriminate this sample from an impulsive noise, we set larger margins for the signal in Step 3. In our examples, we set η = 0.25. In other situations, it could happen that there is no signal within an interval of time, and suddenly a signal sample appears. This sample might be missidentified by the algorithm as an impulse noise. In this situation, if the signal magnitude is larger than (1 + η) max |fρ (n + k)| k ∈ [0, N) for P consecutive samples, then, we go back to Step 1 and recalculate the initial robust STFT for the new sample. To validate the proposed recursive procedure, we compare it with the median-filter procedure in the evaluation of the robust STFT. The same noisy signal is considered as in the previous example. Results of the
13
Fig. 3. Robust STFTs based on the median filter and recursive realization: First row - impulse noise ε = 0.65; Second row - heavy impulse noise ε = 1; First column - median based rSTFT; Second column - recursive calculation; Third and fourth columns - IF estimates obtained by using the TFD depicted in the first and second column.
comparisons are displayed in Fig.3. The rows of this figure correspond to ε = 0.65 and ε = 1, respectively. The first two columns correspond to the median-filter and recursive robust STFTs. The IF estimates based on the peaks of these two robust STFT forms are displayed in the third and fourth columns. Under the same noise conditions, we have also considered a time-varying amplitude signal f(t) = exp(j64πt2 ) exp(−12t2 ). The recursive and median-based robust STFTs are shown in Fig.4. From these examples we can conclude that the robust STFTs, calculated directly and recursively, are very similar. Thus, the robust STFT form computed recursively can be used in the realization of robust TFDs. The overall computational complexity, in the recursive procedure, is decreased about R times, compared to the other methods. In addition, the recursive procedure is appropriate and useful to realize robust TFDs in real time. We should note that some delay may appear, since the calculation of initial STFT, in Step 1 of the algorithm, is done by time-consuming procedures. After that, the calculation can be performed according to the signal sampling period. IV. Statistical Comparison Study In this section, we compare various forms of the robust STFT implemented by using the direct procedure, with the robust STFT implemented by using the recursive procedure. Also, the robust TFDs, realized based on the LAF optimization problem are compared with those realized based the robust STFT. The comparison
14
Fig. 4.
TFDs based on the median filter and recursive realization: (a) Median filter form, ε = 0.65; (b) Recursive
realization, ε = 0.65; (c) Median filter form, ε = 1; (d) Recursive realization, ε = 1.
is given in terms of the mean-square-error (MSE) in the IF estimates. The considered signal was the same as in the previous examples. Signal is considered in the interval T ∈ [−1, 1] sampled with ∆t = 1/512. The Hanning window with Nw = 256 samples is used. In each experiment, the noiseless signal is embedded in the additive noise ν(n) = β(ν1 (n) + jν2 (n)) + ε3 (ν33 (n) + jν43 (n)). We assume that the noise components νi (n), i = 1, 2, 3, 4, are mutually uncorrelated, zero-mean, white, Gaussian with a unit variance. Two series of experiments are performed. In the first one, only the impulsive noise component is considered (i.e., β = 0) with ε ∈ [0, 2]. In the second series, a mixture of Gaussian and impulsive noise is considered, β = 0.4 and ε ∈ [0, 2]. For each noise amount we performed 50 independent trials. The following TFDs are considered in the comparisons: the STFT forms (standard and robust calculated by the iterative, the median-filter, the L-filter (α = 3/8), and the recursive procedures); the WD forms; the SM forms; the LWD forms for L = 2; the PWVD forms. Quadratic and higher order TFDs are calculated based on the LAF (Section II) and based on the STFT (Section III.A, III.C). Some of the most important results are displayed in Fig.5. The MSE of IF estimates is plotted against the impulsive noise amount ε. The first row of the figure shows that the robust STFT forms, including those calculated recursively, perform better than the standard STFT, even for a small amount of noise. The second row of the figure displays the results based on the median STFT form, the median WD and PWVD forms realized by using the LAF and robust STFT. As expected, for small ε the robust WD and PWVD forms outperform the robust STFT. However, for larger values of ε, the robust TFDs based on the LAF are inaccurate, while their counterparts calculated by using the robust STFT outperform the robust STFT for all values of ε. These results clearly confirm that robust TFDs realized by using the robust STFT are not only
15
0 10log10MSE
0 10log10MSE −10
−10
−20
−20 STFTs STFTr STFTm STFTl STFTlrec
−30 −40 −50
0
10
0.5
1
1.5
−40 ε a) 2
0
10
0.5
1
ε b) 2
1.5
10log10MSE
0
−10
−10
−20
−20 STFTm WDm PWVDm rWD rPWVD
−30 −40 −50
−50
10log10MSE
0
STFTs STFTr STFTm STFTl STFTlrec
−30
0
0.5
1
STFTm WDm PWVDm rWD rPWVD
−30 −40
1.5
ε c) 2
−50
0
0.5
1
1.5
ε d) 2
Fig. 5. Statistical analysis of the IF estimator based on the robust TFDs: (a) robust STFT forms and impulse noise; (b) robust STFT forms and mixed Gaussian and impulse noise environment; (c) robust TFDs and impulse noise environment; (d) robust TFDs and mixed Gaussian and impulse noise environment.
accurate but are also more resistant to impulsive noise than those realized by using the LAF. V. Conclusion The robust TFDs realized based on the robust STFT are presented. These TFDs have better T-F resolution than the robust STFT itself. They have stronger resistance to impulsive noise than their counterparts realized by using the LAF based optimization approach. Validity of the presented analysis is demonstrated on several examples. A recursive procedure for the robust STFT computation was also proposed. Results obtained by using this procedure were as accurate as in the other methods of the robust STFT computation. The proposed recursive procedure significantly improves the computational efficiency of the algorithm and offers the possibility to evaluate the robust TFDs in real-time. VI. Appendix Spectral estimation of the complex-valued harmonic f(t) = A exp(jφ(t)) embedded in a white impulse noise ν(t) based on the robust periodogram has been considered in [24]. The robust periodogram is squared
16
magnitude of the robust DFT defined in an implicit manner (4) and determined from the iterative algorithm [24]. Assume the following. 1. The window is bounded, real-valued and even function: w(t) = 0 for |t| ≥ T /2 and w(t) = w∗ (−t). 2. Loss function F () is convex, bounded and twice differentiable function. 3. Loss function F () and noise pdf g() should satisfy the following requirements:
0
Robust Time-Frequency Distributions Based on the Robust Short Time Fourier Transform Igor Djurovi´c1 , LJubiˇsa Stankovi´c1 and Braham Barkat2 1
Elektrotehniˇcki fakultet Podgorica, University of Montenegro, Podgorica, Montenegro. Tel./Fax. +381 81 244 921. E-mail: [email protected], [email protected].
2
Nanyang Technology University, School of EEE, Block S2, Nanyang Avenue, Singapore 639798. E-mail: [email protected]. Abstract Design of time-frequency distributions (TFDs) that are robust to the impulse noise influence is considered. The robust TFDs based on the robust short-time Fourier transform (STFT) are proposed. An efficient procedure to evaluate the robust STFT is given. Robust TFDs based on the robust STFT have better energy concentration around the signal instantaneous frequency (IF) than the robust STFT itself. Also, these TFDs are more resistant to higher impulse noise than the robust TFDs obtained using the local autocorrelation function (LAF) based minimization problem.
I. Introduction An overview of time-frequency distributions (TFDs) that have attracted significant attention in the signal processing community can be found in [1], [2]. Influence of Gaussian noise on TFDs has been extensively studied [3]. The instantaneous frequency (IF), a concept closely related to time-frequency (T-F) analysis, has attracted wide interest. The IF estimation of signals embedded in Gaussian noise has been addressed in several research works, among them [4], [5]. However, only few research works have dealt with other types of noise [6]. In numerous practical situations, environment cannot be modeled with Gaussian noise. Namely, due to natural or man-made phenomenon, the resulting noise could have rare but very strong values. This kind of noise is referred to as impulse noise. More details on the impulse noise and impulse noise modeling can be
2
found in [7]. TFDs are sensitive to this kind of noise. An intensive research effort has been put into design of TFDs that are robust to such a noise [8]-[10]. These TFDs are based on the robust statistics concept, introduced by Huber [11]. This concept has been successfully used in various signal processing applications. Katkovnik introduced this concept into spectral analysis by developing the robust M-Fourier transform (MFT) [8]. In the sequel, transforms that are defined on the basis of the robust statistics concept will be referred to as ‘robust’. The robust M-FT has been defined in an implicit form. An iterative procedure is required for its evaluation. In order to avoid iterative procedures, the marginal-median filter form of the robust FT was proposed in [9]. Both forms of the robust FTs yield accurate spectra estimates for stationary signals embedded in impulse noises. They are almost as accurate as the standard FT for Gaussian noise environment. Generalization to the T-F analysis of non-stationary signals corrupted by impulse noise has been proposed in [9], [10]. The robust Wigner distribution (WD) can be evaluated as a median of the modulated local autocorrelation function (LAF) samples. Note that all minimization problems based on the Hermitian LAF yield real-valued robust TFDs [9]. The L-filter form of robust TFDs is proposed for handling signals corrupted by a mixture of Gaussian and impulse noise in [10]. Higher-order TFDs yield a better energy concentration around their IF than the STFT in the case of noiseless signals. In this paper, we will show that the same fact holds for the robust TFDs. In addition, IF estimates based on peaks of the higher-order robust TFDs have smaller bias as compared to those based on the robust STFT peaks. However, higher-order TFDs are more sensitive to the impulsive noise than the robust STFT. To mitigate this, we propose a new form of the robust TFDs based on the robust STFT. In this way, a significant improvement, comparing to the results obtained by using the robust TFDs based on the LAF, is achieved. In addition, we propose a recursive technique to evaluate the robust STFT. As a consequence, the robust TFDs based on the robust STFT can be efficiently calculated in real-time. After reviewing the existing robust TFDs in Section II, we introduce robust TFDs based on the robust STFT, as well as the recursive procedure to evaluate them, in Section III. A comparison study of various robust TFDs, using IF estimation, is presented in Section IV.
3
II. A Review of the Robust FT and TFDs A. Robust FT Consider discretized version of the useful signal f(n) that is embedded in a white noise ν(n), x(n) = f(n) + ν(n). The standard FT of noisy signal x(n),
X(ω) =
1 N
XN−1 n=0
x(n)e−jωn = mean{x(n)e−jωn | n ∈ [0, N − 1]},
(1)
can be defined as a solution of the optimization problem1 :
X(ω) = arg min J(m; ω) m
(2)
where the penalty function can be written as
J(m; ω) =
XN−1 n=0
F (x(n)e−jωn − m),
(3)
for the loss-function F (e) = |e|2 . The standard DFT, X(ω), is the ML estimate of the useful signal spectra F (ω) for a Gaussian noise environment. It is sensitive to the impulse noise. For a known probability density function (pdf) of the noise the ML estimate of the FT can be determined. The ML estimates are sensitive to variations of the noise pdf. To avoid this difficulty, Huber introduced the M-estimation concept whose main idea is as follows. A class of noises is considered. The ‘worst’ noise is defined as the one whose pdf has the heaviest tail. The ML estimate for this noise is the M-estimate for the used class. The M-estimate produces reasonable accuracy for the entire class. The Laplacian noise was found as the worst one for many impulse noise classes considered in practice. The loss-function F (e) = |e|, that produces the ML estimate for the Laplace noise, results in robust behavior for almost all other impulse disturbances. The robust M-FT in implicit form, obtained as the solution of the 1
Minimizing the penalty function with respect to the complex-valued argument is done in the following manner. Penalty function
J (m) is considered as the 2D function of real-valued arguments J(mr , mi ; ω) where mr and mi are real and imaginary part of m. The searched solution is formed as X(ω) = R(ω) + jI(ω), where (R(ω), I(ω)) are solution of equation (R(ω), I(ω)) = arg J (mr , mi ; ω).
min
(mr ,mi )
4
minimization problem (2)-(3) for F (e) = |e|, is given by [8] XR (ω) =
PN−1
x(n)e−jωn n=0 |x(n)e−jωn −XR (ω)| PN−1 1 n=0 |x(n)e−jωn −XR (ω)|
.
(4)
To avoid computational problems associated with iterative procedures, required for solving (4), an alternative form of the robust M-FT, based on the loss function F (e) = | Re(e)| + | Im(e)|, was proposed in [9]. This alternative form, called the marginal median form of the M-FT, is given by [9] n
o
n
o
XM (ω) = median Re{x(n)e−jωn }| n ∈ [0, N − 1] + j median Im{x(n)e−jωn }| n ∈ [0, N − 1] .
(5)
Transforms (4) and (5) behave similarly. Their accuracy is high for impulse noise. For Gaussian noise they are slightly worse than the standard FT. The median-filter form does not suffer from convergence problem that may exist in (4). B. Robust TFDs Robust Short Time Fourier Transform: The robust STFT forms, which can be seen as an extension of the robust FT, are defined by using the optimization problem ST F T (n, ω) = arg min J(m; n, ω) = arg min m
m
XN−1 k=0
F (x(n + k)e−jωk − m).
(6)
For the loss functions F (e) = |e|2 , F (e) = |e| and F (e) = | Re(e)| + | Im(e)|, yield the standard STFT (ST F TS (n, ω)), the robust STFT in the implicit form (ST F TR (n, ω)), and the robust STFT in the marginalmedian form (ST F TM (n, ω)), respectively. These transforms are given by ST F TS (n, ω) = ST F TR (n, ω) =
XN−1
x(n + k)e−jωk k=0 PN−1 x(n+k)e−jωk k=0 |x(n+k)e−jωk −ST F TR (n,ω)| PN−1 1 k=0 |x(n+k)e−jωk −ST F TR (n,ω)| 1 N
(7) (8)
ST F TM (n, ω) = median{Re{x(n + k)e−jωk }|k ∈ [0, N − 1]} +jmedian{Im{x(n + k)e−jωk }|k ∈ [0, N − 1]}.
(9)
Robust TFDs Based on the LAF: The robust WD was introduced in [9] as a solution of the minimization problem W D(n, ω) = arg min J(m; n, ω) = arg min m
m
XN−1 k=0
F (rxx (n, k)e−2jωk − m)
(10)
5
where the LAF is given by rxx (n, k) = x(n + k)x∗ (n − k). For F (e) = |e|2 , we obtain the standard WD (W DS (n, ω)) and for F (e) = |e|, we obtain the robust WD (W DR (n, ω)), i.e., W DS (n, ω) =
1 N
XN−1 k=0
rxx (n, k)e−2jωk = mean{rxx (n, k)e−2jωk | k ∈ [0, N − 1]},
W DR (n, ω) = median{Re{rxx (n, k)e−2jωk }| k ∈ [0, N − 1]}.
(11) (12)
Here, the loss function F (e) = | Re(e)| + | Im(e)| yields the same solution as the one obtained for F (e) = |e|. The LAFs, rxx (n, k), of some common TFDs are given in Table I [12]-[15], where the acronyms CD, LWD and PWVD stand for the Cohen class of distributions, the L-Wigner distribution and the polynomial Wigner-Ville ∗ (n, −k), the distribution, respectively. It can be shown that for any Hermitian LAF, i.e., rxx (n, k) = rxx
optimization problem (10) results in a real solution for any even loss function. Using the appropriate LAF expression, we define the standard TFD (T F (n, ω)) as the solution of the above optimization problem when F (e) = |e|2 . Similarly, we define the robust TFD (T FR (n, ω)) as the solution obtained for F (e) = |e|. These two quantities are equal to n
T F (n, ω) = mean rxx (n, k)e−2jωk | k ∈ [0, N − 1] n
o
(13) o
T FR (n, ω) = median Re{rxx (n, k)e−2jωk }| k ∈ [0, N − 1] .
Representation
LAF rxx (n, k)
WD
x(n + k)x∗ (n − k)
CD1
cT (n, k) ∗n x(n + k)x∗ (n − k)
LWD
xL (n + k/L)xL∗ (n − k/L)
PWVD2 q = 4
x2 (n + 0.675k)x2∗ (n − 0.675k)x∗ (n + 0.85k)x(n − 0.85k)
(14)
TABLE I 1
cT (n, k) real and symmetric kernel function in the time-lag domain.
2
the PWVD form used in the
paper.
L-filter Forms: The L-filters are effective in filtering stationary signals embedded in a mixture of Gaussian
6
and impulsive noise. They are also successfully applied to non-stationary signals [10] and for design of the robust filters of high-pass data [16]. The L-form of the STFT is given as
ST F TL (n, ω) =
XN−1 l=0
al [r(n, ω; l) + j · i(n, ω; l)],
(15)
where r(n, ω; l) and i(n, ω; l) are elements of the sets R(n, ω) = {Re{x(n + k)e−jωk }| k ∈ [0, N − 1]} and I(n, ω) = {Im{x(n + k)e−jωk }| k ∈ [0, N − 1]}, i.e., r(n, ω; l) ∈ R(n, ω), and i(n, ω; l) ∈ I(n, ω), sorted into non-decreasing order as r(n, ω; l) ≤ r(n, ω; l + 1), i(n, ω; l) ≤ i(n, ω; l + 1). The L-filter forms of the TFDs given in Table I, can be written as
T FL (n, ω) =
XN−1 l=0
al r(n, ω; l),
(16)
where r(n, ω; l) are sorted elements of the set R(n, ω) = {Re{rxx (n, k)e−2jωk }|k ∈ [0, N − 1]}. The filter coefficients al , l = 0, 1, ..., N − 1, are usually determined to satisfy the constraints
PN−1 l=0
al = 1 and al =
aN−l−1 . We will use the α-trimmed form of the L-filter whereby the coefficients are given as: al = 1/[N − α(N − 2)], l ∈ [(N − 2)α, (−N + 2)α + N − 1], α ∈ [0, 0.5], and al = 0 elsewhere. Note that, for α = 0 and α = 0.5, the standard and the median forms follow, respectively. III. Robust TFDs Based on the Robust STFT Realization of standard TFDs, based on the standard STFT, is a well known and described in numerous research papers [13]-[17]. This realization technique is used for many reasons, such as: computational efficiency, to avoid aliasing effects in TFDs, or to produce highly concentrated TFDs with reduced interferences. Since the STFT can be calculated recursively, it means that different TFDs can be realized in real-time. In addition to the just stated benefits, our main objective in this section is to introduce robust TFDs based on the robust STFT, that have stronger resistance to the impulse noise than the previously reviewed robust TFDs. To this aim, a recursive realization of the robust STFT is proposed in this section. As we have previously seen, the robust STFT can be realized by using expressions (8), (9) or (15). In the sequel, we can use any of these forms to realize the robust TFDs. These robust TFDs will be more resilient to impulsive noise than those realized by using the minimization problems discussed earlier. The robust STFT used to realize the robust TFDs will
7
be denoted by ST F Tρ (t, ω), where ρ can be any of the letters M, R or L, used to refer to the robust STFT forms in (8), (9) and (15). The notation of the various robust TFDs, obtained by using the robust STFT, will be preceded by r. A. Proposed Transforms In the light of what we discussed above, and using the same technique to realize standard TFDs based on the standard STFT, we itemize below the different forms of the robust TFDs based on the robust STFT. •
The windowed robust STFT can be calculated by rST F T (n, ω) = ST F Tρ (n, ω) ∗ω W (ω), where W (ω) is the FT of the desired window function and ∗ω represents the convolution in frequency. For instance, for a Hanning window (discrete domain), this convolution reduces to rST F T (n, k∆ω) =
P1
i=−1 ci ST F Tρ (n, (k+
i)∆ω), where c0 = 1/2, c±1 = 1/4, and ∆ω is the frequency step. Similar expressions can be derived for other window types. •
The robust forms of the WD based on the robust STFT are determined as [12] 1 rW D(n, ω) = π
Z
θ
ST F Tρ (n, ω + θ)ST F Tρ∗ (n, ω − θ)dθ.
(17)
The pseudo form of the rWD can be determined in the same manner as the windowed form of ST F Tρ (n, ω). •
The S-method (SM) [12] combines favorable properties of the STFT (non-aliasing for the Nyquist rate sampled signal and cross-terms reduction) and the WD (high concentration). We define its robust counterpart as rSM (n, k) =
1 π
Z
θ
P (θ) ST F Tρ (n, ω + θ) ST F Tρ∗ (n, ω − θ)dθ,
(18)
where P (θ) is the window function in the frequency domain. Wider window means better T-F concentration, while narrower window means better cross-terms interference suppression. Full details on the window length determination can be found in [12]. •
Realization of the quadratic Cohen’s class of TFDs by using the standard STFT was introduced in [17]. We define their corresponding robust Cohen distributions as rCD0 (n, k∆ω) =
X
i
λi |ST F Tρ (n, k∆ω) ∗ω Qi (k∆ω)|2 .
(19)
8
Here, λi represents the eigenvalues of the rotated kernel function in the time-lag domain of the particular member of the CD class, while Qi (k∆ω) represents the FT of the corresponding eigenvectors. Details about this decomposition and its associated quantities can be found in [17]. •
The LWD, proposed in [13], is an appropriate tool for nonlinear FM signals. This distribution can be evaluated recursively by using the WD or the SM. In a similar way, we propose to evaluate its robust counterpart based on the rWD or the rSM defined in (17) and (18). We define the robust LWD by using the recursive expression rLW DL (n, ω) =
Z
1 π
θ
P (θ)rLW DL/2 (n, ω + θ)rLW DL/2 (n, ω − θ)dθ,
(20)
where rLW D1 (n, ω) = rW D(n, ω) or rLW D1 (n, ω) = rSM(n, ω). •
The PWVD is a T-F distribution proposed to deal with polynomial FM signals. For multicomponent signals, this distribution suffers from the presence of cross-terms. To mitigate this problem, an implementation procedure for the fourth order PWVD was proposed in [14]. Its robust counterpart can be defined as rP W V D(n, k∆ω) =
X
l
P (l∆ω)rLW D2 (n, (k + l)∆ω)rW D(n, (k + [l/A])∆ω),
(21)
where A = 0.85/1.35 and [·] denotes rounding to the closest integer value. B. Example In this subsection, we illustrate the reason for introducing robust TFDs based on the robust STFT. Let us consider an FM signal given by n
o
f(t) = exp{jφ(t)} = exp j112π[t arctan(2πt) − ln(1 + 4π2 t2 )/4π] , for −1/2 ≤ t ≤ 1/2, sampled at ∆t = 1/512. Performance of the TFD used in the analysis will be measured based on its accuracy in estimating the signal’s IF given by2 ω(t) = φ0 (t) = 112π arctan(2πt). 2
(22)
Relation (22) represents semi-intuitive definition of the instantaneous frequency for complex-valued signals. We assume that
the complex-valued signal is produced as an analytical extension of the corresponding real-valued signal. The analytical extension is valid mapping to complex domain signals that satisfy conditions defined by Vakman [18]. In order to meet these requirements
9
Peak of the TFD is used as an IF estimator3 ,
ω ˆ (t) = arg max T F (t, ω). ω
(23)
The impulsive noise is added to the noiseless signal, x(t) = f(t) + ε3 (ν13 (t) + jν23 (t)), where νi (t), i = 1, 2, are zero-mean Gaussian noises whose variance is equal to 1. We assume ν1 (t) and ν2 (t) to be uncorrelated. This impulsive noise model was proposed in [8]. We compared the robust TFDs realized by using the robust STFT with those realized by using the LAF optimization problem. Some of these results are depicted in Figs.1,2, for ε = 0.65 and ε = 1, respectively. The following TFDs are shown: standard and median-filter based STFTs (the squared magnitude of the robust STFT is depicted, i.e., robust spectrogram), rWDs and rPWVDs calculated based on the LAF and the rSTFT. It can be seen that in a moderate impulse noise environment (Fig.1) both rWD forms are of the same accuracy. However, the rPWVD produced by using the LAF is inaccurate. In a heavier noise (Fig.2) both robust TFDs calculated by using the LAF are inaccurate, while the rTFDs calculated by using the rSTFT in initial stage are accurate, indicating that these robust TFDs are more resistant to the high impulsive noise. To validate and quantify these results, a statistical performance analysis will be presented in the next section. An explanation of the above behavior can be given as follows. The performance of all robust techniques based on the robust statistics concept degrade as the noise impulse rate increases [11]. For instance, the number of noise impulses in the LAF expression, used for the robust WD, is approximately double of that in the signal used for the robust STFT. For the robust PWVD, the number is much higher. Thus, this representation yields the worst results when the noise we assume that the amplitude in our signal is slow varying comparing to the signal phase |A0 (t)| ¿ |φ0 (t)|. Other forms of complex-domain signal representations are studied in [19] and references therein. 3
TFDs are commonly used as the IF estimators. The simplest used approach is based on the maxima of the T-F representation.
Namely, the standard T-F representations are usually highly concentrated around the IF of the FM signals. The robust TFDs used as the IF estimators should meet the same condition. Katkovnik studied the robust periodogram as an IF estimator [24]. He developed asymptotic expressions for bias and variance for this transform. He has also shown that the robust periodogram has the same bias as the standard one. It can be proven that the same conclusion holds in the case of robust TFDs. This topic is briefly addressed in the Appendix. It can be concluded that the robust TFDs are also highly concentrated along the IF, and that they can be good IF estimators. Detailed numerical study of the robust WD as an IF estimator has been done in [20].
10
Fig. 1. TFDs for moderate impulse noise environment ε = 0.65: (a) Standard STFT; (b) Median-filter based rSTFT; (e) rWD calculated by using the median filter and LAF; (f) rWD calculated by using median filter form of STFT; (i) rPWVD calculated by using median filter and LAF; (j) rPWVD calculated by using median filter form of STFT. Third and fourth columns - IF estimates obtained by using the TFD depicted in the first and second column.
increases. Therefore, we propose elimination of the impulse noise influence by using the rSTFT in initial stage, and after that applying well known procedures for realization of the higher order TFDs based on the STFT. C. Recursive Procedure for the rSTFT Realization Evaluation of various forms of the robust STFT requires computationally demanding iterative or sorting procedures for each point in the T-F domain. To avoid this situation, a recursive procedure to evaluate the robust STFT is proposed [21]. Consequently, the different robust TFDs discussed previously can be implemented in real-time. Let us first recall the recursive procedure to evaluate the standard STFT. A (discrete domain) recursive expression of the standard STFT is given by: ST F TS (n + 1, k∆ω) = [ N1 x(n + N) − N1 x(n) + ST F TS (n, k∆ω)]ej2πk/N . After a number of recursions, and due to accumulation of the roundoff error in the registers, the calculation of the STFT by using the FFT algorithm should be repeated. This
11
Fig. 2. TFDs for high impulse noise environment ε = 1: (a) Standard STFT; (b) Median-filter based rSTFT; (e) rWD calculated by using median filter and LAF; (f) rWD calculated by using median filter form of STFT; (i) rPWVD calculated by using median filter and LAF; (j) rPWVD calculated by using median filter form of STFT. Third and fourth columns - IF estimates obtained by using the TFD depicted in the first and second column.
recursive procedure can be used for a large number of consecutive samples with an acceptable small error. For the robust STFT, which can be evaluated by using any of the forms discussed earlier, the recursive expression cannot be directly used because, in this case, the sample x(n + N) may be heavily corrupted by the impulse noise. To remedy this situation, we need to estimate the noiseless signal samples from the noise observation. Common tool for denoising in the impulse noise environment are median and other related filters [22]. These filters are lowpass and they are not suitable for signals with high spectral content. To this aim robust filters based on the median and myriad filters, admitting negative weights, are proposed in [23]. Here, de-noising is performed on the robust filter designed in the spectral domain [16]. The de-noised (enhanced) signal can be calculated as [16] fρ (n + k) =
N 2π
Z
π
−π
ST F Tρ (n, ω)ejωk dω, k = 0, 1, ..., N − 1
(24)
12
1. Evaluate, for the first sample, the initial robust STFT ST F Tρ (n, ω) by using either (8), (9) or (15). 2. Determine the inverse function fρ (n + k) using (24). Set r = 0 and p = 0. 3. If |x(n + N)| ≤ (1 + η) max |fρ (n + k)| k ∈ [0, N) is satisfied, then, set fρ (n + N) = x(n + N) and p = 0. Otherwise, set fρ (n + N) = fρ (n + N − 1) and p = p + 1. 4. Evaluate the robust STFT at the new time instant n + 1 as ST F Tρ (n + 1, k∆ω) = [ N1 fρ (n + N) −
1 j2πk/N N fρ (n) + ST F Tρ (n, k∆ω)]e
and set r = r + 1 and n = n + 1.
5. If r ≤ R or p ≤ P then and go to Step 3, else go to Step 1. TABLE II Recursive procedure for the evaluation of the robust short-time Fourier transform.
where it is assumed that this inversion significantly reduces the impulsive noise influence. Now, if the magnitude of x(n + N) is smaller than, or close to, the maximum of |fρ (n + k)| (i.e., |x(n + N)| ≤ max |fρ (n + k)| or |x(n + N)| ≈ max |fρ (n + k)|), then we can assume that x(n + N) is not heavily corrupted by the impulsive noise. Otherwise, we assume that it is corrupted by the noise, and we should take the previous known sample fρ (n + N − 1) as an estimate of fρ (n + N). The recursive procedure is given in Table II. Note that R, the maximal number of recurrences used to compute ST F Tρ (n, ω), is now smaller than in the case of the standard STFT due to the constraint imposed by Step 3. In the previous examples, we set R = N/4. In some situations, the value of |x(n + N)| could be larger than magnitude of the samples in the interval [n, n + N − 1]. To discriminate this sample from an impulsive noise, we set larger margins for the signal in Step 3. In our examples, we set η = 0.25. In other situations, it could happen that there is no signal within an interval of time, and suddenly a signal sample appears. This sample might be missidentified by the algorithm as an impulse noise. In this situation, if the signal magnitude is larger than (1 + η) max |fρ (n + k)| k ∈ [0, N) for P consecutive samples, then, we go back to Step 1 and recalculate the initial robust STFT for the new sample. To validate the proposed recursive procedure, we compare it with the median-filter procedure in the evaluation of the robust STFT. The same noisy signal is considered as in the previous example. Results of the
13
Fig. 3. Robust STFTs based on the median filter and recursive realization: First row - impulse noise ε = 0.65; Second row - heavy impulse noise ε = 1; First column - median based rSTFT; Second column - recursive calculation; Third and fourth columns - IF estimates obtained by using the TFD depicted in the first and second column.
comparisons are displayed in Fig.3. The rows of this figure correspond to ε = 0.65 and ε = 1, respectively. The first two columns correspond to the median-filter and recursive robust STFTs. The IF estimates based on the peaks of these two robust STFT forms are displayed in the third and fourth columns. Under the same noise conditions, we have also considered a time-varying amplitude signal f(t) = exp(j64πt2 ) exp(−12t2 ). The recursive and median-based robust STFTs are shown in Fig.4. From these examples we can conclude that the robust STFTs, calculated directly and recursively, are very similar. Thus, the robust STFT form computed recursively can be used in the realization of robust TFDs. The overall computational complexity, in the recursive procedure, is decreased about R times, compared to the other methods. In addition, the recursive procedure is appropriate and useful to realize robust TFDs in real time. We should note that some delay may appear, since the calculation of initial STFT, in Step 1 of the algorithm, is done by time-consuming procedures. After that, the calculation can be performed according to the signal sampling period. IV. Statistical Comparison Study In this section, we compare various forms of the robust STFT implemented by using the direct procedure, with the robust STFT implemented by using the recursive procedure. Also, the robust TFDs, realized based on the LAF optimization problem are compared with those realized based the robust STFT. The comparison
14
Fig. 4.
TFDs based on the median filter and recursive realization: (a) Median filter form, ε = 0.65; (b) Recursive
realization, ε = 0.65; (c) Median filter form, ε = 1; (d) Recursive realization, ε = 1.
is given in terms of the mean-square-error (MSE) in the IF estimates. The considered signal was the same as in the previous examples. Signal is considered in the interval T ∈ [−1, 1] sampled with ∆t = 1/512. The Hanning window with Nw = 256 samples is used. In each experiment, the noiseless signal is embedded in the additive noise ν(n) = β(ν1 (n) + jν2 (n)) + ε3 (ν33 (n) + jν43 (n)). We assume that the noise components νi (n), i = 1, 2, 3, 4, are mutually uncorrelated, zero-mean, white, Gaussian with a unit variance. Two series of experiments are performed. In the first one, only the impulsive noise component is considered (i.e., β = 0) with ε ∈ [0, 2]. In the second series, a mixture of Gaussian and impulsive noise is considered, β = 0.4 and ε ∈ [0, 2]. For each noise amount we performed 50 independent trials. The following TFDs are considered in the comparisons: the STFT forms (standard and robust calculated by the iterative, the median-filter, the L-filter (α = 3/8), and the recursive procedures); the WD forms; the SM forms; the LWD forms for L = 2; the PWVD forms. Quadratic and higher order TFDs are calculated based on the LAF (Section II) and based on the STFT (Section III.A, III.C). Some of the most important results are displayed in Fig.5. The MSE of IF estimates is plotted against the impulsive noise amount ε. The first row of the figure shows that the robust STFT forms, including those calculated recursively, perform better than the standard STFT, even for a small amount of noise. The second row of the figure displays the results based on the median STFT form, the median WD and PWVD forms realized by using the LAF and robust STFT. As expected, for small ε the robust WD and PWVD forms outperform the robust STFT. However, for larger values of ε, the robust TFDs based on the LAF are inaccurate, while their counterparts calculated by using the robust STFT outperform the robust STFT for all values of ε. These results clearly confirm that robust TFDs realized by using the robust STFT are not only
15
0 10log10MSE
0 10log10MSE −10
−10
−20
−20 STFTs STFTr STFTm STFTl STFTlrec
−30 −40 −50
0
10
0.5
1
1.5
−40 ε a) 2
0
10
0.5
1
ε b) 2
1.5
10log10MSE
0
−10
−10
−20
−20 STFTm WDm PWVDm rWD rPWVD
−30 −40 −50
−50
10log10MSE
0
STFTs STFTr STFTm STFTl STFTlrec
−30
0
0.5
1
STFTm WDm PWVDm rWD rPWVD
−30 −40
1.5
ε c) 2
−50
0
0.5
1
1.5
ε d) 2
Fig. 5. Statistical analysis of the IF estimator based on the robust TFDs: (a) robust STFT forms and impulse noise; (b) robust STFT forms and mixed Gaussian and impulse noise environment; (c) robust TFDs and impulse noise environment; (d) robust TFDs and mixed Gaussian and impulse noise environment.
accurate but are also more resistant to impulsive noise than those realized by using the LAF. V. Conclusion The robust TFDs realized based on the robust STFT are presented. These TFDs have better T-F resolution than the robust STFT itself. They have stronger resistance to impulsive noise than their counterparts realized by using the LAF based optimization approach. Validity of the presented analysis is demonstrated on several examples. A recursive procedure for the robust STFT computation was also proposed. Results obtained by using this procedure were as accurate as in the other methods of the robust STFT computation. The proposed recursive procedure significantly improves the computational efficiency of the algorithm and offers the possibility to evaluate the robust TFDs in real-time. VI. Appendix Spectral estimation of the complex-valued harmonic f(t) = A exp(jφ(t)) embedded in a white impulse noise ν(t) based on the robust periodogram has been considered in [24]. The robust periodogram is squared
16
magnitude of the robust DFT defined in an implicit manner (4) and determined from the iterative algorithm [24]. Assume the following. 1. The window is bounded, real-valued and even function: w(t) = 0 for |t| ≥ T /2 and w(t) = w∗ (−t). 2. Loss function F () is convex, bounded and twice differentiable function. 3. Loss function F () and noise pdf g() should satisfy the following requirements:
0
