Sparse TimeFrequency decomposition for multiple signals with same frequencies
Abstract
In this paper, we consider multiple signals sharing same instantaneous frequencies. This kind of data is very common in scientific and engineering problems. To take advantage of this special structure, we modify our datadriven timefrequency analysis by updating the instantaneous frequencies simultaneously. Moreover, based on the simultaneously sparsity approximation and fast Fourier transform, some efficient algorithms is developed. Since the information of multiple signals is used, this method is very robust to the perturbation of noise. And it is applicable to the general nonperiodic signals even with missing samples or outliers. Several synthetic and real signals are used to test this method. The performances of this method are very promising.
1 Introduction
Data is an important bridge connecting the human being and the natural world. Many important information of the real world is encoded in the data. In many applications, frequencies of the signal are usually very useful to reveal the underlying physical mechanism. Therefore, in the past several decades, many researchers devoted to find efficient time frequency analysis methods and many powerful methods have been developed, including the windowed Fourier transform [17], the wavelet transform [5, 17], the WignerVille distribution [9], etc. Recently, an adaptive time frequency analysis method, the Empirical Mode Decomposition (EMD) method [15, 20] was developed which provides an efficient adaptive method to extract frequency information from complicate multiscale data. Inspired by EMD, recently several new time frequency analysis methods have been proposed see e.g. the synchrosqueezed wavelet transform [6], the Empirical wavelet transform [8], the variational mode decomposition [10].
In the last few years, inspired by the EMD method and compressive sensing [3, 2, 7], Hou and Shi proposed a datadriven timefrequency analysis method based on the sparsest timefrequency representation [12]. In this method, the signal is assumed to fit following model:
(1) 
where are smooth functions, . We assume that and are ”less oscillatory” than . After some scaling, we can always make the signal lies in the time interval . So, here we assume the time span of the interval is . We borrow the terminology in EMD method and also call as the Intrinsic Mode Functions (IMFs) [15]. This model is also known as Adaptive Harmonic model (AHM) and is widely used in the timefrequency analysis literatures [1, 6].
The AHM model seems to be natural and simple In this model, the main difficulty is how to compute the decomposition. In the AHM model, the number of degrees of freedom is much larger than the given signal. This makes the possible decomposition is not unique. One essential problem is how to set up a criterion to pick up the ”best” decomposition and this criterion should be practical in computation. Inspired by the compressive sensing, we proposed to decompose the signal by looking for the sparsest decomposition. And the sparsest decomposition is obtained by solving a nonlinear optimization problem formulated as following:
(2) 
where is the dictionary consist of all IMFs which makes the decomposition satisfies the AHM model (see [12] for precise definition of the dictionary).
Two kinds of algorithms were proposed to solve the optimization problem (2). The first one is based on matching pursuit [16, 12] and the other one is based basis pursuit [4, 13]. The convergence of the algorithm based on matching pursuit was also analyzed under the assumption of certain scale separation property [14].
In the previous works, we consider the case that only one signal is given. However, in many application, we can obtain many different signals and these signals have same frequency structure. For instance, in the monitoring of buildings, usually many sensor are put in different locations of the same building to measure the vibration. From these sensors, many measurements of the vibration are obtained. Since the signals from different sensors measure the same building, they should have same instantaneous frequencies which associate with the natrual frequencies of the building. If we analyze the signals from different sensors individually, the structure that these signals share the same instantaneous frequencies is wasted. Intuitively, by taking advantage of this structure, we may get more robust and efficient methods. In this paper, we will propose several methods by exploiting this special frequency structure.
In this paper, we consider multiple signals sharing same instantaneous frequencies. Naturally, we modify the adaptive harmonic model (1) to deal with this kind of signals.
(3) 
where is the number of signals, is the th signal. For each fixed, it is the AHM model (1), i.e. are smooth functions, . We assume that and are ”less oscillatory” than . The main feature of this model is that different signals have same phase functions . We call (3) Multiple Adaptive Harmonic model (MAHM).
Inspired by our previous work on single signal, we also look for the sparsest decomposition which satisfy the MAHM model (3) by solving following optimization problem:
(4)  
The dictionary we used in this paper is same as that in [12]. It can be written as
(5) 
is the collection of all the functions ”less oscillatory” than . In general, it is most effective to construct as an overcomplete Fourier basis given below:
(6) 
where , is the largest integer less than , and is a parameter to control the smoothness of . In our computations, we typically choose .
In the rest of the paper, we will introduce several algorithms to approximately solve above optimization problem (4). First, we give a generic algorithm based on matching pursuit and nonlinear least squares in Section 2. This algorithm can be accelerated by fast Fourier transform (FFT) if the signal is periodic. This will be reported in Section 3. For general nonperiodic signals, in Section 4, we also develop an efficient algorithm based on group sparsity and the algorithm in Section 2. This algorithm can be also generalized to deal with the signals with outliers or missing samples (Section 5). In Section 6, we present several numerical results including both synthetic and real signal to demonstrate the performance of our method. At the end, some conclusion remarks are made in Section 7.
2 Method based on matching pursuit
Following the same idea as that in [12], we proposed a method based on mathcing pursuit to get the sparsest decomposition, see Algorithm 1.
(7)  
To solve the nonlinear leastsquares problem (7) in Algorithm 1, we use the wellknown GaussNewton type iteration. This algorithm is also very similar as that in [12]. The only difference is the update of the phase function. In our model, different signals share the same phase function. So in the algorithm, we only update one common phase function by taking some average among different signals.
(8)  
(9) 
By integrating Algorithm 1 and 2 together, we can get a complete algorithm to compute the sparsest decomposition. The key step and also most expensive step is to solve the leastsquares problem (8). It does not need too much time to solve (8) once, however, we have to solve this leastsquares problem many times to get the final decomposition. This makes that the algorithm is not very practical in applications. In next two sections, we will propose some acceleration algorithms.
3 Periodic signals
If the signals are periodic, we can use a standard Fourier basis to construct the space instead of the overcomplete Fourier basis given in (6).
(10) 
where is a parameter to control the smoothness of functions in and is a positive integer.
In this case, we developed an efficient algorithm based fast Fourier transform to solve the leastsquares problem (8) in [12]. To make this paper selfcontain, we also include the algorithm here.
Suppose that the signal is measured over a uniform grid . Here, we assume that the sample points are fine enough such that can be interpolated to any grid with small error. Let be the normalized phase function and which is an integer.
The FFTbased algorithm to approximately solve (8) is given below:

