main.Rpres

Fourier and wavelet transforms
========================================================
author:Andres Garcia Saravia Ortiz de Montellano
date: 23.11.2015

***
<br/><br/><br/><br/><br/><br/><br/>
![](sage_logo.png)


Representation of functions
========================================================
incremental: true

Let $f:\mathbb{R}\rightarrow\mathbb{C}$

Maclaurin representation:

$$f(x) = a_0 + a_1 x + a_2 x^2 +\cdots $$

Fourier representation:

$$f(x) = a_0 + a_1 \cos{(w_1 x)} + \cdots + b_1 \sin{(w_1 x)} + \cdots$$

equivalently,

$$f(x) = a_0 + a_1 e^{i w_1 x} + a_2 e^{i w_2 x} + \cdots$$


Fourier representations
========================================================
incremental: true

The Fourier representation looks different depending on the domain of the function $f:A\rightarrow\mathbb{C}$, with $A$ being:
<br/><br/>

$\mathbb{R}$: Real numbers

$\mathbb{T}_p$: Circle of length $p$

$\mathbb{Z}$: Integer numbers

$\mathbb{P}_N$: Polygon of $N$ sides

Continuous Fourier Transform (CFT)
========================================================

For functions defined on the real line $\mathbb{R}$. 

$f:\mathbb{R}\rightarrow\mathbb{C}$

$$
  \begin{aligned}
    f(x) &= \int_{-\infty}^\infty F(s) e^{2\pi i s x} ds \\
    F(s) &= \int_{-\infty}^\infty f(x) e^{-2\pi i s x} dx \\
  \end{aligned}
$$

Fourier Series
========================================================

For functions defined on the integers on an interval $\mathbb{T}_p=[0,p)$. 

$g:\mathbb{T}_p\rightarrow\mathbb{C}$ with $g(p) = g(0)$

$$
  \begin{aligned}
    g(x) &= \sum_{k=-\infty}^\infty G(k) e^{2\pi ikx/p} \\
    G(k) &= \frac{1}{p}\int_{0}^p g(x) e^{-2\pi ikx/p} dx \\
  \end{aligned}
$$


Discrete Time Fourier Transform (DTFT)
========================================================

For functions defined on the integers $\mathbb{Z}$. 

$\phi:\mathbb{Z}\rightarrow\mathbb{C}$

$$
  \begin{aligned}
    \phi(n) &= \int_{0}^p \Phi(s) e^{2\pi isn/p} ds \\
    \Phi(s) &= \frac{1}{p}\sum_{n=-\infty}^\infty \phi(n) e^{-2\pi isn/p} \\
  \end{aligned}
$$


Discrete Fourier Transform (DFT)
========================================================

For functions defined on a *polygon* with $N$ vertices $\mathbb{P}_N = \{0,1,2,\ldots, N-1\}$. 

$\gamma:\mathbb{P}_N\rightarrow\mathbb{C}$ with $\gamma(N) = \gamma(0)$

$$
  \begin{aligned}
    \gamma(n) &= \sum_{k=0}^{N-1} \Gamma(k) e^{2\pi ikn/N} \\
    \Gamma(k) &= \frac{1}{N}\sum_{n=0}^{N-1} \gamma(n) e^{-2\pi ikn/N} \\
  \end{aligned}
$$

Parseval's identities
========================================================

$$
  \begin{aligned}
    \int_{-\infty}^\infty f(x)\overline{g(x)}dx 
      &= \int_{-\infty}^\infty F(s)\overline{G(x)}ds,\ \ \mathbb{R} \\
    \int_{0}^p f(x)\overline{g(x)}dx 
      &= p \sum_{k=-\infty}^\infty F(k)\overline{G(k)},\ \ \mathbb{T}_p \\
    \sum_{n=-\infty}^\infty f(n)\overline{g(n)}
      &= p \int_{0}^p F(s)\overline{G(x)}ds,\ \ \mathbb{Z} \\
    \sum_{n=0}^{N-1} f(n)\overline{g(n)} 
      &= N \sum_{k=0}^{N-1} F(k)\overline{G(k)},\ \ \mathbb{P}_N
  \end{aligned}
$$


Discretization and periodization
========================================================

A continuous function $f$ can be made discrete $\phi$ by *$h$-sampling*

$$\phi(n) = f(nh), n\in\mathbb{Z},\ h>0 $$

From a function $f$, under some conditions,  we can construct a periodic function $g$ by *$p$-summation*

$$g(x) = \sum_{m=-\infty}^\infty f(x-mp) $$


Poisson relations
========================================================
incremental: true

$p$-summation

$$g(x) = \sum_{m=-\infty}^\infty f(x-mp) 
\implies G(k) = \frac{1}{p}F\left(\frac{k}{p}\right)$$

