```python
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

import IPython
import scipy.io.wavfile as wav
import scipy.signal as ss
```

# Fonctions utilisees


```python
def synthetize(siglen, deltas, freqs, amps, phases):
    ts = np.arange(siglen).reshape(1, -1)
    pdeltas = np.minimum(deltas, 0)

    logzs = pdeltas + 1j * 2 * np.pi * freqs
    alphas = amps * np.exp(1j * phases)

    alphas = alphas.reshape(-1, 1).dot(np.ones((1, siglen)))
    zst = logzs.reshape(-1, 1).dot(ts)
    return (alphas * np.exp(zst)).T.sum(1)


def hankel(x):
    """ From the scipy.linalg libs """
    c = np.asarray(x).ravel()
    r = np.zeros_like(c)
    
    vals = np.concatenate((c, r[1:]))
    a, b = np.ogrid[0:len(c), 0:len(r)]

    indx = a + b
    return vals[indx]


def pinv(m):
    mH = m.T.conj()
    return np.linalg.inv(mH.dot(m)).dot(mH)


def lsquares(signal, deltas, freqs):
    alphas = pinv(mV).dot(signal)
    amps = np.absolute(alphas)
    phis = np.angle(alphas)
    return amps, phis


def music(signal, n, k):
    freqs = np.linspace(0, 1, 100)
    deltas = np.linspace(-.1, .1, 100)

    dd, ff = np.meshgrid(deltas, freqs)
    # zs = np.exp(deltas + 1j * 2 * np.pi * freqs)
    zz = np.exp(dd + 1j * 2 * np.pi * ff)

    vnz = np.array([np.vander(zs, 32, increasing=True).T for zs in zz])

    pseudospectr = 1 / (np.abs((mWorth.T.conj().dot(vnz))) ** 2).sum(0)
    return dd, ff, np.log(pseudospectr)
```

## Exponential Sinusoidal Model

Petit rappel du model:

$$ 
\begin{align}
s(t) & = \sum_k{a_k ~ e^{\delta_kt} ~ e^{i(2\pi f_k t + \phi_k)}} \\
& = \sum_k{\alpha_k ~ z_k^t},\quad  \alpha_k = a_k ~ e^{i\phi_k},~ z_k = e^{\delta_k + i2\pi f_k}\\
\end{align}
$$

---
## Synthese de signal


```python
siglen = 63
deltas = np.array([0, -.05])
freqs = np.array([1/4, 1/4 + 1 / siglen])
amps = np.array([1, 10])
# phases = np.random.rand(2)
phases = np.zeros(2)

synFrames = synthetize(siglen, deltas, freqs, amps, phases)
print(synFrames.shape)

plt.plot(np.real(synFrames))
plt.title("Synthesized signal")
plt.show()
```


```python
f, ps = ss.periodogram(np.real(synFrames), nfft=1024)
plt.plot(f, ps, label=1024)
f, ps = ss.periodogram(np.real(synFrames), nfft=63)
plt.plot(f, ps, label=63)
plt.title("periodogram for different values of nfft")
plt.legend()
plt.show()
```

---
## High-Resolution methods
### Analyse spectrale


```python
n, k = 32, 2

mX = hankel(synFrames)[:n, :(len(synFrames) - n + 1)]

mERxx = 1 / 126 * mX.dot(mX.T.conj())
u, s, v = np.linalg.svd(mERxx, full_matrices=True)
```


```python
plt.plot(s[:10])
plt.title("first 10 eigenvalues of $\hat R_{xx} $ for our synthesized signal")
plt.show()
```

---
### Methode ESPRIT 
On constate bien la fidelite de la reconstruction sur les deux figures qui suivent : les deux courbes sont parfaitement confondues, les valeurs approximees a un tres haut degre de precision.


```python
mW = u[:, :k]
mWorth = u[:, k:]
mWup = mW[1:, :]
mWdn = mW[:-1, :]

mWdnInv = pinv(mWdn)
mPhi = mWdnInv.dot(mWup)
```


```python
s = np.linalg.eigvals(mPhi)  # linalg.svd return real values
s = sorted(s, key=lambda x: np.abs(x), reverse=True)
recDeltas = np.log(np.absolute(s))
recFreqs = np.angle(s) / 2 / np.pi
poles = np.exp(deltas + 1j * 2 * np.pi * freqs)
```


```python
mV = np.vander(poles, siglen, increasing=True).T
```


```python
recAmps, recPhases = lsquares(synFrames, deltas, freqs)
recSynFrames = synthetize(siglen, recDeltas, recFreqs, recAmps, recPhases)

print(recAmps, amps)
print(recPhases, phases)
print(recDeltas, deltas)
print(recFreqs, freqs)

plt.plot(np.real(recSynFrames), label="rec")
plt.plot(np.real(synFrames), label="orig")
plt.legend()
plt.show()
```


```python
f, ps = ss.periodogram(np.real(synFrames), nfft=1024)
plt.plot(f, ps, label="orig")
f, ps = ss.periodogram(np.real(recSynFrames), nfft=1024)
plt.plot(f, ps, label="rec")
plt.title("periodograms for nfft=1024")
plt.legend()
plt.show()
```

---
### MUSIC method


```python
%matplotlib inline
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.set_xlabel("$\delta$")
ax.set_ylabel("frequences")
ax.plot_surface(*music(synFrames, n, k), cmap="coolwarm")
plt.show()
```


```python
fn = "ClocheA.WAV"
rate, frames = wav.read(fn)
print(rate, frames.shape)

plt.plot(np.arange(len(frames)), frames)
plt.title("amplitudes")
plt.show()
IPython.display.Audio(frames, rate=rate)
```

## Signaux audio


```python
plt.plot(*ss.periodogram(frames))
plt.title("periodogram")
plt.show()
```


```python
n, k, l = 512, 54, 1024
```

aad3b36507e56dbec00ec5247aa4d880864599b6
ipynb
Jupyter Notebook
tp3/main.ipynb
alexis-thual/audio-signal-analysis
226e654d7482869c0976199babebac1c67c3803a
MIT
tp3/main.ipynb
alexis-thual/audio-signal-analysis
226e654d7482869c0976199babebac1c67c3803a
MIT
tp3/main.ipynb
alexis-thual/audio-signal-analysis
226e654d7482869c0976199babebac1c67c3803a
MIT
Qwen/Qwen-72B
1. YES
2. YES