The Fourier Decomposition Method for nonlinear and nonstationary time series analysis

# The Fourier Decomposition Method for nonlinear and nonstationary time series analysis

## Abstract

Since many decades, there is a general perception in literature that the Fourier methods are not suitable for the analysis of nonlinear and nonstationary data. In this paper, we propose a Fourier Decomposition Method (FDM) and demonstrate its efficacy for the analysis of nonlinear (i.e. data generated by nonlinear systems) and nonstationary time series. The proposed FDM decomposes any data into a small number of ‘Fourier intrinsic band functions’ (FIBFs). The FDM presents a generalized Fourier expansion with variable amplitudes and frequencies of a time series by the Fourier method itself. We propose an idea of zero-phase filter bank based multivariate FDM (MFDM) algorithm, for the analysis of multivariate nonlinear and nonstationary time series, from the FDM. We also present an algorithm to obtain cutoff frequencies for MFDM. The MFDM algorithm is generating finite number of band limited multivariate FIBFs (MFIBFs). The MFDM preserves some intrinsic physical properties of the multivariate data, such as scale alignment, trend and instantaneous frequency. The proposed methods produce the results in a time-frequency-energy distribution that reveal the intrinsic structures of a data. Simulations have been carried out and comparison is made with the Empirical Mode Decomposition (EMD) methods in the analysis of various simulated as well as real life time series, and results show that the proposed methods are powerful tools for analyzing and obtaining the time-frequency-energy representation of any data.

{IEEEkeywords}

The Fourier decomposition method (FDM); Fourier intrinsic band functions (FIBFs) and analytic FIBFs (AFIBFs); zero-phase filter bank based multivariate FDM (MFDM); Empirical Mode Decomposition (EMD). \IEEEpeerreviewmaketitle

## 1 Introduction

\IEEEPARstart

The time-frequency representation (TFR) is a well established powerful tool for the analysis of time series signals. There exist many types of time-frequency (TF) analysis methods, e.g. linear (the short-time Fourier transform), quadratic (the Wigner-Ville distribution) and Wavelet transforms. The TFR is achieved by formulation often referred as time-frequency distribution (TFD) and provides insight into the complex structure of a signal consisting of several components. These approaches have many useful applications, however the analysis of nonstationary signals are not well presented by these methods.

Recently developed Empirical Mode Decomposition (EMD) [1] has provided a general method for examining the TFD, and has been applied to all kinds of data. The EMD is an adaptive signal analysis algorithm for the analysis of nonstationary and nonlinear signals (i.e. signals generated from nonlinear systems). The EMD has become an established method for the signal analysis in various applications, e.g. medical studies [2, 3, 4, 5], meteorology [1], geophysical studies [6] and image analysis[7]. The EMD decomposes any given data into a finite number of narrow band intrinsic mode functions (IMFs) which are derived directly from the data, whereas other signal decomposition techniques (like the Fourier, Wavelets, etc.) incorporate predefined fixed basis for signal modeling and analysis. The Ensemble EMD (EEMD) is a noise-assisted data analysis method developed in [8] to overcome the timescale separation problem of EMD. The Multivariate EMD (MEMD) developed in [9] is a generalization of the EMD for multichannel data analysis. The Compact EMD (CEMD) algorithm is proposed in [10] to reduce mode mixing, end effect, and detrend uncertainty present in EMD and to reduce computation complexity of EEMD as well. The IMFs, generated by EMD, are dependent on distribution of local extrema of signal and the type of spline used for upper and lower envelope interpolation. The traditional EMD uses cubic spline for upper and lower envelope interpolation. The EMD algorithm, proposed in [11] to reduce mode mixing and detrend uncertainty, uses nonpolynomial cubic spline interpolation to obtain upper and lower envelopes, and have shown in [12] that it improves orthogonality among IMFs.

The energy preserving property is important for any kind of transformation, and it is obtained by the orthogonal decomposition of signal in various transforms like the Fourier, Wavelet, Fourier-Bessel, etc. The energy preserving property is especially important for the accurate and faithful analysis of three dimensional time-frequency-energy distribution of a signal. The EMD algorithms, proposed in [13], ensure orthogonality or energy preserving property or both in decomposition of signal into IMFs and refereed to as energy preserving EMD (EPEMD).

