\documentclass[12pt]{article}
\usepackage{wsh}
\usepackage[pdftex]{graphicx} 

\newcommand{\sbar}{\bar{s}}
\newcommand{\tbar}{\bar{t}}
\newcommand{\half}{\frac{1}{2}}

\begin{document}  
\bibliographystyle{plain}

\markright{Notes on continuous wavelet transforms --- W.S. Harlan} 
\title{Notes on continuous wavelet transforms} 

\author{William S. Harlan}

\date{Originally November 1996 -- updated October 27, 2017}

\maketitle

\section {Introduction}

First I will provide my own favorite definition of a continuous
wavelet transform, which I wrote down in 1986 to rederive the results in 
Goupillaud et al \cite{GoupillaudOne}.  (At the time, I was working
down the hall from Pierre Goupillaud.)
This continuous transform seemed much better suited 
to geophysical and acoustic applications than certain discrete 
transforms popular at the time, particularly Daubechies wavelets.  
Many wavelets have a broader spectrum than necessary, more suitable for compressing
photographic images with sharp edges.   A continuous transform
can use simple Gaussian tapered monochromatic waves (one or multi-dimensional)
which maximizes the locality in both time/space and frequency/wavenumber.  

First, I will give the definition of what constitutes a
valid ``wavelet'' for the transform.   Second, the inverse
of a wavelet transform will be derived to demonstrate
why the transform works.  Next, I will show how a wavelet-transform
of a stretched function can be written as a simple function
of the wavelet-transform of the original unstretched function.
Such stretching is very useful in pseudo-spectral solutions of differential
equations with spatially varying coefficients.

Next I dwell on the original Morlet wavelet, a Gaussian-tapered sinusoid,
to show how resolution can be managed explicitly in both time and frequency
domains.

(This is an expansion of a much older document \cite{harlan-wavelets} that stressed the
application of wavelet transforms to seismic imaging.  I have removed the geophysical
details and have extended the discussion.)

\section {Legal wavelets}

I will closely imitate the continuous transform originally
provided by Goupillaud, Grossman, and Morlet 
\cite{GoupillaudOne}.  They provide forward
and inverse transforms for related transforms, but
without derivation.

Use the following Fourier convention to link a 
wavelet $w(t)$ and its Fourier transform $\tilde w (s)$:
\begin{eqnarray}
\label{eq:fourier}
\tilde w (s) &=& \int e^{-i 2 \pi s t }  w(t) dt {\rm ~and} \\
w(t) &=& \int e^{i 2 \pi s t }  \tilde w (s) ds . 
\end{eqnarray}
(All integrals with unlabeled extremes are assumed to be
over the entire real line.)

The following integral must exist for a wavelet to be suitable for a continuous transform
over all frequencies.
\begin{equation}
C_w = \int { 1 \over |s| } \tilde w (s) ds .
\end{equation}
Thus, the spectrum of a valid wavelet must approach a zero
value at zero frequency: $\tilde w (s) / |s|^\epsilon $ $\rightarrow 0 $
as $s \rightarrow 0$ for some $\epsilon > 0 $.

In fact many useful wavelets (such as the Morlet wavelet) will be non-zero
at zero frequency and will not be valid for an infinite range
of frequencies.   Discrete transforms will prefer damped least-squares
optimizations rather than simple integrals for transformations 
and will not require an analytic constant for normalization.

This wavelet need not correspond to any specific feature in our data.
Rather, when convolved with data, the wavelet should
suppress all but a band of frequencies.  Because of the
uncertainty principle, I must balance the narrowness
of the bandwidth with the narrowness in time.

I recommend a Gaussian-tapered sinusoid, such as
$w(t) = {\rm cos}(2\pi t) (e^{- \pi t^2 / c^2} )$. 
Such a wavelet is reasonably compact in both
the time and frequency domain. (The constant $c$
adjusts the relative compactness, within
the limits of the uncertainty principle.)

