% Copyright INRIA
%\documentstyle{book}
%\begin{document}
\def\verbatok#1{\expandafter \begin{verbatim}
\input{#1}}
\chapter{Time-Frequency representations of signals}
\label{wigner}
\index{time-frequency}
Numerous tools such as DFT, periodogram, maximum entropy method, have been
presented in previous chapters for the spectral\index{spectral}
analysis of stationary
processes. But, in many practical situations (speech, acoustics,
biomedicine applications, ...), this assumption of stationarity fails to
be true. Specific tools are then to be used when the spectral content of a
signal under consideration is time dependent. Two such tools will be
presented in this chapter: The $Wigner$\index{wigner}-$Ville$ representation
and
the classical $Short$-$Time$ $periodogram$ \index{short-time periodogram},
which are particuliar cases of a more general class of spectral estimators\cite{flandrin}. Nevertheless, due to the superiority of the so-called
$wigner$ $spectrum$, no numerical computation will be done for the
short-time periodogram.
\subsection{The Wigner distribution}
Let be given a discrete-time signal $f(n)$ and is Fourier transform\index{Fourier transform}:
\begin{eqnarray}
F(\nu)=\sum_{n\in {\Bbb Z}}f(n)e^{-jn\nu}
\label{e.wigner.1}
\end{eqnarray}
The Wigner distribution of $f(n)$ is given by:
\begin{eqnarray}
W_f(n,\nu)=2\sum_{k\in {\Bbb Z}}e^{-j2k\nu}f(n+k)f^{*}(n-k)
\label{e.wigner.2}
\end{eqnarray}
Useful in a time-frequency representation framework,
a similar expression may be written for $F(\nu)$:
\begin{eqnarray}
W_F(\nu,n)=\frac{1}{\pi}\int_{-\pi}^{\pi}e^{j2n\zeta}F(\nu+\zeta)F^{*}(\nu-\zeta)d\zeta
\label{e.wigner.3}
\end{eqnarray}
so that:
\begin{eqnarray}
W_F(\nu,n)=W_f(n,\nu)
\label{e.wigner.4}
\end{eqnarray}
illustrating the symmetry between the definitions in both time
and frequency domains.
Among the properties listed in \cite{claas} are the
$\pi$-periodicity and the hermitic property, the last leading to the fact
that the Wigner distribution of real signals is real. One more important
result is the following:
\begin{eqnarray}
\frac{1}{2\pi}\int_{-\pi/2}^{\pi/2}W_f(n,\nu)d\nu=|f(n)|^2
\label{e.wigner.5}
\end{eqnarray}
which means that the integral over a period of the Wigner distribution
is equal to the instantaneous signal power.
\subsection{Time-frequency spectral estimation}
All the properties stated in \cite {claas} for deterministic signals with
finite energy are still valid for random processes with harmonizable
covariance, that is
those for which the covariance function $K(s,t)$ has the decomposition:
\begin{eqnarray}
K(s,t)&=&E(X(s)X^{*}(t))\nonumber \\
&=&\frac{1}{4\pi^2}\int_{-\pi}^{\pi}\int_{-\pi}^{\pi}e^{i(\lambda s-\mu t)}f(\lambda , \mu)d\lambda d\mu
\label{e.wigner.6}
\end{eqnarray}
The same way an harmonic analysis of stationary random processes with the Fourier transform leads to the classical spectrum,
a mixed time-frequency analysis can be done with the Wigner distribution,
leading to a ``time-frequency spectrum'', which reduces to the first with the
stationary assumption.
In \cite{flandrin}, a general class of spectral estimators is considered:
\begin{eqnarray}
\hat{W}(n,\nu;\phi)=\frac{1}{\pi}\sum_{k\in {\Bbb Z}}\sum_{m\in {\Bbb Z}}\int_{-\pi}^{\pi}e^{ipm}\phi(p,2k)X(n+m+k)X^{*}(n+m-k)e^{-i2\nu k}dp
\label{e.wigner.7}
\end{eqnarray}
where $\phi(p,2k)$ is the Fourier transform of a data window caracterizing
the weighting on the products. A first choice
of $\phi(p,2k)$ is:
\begin{eqnarray}
\phi_{STP}(p,2k)=\frac{1}{2N-1}\sum_{j\in {\Bbb Z}}h_N(j+k)h_N^{*}(j-k)e^{ipj}
\label{e.wigner.8}
\end{eqnarray}
leading to the well-known $short$ $time$ $periodogram$:
\begin{eqnarray}
STP(n,\nu)=\frac{1}{2N-1}|\sum_{j\in {\Bbb Z}}X(k)h_N(k-n)e^{-i\nu n}|^2
\label{e.wigner.9}
\end{eqnarray}
which actually is the classical periodogram of the signal to which has
been applied a sliding window; smoothed versions may be obtained for this
estimator\cite{flandrin}
If now $\phi(p,2k)$ is chosen to be:
\begin{eqnarray}
\phi_{SPW}(p,2k)=|h_N(k)|^2\sum_{m\in {\Bbb Z}}g_M(m)e^{-ipm}
\label{e.wigner.10}
\end{eqnarray}
equation (\ref{e.wigner.6}) particularizes in the $smoothed$
$pseudo$-$wigner$ spectral estimator:
\begin{eqnarray}
PW(n,\nu)=2\sum_{k\in {\Bbb Z}}e^{-i2\nu k}|h_N(k)|^2\sum_{m\in {\Bbb Z}}g_M(m)X(n+m-k)X^{*}(n+m-k)
\label{e.wigner.11}
\end{eqnarray}
where $h_N(k)$ and $g_M(m)$ are windows with respective length
$2N-1$ and $2M-1$.
One major advantage in chosing
$\phi_{SPW}(n,2k)$ is that it is a separable function, so that independent
smoothing can be applied in the time and frequency directions, as opposed
to $\phi_{STP}(n,2k)$ which is governed by ``uncertainty relations''
between the time and frequency weightings. Thence,
the bias, known to be introduced by weighting functions in spectral
estimation, can be controlled separately for the pseudo-wigner estimator.
Morever, in the case of unsmoothed pseudo-wigner ($M=1$), the bias in the
time direction vanishes while always present for the short time periodogram.
Now, we can compute the wigner spectrum estimator:
let $x(n)$ denote the analytical signal of the sampled realization of
the process $X(t)$. Setting $M=1$ (unsmoothed estimator)
and $\nu_l\stackrel{\mbox{def}}{=}\pi(l/N)$,
we can write with the help of (\ref{e.wigner.10}):
\begin{eqnarray}
PW(n,\nu_l)&=&2\sum_{k=-N+1}^{N-1}e^{-i2k\pi(l/N)}|h_N(k)|^2x(n+k)x^{*}(n-k) \\
&=&2(2Re\{\sum_{k=0}^{N-1}e^{-i2k\pi(l/N)}|h_N(k)|^2x(n+k)x^{*}(n-k)\}-|x(n)|^2) \nonumber
\label{e.wigner.12}
\end{eqnarray}
and this expression may be easily computed via the FFT with the Scilab
function {\tt wigner}.
\paragraph{Example 1}: The following example is taken from \cite{flandrin}: the signal is a finite
duration sinusoid modulated by a parabola:
\[
S(t)=\left\{ \begin{array}{ll}
p(t)\sin(\frac{2\pi}{16}t+u(t-488)\pi) & 408 < t < 568 \\ \\
0 & 0 \le t \le 408 \; ; \; 567 \le t \le 951
\end{array}
\right.
\]
$p(t)$ is the parabola taking its maximum in $t=488$, $u(t)$ is the
unit step function, $h_N$ is the $64$-point rectangular window; the time
and frequency increments are respectively equal to $12$ and $\pi/128$; $M$
has been set to one (unsmoothed estimator). The signal generation and wigner
spectrum computation are then as follows:
\verbatok{\Diary wigner1.dia}
\end{verbatim}
%%%%%%%%%%%
\input{\Figdir wigner1.tex}
\dessin{{\tt exec('wigner1.code')} Wigner analysis. Sinusoid modulated by a parabola}{wigner.fig1}
%%%%%%%%%%
%\paragraph{Example 2}: For the purpose of easy visualization, this example
%presents the case of a linearly sweeped frequency sinusoid: the frequency
%is swept from $50$ $Hz$ to $55$ $Hz$ with a sampling frequency $200$ $Hz$:
%see figure \ref{wigner.fig2}.
%The results of a Wigner analysis appear on figure\ref{wigner.fig3}.
%%%%%%%%%%%
%\input{\Figdir wigner2.tex}
%\dessin{{\tt exec('wigner2.code')} Wigner analysis. Linear FM signal}{wigner.fi%g2}
%%%%%%%%%%
%\end{document}