In spite of considerable success of EMD, all of the EMD algorithms are based on empirical, heuristic and ad hoc procedure that make them hard to analyze mathematically, and EMD may suffer from mode mixing, detrend uncertainty, aliasing and end effect artefacts [14]. There is a lack of mathematical understanding of the EMD algorithms, e.g. IMFs dependence on the number of sifting and the stopping criteria, convergence property and stability to noise perturbation. In spite of all these limitations, EMD is the widely used nonstationary data analysis method. Therefore, in this paper, EMD is used as a reference to establish the validity, reliability and calibration of the proposed method.

Since many decades, there is an understanding in the literature (e.g. [1, 14, 24]) that the Fourier methods are, directly, not suitable for nonlinear and nonstationary data analysis, and various reasons (e.g. linearity, periodicity or stationarity) are provided to support it. The Fourier transform is valid under very general Dirichlet conditions (i.e. signal is absolutely integrable with finite number of maxima and minima, and finite number of finite discontinuities in any finite interval) thus include nonlinear and nonstationary signals as well. Therefore, in this study, we explore and provide algorithms to analyze nonlinear and nonstationary data by the Fourier method termed as the Fourier Decomposition Method (FDM), which generates small number of Fourier intrinsic band functions (FIBFs). It is already well established that the Fourier method is best tool for spectrum analysis and, in this study, we show that the Fourier method is also a best tool for time-frequency analysis and processing of any signal. The power of the Fourier transform can also be realized from the fact that the analytic representation and, hence, the Hilbert transform of a signal are, inherently, present in the Fourier transform.

In this paper, we also propose a method, which captures the features of the MEMD, using a zero-phase filter bank (ZPFB) approach to construct the multivariate FIBFs and residue components. This multivariate FDM (MFDM) algorithm generates matched multivariate FIBFs and residue through zero-phase filtering. Thus, we propose an adaptive, data-driven, ZPFB based time-frequency analysis method.

For the adaptive data analysis approach, the most difficult challenge has been to establish a general adaptive decomposition method of analysis without a priori basis. In this study, we propose the FDM and MFDM general adaptive data analysis methods that are inspired by the EMD methods and their filter bank properties [15, 16]. This paper is organized as follows: In section II the EMD algorithm is briefly presented. We propose the Fourier decomposition method (FDM) in section III. We propose the ZPFB based multivariate FDM algorithm in section IV. Simulation results are presented in section V. Finally conclusions are presented in section VI.

## 2 Brief overview of the EMD algorithm

There are various tools for nonstationary data processing such as the spectrogram; the wavelet analysis; the Wigner-Ville distribution; evolutionary spectrum [17]; the empirical orthogonal function expansion (EOF) (or principal component analysis or singular value decomposition); Synchrosqueezed wavelet transforms [18]; the EMD, etc.

The EMD is a well established signal decomposition method that decomposes nonlinear and nonstationary data into a set of finite band-limited IMFs and residue through the sifting process. The decomposed signal is expressed as the sum of IMF components plus the final residue as

 x(t)=ℓ∑i=1yi(t)+rℓ(t)=ℓ+1∑i=1yi(t), (1)

where is the IMF and is final residue. The IMFs admit amplitude-frequency modulated (AM-FM) representation [18] (i.e. , with ) and well-behaved Hilbert transforms [1]. For any IMF , its Hilbert transform is defined as convolution of and , i.e. and the Hilbert transform emphasizes the local properties of . An analytic signal can be represented by , where and are instantaneous amplitude and phase of . The instantaneous frequency (IF) of is defined as: . The physical meaning of IF constrains that must be a mono-component function of time. The Bedrosian and Nuttall theorems [19, 20] further impose non-overlapping spectra constraints on the pair [] of a signal .

