Discrete Time-Frequency Signal Analysis and Processing Techniques for Non-Stationary Signals

Show more

1. Introduction

In nature, most of the signals are non-stationary and time-varying signals. Further, the classical and modern methods are widely used to process the stationary signals in which they transform the signals from time-domain to frequency-domain and vice versa. The stationary signals do not change in their statistical properties over the length of the analysis time. Many signals of biological origin are varying in a random manner called non-stationary signals and are changing their properties over the length of the analysis time. The basic idea of time-frequency analysis is to design a joint function, which can describe the characteristics of signals on a time-frequency plan. Time-frequency transforms map a one-dimensional function of time x(t) into a two-dimensional function of time and frequency x(t, f) [1] .

In order to process such non-stationary signals, time-frequency analysis and processing methods are required. Generally, they fall into one of the two categories of time-frequency distributions (TFDs), the linear time-frequency distributions and the quadratic time-frequency distributions (QTFDs). The TFDs give useful information about frequency changes over time. The signal component could be considered as energy continuity in time without abrupt changes in frequency [2] .

Non-stationary signals comprise of mono component or multi-component. Linear TFDs, such as short-time Fourier transform (STFT), which is often used as a first choice of tool in time-frequency analysis, due to their simplicity in usage, well-established algorithm and analysis technique [3] . In order to get enhanced time-frequency resolution QTFDs have been introduced. QTFD classes are non-linear methods in which Wigner-Ville Distribution (WVD) is the primary distributions of QTFD class, from which so many classes called Cohen’s TFDs, have been introduced for various non-stationary signal-processing applications. Consequently, studies on the TFRs have been applied to analyze, modify and synthesize non-stationary signals or time-varying signals. In this paper, two types of time-frequency representation techniques are considered; Linear Time frequency distribution and quadratic time frequency distribution and their principle properties are investigated. The realization of this distribution for hardware and software platforms requires a discrete version. As a result, algorithms were developed for discrete time-frequency STFT and WVD techniques and were tested on non-stationary signals for joint time-frequency analysis.

2. Short-Time Fourier Transformation

STFT is one of the linear time-frequency representations based on the straightforward approach of slicing the waveform of interest into a number of short segments and performing the analysis on each of these segments, using standard Fourier transform. A window function is applied to segment the data, which effectively isolates the segment from the overall signal data, since the segment within the window is assumed as stationary and provides time localization. Then, Fourier Transform is applied to the windowed data and the spectrum or spectrogram could be calculated from the estimated Fourier coefficients.

The STFT of the signal x(t) is given by [4]

$X\left(t,f\right)={\displaystyle {\int}_{t-\tau /2}^{t+\tau /2}x\left(\tau \right)w\left(\tau -t\right){\text{e}}^{-j2\text{\pi}f\tau}\text{d}\tau}$ (1)

where $w\left(t-\tau \right)$ is a window function and $\tau $ is the variable that slides the window across the signal, x(t).

The discrete version of STFT of the signal x(n) is given by

$X\left(m,k\right)={\displaystyle {\sum}_{n=1}^{N}x\left(n\right)\left[w\left(n-k\right){\text{e}}^{-j\omega kn/N}\right]}$ (2)

where n is the time index, k is the frequency index and $w\left(n-k\right)$ is the analysis window that selectively determines the portion of x(n) for analysis. X(m, k) can be expressed as convolution of the signal $x\left(n\right){\text{e}}^{-j\omega kn/N}$ with the window function $w\left(n-k\right)$ . The spectrogram is the square of the magnitude of the STFT obtained in (2)

$PSD\left(t,\omega \right)={\left|X\left(m,k\right)\right|}^{2}$ (3)

Upon selection of discrete STFT, the next step is to select an appropriate window and its size where two closest sinusoids can be distinguished using Equation (3). However, non-stationary signals may involve a large number of sinusoids in close proximity. This results in a very small Δf and consequently a large window is required. This makes the STFT very similar to the Fourier transform and will hamper temporal resolution. In order to select an appropriate window size a novel empirical model is proposed in [5] [6] , which adaptively selects a window size and is given by