A valid wavelet can be used in an unfamiliar construction
of a delta function.
\begin{equation}
\label{eq:delta}
\int w(ut) du = C_w \cdot \delta (t) .
\end{equation}
In effect, this integral says that many wavelets, stretched
uniformly over all scales, will sum destructively everywhere but
at $t=0$.
\begin{equation}
{\bf Proof:~}  \int w(ut) du 
= \int \left \{ \int  {1 \over |u|} \tilde w ( { s \over u } ) 
e^{i 2 \pi s t } ds \right \} du \ldots
\end{equation}
Change variables from $(s,u)$ to $(s'=s,u'=s/u)$.
The Jacobian is \\
$\partial (s,u) / \partial (s',u') = |s'/u'^2|$.
Thus,
\begin{eqnarray} 
\int w(ut) du 
&=& \int \int  \left | {u' \over s'} \right | \tilde w ( u' ) e^{i 2 \pi s' t } 
\left | { s' \over u'^2 } \right | ds' du' \\
&=& \int \int  {1 \over |u'|} \tilde w ( u' ) e^{i 2 \pi s' t } ds' du' \\
&=& \left \{ \int {1 \over |u'|} \tilde w ( u' ) du'  \right \}  \cdot 
\left \{ \int e^{i 2 \pi s' t } ds' \right \} \\
&=& C_w \cdot \delta (t) . 
\end{eqnarray}
Here, I use the familiar Fourier definition of a delta function,
with the same restrictions.

The constraint of symmetry will also simplify the definition
of a wavelet transform in the following section, but
this assumption is optional:
\begin{equation}
w(t) = w(-t) .
\end{equation}

\section {Forward and inverse transforms}

Define a wavelet transform $F(u,a)$ of an $L^2$ function $f(t)$ by 
the following 
\begin{equation}
\label{eq:defn}
F(u,a) = \int w[u(a - t)] f(t) dt .
\end{equation}
This transform decomposes the function as a function of
position $a$ and local frequency $u$.  When calculated
discretely, the sampling need not be uniform over $a$ and $u$,
but the weights should reflect the above integration.

Goupillaud et all 
\cite{GoupillaudOne} prefer $2^u$ or $u^{-1}$ instead of $u$ 
as the stretch factor.
If you prefer an asymmetric wavelet, then you must use
two transforms---the above, and the integral with the time-reversed
wavelet.)

I find the following inverse, which simply scales and sums the transform
with uniform weight over local frequency $u$:
\begin{equation}
\label{eq:inverse}
f(t) = C_w^{-1} \int F(u, a=t) du .
\end{equation}
\begin{eqnarray}
{\bf Proof:~} 
C_w^{-1} \int F(u,a=t) du 
&=& C_w^{-1} \int \left \{ \int w[u(t - t' )] f(t') dt' \right \} du  \\
&=& C_w^{-1} \int \left \{ \int w[u(t - t' )] du \right \} f(t') dt' \\
&=& C_w^{-1} \int C_w \delta(t-t') f(t') dt' \\
&=& f(t). 
\end{eqnarray}
Here I make use of equation~(\ref{eq:delta}).

Most all these equations were implicit in Goupillaud et al 
\cite{GoupillaudOne},
although I personally found it nontrivial to fill in the 
missing derivation of the results.  I was also happy 
to stumble on equation (\ref{eq:delta}), which continues
to fascinate me.

\section {Stretching}

One application that frequently appears is a differential stretching of a function.

We will write
\begin{equation}
g(\tau) = f[t(\tau)].
\end{equation}
The stretching function $t(\tau)$ gives the input coordinate $t$ 
as a monotonically increasing
function of the output coordinate $\tau$.  
Moreover, I will assume that the stretching
function can be well approximated locally by a straight line.
\begin{eqnarray}
\label{eq:approx}
t(\tau) &\approx& t(\tau_0 ) + (\tau - \tau_0 ) 
\cdot {dt \over d\tau} ( \tau_0 ) ~~~~{\rm and }\\
\tau(t) &\approx& \tau_0 + [ t - t (\tau_0 ) ] 
\cdot \left [ {dt \over d\tau} ( \tau_0 ) \right ]^{-1} . 
\end{eqnarray}
Define forward and inverse wavelet transforms 
of $g(\tau)$ with the same wavelets as in definition~(\ref{eq:defn})
and inverse~(\ref{eq:inverse}).
\begin{equation}
G(v, b ) = \int w [ v ( b  - \tau ) ]  g (\tau) d\tau  ~~~~ {\rm and}  
\end{equation}
\begin{equation}
g(\tau) = C_w^{-1} \int G(v, b  = \tau ) dv.
\end{equation}
We would like to be able to express the transform $G(v, b)$
as a simple function of $F(u,a)$.  The approximation~(\ref{eq:approx}) 
allows us to write
\begin{equation}
\label{eq:stretch}
G(v,b) \approx 
\left | {dt \over d\tau} ( b ) \right |^{-1} 
 ~ F \left \{ u=v \cdot  
\left [ {dt \over d\tau} ( b ) \right ]^{-1} , ~ a = t(b) \right \} .
\end{equation}
We can perform an equivalent stretch on the transformed
input by stretching the position $a$, and by scaling
the local frequency $u$.  Both operations require a simple 
two-dimensional mapping of the transformed function.
The amplitude must be scaled as well.
\begin{eqnarray} 
{\bf Proof: ~} 
G ( v, b) &=& \int w[v(b-\tau)] \cdot f[t(\tau)] d \tau \\
&\approx& \int w [ v(b-\tau)] \cdot f \left [ t(b) + (\tau-b) 
\cdot {dt \over d\tau} ( b ) \right ]  d\tau \ldots 
\end{eqnarray}
Substitute the variable of integration:
\begin{eqnarray}
t' & \equiv& t(b) + (\tau-b) \cdot  {dt \over d\tau} ( b ) 
 ~~~~ {\rm and} \\
\tau &=& b + [ t' - t(b) ] 
\cdot \left [ {dt \over d\tau} ( b ) \right ]^{-1} .
\end{eqnarray}
Thus,
\begin{eqnarray}
G ( v, b) &=& \int w \left \{ v \cdot [t(b) - t'] \cdot
\left [ {dt \over d\tau} ( b ) \right ]^{-1} \right \} \cdot f(t') 
\cdot \left | {dt \over d\tau} ( b ) \right |^{-1} ~ dt'  \\
&=& \left | {dt \over d\tau} ( b ) \right |^{-1} 
\int w \left \{ v 
\left [ {dt \over d\tau} ( b ) \right ]^{-1} 
\cdot [t(b) - t'] 
\right \}  \cdot f(t') d t'  \\
&=& \left | {dt \over d\tau} ( b ) \right |^{-1} 
 ~ F \left \{ u=v \cdot  
 \left [ {dt \over d\tau} ( b ) \right ]^{-1} , ~ a = t(b) \right \} . 
\end{eqnarray}
The symmetry of $w(t)$ was not necessary here.
     
Equation~(\ref{eq:stretch}) is the result which allows us to perform
the stretch directly on the wavelet-transformed function.

\section{Morlet wavelets}

A preferred standard wavelet is the Morlet wavelet
\cite{GoupillaudOne}, constructed from  a symmetric cosine tapered by a Gaussian.
\begin{eqnarray}
m(t) &\equiv& m_g(t) m_c(t) , \\
\mbox{ where } 
m_g(t) &\equiv& e^{ - \pi (\sbar t / h )^2 } \\
\mbox{ and }
m_c(t) &\equiv& \cos ( 2 \pi \sbar t) ,  \\
\mbox{ so }
m(t) &=& e^{ - \pi (\sbar t / h )^2 } \cos ( 2 \pi \sbar t) .
\end{eqnarray}
The constants $\sbar$ and $h$ do not vary during a wavelet transform.
If we take the dimension $t$ as units of time in seconds,
then $\sbar$ is the central frequency of the wavelet in hertz.  
The ``half-width'' $h$ controls the relative width of the taper in wavelengths.

The cosine has a base frequency of $\sbar$ and a wavelength of $1/\sbar$.
The inclusion of a variable $\sbar$ for central frequency is just a convenience for 
managing units;  $\sbar$ can be set to a constant without loss of generality.
When we scan frequencies with a wavelet transform~ (\ref{eq:defn}), 
we can set $\sbar$ to 1 and substitute $u$ as frequency.

The envelope has a peak value of 1 at $t=0$ and drops to approximately half of this value
at a separation of $h$ wavelengths.
\begin{eqnarray}
m_g(0) &=& 1 ; \\
m_g[h/(2 \sbar)] &=& m_g[-h/(2 \sbar)] \\
&=& \exp (- \pi / 4) \approx 0.456 .
\end{eqnarray}
And so we call $h$ the half-width, in wavelengths, of the Gaussian taper.

The integral of the Gaussian envelope in time $t$ is
\begin{eqnarray}
\int m_g(t) dt = \int e^{- \pi (\sbar t / h )^2 } dt = \frac{h}{\sbar} ,
\end{eqnarray}
equal to the half-width in units of time.

As an example, figure~\ref{fig:mt} shows one Morlet wavelet. 
We assume a time sampling interval of $1$ for a maximum Nyquist frequency of $0.5$.
We choose a central frequency of $\sbar = 0.25$, and a halfwidth in wavelengths of $h=4$.   
We can see that one wavelength is $1/\sbar = 4$ time samples in length.
The positive lobes at times of $t=\pm 8$ (third most positive values)
are half the height of the peak.  These half-strength peaks are separated by 
four wavelengths, as appropriate for $h=4$.  The two peaks are also 16 time samples apart, as appropriate
for four wavelengths of 4 samples each ($h/\sbar = 4/0.25 = 16$).

\begin{figure}[t]
\begin{center}
\includegraphics[width=130mm]{fig_mt}
\end{center}
\caption{ Morlet wavelet with central frequency of $\sbar = 0.25$ and half-width $h=4$.}
\label{fig:mt}
\end{figure}

We can take the Fourier transform~(\ref{eq:fourier}) of the Gaussian envelope as
\begin{eqnarray}
\tilde m_g(s) &\equiv& \int e^{-i 2 \pi s t } m_g(t) dt =
\int e^{-i 2 \pi s t } e^{ - \pi (\sbar t / h )^2 }  dt \\
\label{eq:fe}
&=& \frac{h}{\sbar} e^{ - \pi ( h s /\sbar )^2 } .
\end{eqnarray}
The equivalent half-width of this Gaussian is $\sbar / h$,
approximately spanning the range of values more than half of the peak value.
This half-width $\sbar / h$ in units of frequency is just the reciprocal of the
half-width $h/\sbar$ in units of time.


Notice that the frequency half-width is directly proportional to the reference frequency $\sbar$.
As we scan frequencies with a wavelet transform, we will see lower frequencies resolved
with narrower ranges of frequencies.  The resolution is constant in the log of frequency,
which is why Morlet recommends sampling frequencies geometrically in octaves rather than linearly
in frequency.  When we construct a discrete basis of these wavelets for a discrete transform,
there is no need to sample frequencies much more densely than this half-width.

The Fourier transform of the cosine is
\begin{eqnarray}
\tilde m_c(s) &\equiv& \int e^{-i 2 \pi s t } m_c(t)  dt
= \int e^{-i 2 \pi s t } \cos (2 \pi \sbar t) dt \\
&=& \half \delta (s+\sbar)  + \half \delta(s - \sbar).
\end{eqnarray}
The Fourier transform of the Morlet wavelet is the convolution of these two functions, for
\begin{eqnarray}
\tilde m(s) &\equiv&  
\int e^{-i 2 \pi s t }  m (t) dt =
\int e^{-i 2 \pi s t }  m_g (t) m_c (t) dt \\
&=& \tilde m_g(s) * \tilde m_c(s) = \int \tilde m_g(s') \tilde m_c(s -s') d s' \\
\label{eq:fm}
&=& 
\frac{h}{2 \sbar} e^{- \pi [ h ( s + \sbar) /\sbar ] ^2 }  +
\frac{h}{2 \sbar} e^{- \pi [ h ( s - \sbar) /\sbar ] ^2 }.
\end{eqnarray}
This gives us two Gaussians centered on the central frequency and mirrored on each side of zero.

Figure~\ref{fig:mf} shows the Fourier transform of the Morlet wavelet in figure~\ref{fig:mt},
also with $\sbar = 0.25$ and $h=4$.

\begin{figure}[t]
\begin{center}
\includegraphics[width=130mm]{fig_mf}
\end{center}
\caption{ Fourier transform of Morlet wavelet in figure~\ref{fig:mt} with central frequency of $\sbar = 0.25$ and half-width $h=4$.}
\label{fig:mf}
\end{figure}

The resolution we care about is the narrowness of the spectrum about the central frequency $\sbar$.
If we measure the standard deviation of one side of the Fourier transform~(\ref{eq:fm}),
we are effectively measuring the width of the Fourier transformed envelope~(\ref{eq:fe}):

For quantifying resolution tradeoffs, a useful measure is the normalized second moment, or dispersion, in each domain.
In time, the integral of the squared Gaussian taper is 
\begin{eqnarray}
\int | m_g(t) |^2  dt = 
\int e^{- 2 \pi (\sbar t / h )^2 } dt = \frac{1}{\sqrt{2}} \frac{h}{\sbar} ~,
\end{eqnarray}
and the second moment is
\begin{eqnarray}
\int | m_g(t) |^2  t^2 dt = 
\int e^ { - 2 \pi (\sbar t / h )^2 } t^2 dt = \frac{1}{4 \sqrt{2} \pi}\frac{h^3}{\sbar^3} ~.
\end{eqnarray}
So the normalized second moment is
\begin{eqnarray}
\Delta ^ 2 t  &=& 
{ \int | m_g(t) |^2  t^2 dt \over
\int | m_g(t) |^2  dt  }
= \frac{1}{4 \pi} \frac{h^2}{\sbar^2} ;  \\
\Delta t  &=& \frac{1}{2 \sqrt{\pi}} \frac{h}{\sbar} ~.
\end{eqnarray}

Similarly we can measure second moments in the frequency domain.
\begin{eqnarray}
\Delta ^ 2 s  &=&
{\int | \tilde m_g(s) |^2  s^2 ds \over
\int | \tilde m_g(s) |^2  ds }  =
\frac{1}{4 \pi}\frac{\sbar^2}{h^2} ;  \\
\Delta s  &=& \frac{1}{2 \sqrt{\pi}}\frac{\sbar}{h} ~.
\end{eqnarray}

The uncertainty relation for resolution in time and frequency simplifies to
\begin{eqnarray}
\Delta s \Delta t = \frac{1}{4 \pi}.
\end{eqnarray}
This inverse relationship is independent of both the half-width $h$ and the central frequency $\sbar$.
The half-width in wavelengths $h$ controls the width of resolution in each domain,
but cannot change this product of two widths.
Improving the resolution in one domain necessarily reduces the resolution in the other.
If we double the width of the envelope in the time domain, then we necessarily halve
the width in the frequency domain.

\bibliography{wsh}
\end{document}
