Novel Fourier quadrature transforms and analytic signal representations for nonlinear and non-stationary time-series analysis

The Hilbert transform (HT) and associated Gabor analytic signal (GAS) representation are well known and widely used mathematical formulations for modelling and analysis of signals in various applications. In this study, like the HT, to obtain quadrature component of a signal, we propose novel discrete Fourier cosine quadrature transforms (FCQTs) and discrete Fourier sine quadrature transforms (FSQTs), designated as Fourier quadrature transforms (FQTs). Using these FQTs, we propose 16 Fourier quadrature analytic signal (FQAS) representations with following properties: (1) real part of eight FQAS representations is the original signal, and imaginary part of each representation is FCQT of real part; (2) imaginary part of eight FQAS representations is the original signal, and real part of each representation is FSQT of imaginary part; (3) like the GAS, Fourier spectrum of all FQAS representations has only positive frequencies; however, unlike the GAS, real and imaginary parts of FQAS representations are not orthogonal. The Fourier decomposition method (FDM) is an adaptive data analysis approach to decompose a signal into a set Fourier intrinsic band functions. This study also proposes new formulations of the FDM using discrete cosine transform with GAS and FQAS representations, and demonstrates its efficacy for improved time-frequency-energy representation and analysis of many real-life nonlinear and non-stationary signals.

Novel Fourier quadrature transforms and analytic signal representations for nonlinear and non-stationary time-series analysis Pushpendra Singh School of Engineering and Applied Sciences, Bennett University, Greater Noida, India PS, 0000-0001-5615-519X The Hilbert transform (HT) and associated Gabor analytic signal (GAS) representation are well known and widely used mathematical formulations for modelling and analysis of signals in various applications. In this study, like the HT, to obtain quadrature component of a signal, we propose novel discrete Fourier cosine quadrature transforms (FCQTs) and discrete Fourier sine quadrature transforms (FSQTs), designated as Fourier quadrature transforms (FQTs). Using these FQTs, we propose 16 Fourier quadrature analytic signal (FQAS) representations with following properties: (1) real part of eight FQAS representations is the original signal, and imaginary part of each representation is FCQT of real part; (2) imaginary part of eight FQAS representations is the original signal, and real part of each representation is FSQT of imaginary part; (3) like the GAS, Fourier spectrum of all FQAS representations has only positive frequencies; however, unlike the GAS, real and imaginary parts of FQAS representations are not orthogonal. The Fourier decomposition method (FDM) is an adaptive data analysis approach to decompose a signal into a set Fourier intrinsic band functions. This study also proposes new formulations of the FDM using discrete cosine transform with GAS and FQAS representations, and demonstrates its efficacy for improved time-frequency-energy representation and analysis of many real-life nonlinear and non-stationary signals.

Introduction
The Fourier theory is the most important mathematical tool for analysis and modelling physical phenomena and engineering systems. It has been used to obtain solutions in almost all fields of science and engineering problems. It is the fundamentals of a representations. The FDM is computationally efficient due to FFT implementation, produces desired number of FIBFs with required cut-off frequencies. These attributes of the FDM are useful in many applications.
Thus, in this study, we present FQTs as effective alternatives to the HT, and FSAS representations as alternatives to the GAS representation for nonlinear and non-stationary time-series analysis. Generally, we acquire a continuous time (CT) signal and convert it to discrete-time (DT) signal for efficient compression, storage, transmission, representation and analysis. So, we consider and discuss only DT representations in detail, and discussion of its CT counterpart is limited, but can be easily obtained by analogy to the DT part. This study is organized as follows: a brief overview of the analytic signal representation and the FDM is presented in §2. FQTs, FSAS representations and new formulations of the FDM are presented in §3. Simulation results and discussions are presented in §4. Section 5 presents conclusion of the work.

