% 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}