All IMFs must satisfy two basic conditions: (1) In the complete range of time series, the number of extrema (i.e. maxima and minima) and the number of zero crossings are equal or differ at most by one. (2) At any point of time in the complete range of time series, the average of the values of upper and lower envelopes, obtained by the interpolation of local maxima and the local minima, is zero. The first condition ensure that IMFs are narrow band signals and the second condition is necessary to ensure that the IF does not have redundant fluctuations because of asymmetric waveforms [1].

## 3 The Fourier Decomposition Method

We propose a class of functions, termed as the Fourier intrinsic band functions (FIBFs), belonging to , here with the following formal definition.

###### Definition 1

The Fourier intrinsic band functions (FIBFs), , are functions that satisfy the following conditions:
(1) The FIBFs are zero mean functions, i.e. .
(2) The FIBFs are orthogonal functions, i.e. , for .
(3) The FIBFs provide analytic FIBFs (AFIBFs) with instantaneous frequency (IF) and amplitude always greater than zero, i.e. , with , .

Thus, the AFIBFs are monocomponent signals and, physically, the IF has meaning only for monocomponent signals, i.e., signal has only one frequency or a narrow range of frequencies varying as a function of time [24]. Thus, the FIBF is sum of zero mean sinusoidal functions of consecutive frequency band.

The main objective of this study is to obtain unique representation of multicomponent signal as a sum of constant and monocomponent signals, i.e. signals which can be represented by the following model [24]:

 x(t)=M∑i=1yi(t)+n(t), (2)

where is a noise representing any residue (constant or trend) components and the are single component nonstationary signals, which in our proposed framework would be the FIBFs, defined above.

The necessary conditions [1], for a basis to represent a nonlinear and nonstationary time series, are completeness, orthogonality, locality, and adaptiveness. The FIBFs, intrinsically, follow all the necessary conditions by virtue of the decomposition.

The available data are usually of finite duration, nonstationary and generated from the systems that are generally nonlinear. Let be a time limited real valued signal which follows the Dirichlet conditions. We construct the periodic signal as such that , where , for and zero otherwise. The Fourier series expansion of is given by

 x\tiny T0(t)=a0+∞∑k=1[akcos(kω0t)+bksin(kω0t)], (3)

where frequency (rad/s) , , and . We write (3) as

 x\tiny T0(t)=a0+12∞∑k=1[ckexp(jkω0t)+c∗kexp(−jkω0t)], (4)

where and . From (4), it is clear that

 x\tiny T0(t)=a0+Re{z\tiny T0(t)}, (5)

where analytic function

 z\tiny T0(t)≜∞∑k=1ckexp(jkω0t) (6)

is complex conjugate of and denotes the real part of . We write as

 z\tiny T0(t)=M∑i=1ai(t)exp(jϕi(t)), (7)

where, in forward search (i.e. low to high frequency scan) of AFIBFs, , , , , i.e.

 ai(t)exp(jϕi(t))=Ni∑k=Ni−1+1ckexp(jkω0t), (8)

with and . The FIBFs are the real part of AFIBFs presented in Eq. (8). In order to obtain minimum number of AFIBFs in low to high frequency scan (LTH-FS), for each , start with and append more term till we reach the maximum value of such that and

 ai(t),ωi(t)=dϕi(t)dt≥0,∀t. (9)

It is easy to observe that such a decomposition is always possible.

Similarly, in reverse search (i.e. form high to low frequency scan (HTL-FS)) of AFIBFs, we obtain , , , , and the lower and upper limits of sum in Eq. (8) would change to , respectively, with , . Here, we start with , decrease and select minimum value of such that and Eq. (9) is satisfied for .

Observe that  (7) has precisely the form that in the literature [1] is termed as a generalized Fourier expansion. Moreover, it is worth noting that this representation is complete, orthogonal, local, adaptive and purely Fourier based. Thus, we have obtained a generalized Fourier expansion of a time series in Eq. (7) by the Fourier method itself. The variable amplitude and the IF have improved the efficiency of the expansion by expanding the signal into finite number of analytic FIBFs, in (7), and enabled the expansion to accommodate nonstationary data. Thus, we have obtained a variable amplitude and frequency representation, whereas, the classical Fourier expansion provides the constant amplitude and fixed-frequency representation.