$W=\frac{3{B}_{s}{f}_{s}}{\mu}$ (4)

where f_{s} is the sampling frequency and μ = 386.3 for
$\Delta f=\frac{\mu}{3}$ . For rectangular

window, B_{s} = 2, Hanning/Hamming window B_{s} = 4 and for Blackman window B_{s} = 6.

3. Wigner and Wigner-Ville Distributions

All Quadratic Time-Frequency representations should satisfy the time and frequency shift invariance belong to general class of distributions introduced by Cohen and are given by the following expression [7]

$w\left(t,f\right)=\frac{1}{\text{2\pi}}{\displaystyle \iiint {\text{e}}^{-j\theta t}{\text{e}}^{-j\tau \omega}{\text{e}}^{-j\theta u}\varnothing \left(\theta ,\tau \right)x\left(u+\frac{\tau}{2}\right){x}^{*}\left(u-\frac{\tau}{2}\right)\text{d}u\text{d}\tau \text{d}\theta}$ (5)

where x(u) is the time signal, x*(u) is its complex conjugate and $\varnothing \left(\theta ,\tau \right)$ is an arbitrary function called the kernel. By choosing different kernels, different distributions are obtained. Wigner distribution is obtained by taking $\varnothing \left(\theta ,\tau \right)=1$ . Here, the range of all integrations is from −∞ to ∞.

A real valued signal x(t) is used in WDF, which has positive and negative frequency components and introduced aliasing or cross-terms between positive and negative frequencies in time-frequency domain.

Wigner-Ville Distribution

A simple approach to avoid aliasing is to use an analytic signal before computing the WDF. Ville (1948) proposed the use of the analytic signal in time-frequency representations of a real signal. An analytic signal is a complex signal that contains both real and imaginary components. The advantage of using analytical signal is that in the frequency domain the amplitude of negative frequency components are zero. The imaginary part is obtained by Hilbert transform. The analytic signal may be expressed by, [8] [9] ,

$z\left(t\right)=x\left(t\right)+jH\left[x\left(t\right)\right]$ (6)

where H[x(t)] is the Hilbert transform, which is generated by the convolution of the impulse response h(t) of 90˚ phase shift as follows

$H\left[x\left(t\right)\right]=x\left(t\right)\ast h\left(t\right)$ (7)

