```python
import numpy as np
import scipy.signal as sl
import scipy.linalg as la
import pandas as pd
from autocorr import acf

from scipy.fftpack import fft, ifft, dct, idct

from matplotlib import pyplot as plt
from matplotlib.colors import to_rgba

%matplotlib inline

import warnings
warnings.simplefilter(action='ignore', category=FutureWarning)
```


```python
def figure(figsize=(10, 5)):
    plt.figure(figsize=figsize)
```


```python
c1 = '#2b2a2b'
c2 = '#863234'
c3 = '#9f847b'
```

# Correlation Function

\begin{align}
     C_{AB}(\tau) \equiv \langle A B \rangle(\tau) 
         = \lim_{t_0 \rightarrow \infty} ~ \frac{1}{t_0} \int_0^{t_0} \mathrm d t ~ 
             A (t + \tau) ~ B (t)
\end{align}


```python
def correlate(f1, f2, window=True):
    t0 = min(len(f1), len(f2))
    
    corr = np.zeros(t0)
    
    for t in range(t0):
        corr[t] = f1[t : t0] @ f2[: t0 - t] 
        # normalize
        corr[t] /= (t0 - t)
        
    if window:
        # apply window function
        wd = np.hanning(2 * t0)
        corr *= wd[t0:]
        
    return corr
```


```python
# Reference data: fast + slow oscillation

t_max = 20
dt = 0.01
w1 = 1 / 1
w2 = 1/ 13

A1 = 1
A2 = 0.33

t = np.arange(0, t_max+dt, dt)
N = len(t)
```


```python
# noise
rand = np.random.RandomState(0)

noise = (1 * rand.randn(len(t)) * ( 1 - np.exp(- t / 10)))
```


```python
f = A1 * np.cos(2*np.pi * w1 * t) + A2 * np.cos(2*np.pi * w2 * t) + noise

figure()
plt.plot(t, f)
```


```python
# mean and variance
print(f'{f.mean():.5f} +- {f.var():.5f}')
```

    -0.01435 +- 0.90553



```python
C = correlate(f, f, window=False)
C2 = np.correlate(f, f, 'full')[len(f)-1:] / len(f)

C2s = sl.correlate(f, f)[len(f)-1:]
# Normalize
C2s /= np.arange(len(f), 0, -1)
# window
# C2s *= np.hanning(2 * len(f))[len(f):]
C2s *= sl.windows.hann(2 * len(f))[len(f):]

# Yair's correlation function (includes proper normalization and windowing and)
C3 = acf(f[:, None], lag=len(f))

fig, (ax1, ax2, ax3, ax4, ax5) = plt.subplots(nrows=5, sharex=True)

ax1.plot(t, C)
ax2.plot(t, C2)
ax3.plot(t, C2s)
ax4.plot(t, C3)
ax5.plot(t, (C - C2))
ax5.plot(t, (C - C2s))
ax5.set_ylim(0.1 * np.array([-1, 1]))
```


```python
np.arange(len(f), 0, -1)
```




    array([2002, 2001, 2000, ...,    3,    2,    1])




```python
ref_C = correlate(f, f)

for c in (C2, C2s, C3):
    print(la.norm(ref_C - c) < 1e-12)
```

    False
    True
    True


# Spectral Density / Power Spectrum

\begin{align}
     S(\omega) 
         = \left\vert ~
         \int_{-\infty}^{\infty} \mathrm d t ~ f(t) ~ \mathrm e^{- \mathrm i \omega t} 
         ~ \right\vert^2
\end{align}


```python
k = np.linspace(0, 1 / dt, N)

S = abs(fft(f))**2
```


```python
tm = 200
plt.plot(k[:tm], S[:tm], color=c1)
plt.axvline(w1, color=c2, alpha=0.75)
plt.axvline(w2, color=c2, alpha=0.75)
```

# Wiener Khintchine Theorem