Interpolate from in the physical space to a uniform mesh in the coordinate to get and compute the Fourier transform :
(11) where are uniformly distributed in the coordinate,i.e. . And the Fourier transform of is given as follows
(12) where .

Apply a cutoff function to the Fourier Transform of to compute and on the mesh of the coordinate, denoted by and :
(13) (14) is the inverse Fourier transform defined in the coordinate:
(15) 
Interpolate and from the uniform mesh in the coordinate back to the physical grid points :
(16) (17)
The lowpass filter in the second step is determined by the choice of .
In this paper, we choose the following lowpass filter to define :
(18) 
By incorporating the FFTbased solver in Algorithm 2, we get an FFTbased iterative algorithm which is summarized in Algorithm 3.
(19) 
Remark 3.1.
We remark that the projection operator is in fact a lowpass filter in the space. For nonperiodic data, we apply a mirror extension to before we apply the lowpass filter.
GaussNewton type iteration is sensitive to the initial guess. In order to abate the dependence on the initial guess, we gradually increase the value of to improve the approximation to the phase function so that it converges to the correct value. The detailed explanation can be found in [12].
4 nonperiodic signals
In most of the applications, the signals are not periodic. To apply the FFTbased algorithm in previous section, we have to do periodic extension for general nonperiodic signals. In this paper, the signals are assumed to satisfy the MAHM model (3). One consequence of the MAHM model is that the phase function can be used as a coordinate and the signals are sparse over the Fourier basis in coordinate. Then, one natural way to do periodic extension is to look for the sparsest representation of the signal over the overcomplete Fourier basis defined in coordinate.
More precisely, the extension is obtained by solving an minimization problem.
(20) 
is the sample of the signal at , is the number of sample points. may not be uniformly distributed, however, we assume that the sample points are fine enough such that can be interpolated over any other grids without loss of accuracy. is a matrix consisting of basis functions, is the number of Fourier modes. . is a positive integer determined by the length of the period of the extended signal which will be discussed later.
To get the periodic extension of signals, one way is to solve (20) times independently. However, these signals are not independent. They share the same phase functions. To take advantage of this structure, inspired by the methods for simultaneous sparsity approximation [18, 19], we propose to solve following optimization problem to get the periodic extension of signals simultaneously.
(21) 
where is a matrix, , is a matrix and
(22) 
The problem left is how to solve the optimization problem (21) efficiently. Since the sample points are fine enough, we can assume that the sample points are uniformly distributed in coordinate, i.e., and . Otherwise, we could use some interpolation method (for instance cubic spline interpolation) to get the signals over the uniformly distributed sample points. Now, we define a matrix, , consisting of complete Fourier basis,
(23) 
represents the element of at th row and th column. Using the property of the Fourier basis, we know that is an orthonormal matrix which will give us lots of convenience in deriving the fast algorithm.
For the convenience of notation, denote be the set of index of ,
(24) 
Then, for any , is a sample point. We also define an extension of by zero padding,
(25) 
where .
First, we remove the constraint in (21) by introducing a penalty term,
(26) 
Here is a parameter of the penalty. We will let goes to infinity later. Then the solution of (26) will converge to the solution of (21).
Let , and . Using the Augmented Lagrange Multiplier method (ALM) to solve the unconstrained optimization problem (26), we get following algorithm. Let , and repeat following two steps until converge
Here is a parameter. In above ALM iteration, the main computational load is to solve the optimization problem in the first step. Notice that the objective functional depends on two terms, and . It is natural to minimize the functional alternately by optimizing one of these two terms when the other one is fixed like that in the split Bregman iteration [11]. Using this idea, we get following iterative algorithm to solve (26).
Let , and repeat
The good news is that the optimization problems in the first and second step can be solved explicitly and we can use fast Fourier transform to accelerate the computation. First, let us see the first optimization problem. Using the fact that is an orthonormal matrix, we have
(27) 
where is the conjugate transpose of . Here is a shrinkage operator which is defined as following. For any
(28) 
If is a matrix and is its th row, then is a matrix of the same size as and the th row is defined to be . Notice that is the matrix consisting of Fourier basis in coordinate, so the matrixmatrix multiplication can be evaluated by applying discrete Fourier transform of each column of the second matrix which can be accelerated by FFT.
In the second optimization problem, recall that is a penalty parameter in (26) which is the larger the better. Now, let goes to infinity, then we can solve
(29) 
Summarizing above derivation, we get a fast algorithm to solve (21), see Algorithm 4.
By combining Algorithm 4 and the algorithm in previous section, we can get a complete algorithm to deal with nonperiodic data. Before giving the complete algorithm, there is still one issue we need to address. In the derivation of Algorithm 4, we assume that the sample points are uniformly distributed in coordinate with given . For general signal which may not be uniformly sampled in coordinate, we need to do interpolation before applying Algorithm 4.
In order to do the interpolation, first, we need to give the interpolation points which are uniformly distributed in coordinate. In this paper, the interpolation points are chosen to be , . We set and . This choice is corresponding to the 2fold overcomplete Fourier basis used in (6). where is the number of sample points of the original signal. The original signals are interpolated over the interpolation points by cubic spline interpolation.
Now, we have an algorithm for nonperiodic signal summarized in Algorithm 5.
(30) 
5 Signals with outliers or missing samples
To deal with the signal with outliers, we need to enlarge the dictionary to include all the impulses
In this case, the optimization problem is formulated in the following way,
(31) 
where and are same as those in (21)
Following the similar derivation in Section 4, we obtain Algorithm 6 to solve above optimization problem (31).
By using the same interpolation procedure in Section 4, we can integrate Algorithm 6 with the algorithm in Section 3 to get the method to deal with signals with outliers.
For the signals with missing samples, first, we assign the value of the signal as the average of the signal at the locations where the sample was lost and then treat the missing sample as the outliers.
6 Numerical results
In this section, we use two numerical examples, one is synthetic and one is real data, to demonstrate the performance of the method proposed in this paper. The first example is a simple synthetic signal which is used to test the robustness of the method to the perturbation of white noise.
Example 1: Synthetic signal with white noise
In this example, the signals are generated as following:
(32) 
where are independent white noise with standard deviation . The number of sample points is 512 and the sample points are uniformly distributed over .
The original signals are plotted in the left figure of Fig. 1. As we can see that the noise is so large such that we can not see any pattern of the original clean signal. The recovered instantaneous frequency is shown in the right figure of Fig. 1. If we recover the frequency from each signal separately, since the noise is too large, the frequencies are totally wrong. However, if we use the structure that 10 signals have same instantaneous frequency, the frequency recovered is much better.
Example 2: Cable tension estimate
This example is a real problem about the tension estimate of the cables in bridge. Before demonstrating the numerical results, we give some backgrounds first.
For large span bridges, such as cablestayed bridge and suspension bridge, the cables are a crucial element for overall structural safety. Due to the moving vehicles and other environmental effects, the cable tension forces vary over time. This variation in cable tension forces may cause fatigue damage. Therefore, estimation of the timevarying cable tension forces is important for the maintenance and safety assessment of cablebased bridges.
One often used method to estimate the cable tension force is based on the instantaneous frequency estimate of the cable vibration signal. According to the flat taut string theory that neglects both sagextensibility and bending stiffness, cable tension force, , can be calculated by
(33) 
where is the timevarying th natural frequency in and , are mass density and length of cable.
Also from the flat taut string theory, an important and useful feature of the vibriations of the cable is that the natural frequencies of the higher modes are integer multiples of the fundamental frequency, that is . This feature means that we can combine the information of different modes together to recover the instantaneous frequency. Then the method developed in this paper for multiple signals applies after some minor modifications.
The original experimental signal is given in Fig. 2. Obviously, the signal has some outliers. In the original signal, about 10% of the samples were lost which is demonstrated very clearly in the zoom in picture in Fig. 2. The tension force estimated by our method is given in Fig. 3 and 4. The tension forces obtained from different modes are shown in Fig. 3. If only one mode is used, the estimation of the force is not very accurate. There are many oscillations although the overall picture match with the measured force. If we use 15 modes together to estimate the tension force, the result is much better, Fig. 4.
7 Conclusion
In this paper, we consider multiple signals and these signals share the same instantaneous frequencies. This kind of signals emerge in many scientific and engineering problems. By exploiting the structure of the common frequencies, we develop some algorithms to find the sparsest timefrequency decomposition. These algorithms can be accelerated by FFT, so they are very efficient. The algorithms are also very robust to noise, since the information of multiple signals are used simultaneously.
Acknowledgments. This work was supported by NSF FRG Grant DMS1159138, DMS1318377, an AFOSR MURI Grant FA95500910613 and a DOE grant DEFG0206ER25727. The research of Dr. Z. Shi was supported by a NSFC Grant 11201257.
Footnotes
 The authors would like to thank Prof. Yuequan Bao of the School of Civil Engineering of Harbin Institute of Technology for the experimental data and for permission to use it in this paper.
