next up previous [pdf]

Next: ESTIMATION OF TIME-VARYING AVERAGE Up: Liu etc.: Time-frequency analysis Previous: Introduction

Time-frequency analysis using local attributes

The Fourier transform has found various applications in signal analysis. The classic Fourier transform indicates the presence of different frequencies within the analysis window, but does not show where in that window the particular frequency components reside. Localized frequency information can be obtained by computing the Fourier transform with a temporally shifted window. Such an approach to time-frequency analysis is known as the STFT (Allen, 1977). The window function is commonly parameterized by window size, overlap, and taper. Once the window function has been chosen for the STFT, temporal and spectral resolutions are fixed for the entire time-frequency map.

The S-transform (Stockwell et al., 1996) is similar to the STFT, but with a Gaussian-shaped window whose width scales inversely with frequency. The expression of the S-transform is

\begin{displaymath}
S(\tau,f)=\int_{-\infty}^{\infty}S(t)\left\{\frac{\vert f\v...
...frac{-f^2(\tau-t)^2}{2}\right ] \exp(-2\pi ift)\right\} dt,
\end{displaymath} (1)

where $s(t)$ is a signal and $\tau$ is a parameter which controls the position of the Gaussian window. The S-transform is conceptually a hybrid of the STFT and wavelet analysis, containing elements of both but having its own properties. Like STFT, the S-transform uses a window to localize the complex Fourier sinusoid, but, unlike the STFT and analogously to the wavelet transform, the width of the window scales with frequency.

Consider a signal $f(x)$ on $[0,L]$. The Fourier series, assuming a periodic extension of the boundary conditions, can be expressed as

\begin{displaymath}
f(t)\approx a_{0}+\sum_{k=1}^{\infty}\left[ a_{k} \cos \lef...
...\right)+b_{k} \sin \left( \frac{2k\pi t}{L}\right)\right] ,
\end{displaymath} (2)

