EURASIP Journal on Applied Signal Processing 2005:8, 1261–1268 c ○ 2005 Hindawi Publishing Corporation Efficient Analysis of Time-Varying Multicomponent Signals with Modified LPTFT

This paper presents efficient algorithms for the analysis of nonstationary multicomponent signals based on modified local polynomial time-frequency transform. The signals to be analyzed are divided into a number of segments and the desired parameters for computing modified the local polynomial time-frequency transform in each segment are estimated from polynomial Fourier transform in the frequency domain. Compared to other reported algorithms, the length of overlap between consecutive segments is reduced to minimize the overall computational complexity. The concept of adaptive window lengths is also employed to achieve a better time-frequency resolution for each component. Numerical simulations with synthesized multicomponent signals show that the proposed ones achieve better performance on instantaneous frequency estimation with greatly reduced computational complexity.


INTRODUCTION
Due to their superior performance in dealing with nonstationary signals, time-frequency transforms (TFTs) have found various applications in many areas including communications, multimedia, mechanics, and biology [1]. The most popular and simplest TFT is short-time Fourier transform (STFT) that has been widely used for many practical applications [1,2]. Nevertheless, the STFT suffers from low resolution when the analyzed signal is highly nonstationary. Local polynomial time-frequency transform (LPTFT), referred to as the generalization of STFT, was reported to provide high resolution for nonstationary signals [3,4] with a local polynomial function approximating to the frequency characteristics. Unfortunately, the estimation of a number of extra parameters required by LPTFT computation results in a heavy computational load. This is mainly due to the long overlap between consecutive signal segments for which the estimation process is implemented [4]. In order to reduce the computational complexity, attempts can be made to reduce the length of overlap between the consecutive segments. However, problems of reduced resolution in the time-frequency domain have to be solved by using more effective methods of window length selection. This paper presents analysis algorithms for time-varying multicomponent signals containing white Gaussian and/or impulse noises. Different from previously reported algorithms, the proposed modified local polynomial timefrequency transform (MLPTFT) reduces the overlap length between consecutive segments to minimize the number of segments to be processed. Effective methods of estimating the MLPTFT parameters from each signal segment are presented. Deterioration of resolution due to the reduction of overlap length is avoided by using adaptive window lengths. This paper is organized as follows. Section 2 provides a brief review of the LPTFT and introduces the MLPTFTs for the analysis of multicomponent signals containing different noises. Sections 3 and 4 present the details of parameter estimation and window length selection. Simulation results are given in Section 5 to show the effectiveness of the proposed algorithms.

Segmentation
The signal with noise to be analyzed is defined as  where w(t) represents white Gaussian noise or impulse noise and s(t) contains mono-or multi-nonstationary components in the time-frequency domain. In the paper, the impulse noise either belongs to Middleton class A [5] or is defined as w(t) = w 3 1 (t) + jw 3 2 (t) [4], where w 1 (t) and w 2 (t) represent mutually independent white Gaussian noises with unit variances. It is also assumed that the sampling frequency of the discrete data is normalized to be one Hz and parameter t takes integer values. The input signal x(t) is divided into many small segments with a window function h(τ) in the time domain. The jth signal segment is defined as where x is the function to return the largest integer that is equal to or smaller than x, N is the length of signal x(t), Q, which is assumed to be an odd number without loss of generality, is the length of the window h(τ) or equivalently the length of the signal segment, and α represents the length of the overlap between the consecutive signal segments. Figure 1 shows examples for α =0, Q−1, and (Q−1)/2 with Q = 5. Heavy computational complexity is needed for estimating the extra parameters required by LPTFT computation if the overlap length is large because the number of signal segments to be processed is accordingly increased.