For each FIBFs, the amplitude and IF are functions of time, therefore, we define the three dimensional time-frequency distribution of amplitude as the Fourier-Hilbert spectrum (FHS) . The marginal Hilbert spectrum (MHS), derived from Hilbert spectrum, is defined in [1]. Similarly, here we define the marginal Fourier-Hilbert spectrum (MFHS) from the FHS as follows:

 h(f)=∫T00H(f,t)dt. (10)

The marginal Fourier-Hilbert spectrum offers a measure of total amplitude (or energy) contribution from each value of frequency in a probabilistic sense. The frequency in either or has a different meaning from the Fourier spectral analysis [1]. The presence of energy at each frequency in MFHS means that, in the total duration of the signal, there is a higher likelihood for such a wave (FIBF) to have appeared locally. The frequency in the MFHS indicates only the likelihood that an oscillation with such a frequency exists. The exact occurrence time of that oscillation is given in the full Fourier-Hilbert spectrum. We can also define the instantaneous energy density, which can be used to measure the fluctuation of energy with time, as

 E(t)=∫fM0H2(f,t)df, (11)

where is a maximum frequency of signal. From Eq. (3), we obtain the energy of signal (or power of signal ) by the Parseval’s theorem as and from Eq. (6) energy of the analytic signal (or power of signal ) as , therefore, relation between and is given by

 Ex=a20+Ez2. (12)

Hence, energy of zero mean signal is half of the energy of its analytic signal.

Since in practice the continuous time signals are, generally, discretized for further processing by a computing device, so we present the FDM for discrete signal. Let, , be a discrete signal of length . Using the discrete Fourier transform (DFT), we can write as

 x[n]=N−1∑k=0X[k]exp(j2πknN), (13)

where is the DFT of signal . Let be an even number (we can proceed in the similar fashion when is an odd number), then and are real numbers; and we can write as

 x[n]=X[0]+N2−1∑k=1X[k]exp(j2πknN)+X[N2]exp(jπn)+N−1∑k=N2+1X[k]exp(j2πknN). (14)

Since is real, therefore, is complex conjugate of and we can write (14) as

 x[n]=X[0]+2Re{z1[n]}+X[N2](−1)n, (15)

where denote the real part of . Now, we write analytic signal as

 N2−1∑k=1X[k]exp(j2πknN)=M∑i=1ai[n]exp(jϕi[n]), (16)