where $a_{k}$ and $b_{K}$ are the series coefficients. The relationship between $k$ and frequency $f$ is $k=Lf$. In the case of the discrete Fourier transform, frequency is finite. Letting $ \mathbf{\Psi}_{k}(t)$ represent the Fourier basis
\begin{displaymath}
\mathbf{\Psi}_{k}(t) = \left[\begin{array}{ll} \mathbf{\Psi...
... \sin \left( \frac{2k \pi t}{L}\right) \end{array} \right],
\end{displaymath} (3)

and $\mathbf{C}_k$ represent the series coefficients
\begin{displaymath}
\mathbf{C}_{k} = \left[a_{k} \quad b_{k}\right],
\end{displaymath} (4)

where $b_{0}=0$, equation 1 can be written as
\begin{displaymath}
f(t)=\sum_{k=0}^{\infty}\mathbf{C}_{k}\mathbf{\Psi}_{k}(t).
\end{displaymath} (5)

If frequency is finite, the range of $k$ becomes $[0,N]$, $N = k_{max} = Lf_{max}$. In linear notation, $\mathbf{C}_k$ can be obtained by solving the least-squares problem
\begin{displaymath}
\min_{\mathbf{C}_{k}}\left\vert \left\vert f(t)-\sum_{k=0}^...
...{C}_k \mathbf{\Psi}_{k}(t)\right\vert \right\vert _{2}^{2},
\end{displaymath} (6)

where $\left\vert \left\vert \right\vert \right\vert$denotes the squared $L-2$ norm of a function. By allowing coefficients $\mathbf{C}_k$ to change with time $t$, we can define the timevarying coefficients $\mathbf{C}_{k}(t)$ via the following least-squares problem
\begin{displaymath}
\min_{\mathbf{C}_{k}}\sum_{k=0}^{N}\left\vert \left\vert f(...
...{C}_{k}\mathbf{\Psi}_{k}(t)\right\vert\right\vert _{2}^{2}.
\end{displaymath} (7)

The Fourier coefficient $\mathbf{C}_{k}(t)$ in equation 7 is a function of time $t$ and frequency $f = k∕L$ and $\mathbf{C}_{k}(t)=\left[a_{k} \quad b_{k}\right]$. The numerical support of frequency f can be between zero and the Nyquist frequency (Cohen, 1995), and the interval of frequency can be $\Delta f = 1∕L $. In practical applications, the range of frequencies can also be assigned by the user. In matrix notation, equation 7 can be written as
\begin{displaymath}
\left[f(t) \quad f(t) ... \quad f(t)\right]^{T} \approx \le...
...\mathbf{C}_{1}(t)\quad...\quad\mathbf{C}_{N}(t)\right]^{T},
\end{displaymath} (8)

where $D\{...\}$ denotes a diagonal matrix which is composed from the elements of $ \mathbf{\Psi}_{k}(t)$.

The problem of minimization in equation 7 is mathematically ill-posed because it is severely underconstrained: There are more unknown variables than constraints. To solve this ill-posed problem, we force the coefficients $\mathbf{C}_{k}(t)$ to have a desired behavior, such as smoothness. With the addition of a regularization term, equation 7 becomes

\begin{displaymath}
\min_{\mathbf{C}_{k}(t)}\sum_{k=0}^{N}\left\vert\left\vert ...
...t\vert\right\vert _{2}^{2}+R\left[\mathbf{C}_{k}(t)\right],
\end{displaymath} (9)

where $R$ denotes the regularization operator. If we use classic Tikhonov regularization (Tikhonov, 1963), equation 9 can be written as
\begin{displaymath}
\min_{\mathbf{C}_{k}(t)}\sum_{k=0}^{N}\left\vert\left\vert ...
...ft[\mathbf{C}_{k}(t)\right]\right\vert\right\vert _{2}^{2},
\end{displaymath} (10)

where $\mathbf{D}$ is the Tikhonov regularization operator (roughening operator) and $\varepsilon $ is a scalar regularization parameter.

Shaping regularization (Fomel, 2007b) provides a particularly convenient method of enforcing smoothness in iterative optimization schemes. In the appendix, we review the method of shaping regularization in detail. Fomel (2009) used shaping regularization to constrain the problem of nonstationary regression. In this paper, we use shaping regularization, analogous to the problem of nonstationary regression, to constrain estimated coefficients. We choose our shaping operator to be Gaussian smoothing with an adjustable radius. In shaping regularization, the radius of the Gaussian smoothing operator controls the smoothness of coefficients $\mathbf{C}_{k}(t)$.

Once we obtain time-varying Fourier coefficients $\mathbf{C}_{k}(t)=\left[a_{k}(t) \quad b_{k}(t)\right]$, the time-frequency map is defined as

\begin{displaymath}
F(t,f=k/L)=\sqrt{a_{k}^{2}(t)+b_{k}^{2}(t)}.
\end{displaymath} (11)

It is also possible to do an invertible time-frequency transform, as shown by Liu and Liu and Fomel (2010).

Consider a simple signal (Figure 1a), which includes three monophonic components at 10, 20, and 30 Hz and two broad-band spikes at 2 and 2.3 s. Time-frequency maps generated by the S-transform and our method are shown, respectively, in Figure 1b and 1c. Frequency components appear to be represented well from low to high frequencies in the proposed method. In comparison with the S-transform, the proposed method provides superior resolution for monophonic waves. At the lower frequencies, the S-transform is as good as the proposed method on the spectral resolution for monophonic waves. Note that the S-transform represents spikes better at high frequencies. However, because the width of analysis window is large at low frequencies, the S-transform provides poor time resolution at low frequencies for spikes. When we use the shaping operator as a regularization term, there is an edge effect from smoothing, which appears near 0 s and 4 s (Figure 1c). Figure 1d displays the time-frequency map using a different parameter (30-point smoothing radius) to demonstrate adjustable time-frequency representation of the proposed method.

s-0 st-0 proj-0 proj-0-30
s-0,st-0,proj-0,proj-0-30
Figure 1.
(a) Synthetic signal with three constant frequency components and two spikes. (b) Time-frequency map of the S-transform. (c) Time-frequency map of the proposed method with smoothing radius of 15 points. (d) Time-frequency map of the proposed method with smoothing radius of 30 points.
[pdf] [pdf] [pdf] [pdf] [png] [png] [png] [png] [scons]

To show the effect of varying frequencies, we provide a composite chirp signal (Figure 2a), which includes two parabolic frequencies, each having a constant amplitude. Figure 2b and 2c shows the results of the S-transform and the proposed method with 10-point smoothing radius, respectively. The two frequency components are detected with higher time and frequency resolution by our method. The S-transform has high spectral resolution near low frequencies, but loses resolution at high frequencies. Resolution of the time-frequency map deteriorates in both methods when the curvature of the time-frequency curve becomes large (as indicated by the arrow).

s-1 st-1 proj-1
s-1,st-1,proj-1
Figure 2.
(a) Synthetic chirp signal with two parabolic frequency changes. (b) Time-frequency map of the S-transform. (c) Time-frequency map of the proposed method.
[pdf] [pdf] [pdf] [png] [png] [png] [scons]

In the proposed method, the smoothing radius used in shaping regularization is a parameter which controls the smoothness of the model. It is different from the method of the STFT and the S-transform, in which division into local windows is applied to the data to localize frequency content in time. In contrast to the wavelet-based methods (Chakraborty and Okaya, 1995; Wang, 2007; Sinha et al., 2005), our method, as a straightforward extension of the classic Fourier analysis, is different because of the choice of basis functions. The wavelet-based methods need to decide the wavelets, which can be regarded as patterns of seismic data. The closer the selected wavelets are to the true pattern of the input data, the better the result of the time-frequency analysis. In this paper, we use more general Fourier basis functions to compute the time-frequency map.


next up previous [pdf]

Next: ESTIMATION OF TIME-VARYING AVERAGE Up: Liu etc.: Time-frequency analysis Previous: Introduction

2013-07-26