The MLPTFT
The local polynomial time-frequency transform (LPTFT) of x(t) is defined as [3] LPTFT where h(τ) is the window function with length Q, and l(t) = [l 1 (t), . . . , l M−1 (t)] are the parameters related to the derivatives of the instantaneous frequency of x(t) [3]. The LPTFT is based on the idea of fitting an (M − 1)th-order polynomial function approximation of the frequency of x j defined in (2) with α = Q − 1 to determine the nonparametric characteristic of the signal [3]. In addition to the calculation of (3), other processing costs for the LPTFT are for the estimation of both the time-varying parameter l(t) and window length Q. It has been shown in [3] that the LPTFT can yield high resolution of the time-varying frequency provided that these parameters are accurately estimated and properly updated. The LPTFT cannot be directly used for signals containing multiple components because individual signal components have their own parameter l(t) and window length Q. We define the MLPTFT p for signals containing p components with sets of parameters L(t) : {l i (t); 1 ≤ i ≤ p} and window length Q : where is the scaling factor keeping the signal energy unchanged and · 2 is the second-norm operation in terms of τ. Now we consider the performance of the MLPTFT p for signals containing p nonparallel components with Mthorder polynomial phase defined as where A i is the amplitude of the ith component x i (t). It was shown that the optimal {l i (t)} for the ith component is given as [3] where where * is the convolution operator and FT τ represents the Fourier transform in terms of τ. The first term in (8) is the useful autoterm and the second one contains the undesirable cross-terms. Generally, the cross-terms have much smaller magnitudes, compared with that of the autoterm because the phase factor of e − j2π M m=2 [ki,m−kq,m](τ m /m!) is spread in the frequency domain [6]. The MLPTFT p generally has fewer crossterms than the bilinear TFT [1] and can be approximately viewed as the sum of the optimal LPTFTs of each component. Furthermore, the constant scaling factor a for each window in (4) or (8) keeps the signal energy ratio between components approximately unchanged after the MLPTFT p computation because the influence of the cross-terms is trivial. It is worth mentioning that a similar modified form of the LPTFT, which also uses summation of the LPTFTs with several parameters suitable for each component, was proposed in [4]. However, the estimation of L(t) is based on maximizing the concentration measure [4], which is different from our proposed MLPTFT p .

Robust modified LPTFT
The MLPTFT p is only suitable for signals containing Gaussian noise and achieves poor performance for signals containing impulse noise. It is known that for stationary signals, robust FT (RFT) [7,8] has been developed to deal with impulse noises. Similarly, the robust MLPTFT p (RMLPTFT p ) is defined for nonstationary signals with impulse noises. The MLPTFT p defined in (4) is also expressed as The RFT [7,8] can be conveniently used to replace the FT τ in (9) to define RMLPTFT p as follows. We use the suboptimal marginal-median form of the RFT [8] and replace FT τ with median operation as where median τ refers to selecting the median value with regard to τ on the real and imaginary parts, respectively. With the MLPTFT p and RMLPTFT p in (9) and (10) introduced for p-component nonstationary signals with different kinds of noises, major steps for the estimation of the parameter L(t) and window length Q are presented in detail in the following two sections, respectively.

ESTIMATION OF L(T)
This section considers methods estimating L(t) from the segments achieved by (2) to minimize the computational complexity without deteriorating the smoothness of the spectrum. As mentioned previously, the signal is generally divided into many segments and the L(t) estimation is based on the idea of finding an (M − 1)th-order polynomial function to approximate the frequency characteristics of each signal segment. In the previously reported methods, the overlap factor α equals Q − 1 which means that there are N segments of length Q for an N-point input sequence. In general, several MLPTFT p s with different sets of parameters are computed for each signal segment. For segment x j , for example, the L( j) that yields the maximum value [2] or values larger than a threshold [4] is selected. Because two consecutive signal segments overlap heavily, that is, two adjacent signal segments differ by only one data point, this method requires a large computational load [4].
To reduce the overall computational load, it is necessary to minimize the length of overlap between consecutive segments, such as α < Q − 1. For segment x j , the set of parameters L( j(Q − α + τ)) within the duration −(Q − α − 1)/2 ≤ τ ≤ (Q − α − 1)/2 is estimated simultaneously. As shown in Figure 1, the parameters for the shaded data intervals are estimated from the corresponding signal segment. For example, L(2), L(3), and L(4) are estimated from segment x 2 when α = 2 in Figure 1b. In this way, only N/(Q − α) instead of N signal segments are processed to acquire L(t) at all time instants. Generally, α controls the tradeoff between the computational load and the smoothness of the spectrum. When α = 0, there is no overlap and only N/Q segments are processed, which reduces the computational load Q times compared with that with α = Q − 1 in the previously reported method. In general, the MLPTFT p with L(t) estimated with α = 0 yields satisfied performance to achieve a good polynomial function approximation to the frequency components if the window length Q is small enough, which is further illustrated in the first experiment of Section 5. Otherwise, this arrangement may cause problems of unsmoothness of the frequency components when consecutive segments are connected because longer windows are prone to larger differences between the estimated parameters for the consecutive segments. Under this circumstance, the overlapping factor α is required to be increased. According to the types of noises encountered in practice, the following methods of parameter estimation are presented for each case.

