Sparse Time-Frequency decomposition for multiple signals with same frequencies

Sparse Time-Frequency 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 data-driven time-frequency 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 Wigner-Ville 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 data-driven time-frequency analysis method based on the sparsest time-frequency 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 time-frequency 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.

0:  Signals , .
0:  Phase functions and the corresponding envelopes: .
1:  Set and .
2:  while   do
3:     Solve the following nonlinear least-square problem:
(7)
4:     Update the residual
5:     Set .
6:  end while
Algorithm 1 Method based on matching pursuit

To solve the nonlinear least-squares problem (7) in Algorithm 1, we use the well-known Gauss-Newton 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.

0:  Initial guess of the phase function .
0:  Phase functions and the corresponding envelopes: .
1:  while   do
2:     Solve the following least-square problem for each signal , :
(8)
3:     Calculate averaged updation of frequency:
(9)
where and denote the derivative of with respect to .
4:     Update
where is chosen to make sure that is monotonically increasing:
5:  end while
Algorithm 2 Gauss-Newton iteration to solve the nonlinear least-squres problem

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 least-squares problem (8). It does not need too much time to solve (8) once, however, we have to solve this least-squares 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 over-complete 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 least-squares problem (8) in [12]. To make this paper self-contain, 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 FFT-based 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 low-pass filter in the second step is determined by the choice of .

In this paper, we choose the following low-pass filter to define :

(18)

By incorporating the FFT-based solver in Algorithm 2, we get an FFT-based iterative algorithm which is summarized in Algorithm 3.

0:  Initial guess of the phase function .
0:  Phase functions and the corresponding envelopes: .
1:  while  do
2:     while   do
3:        Interpolate to a uniform mesh in the -coordinate to get and compute the Fourier transform .
4:        Apply a cutoff function to the Fourier Transform of to compute and on the mesh of the -coordinate, denoted by and .
5:        Interpolate and back to the uniform mesh of .
6:        Calculate averaged updation of frequency:
(19)
where and denote the derivative of with respect to .
7:        Update
where is chosen to make sure that is monotonically increasing:
and is the projection operator to the space .
8:     end while
9:     
10:  end while
Algorithm 3 (FFT-based algorithm to solve the nonlinear least-squres problem)
Remark 3.1.

We remark that the projection operator is in fact a low-pass filter in the -space. For non-periodic data, we apply a mirror extension to before we apply the low-pass filter.

Gauss-Newton 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 FFT-based 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 over-complete 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 matrix-matrix 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.

0:  .
0:  Fourier coefficients on over-complete Fourier basis .
1:  while   do
2:     , where means apply Fourier transform on each column.
3:     
4:     
5:  end while
6:  .
Algorithm 4 (Fourier extension by group sparsity)

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 2-fold over-complete 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.

0:  Initial guess of the phase function .
0:  Phase functions and the corresponding envelopes: .
1:  while  do
2:     while   do
3:        Interpolate to a uniform mesh in the -coordinate to get .
4:        Using Algorithm 4 to get the Fourier coefficients of on the over-complete Fourier basis.
5:        Apply a cutoff function to the Fourier Transform of to compute and on the mesh of the -coordinate, denoted by and .
6:        Interpolate and back to the uniform mesh of .
7:        Calculate averaged update of frequency:
(30)
where and denote the derivative of with respect to .
8:        Update
where is chosen to make sure that is monotonically increasing:
and is the projection operator to the space .
9:     end while
10:     
11:  end while
Algorithm 5 (Algorithm for nonperiodic signal)

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).

0:  .
0:  Fourier coefficients on over-complete Fourier basis , estimate of outliers .
1:  while   do
2:     , where means apply Fourier transform on each column.
3:     
4:     , where is a shrinkage operator. For any ,
For any matrix , .
5:     
6:  end while
7:  , , is defined in (24).
Algorithm 6 (Algorithm for signals with outliers)

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.

Figure 1: Left: 10 measurements generated by . Bold black: signal without noise; thin lines: 10 measurements with larger noise; Right: instantaneous frequency recoverd from signal with 10 different measurements with larger noise. Bold blue: exact frequency; bold black: recovered frequency from 10 signals; thin lines: frequency obtained from 10 signals separately.