$h\left(t\right)=\{\begin{array}{l}2\frac{{\mathrm{sin}}^{2}\left(\text{\pi}\frac{t}{2}\right)}{\text{\pi}t},\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}t\ne 0\\ 0,\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.05em}}\text{\hspace{0.05em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}t=0\end{array}$

The discrete form of the equation is given by,

$H\left[x\left(n\right)\right]={\displaystyle {\sum}_{k=-\infty}^{\infty}h\left(n-k\right)x\left(k\right)}$ (8)

Substituting the kernel $\varnothing \left(\theta ,\tau \right)=1$ in Equation (5), the continuous time WVD is obtained for continuous time signal

$Wx\left(t,f\right)={\displaystyle {\int}_{-\infty}^{\infty}z\left(t+\frac{\tau}{2}\right){z}^{*}\left(t-\frac{\tau}{2}\right){\text{e}}^{-j2\text{\pi}f\tau}\text{d}\tau}$ (9)

where t is time domain variable, f is frequency domain variable and z(t) is analytical associate of the real signal x(t) obtained from Hilbert Transform. The Wigner-Ville Distribution (WVD) is the most powerful and fundamental time frequency representation [10] . The superior properties of the WVD over the STFT technique make it ideal for signal processing in such diverse fields as radar, sonar, speech, seismic and biomedical analysis [11] [12] . For these applications, there is a need of a flexible Wigner-Ville Distribution for non-stationary signal analysis.

The Discrete version of WVD of the signal x(n) is given by [13] [14] .

$W\left(n,m\right)=2{\displaystyle {\sum}_{k=-\infty}^{\infty}{\text{e}}^{-\frac{2\text{\pi}mn}{N}}z\left(n+k\right){z}^{*}\left(n-k\right)}$ (10)

$w\left(n,m\right)={\displaystyle {\sum}_{m=-\infty}^{\infty}{\text{e}}^{-\frac{2\text{\pi}nm}{N}}{R}_{xx}\left(n,k\right)}=FF{T}_{k}\left[{R}_{xx}\left(n,k\right)\right]$ (11)

where t = nTs and f = m/(NTs).

The WVD uses a variation of autocorrelation, where time remains in the result. This is achieved by comparing the waveform with itself for all possible lags, i.e., the comparison is done for all possible values of time. This comparison gives rise to the defining equation called instantaneous auto-correlation function for continuous time signal

${R}_{xx}\left(t,\tau \right)=z\left(t+\frac{\tau}{2}\right){z}^{*}\left(t-\frac{\tau}{2}\right)$ (12)

Its discrete version is

${R}_{xx}\left(n,k\right)=z\left(k+n\right){z}^{*}\left(k-n\right)$ (13)

where $\tau $ and n are the time lags as in autocorrelation, and * represents the complex conjugate of the signal z. The instantaneous autocorrelation function retains both lag and time. Some important properties examined shows that WVD is always real and satisfies time and frequency marginal characteristics.

4. STFT Algorithm Implementation

Fast Fourier Transform (FFT) is applied using straight forward approach using the separate function [B,t,f] coded in MATLAB. Here, B is a complex matrix containing the magnitude and phase of the STFT frequency spectrum with the rows encoding the time-axis and the column representing the frequency-axis and t and f are optional argument vectors that can be helpful in plotting.

5. DWVD Algorithm Implementation

The following are the steps involved to develop the algorithm:

Step 1: Convert the real signal into analytical signal using Hilbert transform.

Step 2: Compute the WVD using a separate function. The input function has arguments x and f_{s}.

Step 3: Compute the instantaneous autocorrelation using loop to construct an array.

Step 4: Find the WVD magnitude using FFT.

6. Spectrum Analysis of the Proposed Time-Frequency Distributions

1) Two Sequential Sinusoid

The proposed STFT and WVD techniques are tested over different inputs such as a two sequential sinusoid of 10 Hz and 50 Hz. The sinusoid is preceded and followed by a time gap of 0.5 sec. The simulated signals are shown in Figure 1.

The STFT magnitude spectrum and contour plot shown in Figure 2 and Figure 3 produces a time-frequency plot with the step change in frequency approximately at the correct time, but the step change of time and frequency are not defined very preciously.

The lack of finite support in either time or frequency is evident from the appearance of energy slightly before 0.5 sec and after 0.5 sec and energy other than 10 and 50 Hz as shown in Figure 3. Further, when the window length is increased, the frequency resolution increases but there is a decreases in the time resolution.

The DWVD magnitude spectrum and contour plot shown in Figure 4 and Figure 5 respectively, which generates a time-frequency plot with the step change in frequency approximately at the correct time but, the step change in time and frequency is defined very accurately. From this, it is inferred that the DWVD preserves the time and frequency properties, marginally. Further, it has finite support either in time or in frequency, which is obvious from the non-appearance of energy slightly before 0.5 sec and after 0.5 sec and energy other than 10 and 50 Hz. This is because the WVD not uses the window or

Figure 1. Two sequential sinusoids of 10 Hz and 50 Hz with a time gap of 0.5 sec.

Figure 2. Magnitude spectrum of two sinusoids using STFT.

Figure 3. Contour plot of two sinusoids using STFT.

kernel function $\varnothing \left(\theta ,\tau \right)=1$ , which introduces the cross term due to instantaneous autocorrelation.

2) Chirp Signal

A sinusoid that has increases in frequency over time is called a chirp signal. This signal can be generated by multiplying the argument of a sine function by a

