<jupyter_start><jupyter_text>Fonctions utilisees<jupyter_code>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)<jupyter_output><empty_output><jupyter_text>Exponential Sinusoidal ModelPetit 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<jupyter_code>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()
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()<jupyter_output><empty_output><jupyter_text>--- High-Resolution methods Analyse spectrale<jupyter_code>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)
plt.plot(s[:10])
plt.title("first 10 eigenvalues of $\hat R_{xx} $ for our synthesized signal")
plt.show()<jupyter_output><empty_output><jupyter_text>--- 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.<jupyter_code>mW = u[:, :k]
mWorth = u[:, k:]
mWup = mW[1:, :]
mWdn = mW[:-1, :]

mWdnInv = pinv(mWdn)
mPhi = mWdnInv.dot(mWup)
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)
mV = np.vander(poles, siglen, increasing=True).T
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()
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()<jupyter_output><empty_output><jupyter_text>--- MUSIC method<jupyter_code>%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()
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)<jupyter_output>8000 (95997,)<jupyter_text>Signaux audio<jupyter_code>plt.plot(*ss.periodogram(frames))
plt.title("periodogram")
plt.show()
n, k, l = 512, 54, 1024<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1802247191011236, 'fraction_numerical': 0.025617977528089888, 'mean_word_length': 4.591708542713568, 'pattern_counts': {'":': 0, '<': 24, '<?xml version=': 0, '>': 24, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '6', 'dir': 'jupyter-structured-clean-dedup', 'id': '386952', 'license': 'MIT', 'n_tokens_mistral': 2091, 'n_tokens_neox': 1921, 'n_words': 379, 'path': 'tp3/main.ipynb', 'repo_name': 'alexis-thual/audio-signal-analysis', 'size': '4450'}