where, in forward search (i.e. low to high frequency scan) of AFIBFs, we obtain , , , , i.e.

 ai[n]exp(jϕi[n]=Ni∑k=Ni−1+1X[k]exp(j2πknN), (17)

with and . In order to obtain minimum number of FIBFs in LTH-FS, for each , we scan from () to , obtain maximum value of such that and phase is a monotonically increasing function, i.e. or

 ωi[n]=(ϕi[n+1]−ϕi[n−1]2)≥0,∀n (18)

and for and . Observe that such a decomposition always exists.

Similarly, in HTL-FS for FIBFs, we obtain , , , , and the lower and upper limits of sum in Eq. (17) will change to , respectively, with , . Here, for each , we scan from () to , obtain the minimum value of such that and phase is a monotonically increasing function.

Thus, FDM provides two views, low to high frequency and high to low frequency view, of the signal and generate two set of time-frequency-energy distribution. Depending on the signal, both view may be same or sometimes they reveal two different types of features of the signal. The FDM is summarized in Algorithms A and B. The FDM with Fourier transform (FT) and discrete time Fourier transform (DTFT) is summarized in appendix.

Algorithm A: The FDM algorithm (LTH-FS) to obtain AFIBFs, for with and
Obtain .
Set , obtain maximum value of such that and phase of is a monotonically increasing function, i.e. , .
Algorithm B: The FDM algorithm (HTL-FS) to obtain AFIBFs, for with and
Obtain .
Set , obtain minimum value of such that and phase of is a monotonically increasing function, i.e. , .

## 4 Multivariate Fourier Decomposition Method

From (8) and (17), we observe that the operation that generates the FIBFs is nothing but the Fourier based zero-phase filtering (ZPF). This is another motivation, in addition to the FB properties of IMFs, to use the Fourier or other methods of zero-phase filtering to decompose any data into a set of FIBFs. The ZPF of a real valued signal by zero-phase filter () can be obtained by two methods: (1) Convolution method, i.e. , where and , where is a real sequence. (2) The Fourier method, i.e. set at desired frequency band and otherwise, obtain output by the inverse DFT, i.e. , where .

We use ZPF that does not shift the essential features of the signal, and propose a multivariate FDM (MFDM) algorithm to generate multivariate FIBFs (MFIBFs) and residue as follows: Apply zero-phase high pass filtering (ZP-HPF) with cutoff frequency to each of the components of the P-variate (P-channel) time series and obtain first set of MFIBF . The first set of residue is obtained as follows:

 rp1(t)=xp(t)−yp1(t),p=1,2,⋯,P. (19)

Apply ZP-HPF with cutoff frequency to set of residue and obtain second set of MFIBF . The second set of residue is obtained as

 rp2(t)=rp1(t)−yp2(t)p=1,2,⋯,P. (20)

We can repeat this ZP-HPF procedure times and obtain final set of MFIBF and residue (with cutoff frequency )

 rpℓ(t)=rp(ℓ−1)(t)−ypℓ(t)p=1,2,⋯,P. (21)

Through the addition of (19), (20) and (21) we obtain expression, similar to (1), for P-variate time series as

 xp(t)=ℓ∑i=1ypi(t)+rpℓ(t),p=1,2,⋯,P. (22)

When we use the Fourier based zero-phase filtering, as in Eq. (8), to obtain MFIBFs, first two conditions of FIBFs are fully satisfied and the third one is approximately satisfied (i.e. satisfied in all practical sense), obviously, it can not be guaranteed simultaneously to all P-channel data. This is similar to MEMD algorithm problem where in derivation of multivariate IMFs, first condition of IMF is not imposed [9].

The question is, how to obtain cutoff frequencies (CFs) corresponding to zero-phase high pass filters ? There is lot of flexibility and are various ways to select CFs, e.g. dyadic (i.e. , where is the maximum frequency of a signal and for the sampled signal, maximum frequency is () half of the sampling frequency), non-dyadic, uniform and non-uniform CFs. We can take the Fourier transform of signal to obtain its spectrum details and make strategy to decide CFs.

For narrowband signal, we define ratio of center frequency () to bandwidth (BW) as

 m=fCi/(fHi−fLi),fCi=(fHi+fLi)/2, (23)

where is the highest frequency and is the lowest frequency of th band of a filter bank. From (23) we obtain

 fLi=[(2m−1)/(2m+1)]fHi,m>1/2. (24)

From (23) and (24), we observe that the ratios, for the consecutive th and th bands, of center frequencies , CFs and BWs can be taken as a constant, i.e.

 fCi/fCi+1=fci/fci+1=(fHi−fLi)(fHi+1−fLi+1)=l. (25)

From (23), (24) and (25), we obtain or with , and as [, ], [, ]. Here, we have liberty to select any suitable value of or , and greater the value of (or lesser the value of ) narrower the band, whereas in the case of dyadic FB and are fixed values. If required, we can vary the value of (or ) for each band rather than taking the fixed value. Thus, we here propose the compact and elegant way to decide CFs as summarized in Algorithm C. Algorithm C: An algorithm to obtain cutoff frequencies . Select suitable value of and set . Set . Set . Repeat step 2 to 3 for .

In MFDM, we can use zero-phase low pass filtering (ZP-LPF) in place of ZP-HPF to decompose signal in order of residue to first MFIBFs, i.e. . We use zero-phase filtering as it preserves salient features (e.g. maxima, minima, etc.) in the filtered time waveform exactly at the time where those features occur in the unfiltered waveform, whereas conventional (non zero-phase) filtering shifts the features in the signal and hence cannot be used. The zero-phase filtering of time series can be obtained through the finite impulse response (FIR) or infinite impulse response (IIR) filters.

Similar to the MEMD and noise-assisted MEMD (NA-MEMD) [16], this MFDM algorithm produces the equal number of scale-aligned MFIBFs for all channels and preserving joint channel properties that make it suitable for direct multichannel modelling. The FDM does not suffer from mode mixing, detrend uncertainty and end effect artefacts as extraction of FIBFs does not depend on distribution of local extrema across the range of signal.

Table 1 presents comparisons among Fourier, Wavelet, EMD [27] and FDM methods in data analysis.

## 5 Simulation results

The online available MATLAB software of MEMD [21], EMD and EEMD [22] have been used in simulation results.

### 5.1 Multivariate data decomposition

We used quadri-variate time series signal, which is summation of sinusoids (with combination of frequencies ) and Gaussian white noise of zero mean and standard deviation of 0.2., i.e. (for ). A sinusoid is present to first, second and fourth channels; a sinusoid is present to first, second and third channels; a sinusoid is common to all channels; a sinusoid is present to first, third and fourth channels. On the same machine, computation time for MFDM is 0.45 sec. and for MEMD is 69.5 sec. in this simulation. The MFDM algorithm, similar to MEMD, generating perfectly aligned intrinsic bands, as shown in Figure 1, in all the four channels, whereas, MFDM is computationally more efficient.

### 5.2 Intermittency and mode mixing

The intermittency, in the time series, is a main cause of mode mixing [e.g. signal in Figure 2(a)] and mode splitting in EMD algorithm. These issues are mitigated by EEMD and NA-MEMD. The decompositions of a signal through FDM algorithm is shown in Figure 2(b) without end effect artefacts. The MFDM is able to localize the mono-component sinusoid within a single FIBF and outperforming EMD 2(c) and EEMD 2(d). On the same machine, computation time for MFDM, EMD and EEMD are 0.56 sec., 0.21 sec. and 77.18 sec., respectively. The ensemble size for EEMD was N = 500 with the 16.94 dB signal-to-noise power ratio.

### 5.3 Time-Frequency-Energy Analysis

Figure 3 shows time-frequency-energy (TFE) estimates for a nonstationary signal mixture of a linear chirp and frequency modulated (FM) sinusoid, obtained using the FDM and EMD. There is a enhanced TFE tracking when using FDM with low to high frequency scan and other one (HTL-FS) is similar to plot obtained by EMD algorithm.

### 5.4 Intrawave frequency modulation

First, we decompose the following signal that has intrawave frequency modulation and it is considered challenging because the instantaneous frequency itself has very high frequency modulation [23]

 x(t)=11.2+cos(2πt)+cos(32πt+0.2cos(64πt))1.5+sin(2πt) (26)

Second, we consider a model wave

 x(t)=cos(ωt+ϵsinωt) (27)

that satisfies the following highly nonlinear differential equation [1]

 d2x(t)dt2+[ω+ϵωcos(ωt)]2x(t)−[ϵω2sin(ωt)]√1−x2(t)=0 (28)

with and . We demonstrate that, Figure 5, 6 and 7, our method well applies to these challenging cases with good accuracy. These examples clearly demonstrate that the FDM can indeed analyze nonlinear signals and it is a nonlinear decomposition method.

### 5.5 Analysis of a white Gaussian noise

Figure 8 shows the TFE analysis of a white Gaussian noise (with zero mean, unit variance, 1024 samples and sampling frequency Hz) obtained from the FDM and EMD algorithm. Clearly, both LTH-FS and HTL-FS views of TFE is similar and complete data is decomposed in FIBFs. Figure 9 shows the power spectral density (PSD) plot of same white Gaussian noise with the FDM and EMD algorithm. The FDM has dived the complete data in narrowband and orthogonal FIBFs. Both LTH-FS and HTL-FS views of PSD looking similar but FIBFs have different frequency band, e.g. approximately 40 Hz is cutoff frequency for one of the band in PSD (LTH-FS), whereas, it is mid frequency of the one band in other PSD (HTL-FS) view. There are enhanced TFE and PSD tracking when using FDM.

### 5.6 TFE Analysis of unit sample sequence

The unit sample sequence is defined as at and zero otherwise. By using the relation , we obtain the analytic representation of as

 z[n] =sin(π(n−n0))+j[1−cos(π(n−n0))]π(n−n0) (29) z[n] =a[n]exp(jϕ[n]),

where real part of is , , and hence which corresponds to half of the Nyquist frequency, i.e. Hz. Figure 10 shows the plots of real, imaginary part and absolute value of with , sampling frequency Hz, length . Theoretically, this clearly indicate that most of the energy of signal is concentrated at time sec. () and frequency Hz.

Figure 12 shows the plots of FIBFs obtained by the FDM and Figure 11 shows the TFE analysis of unit sample sequence from the FDM, EEMD and continuous wavelet transform (CWT) methods. This clearly indicate that the TFE plot obtained by the FDM method is same as theoretical estimation, whereas there is energy spread over a range of frequencies and lack of accuracy in TFE plot by the EEMD and CWT methods.

The uncertainty principle is a consequence of the Fourier transform (or any other type of integral transform) pair in the frequency definition. Therefore, its limitation could only be applied to such integral transforms, in which time would be ‘integrated out’ or smeared over the integral time limit. Consequently, if we avoid an integral transform in the frequency computation, time-frequency analysis of signal would not be bounded by the uncertainty principle [26]. The IF is defined through differentiation of phase rather than integration and, hence, overcome the restriction of the uncertainty principle. One can obtain arbitrary precision on time and frequency resolution subject only to the sampling rate in case of discrete time signal. The uncertainty principle in signal analysis state that the finer the time resolution one wants, the cruder the resulting frequency resolution would be. However, this example clearly demonstrate that the FHS by FDM is indeed not limited by uncertainty principle and signal can be highly concentrated in both time and frequency domain.

### 5.7 Application to a earthquake signal analysis

Earthquake time series signal is nonlinear and nonstationary data. The Elcentro Earthquake data (sampled at ) has been taken from [25] and is shown in Figure 13. The critical frequency range that matter in the structural design is less than , and from the Fourier based power spectral density (PSD), marginal spectrum by FDM and marginal spectrum by EMD in Figure 13 show that almost all the energy in this data is within . The instantaneous energy fluctuations by the FDM and EMD methods, as shown in Figure 14, are similar in nature. The TFE distribution by the FDM, EMD and CWT methods are shown in Figure 15, and all three methods indicate maximum energy concentration around and 2 second. There is enhanced TFE tracking by FDM methods as it provide better details of how the different waves arrive from the epical center to the recording station, e.g. the compression waves of small amplitude but higher frequency range of to , the shear and surface waves of strongest amplitude and lower frequency range of below which does most of the damage, and other body shear waves which are present over the full duration of the data span.

### 5.8 Application to a speech signal analysis

To illustrate the advantages of FDM in speech signal analysis, we decompose a quasiperiodic voice signal (small-cap I vowel). In Figure 16, the FIBFs to seem to extract most of the noise, while to seem to express most of the signal information at high frequencies. Moreover, the FIBF captures accurately , the fundamental frequency (inverse of glottal pulse or pitch period length) of the voice signal. In contrast, EMD (Figure 17 (a)) presents mode mixing between modes two and three, between modes three and four, and between modes four and five. Also, the EMD is not able to catch in any mode. In Figure 17 (b), the EEMD’s is able to catch in mode , although it can be observed that positive minima and negative maxima in the mode and , there is clear violation to IMF condition (1). Besides this problem, the major drawbacks in this case are huge computation complexity for this ensemble size and the reconstruction error that may be significant. The computation complexity can be reduced by selecting the smaller ensemble size but that increases the reconstruction error. Figure 18 shows the plots of TFE obtained by the FDM, EMD and EEMD algorithms for same data (vowel ‘small-cap I’). Clearly, there is enhanced TFE tracking when using FDM. The results in this application as well are in clear favor of the Fourier method proposed in this paper.

## 6 Conclusion

In this paper, we have proposed: (1) The Fourier Decomposition Method (FDM), for nonlinear and nonstationary time series analysis, which decomposes any data into a small number of ‘Fourier intrinsic band functions’ (FIBFs). The FDM is a generalized Fourier expansion with variable amplitudes and frequencies of a time series by the Fourier method itself. (2) The zero-phase filter bank based multivariate FDM (MFDM) algorithm, for the analysis of multivariate nonlinear and nonstationary time series, which is generating finite number of band limited multivariate FIBFs (MFIBFs). (3) An algorithm to obtain cutoff frequencies required in MFDM algorithm for zero-phase high or low pass filtering of multivariate signals.

The fundamental and conceptual contributions of the this study are the Fourier based decomposition method (i.e. FDM) and the introduction of the FIBFs. The FIBFs form the basis of the decomposition that are complete, orthogonal, local and adaptive. The instantaneous frequencies of the FIBFs produce a time-frequency-energy distribution of any signal. A time-frequency-energy distribution of a signal is used in various fields of science and engineering for analysis of physical phenomena and engineering systems. The proposed methods produce the final presentation of the results in a time-frequency-energy distribution that reveals the imbedded structures of a signal. Unlike the various EMD algorithms, the FDM and MFDM are mathematically well defined, supported by the well established theories of filter and the Fourier transforms. The FDM and MFDM methods do not suffer from mode mixing, detrend uncertainty and end effect artefacts as extraction of FIBFs does not depend on distribution of local extrema across the range of data. Simulation results demonstrate the efficacy of the proposed methods.

### .1 The FDM for continuous time real function

Let be a non-periodic, real function of time and follow the Dirichlet conditions, then the Fourier transform (FT) of is defined as

 X(f)=∫∞−∞x(t)exp(−j2πft)dt (30)

and inverse Fourier transform is defined as

 x(t)=∫∞−∞X(f)exp(j2πft)df (31)

It is easy to show that . From Eq. (30) and, hence, we rewrite Eq. (31) as

 x(t)=∫∞0[X(f)exp(j2πft)+X∗(f)exp(−j2πft)]df (32)

In this Eq., second term is complex conjugate of first term. As is a real function, we can write

 x(t)=Re{z(t)} (33)

where analytic function and denote real part of function . We write Eq.(33) as

 2∫∞0X(f)exp(j2πft)df=M∑i=1ai(t)exp(jϕi(t)) (34)

where (with )

 ai(t)exp(jϕi(t))=2∫fifi−1X(f)exp(j2πft)df, (35)

for . To obtain minimum number of AFIBFs in low to high frequency scan, for each , start with , increase and select the maximum value of such that and

 ai(t),fi(t)=12πdϕi(t)dt≥0,∀t. (36)

Similarly, in high to low frequency scan, the lower and upper limits of integration in (35) will change to to , respectively, with , and we can obtain minimum number of AFIBFs by selecting the minimum value of such that and Eq. (36) is satisfied.

From Eq. (32) and (33), it is easy to obtain relationship between the energy of original signal and energy of its analytic signal as

 Ex=Ez2, (37)

that is, the energy of analytic signal is twice of the energy of original signal.

The IF characterizes a local frequency behavior as a function of time. In a dual way, the local time behavior as a function of frequency is characterized by the group delay (GD) : . The GD measures the average time arrival of the frequency . In general, the IF and GD define two different curves in the time-frequency plane. Similar to the IF process, for causal signal , we obtain

 Extra open brace or missing close brace (38)

with such that , for .

### .2 The FDM for discrete time real function

Let be a non-periodic, real function of time and follow the Dirichlet conditions, then the discrete time Fourier transform (DTFT) of is defined as

 X(ω)=∞∑−∞x[n]exp(−jωn) (39)

and inverse discrete time Fourier transform (IDTFT) is defined as

 x[n]=12π∫π−πX(ω)exp(jωn)dω (40)

It is easy to show that . From Eq. (39) and, hence, we rewrite Eq. (40) as

 x[n]=12π[∫π0X(ω)exp(jωn)dω+∫π0X∗(ω)exp(−jωn)dω]. (41)

In this Eq., second term is complex conjugate of first term. As is real function, we can write

 x[n]=Re{z[n]} (42)

where analytic signal