(b) Estimation for signals with impulse noise
The robust PFT (RPFT) is derived for the estimation of coefficients from signals with impulse noise. The estimation of phase parameters is the same as that presented in the previous section except that the RPFT uses the RFT [7,8] instead of the FT in the PFT computation. Figure 2 compares the performances achieved by the PFT and RPFT of the signal expressed as x(t) = e − jπ0.002t 2 + e jπ(0.002t 2 +0.3t) with impulse noise w(t) = 0.5[w 3 1 (t) + jw 3 2 (t)]. It is shown that two chirp rates are easily identified from the RPFT shown in Figure 2a rather than from the PFT in Figure 2b.

WINDOW LENGTH ESTIMATION
In the previous section, L(t) is estimated based on the idea of modeling each segment as an Mth-order polynomial phase signal. Therefore, the window length used in the MLPTFT p or RMLPTFT p is the same as the length of the segment. It is known that there is a tradeoff between the window length and the resolution of the MLPTFT p [9]. In general, approximation errors increase with the window length if the order of the MLPTFT p is lower than that of the phase of the signal segment. For a polynomial phase component whose order is not higher than that of MLPTFT p , on the other hand, the MLPTFT p gives a better resolution if a longer window (or segment) is used. For a good compromise, it is always desired that the length of the segment is adaptively matched to the characteristics of the signal components. In our analysis, the initial window length is selected to be small enough to provide acceptable accuracy of the approximation and the actual length of the window is increased according to the properties of consecutive signal segments.
Since we intend to increase the window length used if consecutive segments have the same polyphase model, we assume that two consecutive segments, the jth and ( j + 1)th segments, belong to the same polynomial phase model. If the jth segment has the phase 2π M m=0 k i,m t m , the phase of the ( j + 1)th segment should be 2π M m=0 k i,m (t + (Q − α)) m because the ( j + 1)th segment is delayed by a time interval of the segment overlap compared with that of the jth segment. The difference between the coefficients of the consecutive segments is calculated by where C s m = s!/(m!(m − s)!). For clarity of presentation, we define where b m is the constant coefficient associated with t m term on the right-hand side of (13). We represent the coefficients of the polynomial function estimated from the jth and ( j + 1)th segments with a j,m and a j+1,m , respectively, where 1 ≤ m ≤ M. The difference (a j+1,m − a j,m ) is compared with b m . If each |a j+1,m − a j,m − b m | is smaller than a predefined threshold T m , these two segments have the same polyphase model and the length of the window increases by Q − α.
The final window length is the total length of the consecutive segments that have the same polynomial function model. In general, T m is defined as the summation of two values. The first value is the deviation caused by the estimation of k j,m in the consecutive segments due to the noise influence. This value is decided by the statistical performance of the estimation method introduced in Section 3. The second one is the difference defined by the user. This factor controls the tradeoff between the resolution and distortion from the real time-frequency characteristics of the signal. A larger value is used if the resolution is the main consideration of the applications. In our simulations in the next section, we define the first value as the bias introduced by the gridding operation during the estimation of L(t) and the second value as zero.
From the previously described estimation methods of L(t) and window lengths, the main computational complexity is the computation of several PFTs so that L(t) leading to the maximum peak values can be selected. It can be easily seen that compared with the algorithm reported in [3], the computational complexity for L(t) estimation is significantly reduced. This is because, with the segmentation method shown in Figure 1, the number of segments for an N-point sequence is reduced to be N/(Q − α) in comparison with N segments needed in [3]. It is worth mentioning that the computational complexity can be further reduced if an efficient algorithm, such as a high-order ambiguity function [10] and its variations [11], is used in estimating the polynomial phase parameters instead of PFT. The estimation of window lengths requires overheads for computation of (12) and the costs of comparison with the given threshold is trivial. It should be noted that the presented computation process of MLPTFT p can be extended to deal with the mixture of Gaussian and impulse noise by using L-estimation-based transforms [12] instead of FT in the same way as that MLPTFT p is extended to RMLPTFT p .