Example 2: Cable tension estimate1

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 cable-stayed 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 time-varying cable tension forces is important for the maintenance and safety assessment of cable-based 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 sag-extensibility and bending stiffness, cable tension force, , can be calculated by

(33)

where is the time-varying 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 1-5 modes together to estimate the tension force, the result is much better, Fig. 4.

Figure 2: Left: vibration signal of the cable; Right: zoom in of the signal to demonstrate the missing samples (flat segments).
Figure 3: Cable force estimated by single natural mode 1-5 from up to bottom.
Figure 4: Cable force estimated by using 1-5 natural modes together.

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 time-frequency 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 DMS-1159138, DMS-1318377, an AFOSR MURI Grant FA9550-09-1-0613 and a DOE grant DE-FG02-06ER25727. The research of Dr. Z. Shi was supported by a NSFC Grant 11201257.

Footnotes

  1. 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

  1. 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
  2. E. Cands and T. Tao, Near optimal signal recovery from random projections: Universal encoding strategies?, IEEE Trans. on Information Theory, 52(12), pp. 5406-5425, 2006.
  3. E. Candes, J. Romberg, and T. Tao, Robust uncertainty principles: Exact signal recovery from highly incomplete frequency information, IEEE Trans. Inform. Theory, 52, pp. 489-509, 2006.
  4. S. Chen, D. Donoho and M. Saunders, Atomic decomposition by basis pursuit, SIAM J. Sci. Comput., 20, pp. 33-61, 1998.
  5. I. Daubechies, Ten Lectures on Wavelets, CBMS-NSF Regional Conference Series on Applied Mathematics, Vol. 61, SIAM Publications, 1992.
  6. I. Daubechies, J. Lu and H. Wu, Synchrosqueezed wavelet transforms: an empirical mode decomposition-like tool, Appl. Comp. Harmonic Anal., 30 (2011), pp. 243-261.
  7. D. L. Donoho, Compressed sensing, IEEE Trans. Inform. Theory, 52, pp. 1289-1306, 2006.
  8. K. Dragomiretskiy and D. Zosso, Varitational Mode Decomposition, IEEE Trans. Signal Process, 62, pp. 531-544, 2014.
  9. P. Flandrin, Time-Frequency/Time-Scale Analysis, Academic Press, San Diego, CA, 1999.
  10. J. Gilles, Empirical Wavelet Transform, IEEE Trans. Signal Process, 61, pp. 3999-4010, 2013.
  11. Tom Goldstein and Stanley Osher, The Split Bregman Method for -Regularized Problems, SIAM J. Imaging Sci., 2, pp. 323-343, 2009.
  12. T. Y. Hou and Z. Shi, Data-Driven Time-Frequency Analysis, Appl. Comput. Harmon. Anal., 35, pp. 284-308, 2013.
  13. T. Y. Hou and Z. Shi, Sparse Time-Frequency decomposition by dictionary learning, arXiv cs.IT 1311.1163, 2014.
  14. T. Y. Hou, Z. Shi and P. Tavallali, Convergence of a Data-Driven Time-Frequency Analysis Method, Appl. Comput. Harmon. Anal., 37, pp. 235-270, 2014.
  15. N. E. Huang et al., The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis, Proc. R. Soc. Lond. A, 454 (1998), pp. 903-995.
  16. S. Mallat and Z. Zhang, Matching pursuit with time-frequency dictionaries, IEEE Trans. Signal Process, 41, pp. 3397-3415, 1993.
  17. S. Mallat, A wavelet tour of signal processing: the Sparse way, Academic Press, 2009.
  18. 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. 572-588, 2006.
  19. 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. 589-602, 2006.
  20. Z. Wu and N. E. Huang, Ensemble Empirical Mode Decomposition: a noise-assisted data analysis method, Advances in Adaptive Data Analysis, 1, pp. 1-41, 2009.
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
""
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
   
Add comment
Cancel
Loading ...
130369
This is a comment super asjknd jkasnjk adsnkj
Upvote
Downvote
""
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters
Submit
Cancel

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test
Test description