Figure 4. Magnitude spectrum of two sinusoids using DWVD.

Figure 5. Contour plot of two sinusoids using DWVD.

linearly increasing term. A linearly increasing sine wave that varies between 10 and 200 Hz over a 1 sec period is generated as shown in Figure 6.

The STFT magnitude spectrum and contour plots shown in Figure 7 and Figure 8 respectively exhibit a time-frequency plot with the step change in frequency approximately at the correct time. Further, it shows the chirp signal is linear FM signal, which lead to the appearance of the step changes preciously both in time and frequency axes with no lag in the finite support. Similar to STFT, if there is an increase in window length; DWVD also increases the frequency resolution and decreases the time resolution.

In the DWVD magnitude spectrum and contour plots shown in Figure 9 and Figure 10 respectively, the DWVD produces a time-frequency plot with the step change in frequency approximately at the correct time, which shows that the chirp signal is a linear FM signal and so supports the step changes preciously both in time and frequency domains with finite support. The contour plot shows that the signal has low energy. In both the cases the marginal property of the signals are preserved. Hence, STFT and DWVD are suitable to analyse liner non-stationary signals.

3) ECG Signal

ECG signal is one of the non-linear multicomponent non-stationary signals. The ECG wave form discrete data is imported to MATLAB [15] and is shown in Figure 11.

Figure 6. Chirp signal.

Figure 7. STFT magnitude spectrum of chirp signal.

Figure 8. STFT-contour plot for chirp signal.

Figure 9. DWVD magnitude spectrum of chirp signal.

Table 1 shows the important parameter to design the window function.

From the two STFT spectrums of the Arrhythmia ECG signal, the window width plays a predominant role; since the 128-point window unravels the high frequency components very well than the 32-point, even though there is a compromise in time resolution as shown in Figure 12 and Figure 13. Thus, the

Figure 10. DWVD-contour plot for chirp signal.

Figure 11. Arrhythmia ECG signal with 10 sec time duration.

Figure 12. STFT magnitude spectrum with N = 1024, nfft = 128 and overlap = 64 for the Arrhythmia ECG signal with 10 sec time duration.

Table 1. Window parameter for STFT.

STFT not preserving the marginal property, due to the presence of low energy signals which lead to spectral leakage.

Figure 13. STFT magnitude spectrum with N = 512, nfft = 64 and overlap = 32 for Arrhythmia ECG signal with 10 sec time duration.

Figure 14. DWVD magnitude spectrum with N = 512 for Arrhythmia ECG signal with 10 sec time duration.

Whereas in DWVD magnitude spectrum shown in Figure 14, the actual signals are found to be in short distances so that local oscillation takes place that introduces cross term between the two auto terms. When it is dominant enough, it could not able to provide good time and frequency resolution. From the above discussion, if the signals are multi component non-linear non-stationary signals, the DWVD is not suitable to analyze the signal until the cross term is eliminated.

7. Conclusion

In this work, two time-frequency analysis methods viz., discrete STFT and WVD algorithms were developed and compared with their performance for the purpose of defining and applications of the time-frequency resolution of the non-stationary signals. The performance of these methods was tested in three different non-stationary signals and their merits and demerits were investigated. The results of this study revealed that the time-frequency resolution of the STFT technique is inversely related to the window length. Increasing the window length increases the frequency resolution, but at the cost of reduction in frequency tracking capability. Conversely, WVD has several advantages over the STFT. It reduces the cross terms and sampling frequency by using an analytical signal. A DWVD also maintains some of the properties such as marginal and invariability of the non-stationary signals. It also produces a good spectrum of time-frequency structure. In the DWVD, the kernel $\varnothing \left(\theta ,\tau \right)=1$ introduces cross terms. These cross-terms will be reduced by introducing the window, kernel and adaptive filters, which will make the DWVD a more suitable and powerful tool for non-stationary signal analysis. Since The Wigner Ville distribution preserves all the information, it will be used for two-dimensional signal processing like digital image processing. This work supports the need of using time-frequency distributions when dealing with non-stationary signals.