$p/N$-sampling

$$\phi(n) = f\left(\frac{np}{N}\right)
\implies \Phi(s) = \sum_{m=-\infty}^\infty F\left(s - \frac{mN}{p}\right)$$


Poisson cube
========================================================

![](Poisson-cube.png)


Nyquist-Shannon sampling theorem
========================================================
incremental: true

Assume that $f(t)$ is *$\sigma$-bandlimited*

$$|\omega|>\sigma \implies F(\omega)=0$$

Sample it with an interval $\Delta t$

$$g(n) \equiv f(n\Delta t)$$

Then, $f(t)$ can be uniquely reconstructed only when

$$2\sigma \Delta t <1$$

$$\sigma <  \frac{1}{2\Delta t} \equiv \omega_{Nyq}$$

Example: single sinusoid
========================================================
incremental:true

$$
  \begin{aligned}
    f(t) &= a_1\sin(2\pi \omega_1 t) \\
    F(\omega) 
     &= \int_{-\infty}^{\infty}a_1\sin(2\pi \omega_1 t)e^{-2\pi i \omega t} dt\\
     &= a_1\int_{-\infty}^{\infty} 
      \left(\frac{e^{2\pi i\omega_1 t} - e^{-2\pi i\omega_1 t}}{2i}\right)
      e^{-2\pi i \omega t}dt\\
     &= \frac{a_1}{2i}\left[\delta(\omega-\omega_1) 
        - \delta(\omega+\omega_1)\right]
  \end{aligned}
$$

$$
  P(\omega) \equiv \left|F(\omega)\right|^2= \left(\frac{a_1}{2}\right)^2\left[\delta(\omega-\omega_1) 
        - \delta(\omega+\omega_1)\right]^2$$


Example: single sinusoid
========================================================

```{r, echo=FALSE}
N <- 1024
w1 <- 1
a1 <- 1
x <- seq(from=0, to=5.3,length.out = N)
y <- a1*sin(w1*x*2*pi)
plot(x*w1,y/a1,type="l",
     xlab = expression(paste("w"[1],"t")),
     ylab = expression("f(t)/a"[1]))
```

***
```{r, echo=FALSE}
Dt <- (max(x) - min(x)) / N
Dw <- 1 / (N * Dt)
Y <- fft(y)/N
P <- Mod(Y)^2
freqs <- (seq(from=0, to=N-1) * Dw)
plot(freqs / w1,P/a1^2, type="b", 
     xlim=c(-2,2),
     xlab = expression("w/w"[1]),
     ylab = expression("P/a1"^2))
lines(seq(from=-1, to=-N) * Dw / w1, rev(P/a1^2), type="b")
```


Example: gaussian noise
========================================================

```{r, echo=FALSE}
set.seed(1)
N <- 1024*8
x <- seq(from=0, to=1,length.out = N)
y <- rnorm(N,0,1)
plot(x,y,type="b",
     xlab = "t",
     ylab = "f(t)")
```

***

```{r, echo=FALSE}
Dt <- (max(x) - min(x)) / N
Dw <- 1 / (N * Dt)
Y <- fft(y)/N
P <- Mod(Y)^2
freqs <- (seq(from=0, to=N-1) * Dw)
plot(freqs, P*N,
     type="b", 
     xlim = c(0,N/2),
     xlab = "w",
     ylab = "P")
```


Example: Multiple sinusoids with noise
========================================================

```{r, echo=FALSE}
N <- 1024*4
w1 <- 10
w2 <- 1
a1 <- 2
a2 <- 3

noise <- 10
x <- seq(from=0, to=5,length.out = N)
y0 <- a1*sin(w1*x*2*pi) + a2*cos(w2*x*2*pi)
y  <- y0 + rnorm(N, 0, noise)

par(mfrow=c(2,1))
plot(x,y0,
     type = "l",
     xlab = "t", ylab="f(t)")
plot(x,y,
     type="l",
     xlab = "t", ylab = "f(t)")
par(mfrow=c(1,1))
```

***

```{r, echo=FALSE}
Dt <- (max(x) - min(x)) / N
Dw <- 1 / (N * Dt)
Y <- fft(y)/N
P <- Mod(Y)^2
freqs <- seq(from=0, to=N-1) * Dw
plot(freqs, P, 
     type="b",
     xlab = "w", ylab = "P(w)",
     xlim=c(0,12))
```

Example: Time-varying signal
========================================================

```{r, echo=FALSE}
N <- 1024*4
w1 <- 5
a1 <- 1
mu <- 5
sigma <- 0.5

x <- seq(from=0, to=10,length.out = N)
y <- a1*sin(w1*x*2*pi)
yv <- a1*sin(w1*x*2*pi)*dnorm(x, mu, sigma)
par(mfrow = c(2,1))
plot(x,y,type="l",
     xlab = "t",
     ylab = "f(t)")
plot(x,yv,type="l",
     xlab = "t",
     ylab = "f(t)")
par(mfrow = c(1,1))
```