References
 C. K. Chui and H. N. Mhaskar, Signal decomposition and analysis via extraction of frequencies, Appl. Comput. Harmon. Anal., 2015, http://dx.doi.org/10.1016/j.acha.2015.01.003
 E. Cands and T. Tao, Near optimal signal recovery from random projections: Universal encoding strategies?, IEEE Trans. on Information Theory, 52(12), pp. 54065425, 2006.
 E. Candes, J. Romberg, and T. Tao, Robust uncertainty principles: Exact signal recovery from highly incomplete frequency information, IEEE Trans. Inform. Theory, 52, pp. 489509, 2006.
 S. Chen, D. Donoho and M. Saunders, Atomic decomposition by basis pursuit, SIAM J. Sci. Comput., 20, pp. 3361, 1998.
 I. Daubechies, Ten Lectures on Wavelets, CBMSNSF Regional Conference Series on Applied Mathematics, Vol. 61, SIAM Publications, 1992.
 I. Daubechies, J. Lu and H. Wu, Synchrosqueezed wavelet transforms: an empirical mode decompositionlike tool, Appl. Comp. Harmonic Anal., 30 (2011), pp. 243261.
 D. L. Donoho, Compressed sensing, IEEE Trans. Inform. Theory, 52, pp. 12891306, 2006.
 K. Dragomiretskiy and D. Zosso, Varitational Mode Decomposition, IEEE Trans. Signal Process, 62, pp. 531544, 2014.
 P. Flandrin, TimeFrequency/TimeScale Analysis, Academic Press, San Diego, CA, 1999.
 J. Gilles, Empirical Wavelet Transform, IEEE Trans. Signal Process, 61, pp. 39994010, 2013.
 Tom Goldstein and Stanley Osher, The Split Bregman Method for Regularized Problems, SIAM J. Imaging Sci., 2, pp. 323343, 2009.
 T. Y. Hou and Z. Shi, DataDriven TimeFrequency Analysis, Appl. Comput. Harmon. Anal., 35, pp. 284308, 2013.
 T. Y. Hou and Z. Shi, Sparse TimeFrequency decomposition by dictionary learning, arXiv cs.IT 1311.1163, 2014.
 T. Y. Hou, Z. Shi and P. Tavallali, Convergence of a DataDriven TimeFrequency Analysis Method, Appl. Comput. Harmon. Anal., 37, pp. 235270, 2014.
 N. E. Huang et al., The empirical mode decomposition and the Hilbert spectrum for nonlinear and nonstationary time series analysis, Proc. R. Soc. Lond. A, 454 (1998), pp. 903995.
 S. Mallat and Z. Zhang, Matching pursuit with timefrequency dictionaries, IEEE Trans. Signal Process, 41, pp. 33973415, 1993.
 S. Mallat, A wavelet tour of signal processing: the Sparse way, Academic Press, 2009.
 J. A. Tropp, A. C. Gilbert, and M. J. Strauss, Algorithms for simultaneous sparse approximation. Part I: Greedy pursuit, Signal Processing, special issue ”Sparse approximations in signal and image processing”, 86, pp. 572588, 2006.
 J. A. Tropp, Algorithms for simultaneous sparse approximation. Part II: Convex relaxation, Signal Processing, special issue ”Sparse approximations in signal and image processing”, 86, pp. 589602, 2006.
 Z. Wu and N. E. Huang, Ensemble Empirical Mode Decomposition: a noiseassisted data analysis method, Advances in Adaptive Data Analysis, 1, pp. 141, 2009.