References

[1] Hlawatsch, F. and Boudreaux-Bartels, G.F. (1992) Linear and Quadratic Time-Frequency Signal Representations. IEEE Signal Processing Magazine, 9, 21-67.

https://doi.org/10.1109/79.127284

[2] Damira, M., Victora, S. and Car Zlatana, B. (2015) Optimizing the Reference Signal in the Cross Wigner-Ville Distribution Based Instantaneous Frequency Estimation Method. Procedia Engineering, 100, 1657-1664.

https://doi.org/10.1016/j.proeng.2015.01.540

[3] Boashash, B., Ed. (2003) Time-Frequency Signal and Processing: A Comprehensive Reference. Elsevier Science, Oxford.

[4] Allen, J.B. (1977) Short Term Spectral Analysis, Modification by Discrete Fourier Transform Synthesis. IEEE Transactions on Acoustics, Speech, and Signal Processing, 25, 235-238.

https://doi.org/10.1109/TASSP.1977.1162950

[5] Shibli. N., Khan, O.U. and Tariq, M. (2016) An Efficient Adaptive Window Size Selection Method for Improving Spectrogram Visualization. Computational Intelligence and Neuroscience, 2016, Article ID: 6172453.

[6] Shin, Y.S. and Jeon, J.-J. (1993) Pseudo Wigner-Ville Time-Frequency Distribution and Its Application to Machinery Condition Monitoring. Shock and Vibration, 1, 65-76.

https://doi.org/10.1155/1993/372086

[7] Claasen, T.A.C.M. and Mecklenbrauker, W.F.G. (1980) The Wigner Distribution—A Tool for Time-Frequency Signal Analysis Part Ill: Relations with Other Time-Frequency Signal Transformations. Philips Journal of Research, 35, 372-389.

[8] Marple Jr., S.L. (1999) Computing the Discrete-Time ‘Analytic’ Signal via FFT. IEEE Transactions on Signal Processing, 47, 2600-2603.

https://doi.org/10.1109/78.782222

[9] O’ Toole, J.M., Mesbah, M. and Boashash, B. (2008) A New Discrete Analytic Signal for Reducing Aliasing in the Discrete Wigner-Ville Distribution. IEEE Transactions on Signal Processing, 56, 5427-5434.

https://doi.org/10.1109/TSP.2008.929325

[10] Boashash, B. and Black, P.J. (1987) An Efficient Real -Time Implementation of the Wigner-Ville Distribution. IEEE Transactions on Acoustics, Speech and Signal Processing, ASSP-35, 1611-1618.

[11] Debnath, L. (2002) Recent Developments in the Wigner-Ville Distribution and Time-Frequency Signal Analysis. PINSA, 68, 35-36.

[12] Zhu, Q., Wang, Y.S. and Shen, G.Q. (2012) Research and Comparison of Time-Frequency Techniques for Nonstationary Signals. Journal of Computers, 7, 954-958.

https://doi.org/10.4304/jcp.7.4.954-958

[13] O’Toole, J., Mesbah, M. and Boashash, B. (2005) A Discrete Time and Frequency Wigner Ville Distribution Properties and Implementation,

http://eprints.qut.edu.au

[14] Gulum, T.O., Erdogan, A.Y., Durak-Ata, L., Yildirim, T. and Pace, E. (2013) Elliptic Gaussian Filtering for Time-Frequency Signal Analysis. 2013 IEEE Radar Conference, Ottawa, ON, 29 April-3 May 2013, 1-5.

https://doi.org/10.1109/RADAR.2013.6585973

[15] Mousa, A. and Saleem, R. (2010) The Shortage of WVD in Analyzing Abnormal ECG Signal. IEEE Symposium on Industrial Electronics and Applications (ISIEA 2010), Penang, 3-5 October 2010, 651-654.

https://doi.org/10.1109/ISIEA.2010.5679386