EXPERIMENTAL RESULTS
Two types of signals, which contain a mono-and multicomponent, respectively, are used to test the performance of the proposed algorithms. For simplicity, the second-order MLPTFT p is used in all experiments dealing with the input sequence x(t) with N = 512. The L(t) is estimated from the positions (a i = {a i,1 , . . . , a i,M }, 1 ≤ i ≤ p and M = 2) of the peaks in the PFT or RPFT based on (7). The noise term w(t) is chosen from Gaussian and/or impulse noise.
The first type of signals contains a monocomponent defined as The estimation of instantaneous frequency of x 1 (t) is conducted with Gaussian noise w(t) of different variances. Monte Carlo simulations are performed to obtain the mean square error (MSE) for each estimator. The MSE is defined by is the true instantaneous frequency andf (t) is the estimation of f (t) according to the curve peak positions in the MLPTFT p of x 1 (t). The MSEs are compared for different overlap lengths α = 0, Q/2 , Q/3 , and Q − 1, as shown in Figure 3. A large range of window lengths (6 < Q < 80) have been tried. Only the results using Q = 23, Q = 33, and Q = 43 are shown because the MSEs achieved with these windows are mostly below 10 −2 . It is observed that for high SNRs, smaller window length generally gives smaller MSEs, for example, the curve for Q = 23 gives the best performance when SNR > 7 dB. With low SNRs, that is, SNR ≤ 7 dB, larger window length yields lower MSEs, as seen from Figure 3. This is because that MSEs are mainly influenced by the bias and the variance of the input sequence [3]. When SNR is high, the MSE is mainly affected by the bias which increases with the increase of window length. When SNR is low, the variance of the signal is the dominant factor affecting the MSE. The variances decrease with the increase of window length so that the MSE becomes relatively small. It is worth mentioning that, when SNR is extremely low, for example, below 0 dB, MSEs deteriorate significantly. This is because the windows used in our proposed MLPTFT p are generally with smaller length. The use of narrow windows leads to the increasing influence of noise especially for low SNR. The most important observation made in Figure 3 is that the MSE performances for different overlap lengths are very close to each other regardless of the window lengths, which leads to the conclusion that the decrease of overlap length between segments does not deteriorate the performance of MLPTFT p . In our investigation, several kinds of nonstationary signals, including the fourthorder polynomial phase signals, sinusoidal FM signals, and logarithm FM signals, have been used to reach similar conclusions. Therefore, the reduction of computational complexity due to the reduction of the number of segments to be processed does not obviously deteriorate the MSE performance. We now change w(t) embedded in x 1 (t) to be the impulse noise which belongs to the Middleton class A model [5,13] with the form w(t) = w 1 (t) 10δ(t) and τ i forms Poisson processes. Figure 4 presents the MLPTFT p and RMLPTFT p to clearly show that good performance is achieved by using the RFT in RMLPTFT p instead of the FT in MLPTFT p . We consider the signal contains multiple components, which is defined as x 2 (t) = exp − j 256 π cos √ 2πt 256 + exp − j0.002πt 2 +exp − j0.002πt 2 +0.06πt where w(t) is the Gaussian noise with SNR = 0 dB. This signal is used to test the performance of using adaptive window lengths. Figure 5 shows the MLPTFT p s of x 2 (t) that are computed with and without using adaptive window lengths.
In each segment with α = 0, three peaks are detected from the corresponding PFTs. By comparing the difference between the estimated parameters in the consecutive segments with b m in (13), the algorithm finds that two peaks in the consecutive segments belong to the same components and the window length is enlarged to N to obtain a high resolution. It is shown that the resolution in Figure 5b for the linear component of x 2 (t) is improved significantly and the two parallel chirp components are clearly distinguished in comparison with Figure 5a in which the window length is fixed. Figure 6 shows the comparison between the MLPTFT p and RMLPTFT p when w(t) = 0.9[w 3 1 (t) + jw 3 2 (t)]. The signal s(t) used is the same as in the previous experiment except that the third component is zero. The components smear significantly in Figure 6a because the FT used in the MLPTFT p is not able to substantially suppress the non-Gaussian noise. However, the RMLPTFT p yields much better performance and the two components can be clearly seen in Figure 6b.

CONCLUSION
This paper presents analysis algorithms effectively dealing with time-varying multicomponent signals. In particular, these algorithms allow the reduction of computational complexity by minimizing the length of overlap between consecutive signal segments. Experiments show that by using the proposed algorithms of parameters estimation and adaptive window length, the signals containing both single and multiple components with Gaussian and impulse noises can be more accurately represented in the time-frequency domain.