***

```{r, echo=FALSE}
Dt <- (max(x) - min(x)) / N
Dw <- 1 / (N * Dt)
Y  <- fft(y)/N
Yv <- fft(yv)/N
P  <- Mod(Y)^2
Pv <- Mod(Yv)^2 
freqs <- (seq(from=0, to=N-1) * Dw)
par(mfrow = c(2,1))
plot(freqs / w1,P, type="b", 
     xlim=c(0,2),
     xlab = expression("w/w"[1]),
     ylab = "P")
plot(freqs / w1,Pv, type="b", 
     xlim=c(0,2),
     xlab = expression("w/w"[1]),
     ylab = "P")
par(mfrow = c(1,1))
```

Example: Time-varying signal 2
========================================================

```{r, echo=FALSE}
N <- 1024*4
w1 <- 5
w2 <- 0.7
w3 <- 1
a1 <- 15
a2 <- 3
a3 <- 1
mu <- 5
sigma <- 0.5
noise <- 3

x <- seq(from=0, to=10,length.out = N)

y0 <- a1*sin(w1*x*2*pi)*dnorm(x, mu, sigma)
y1 <- a2*sin(w2*x*2*pi) + a3*sin(w3*x*2*pi)
y <- y0 + y1 + rnorm(n = N, mean = 0, sd = noise)
par(mfrow = c(2,1))
plot(x,y1,type="l", ylim = c(-10, 10),
     xlab = "t",
     ylab = "f(t)")
lines(x,y0, col = "blue", lwd=2)
plot(x,y,type="l",
     xlab = "t",
     ylab = "f(t)")
par(mfrow = c(1,1))
```

***

```{r, echo=FALSE}
Dt <- (max(x) - min(x)) / N
Dw <- 1 / (N * Dt)
Y  <- fft(y)/N
P  <- Mod(Y)^2
freqs <- (seq(from=0, to=N-1) * Dw)
plot(freqs / w1,P, type="b", 
     xlim=c(0,1.2),
     xlab = expression("w/w"[1]),
     ylab = "P")
```


Short-time Fourier Transform (STFT)
========================================================
incremental: true

Make Fourier Transform of *short* segments in the timeseries using a *window* function centered at $\tau$, $w(t-\tau)$

$$\hat{F}(\omega, \tau) = 
\int_{-\infty}^\infty f(t) w(t-\tau)  e^{-2\pi i \omega t} dt$$

**Problem:** Choose an appropriate window width 

Frequency vs. time resolution


Wavelet transforms
========================================================
incremental: true

**Idea:** 

- Narrow windows for large frequencies (good time resolution)
- Wide windows for small frequencies (good frequency resolution)

Choose a *localized wave* (mother wavelet) $\psi(t)$ and define

$$\psi_{a,b}(t)=\frac{1}{\sqrt{a}}\psi\left(\frac{t-b}{a}\right)$$

Wavelet transforms
========================================================
incremental: true

Example: Mexican-hat wavelet

$$\psi(t) = \left(1-t^2\right)e^{-t^2/2}$$
$$\psi_{a,b}(t)=\frac{1}{\sqrt{a}}\psi\left(\frac{t-b}{a}\right)$$

***

```{r,echo=FALSE}
mhw <- function(t, a, b) {
  return((1/sqrt(a)) * (1-((t-b)/a)^2) * exp(-((t-b)/a)^2/2))
}
x <- seq(from=-5, to=5, length.out = 1000)
y0 <- mhw(x,1.2,0)
y1 <- mhw(x,0.2,0)
plot(x,y1, type="l", xlab = "t", ylab = expression(psi))
text(0.6,2, labels = "a=0.2")
lines(x,y0, col="blue")
text(1,0.8, labels = "a=1.2", col="blue")
```


Definition of wavelet transforms
========================================================
incremental: true

Given $f\in L^2(\mathbb{R})$, we define its *continuous wavelet transform* with respect to the wavelet $\psi$ as

$$\mathcal{W}_{\psi}[f] (a,b) = \int_{-\infty}^\infty f(t)\overline{\psi_{a,b}(t)}dt$$

For $\mathcal{W}_\psi[f]$ to be invertible we require 

$$ 0 < C_{\psi}\equiv\int_{-\infty}^\infty
\frac{\left|\hat{\psi}(\omega)\right|^2}{\left|\omega\right|} d\omega
<\infty $$

Example: Time-varying signal CWT
========================================================
incremental: true