A brief overview of the GAS representation, IF and FDM
The GAS representation [2] is a complex-valued function, z[n], that has only positive frequency components in the Fourier spectrum, and it is defined as where the real part of GAS is the original signal and the imaginary part is the HT of the real part, and real and imaginary parts are orthogonal to each other (i.e. inner product hx[n],x[n]i ¼ 0). The HT of a signal is defined asx implemented in real applications, because, its impulse response is unstable (i.e. absolutely not summable as P 1 n¼À1 jh[n]j ¼ 1), non-causal and has infinite time support. Practically, the GAS representation z[n], from a real signal x[n] of length N, is obtained by the inverse discrete Fourier transform (IDFT) as [14,56] , 0 n, k N 2 1, exp (jpn) ¼ (À1) n , and X[N/2] is the highest frequency component of the Fourier spectrum. This is the only practical approach, based on the Fourier theory and being used by Matlab as well, which is available in the literature to obtain the GAS representation that satisfies the following properties: (P1) only positive frequencies are present in the Fourier spectrum, (P2) real part is the original signal ð2:4cÞ   [1] . . . X si [N 2 1]] T are the DCT and DST of ith type, respectively. Thus, we have defined linear transformations of x intox ci , and x intox si with orthogonal transformation matricesS T i C i and C T i S i , respectively. These transformation matrices are orthogonal due to the properties of an orthogonal matrix (i.e. if Q is an orthogonal matrix, then so is Q T and Q T ¼ Q 21 ; if Q 1 and Q 2 are orthogonal matrices, then so is Q 1 Q 2 ). The continuous time FCQT, FSQT and corresponding FSAS representations are defined in appendix A. The proposed two-dimensional FSAS representations are defined in appendix B. Now, we consider the complete process of obtaining FQT, corresponding FSAS and FDM using DCT-2 as follows. The DCT-2 of a sequence, x[n] of length N, is defined as [15] and inverse DCT (IDCT) is obtained by where normalization factors s k ¼ 1= ffiffi ffi 2 p for k ¼ 0, and s k ¼ 1 for k = 0. If consecutive samples of sequence x[n] are correlated, then DCT concentrates energy in a few X c2 [k] and decorrelates them. The DCT basis sequences, cos(pk(2n þ 1)/2N), which are a class of discrete Chebyshev polynomials [15], form an orthogonal set as inner product kcos( pk(2n þ 1)/2N), cos(pm(2n þ 1)/2N)l ¼ 0 for k = m, and The proposed FCQT (3.7) uses the basis sequences, sin(pk(2n þ 1)/2N), whose inner product is defined as and, therefore, the set of these basis vectors form an orthogonal set for 1 k N 2 1. From (3.4) and (3.7), one can easily observe that the FCQT of a constant signal, like the HT, is zero. The vector space theoretic approach and explanation of this proposed transformation is as follows. Let V and W be vector spaces, a function T:V ! W is called a linear transformation if for any vectors v 1 , v 2 [ V and scalar c, Here in this study, V is a vector space spanned by the set of cosine basis vectors, fcos( pk(2n þ 1)/2N)g for 0 k N 2 1, with dimension N, and W is a vector space spanned by the set of sine basis vectors, fsin( pk(2n þ 1)/2N )g for 0 k N 2 1, with dimension (N 2 1). Therefore, these two vector spaces are homomorphic, transformation is linear and non-invertible as a constant vector is mapped to the zero vector. However, for a zero-mean vector x[n] (i.e. X c2 [0] ¼ 0), these two vector spaces are isomorphic, transformation is linear and invertible, and a set of orthogonal cosine basis vectors are mapped to the set of orthogonal sine basis vectors, i.e. fcos(pk(2n þ 1)/2N)g7 !fsin( pk(2n þ 1)/2N)g or T(cos(pk(2n þ 1)/2N)) ¼ (sin( pk(2n þ 1)/2N)) for 1 k N 2 1. Thus, the linear transformation of x[n] intox c2 [n], defined in (3.3) and (3.7) with transformation matrixS T 2 C 2 , is (a) non-invertible if the mean of signal x[n] is non-zero and (b) invertible if the mean of signal x[n] is zero. In practice, we can always remove the mean from x[n] to make it a zero-mean vector that leads the proposed transformation to be isomorphic.
We proposed and derived 16 different types of FSAS representations in (3.3) using eight types of DCTs and eight types of DSTs. However, in this study, we especially consider DCT-2 in detail to obtain the FSAS as , real part ofz c2 [n] is the original signal, and imaginary part ofz c2 [n] is the FQT of real part. It is worthwhile to note that the signal x[n] and its FQTx c2 [n] are not orthogonal, i.e. hx[n],x c2 [n]i = 0. To prove this, we consider the inner product of basis sequences, cos( pk(2n þ 1)/2N) and sin(pm(2n þ 1)/2N), and using trigonometric manipulation [2cos(a)sin . We obtain sum of exponential series as Now, we propose the DCT-based FDM (i.e. DCT-FDM), and devise two new approaches to decompose a signal into a small set of AM -FM signals (i.e. FIBFs) and corresponding analytic representations. In the first approach, we obtain the decomposition of signal x[n] (excluding the DC component) into a set of FIBFs, fx i [n]g M i¼1 , and corresponding AFIBFs, {z c2i [n]} M i¼1 , using the FSAS representation (3.11) as , which can be written as In the second approach, using the DCT-FDM, we decompose the signal x[n] into the same set of FIBFs, fx i [n]g M i¼1 , and apply the HT to obtain corresponding set of AFIBFs, {ẑ i [n]} M i¼1 , using the GAS representation asẑ It is pertinent to note that in the FDM (3.15) all five properties (P1) -(P5) of the GAS are satisfied; however, in the FDM (3.14), which uses FSAS, only the first two properties (P1) and (P2) of the GAS are being satisfied. The block diagrams of the FDM, using DCT-based zero-phase filter-bank to decompose a signal into a set of desired frequency bands, are shown in figure 1 We have used ZPF where frequency response of filter is one in the desired band and zero otherwise. Moreover, one can use any other kind of ZPF (i.e. such as the Gaussian filter to decompose a signal. The FDM advocates to use ZPF because it preserves salient features such as minima and maxima in the filtered waveform exactly at the position where those features occur in the original unfiltered waveform. Based on the requirement and type of application, we can devise methods to select ranges of frequency parameter k in (3.17), e.g. one can use three approaches to divide complete frequency band of a signal under analysis into small number of equal, dyadic and equal-energy frequency bands.
To obtain TFE representation of a signal using (3.13) and (3.15), we write FSAS and GAS in polar representation, for 1 i M, as Algorithms to implement FDM using the DFT and finite impulse response (FIR) filtering are presented in [14,41]. We summarize the steps to implement the proposed FDM, using DCT to decompose a signal x[n] into a set of M-AFIBFs, in algorithm 1, and complete Matlab implementation of the FDM using DFT, DCT and FIR filtering is made publicly available to download at the Dryad Digital Repository at https://doi.org/10.5061/dryad.jc21t36 [58] and at https://www.rese archgate.net/publication/326294577_MATLABCodeOfFDM_DCT_DFT_FIR_FSASJuly2018. From figure 1 as well as algorithm 1, we observe that FDM implementation requires one DCT/DFT and M number of IDCT/IDFT, so computational complexity is (M þ 1) times of FFT algorithm that has complexity of order O(N log 2 N).

Results and discussions
In this section, to demonstrate the efficacy of the proposed methods-DCT-and FSAS-based FDM (DCT-FSAS -FDM), DCT-and GAS-based FDM (DCT -GAS-FDM)-we consider many simulated as well as real-life data, and compare the obtained results with other popular methods such as EMD-, CWT-, DFT-and GAS-based FDM (DFT -GAS-FDM), FIR-and GAS-based FDM (FIR -GAS-FDM). We primarily consider those signals which have been widely used in literature for performance evaluation and results comparison among proposed and other existing methods.

A unit sample sequence analysis
First, we consider an analysis of unit sample (or unit impulse or delta) sequence using the proposed FSAS and compare the results with GAS representation. A unit sample function is defined as d[n 2 n 0 ] ¼ 1 for n ¼ n 0 , and d[n 2 n 0 ] ¼ 0 for n = n 0 . Using the inverse Fourier transform, analytic representation, is d[n 2 n 0 ] ¼ sin(p(n 2 n 0 ))/p(n 2 n 0 ),a[n] ¼ sin ((p=2)(n À n 0 ))=((p=2)(n À n 0 ))) (IA), f[n] ¼ (p/2)(n 2 n 0 ) (IP), and therefore v[n] ¼ p/2 (IF) which corresponds to half of the Nyquist frequency, that is, F s /4 Hz. Theoretically, this representation demonstrates that most of the energy of signal d[n 2 n 0 ] is concentrated at time t ¼ 4.99 s (n 0 ¼ 499) and frequency f ¼ 25 Hz, where sampling frequency, F s ¼ 100 Hz, and length of signal, N ¼ 1000, are considered. The plots of IA, IF and TFE using the GAS and FSAS representations are shown in figure 2a,b, respectively. In the GAS representation, figure 2a, IF is varying between 0 and 50 Hz at both ends of signal and converges to theoretical value only at the position of delta function. On the other hand, the FSAS representation yields better results as it produces correct value of IF, figure 2b, for all the time.

Chirp signal analysis
In this example, we consider an analysis of a non-stationary chirp signal using the proposed FSAS and compare the results with GAS approach, with sampling frequency F s ¼ 1000 Hz, time duration t [ [0, 1) s and frequency f [ [5, 100) Hz. Figure 3a shows chirp signal (5 -100 Hz) (i), proposed FCQT (ii) obtained using (3.7), which is also imaginary part of the proposed FSAS (3.11), and HT (iii) that is imaginary part of the GAS representation (2.3), which has unnecessary distortions at both ends of the signal. The TFE distributions obtained (i) using the FSAS representation is shown in figure 3b and (ii) using the GAS representation is shown in figure 3c, which has lot of energy spreading over wide range of time -frequency plane, at both ends of the signal under analysis.

Speech signal analysis
Estimation of the instantaneous fundamental frequency F 0 is the most important problem in speech processing, because it arises in numerous applications, such as language identification [58], speaker recognition [59], emotion analysis [60], speech compression and voice conversion [61,62]. Figure 4a shows a segment of voiced speech signal (i) from the CMU Arctic database [63] sampled at F s ¼ 32 000 Hz, corresponding electroglottograph (EGG) signal (ii) and differenced EGG (iii) signal. Figure 4b(i) shows magnitude spectrum of speech signal where F 0 is around 132 Hz, and energy of speech signal is concentrated on harmonics of F 0 , and the magnitude spectrum of differenced EGG (DEGG) signal is shown in plot (ii) that also indicate F 0 around 132 Hz. Figure 4c shows a decomposition of speech segment into a set of 10 FIBFs (FIBF1-FIBF10), a high-frequency component (FIBF11) and a low-frequency component (LFC) using the FDM, where FIBF1 captures F 0 , and harmonic components of F 0 are encapsulated in FIBF2-FIBF10. Figure 4d presents TFE representation obtained from the FDM where F 0 and its harmonics are clearly separated in distinct frequency bands.
Jain & Pachori [53] studied event-based method for F 0 estimation from voiced speech based on eigenvalue decomposition of Hankel matrix, and method uses iterative approach to estimate F 0 , which is computationally complex; moreover, it is difficult for method to estimate all harmonics of F 0 accurately. On the other hand, the FDM-based approach presented in this study is able and efficient to implement, achieve and follow the proposed model (2.7) precisely.

Noise removal from ECG signal
ECG signals, records of the electrical activity of the heart, are used to examine the activity of human heart. There are various problems that may arise while recording an ECG signal. It may be distorted due to the presence of various noises such as channel noise, baseline wander (BLW) noise (of generally below 0.5 Hz), power-line interference (PLI) of 50 Hz (or 60 Hz) and physiological artefacts. Owing to these noises, it becomes difficult to diagnose diseases, and thus appropriate treatment may be impacted [64]. The BLW and PLI removal from an ECG signal (obtained from MIT-Arrhythmia database, Fs ¼ 360 Hz) using the proposed method is shown in figure 5: (a) a segment of clean ECG signal (i), ECG signal heavily corrupted (i.e. signal-to-noise ratio (SNR) is low, in the range %218.4 dB) with BLW and PLI (ii) noises, (b) ECG signal after noise removal (i), separated BLW (ii) and PLI (iii). Thus, proposed FDM can be used to remove BLW and PLI noises, and recover ECG signal even in scenarios where SNR is rather poor. To present a numerical comparison, we consider input SNR i and output SNR o , which are defined as

Trend and variability estimation from a nonlinear and non-stationary time series
Estimating trend and performing detrending operations are important steps in numerous applications, e.g. in climatic data analyses, the trend is one of the most critical parameters [65]; in spectral analysis and in estimating the correlation function, it is necessary to remove the trend from the data to obtain meaningful results [65]. Thus, in signal and other data analysis, it is crucial to estimate the trend, detrend the data and obtain variability of the time series at any desired timescales of interest.
Here, we consider the global surface temperature anomaly (GSTA) data (of 148 Years from 1856 to 2003, publicly available at http://rcada.ncu.edu.tw/research1_clip_ex.htm.)) analysis, its trend and variabilities obtained from the proposed method, as shown in figure 6 To demonstrate the efficacy of the proposed FDM to estimate a trend and variability from data at any desired timescale, e.g. we estimated the trends and variabilities of GSTA data in various timescales as shown in figures 6d, 7 and 8. Figure 6d shows GSTA data (blue solid), and trends (red dashed) in multiple of 10 (i.e. 10 -80) years or longer timescale variations of temperature. Figure 7 presents: (a) GSTA data, and its trend in 15, 25 and 53 years or longer timescale and (b) GSTA data, and its trend in 75, 100 and 125 years or longer timescale. The variabilities of the annual GSTA data corresponding to various trends (i.e. data minus corresponding trends) with their timescales are given in figure 8: (a) variabilities corresponding to trends in figure 6d and (b) variabilities corresponding to trends in figure 7. From these figures, it is clear that the variabilities up to 53 years or longer timescale are similar and do not contain any lowfrequency component (i.e. there is no mode-mixing). However, the variabilities from 54 years (experimentally found accurately 54 years as shown in figure 9) or longer timescale are different from up to 53 years and also contain low-frequency components (i.e. there is mode-mixing).
The trend and variability analysis of this climate data is also performed by the EMD algorithm in [65]. The EMD algorithm decomposes the climate data into a set of IMFs which have overlapping frequency spectra and their cut-off frequencies are also not under control, therefore the trend and variability of data at desired timescales with clear demarcation cannot be obtained. For example, we have shown that there is a significant change in the trend and variability of data in just 1 year difference of timescale (observe  figure 9). This kind of fine control along with determination of trend and variability of data in desired timescale with clear demarcation can be obtained by the proposed method, which is not achievable by the EMD and its related algorithms.

Seismic signal analysis
An earthquake time series is nonlinear and non-stationary in nature, and in this example, we consider the Elcentro Earthquake data. The Elcentro Earthquake time series (which is sampled at F s ¼ 50 Hz) has been downloaded from http://www.vibrationdata.com/elcentro.htm and is shown in figure 10a(i). The most critical frequency range that matters in the structural design is less than 10 Hz, and the Fourier power spectral density (PSD), figure 10a(ii), shows that almost all the energy in this earthquake time series is present within 10 Hz.  Gravitational waves (GWs), predicted in 1916 by Albert Einstein, are ripples in the space-time continuum that travel outward from a source at the speed of light, and carry with them information about their source of origin. The GW event GW150914, produced by a binary black hole merger nearly 1.3 billion light years away [50], marks one of the greatest scientific discoveries in the history of human life. In this example, using the proposed method, we consider the noise removal, IF estimation and TFE representation of the GW event data which is publicly available at https://losc.ligo.org/ events/GW150914/ (sampling rate F s ¼ 16384 Hz). The GW signal sweeps upwards in frequency from 35 to 250 Hz, and amplitude strain increases to a peak GW strain of 1.0 Â 10 221 [50]. The noise removal and an accurate IF estimation of the GW data is important because IF is the basis to estimate many parameters such as velocity, separation, luminosity distance, primary mass, secondary mass, chirp mass, total mass and effective spin of binary black hole merger [50]. Figure 11 shows (a) the GW H1 strain data GW150914 (i) observed at the Laser Interferometer Gravitational-Wave Observatory (LIGO) Hanford, which is heavily corrupted with noise, the Fourier spectrum (ii) of the GW data, which is not able to capture the non-stationarity (i.e. upwards sweep in the frequency and amplitude) present in the data, (b) the TFE representation of the GW H1 data without decomposition by proposed method, which shows the signal frequency is increasing with time but having lots of unnecessary fluctuations in frequency due to noise, and (c) decomposition of data into a set of six FIBFs (FIBF1 (25- The further analysis, comparison, residue and TFE estimation of GW event GW150914 H1 strain data (captured at LIGO Hanford) using the proposed FDM are shown in figure 12: (a) GW H1 strain data ((i) red dashed line), proposed reconstruction ((i) blue solid line), and estimated residue component (ii), (b) numerical relativity (NR) data [50] ((i) red dashed line), and proposed reconstruction ((i) blue solid line), and difference between NR and reconstructed data (ii); TFE estimates of the NR data and the reconstructed data are shown in (c) and (d), respectively. The very same analysis of the GW event GW150914 L1 strain data (captured at LIGO Livingston) using the proposed method is shown in figure 13. These examples clearly demonstrate the efficacy of the proposed method for the analysis of real-life non-stationary signals such as speech ( §4.3), ECG ( §4.4), climate ( §4.5), seismic ( §4.6) and gravitational ( §4.7) time-series. This study is aimed to complement the current nonlinear and non-stationary data processing methods with the addition of the FDM, which is based on the DCT, discrete FCQT and zero-phase filter approach using GAS and FSAS representations.

Conclusion
The HT, to obtain quadrature component and GAS representation, is most well known, widely used and a key mathematical representation for modelling and analysis of signals to reveal underlying physical phenomenon in various applications. The fundamental contributions and conceptual innovations of this study are introduction of the new discrete FCQTs and discrete FSQTs, designated as FQTs, as effective alternatives to the HT. Using these FQTs, we proposed Fourier quadrature analytic signal (FQAS) representations, as coherent alternatives to the GAS, with following properties: (1) real part of FQAS is the original signal and imaginary part is the FCQT of the real part, (2) imaginary part of FQAS is the original signal and real part is the FSQT of the imaginary part, (3) like the GAS, Fourier spectrum of the FQAS has only positive frequencies; however, unlike the GAS, the real and imaginary parts of the proposed FQAS representations are not orthogonal to each other. Moreover, continuous-time FQTs, FQAS representations, and two-dimensional extension of these formulations are also presented. The Fourier theory and ZPF-based FDM is an adaptive data analysis approach to decompose a signal into a set of small number of Fourier intrinsic band functions (FIBFs). This study also proposed a new formulation of the FDM using the discrete cosine transform and the discrete cosine quadrature transform with the GAS and FQAS representations, and demonstrated its efficacy for analysis of nonlinear and non-stationary simulated as well as real-lifetime series such as instantaneous fundamental frequency estimation from unit sample sequence, chirp and speech signals, baseline wander and power-line interference removal from ECG data, trend and variability estimation and analysis of climate data, seismic time series analysis, and finally noise removal, IF and TFE estimation from the gravitational wave event GW150914 strain data.
Many time-frequency representation approaches such as the EMD algorithms, FDM and VMD are primarily based on the HT and GAS representation which are well known and widely used in numerous applications. This study proposes 16 FQTs as alternatives of the HT, 16 FQAS representations as alternatives of the GAS, and presents detailed simulation results using DCT-2-based FQT and FQAS, which generates better results (e.g. in terms of accurate instantaneous frequency estimation, higher output SNR, trend and variability estimation in desired timescale) in many but limited applications considered in this work. Moreover, the proposed novel methodology provides alternative ways of time-series representation and analysis, advances the existing FDM, which is useful many applications, and therefore, these transforms and corresponding analytic signal representations are required. Finally, as the DFT and DCTs are widely used tools in numerous different sets of applications, future direction of research would be to explore proposed FQTs and corresponding analytic signal representations in various applications.
Data accessibility. The FDM Matlab code is publicly available for download at the Dryad Digital Repository (doi:10.5061/ dryad.jc21t36) [58] and https://www.researchgate.net/publication/326294577_MATLABCodeOfFDM_DCT_DFT_ FIR_FSASJuly2018. All data used in this study are publicly available and links are included in the references of paper.