\begin{align}
     S(\omega) 
         = \int_{-\infty}^{\infty} \mathrm d \tau ~ 
             C(\tau) ~ \mathrm e^{- \mathrm i \omega \tau}
         = 2 \int_{0}^{\infty} \mathrm d \tau ~ 
             C(\tau) ~ \cos(\omega \tau)
\end{align}


```python
k = np.linspace(0, 1 / 2 / dt, N)
```


```python
SWK = 2 * dct(C2s)
```


```python
tm = 200
plt.plot(k[:tm], SWK[:tm], color=c1)
plt.axvline(w1, color=c2, alpha=0.75)
plt.axvline(w2, color=c2, alpha=0.75)
```

\begin{align}
     C(\tau) 
         = \frac{1}{2 \pi} \int_{-\infty}^{\infty} \mathrm d \omega ~ 
             S(\omega) ~ \cos(\omega \tau)
\end{align}


```python
C3 = idct(S)
C4 = idct(SWK)
```


```python
fig, (ax1, ax2) = plt.subplots(nrows=2, sharex=True)
ax1.plot(C3)
ax2.plot(C4)
```

# Telegraph Equation
\begin{align}
     X(t) = (-1)^{N(t)}~,
\end{align}

where $N(t)$ is the number of events generated by a Poisson process in the time interval $[0, t)$.

A Poisson process can be modelled as and infinitely dense series of coin tosses with a probability

$$ P(\text{Head}) = \lambda \mathrm d t $$
to flip a head.

The probability of tossing $n$ times a head within $[0, t)$ is then given by

$$ P(N=n, t) = \frac{(\lambda t)^n ~ \mathrm{e}^{-\lambda t}}{n!}~.$$

See https://www.probabilitycourse.com/chapter11/11_1_2_basic_concepts_of_the_poisson_process.php


```python
lam = 1

nt = 10000
ts = np.linspace(0, 500, nt)
dt = ts[1] - ts[0]
dt, lam * dt
```




    (0.05000500050005, 0.05000500050005)




```python
sample = np.zeros(nt)

rand = np.random.RandomState(1)
rands = rand.rand(nt)
sample[rands < lam*dt] = 1
sample.sum()
```




    499.0




```python
signal = (-1) ** sample.cumsum()
```


```python
plt.plot(ts, signal)
plt.xlim([0, 20])
```


```python
%%time
C_telegraph = correlate(signal, signal)
```

    CPU times: user 201 ms, sys: 308 µs, total: 202 ms
    Wall time: 50.5 ms



```python
%%time
C_telegraph = np.correlate(signal, signal, mode='full')[nt-1:] / nt
```

    CPU times: user 99.6 ms, sys: 0 ns, total: 99.6 ms
    Wall time: 24.8 ms



```python
%%time
C_telegraph = sl.correlate(signal, signal)[nt-1:] / nt
```

    CPU times: user 8.34 ms, sys: 0 ns, total: 8.34 ms
    Wall time: 1.84 ms



```python
fig, ax = plt.subplots()

ax.plot(ts, np.exp(-  2 * ts * lam), linewidth=3)
ax.plot(ts, C_telegraph, '.')
ax.set_xlim([0, 5 * 1 / lam])
ax.axvline(1 / lam)

ax.legend([r'$\exp(- \lambda \tau)$', r'$C (\tau)$', '1 / $\lambda$'])
ax.set_xlabel(r'$\tau$')

plt.show()
```

6bf28b7ef68b5ec7a8ee149b31f69bf6abd2989d
ipynb
Jupyter Notebook
mathematics/correlation_function/autocorrelation.ipynb
flokno/python_recipes
a09da65528ce3a2f1fd884aea361275b9aab0c15
0BSD
2020-05-06T14:38:16.000Z
2022-01-05T09:43:04.000Z
mathematics/correlation_function/autocorrelation.ipynb
flokno/python_recipes
a09da65528ce3a2f1fd884aea361275b9aab0c15
0BSD
mathematics/correlation_function/autocorrelation.ipynb
flokno/python_recipes
a09da65528ce3a2f1fd884aea361275b9aab0c15
0BSD
Qwen/Qwen-72B
1. YES
2. YES