```{r, echo=FALSE}
N <- 1024*256
w1 <- 5
w2 <- 0.7
w3 <- 1
a1 <- 15
a2 <- 3
a3 <- 1
mu <- 5
sigma <- 0.5
noise <- 3

x <- seq(from=0, to=10,length.out = N)

y0 <- a1*sin(w1*x*2*pi)*dnorm(x, mu, sigma)
y1 <- a2*sin(w2*x*2*pi) + a3*sin(w3*x*2*pi)
y <- y0 + y1 + rnorm(n = N, mean = 0, sd = noise)
par(mfrow = c(2,1))
plot(x,y,type="l",
     xlab = "t",
     ylab = "f(t)")
Dt <- (max(x) - min(x)) / N
Dw <- 1 / (N * Dt)
Y  <- fft(y)/N
P  <- Mod(Y)^2
freqs <- (seq(from=0, to=N-1) * Dw)
plot(freqs / w1,P, type="b", 
     xlim=c(0,1.2),
     xlab = expression("w/w"[1]),
     ylab = "P")
par(mfrow=c(1,1))
```

***

Mexican hat CWT

![](tvs1.png)

Example: Time-varying signal CWT 2
========================================================
incremental: true

Haar CWT

![](tvs2.png)

***

Morlet CWT

![](tvs3.png)


Example: Time-varying signal CWT 3
========================================================
incremental: true

Morlet

![](tvs4.png)

***

Mexican-hat

![](tvs5.png)


Parseval's relation for wavelets
========================================================

$$
\int_{-\infty}^\infty\int_{-\infty}^\infty \mathcal{W}_\psi f (a,b)\overline{\mathcal{W}_\psi g (a,b)} \frac{da db}{a^2} = C_{\psi}\langle f,g\rangle
$$

where

$$
C_\psi \equiv \int_{-\infty}^\infty \frac{|\hat{\psi}(\omega)|^2}{|\omega|}d\omega 
$$


Inverse of a wavelet transform
========================================================
incremental: true

$$f(t) = \frac{1}{C_\psi}\int_{-\infty}^\infty\int_{-\infty}^\infty \mathcal{W}_\psi[f](a,b) \psi_{a,b}(t) \frac{da\ db}{a^2}$$

only when

$$ 0 < C_{\psi}\equiv\int_{-\infty}^\infty
\frac{\left|\hat{\psi}(\omega)\right|^2}{\left|\omega\right|} d\omega
<\infty $$


Discrete wavelet transform
========================================================
incremental: true

Change the continuous version

$$\psi_{a,b}(t) = \frac{1}{\sqrt{a}} \psi\left(\frac{t-b}{a}\right),\ \ a,b\in\mathbb{R},\ a\neq 0 $$

to a discrete version

$$\psi_{m,n}(t)=2^{-m/2} \psi \left( 2^{-m} t - n \right),\ \ \ n,m \in \mathbb{Z}$$

When can we recover $f(t)$ from $\mathcal{W}_\psi[f](m,n)$ ?


Orthonormal basis
========================================================

Find $\psi_{m,n}$ such that they form a complete and orthonormal basis in $L^2(\mathbb{R})$:

$$f(t) = \sum_{m,n=-\infty}^\infty \langle f,\psi_{m,n}\rangle\ \psi_{m,n}(t)$$


Multiresolution analysis (MRA)
========================================================
incremental: true

![](MRA.png)

> *MRA is really an effective mathematical framework for hierarchical decomposition of an image (or signal) into components of different scales (or frequencies).*


Example application
========================================================

**Time-frequency analysis of solar $p$-modes**

*F. Baudin, A. Gabriel, D. Gibert (1994)*

***

![](Baudin.png)


Example application 2
========================================================

**Wavelets: a powerful tool for studying rotation, activity, and pulsation in Kepler and CoRoT stellar light curves**

*J. P. Bravo, S. Roque, R. Estrela, I. C. Leão, and J. R. De Medeiros (2014)*

***

![](kprl-example.png)


Thank you!
========================================================

![](coffee_is_essential.jpg)

***

[https://github.molgen.mpg.de/saravia/wavelets-SAGE](https://github.molgen.mpg.de/saravia/wavelets-SAGE)

saravia@mps.mpg.de

ags3006@gmail.com

Orthogonality relations
========================================================

$$\int_{-\infty}^\infty e^{2\pi i (x-x')}dx' = \delta(x)$$

$$\int_0^p e^{2\pi i x (k-l)/p}dx =
  \begin{cases} 
    p &\mbox{if}\ \ k=l \\
    0 &\mbox{otherwise}
  \end{cases}
$$

$$\sum_{n=0}^{N-1} e^{2\pi i n (k-l)/N} =
  \begin{cases}
    N &\mbox{if}\ \ k = mN\ ,m\in\mathbb{Z} \\
    0 &\mbox{otherwise}
  \end{cases}
$$