Protokol k projektu ISS 2025/26: Music Search

Autor: Ľubor Pačaj (xpacajl00)

1 Úvod

V tomto projekte bolo cieľom vytvoriť jednoduchú aplikáciu na vyhľadávanie hudby (Music Information Retrieval, MIR). Pre vypracovanie som si zvolil Python Jupyter notebook, nakoľko sa tu nachádzala kostra programu a veľa python knižníc, ktoré mi prácu výrazne uľahčili. Využíval som rady zo zadania a nakoniec sa mi podarilo vytvoriť dostatočne funkčnú verziu.

Pomocné funkcie zo zadania:

# potrebne
import os
import re
import soundfile as sf
from IPython.display import Audio
from IPython.display import display
import numpy as np
# odporucane ...
from scipy import signal
from scipy.io import wavfile
from scipy.fft import fft, ifft, fftfreq
import scipy.io
import matplotlib.pyplot as plt
from matplotlib import cm

# nacitanie dat
login = "xpacajl00"

for name in ("known.zip", "valid.zip", login + ".zip"):
  folder = name.replace(".zip", "")
  if not os.path.exists(folder):
    file = "https://www.fit.vut.cz/study/course/ISS/public/proj2025-26/" + name
    !wget -q $file
    !unzip -o -q $name
    !rm $name

# nacitanie znamych dat - funkcia vracia velku maticu so vsetkymi signalmi
def load_data (S, dirname, count, no_samples):
  ii = 0
  for one in np.arange(count):
    S[ii], Fs = sf.read(dirname + "/" + str(one) + ".wav")
    ii = ii+1

Fs = 16000
# nacitanie znamych dat
N_known = 706; duration_known = 10; no_samples_known = Fs * duration_known
known_signals=np.zeros([N_known, no_samples_known]); load_data(known_signals, "known", N_known, no_samples_known)

# nacitanie validacnych dat
N_valid = 50; duration_valid = 5; no_samples_valid = Fs * duration_valid
valid_signals=np.zeros([N_valid, no_samples_valid]); load_data(valid_signals, "valid", N_valid, no_samples_valid)

# nacitanie xlogin dat 
N_eval = 50; duration_eval = 5; no_samples_eval = Fs * duration_eval
eval_signals=np.zeros([N_eval, no_samples_eval]); load_data(eval_signals, login, N_eval, no_samples_eval)

2 Popis riešenia

2.1 Začiatok práce

Pri premýšľaní a analyzovaní zadania projektu som si vypočul niekoľko zvukových nahrávok. Hneď 1. z nahrávok v priečinku s mojím loginom sa mi zapáčila a pri neskoršom hľadaní sa mi podarilo ju nájsť medzi známymi nahrávkami. Použil som ich teda ako testovacie nahrávky pri vykresľovaní a vyhodnocovaní mojich prístupov.

2.2 Vykresľovanie spektrogramov

Na začiatok som vytvoril jednoduchú funkciu ktorá vykreslí spektrogram z nahrávky. Aby som si ušetril prácu s implementáciou využil som funkciu scipy.spectrogram. Táto funkcia vykoná Rýchlu Fourierovu trasnformáciu (FFT) a vráti tri hodnoty (frekvenciu, čas a amplitúdovú zložku).

$$X[k] = \sum_{n=0}^{N-1} x[n] \cdot e^{-j 2\pi \frac{k}{N} n}$$ $$|X[k]| = \sqrt{\text{Re}(X[k])^2 + \text{Im}(X[k])^2}$$

Pri tvorbe spektrogramu som sa rozhodol využiť Hanningovo okno, ktoré podľa mojich testovaní poskytlo najlepšie výsledky. Pri parametroch som musel zvoliť kompromis medzi časovým a frekvenčným rozlíšením. Ako parametre som zvolil nasledovné:

• Dĺžku rámca(nperseg) som zvolil 512. Pri vzrokovacej frekvencii 16kHz zodpovedá 512 vzoriek výslednému frekvenčnému rozlíšeniu 31.25 Hz. [Fs/N = 31.25 Hz]. Toto frekvenčné rozlíšenie mi po odskúšaní ponúka ideálny kompromis pre identifikáciu tónov v čase. Pri dlhších rámcoch sa mi síce frekvenčné detaily zvýšili ale nebolo možné tak presne určiť časovú lokalizáciu, teda spektrogram bol viac rozmazaný.

• Prekyv(noverlap) som vždy volil analogicky ku nperseg ako 50% dĺžky rámca [N/2 = 256], čo je odporúčaná hodnota pri použití Hanningovho okna, aby sa vyladili prechody medzi rámcami.

Nakoniec som ich zlogaritmoval na dB a vykreslil. $$X_{\text{dB}}[k] = 10 \cdot \log_{10} (|X[k]|)$$

Nižšie sú ukážky spektrogramov dvoch nahrávok:

• known_signals[114], ktorá obsahuje originálnu nahrávku z databázy známych nahrávok.
• eval_signals[0], ktorá je 1. z nahrávok na analýzu z priečinka môjho loginu. Jedná sa o 5 sekundovú orezanú a zašumenú verziu known_signals[114].

Ďalej v protokole na nich budem odkazovať týmito názvami.

# Funkcia na ukazanie spektrogramu danej nahravky
def show_spectrogram(arr, idx):
    # Priradi pole a index
    if arr == "known_signals":
        sig = known_signals[idx]
    elif arr == "valid_signals":
        sig = valid_signals[idx]
    elif arr == "eval_signals":
        sig = eval_signals[idx]
    else:
        return
    name = f"{arr}[{idx}]"
    # Vrati magnitudu
    f, t, mag = signal.spectrogram(sig, fs=Fs, window="hann", nperseg=512, noverlap=256)
    # Zlogaritmovanie (na decibely)
    # + 1e-10 konstanta vsade kvoli python vypoctu
    mag_log = 10 * np.log10(mag + 1e-10)
    plt.figure(figsize=(12, 6))
    plt.pcolormesh(t, f, mag_log, shading="gouraud", cmap="turbo")
    plt.ylabel("Frekvencia [Hz]")
    plt.xlabel("Čas [s]")
    plt.title(f"Spektrogram nahrávky {name}")
    plt.colorbar(label="Intenzita [dB]")
    plt.show()

Nahávka known_signals[114]:

# Booty booty nahravka známa ;)
if not os.environ.get('CI'):
  display(Audio(known_signals[114], rate=Fs))
show_spectrogram("known_signals", 114)

Nahávka eval_signals[0]:

# Booty booty nahravka nájdená ;)
if not os.environ.get('CI'):
  display(Audio(eval_signals[0], rate=Fs))
show_spectrogram("eval_signals", 0)

2.3 Úpravy jednotlivých nahrávok

Pri riešení som si nechal poradiť zadaním a rozhodol sa každú nahrávku orezať od zbytočných vysokých a nízkych frekvencií a taktiež normalizovať pre neskoršie porovnanie. To som docielil vo funkcii get_spectrogram():

1 Získanie amplitúdovej zložky (rovnako ako vyššie).

2 Zlogaritmovanie (na decibely rovnako ako vyššie).

3 Orezanie frekvencie ("pásmová propust")
Orezanie frekvencie bolo výhodné pre odstránenie potenciálenho šumu na vysokých a nízkych frekvenciách. Tieto frekvencie neobsahujú relevantné pre celkové porovnávanie hudby. Po odskúšaní niekoľko hodnôt som sa rozhodol orezať frekvencie od 250Hz [index 8] do 3969Hz [index 128] (kvôli priebežnému porovnávaniu som nechal rozsah ako argument funkcie)

4 Normalizácia bola nevyhnutná na správne porovnávanie nahrávok. Pri riešení som využil dva prístupy:

4.1. Globálna (Z-score) normalizácia (v kóde "Z")

Ako protný navrh som zvoli Z-score normalizáciu.
$$ S_{norm}[f, t] = \frac{S_{log}[f, t] - \mu}{\sigma} $$ Vypočíta sa stredná hodnota a smerodajná odchýlka z celej matice spektrogramu. Odčítaním strednej hodnoty od priemeru matice a následným vydelením smerodajnou odchýlkou sa hodnoty viac priblížia nule (zníži sa rozsah na osi intenzity) bez zmeny celkového výzoru spektrogramu. To mi pomohlo zjednotiť rozdiely v hlasitosti medzi nahrávkami pre lepšie porovnávanie medzi sebou.
Tento návrh mi validačnom datasete fungoval obstojne (Top 1 96% Top 5 100%) ale stále sa mi nepodarilo odsiahnuť ideálny výsledok, kvôli šumu v nahrávkach. (inšpirácia [https://www.geeksforgeeks.org/data-analysis/z-score-normalization-definition-and-examples/])

4.2. Normalizácia priemerným spektrom (v kóde "P")

Ako ďalší prístup som sa skúsiť návrh zo zadania, normalizáciu priemerným spektrom. $$ S_{norm}[f, t] = \frac{S_{log}[f, t] - \mu_f}{\sigma_f + \epsilon} $$ Funguje podobne ako Z-score normalizácia ale namiesto vypočítavania strednej hodnoty a smerodajnej odchýlky z celej matice sa vypočítavajú pre každú frekvenciu v čase zvlášť. Týmto sa docieli odstránenie šumu a zároveň sa normalizuje hlasitosť rovnako ako pri Z-score normalizácii.
Touto metódou sa mi celkovo na validačnom datasete zvýšila presnosť na (Top 1 100% Top 5 100%).

Nižšie sa nachádzajú ukážky jednotlivých krokov na spektrogramoch:

# Ziska spektrogram, logaritmuje, oreze a normalizuje
# Oreze podla parametrov od from_freq (Hz) do to_freq (Hz)
# Normalizuje podla typu
def get_spectrogram(sig, from_freq, to_freq, show, type):
    # 1. Vypocet spektrogramu (fft a ponechana iba amplitudova zlozka)
    #   f=freq ; t=čas; mag=amplituda
    #   (f_rozlisenie = 16000(Fs)/512(nperseg) = 31.25 Hz)
    f, t, mag = signal.spectrogram(sig, fs=Fs, window="hann", nperseg=512, noverlap=256)

    # 2. Zlogaritmovanie (na decibely)
    mag_log = 10 * np.log10(mag + 1e-10)  # + 1e-10
    # Ukazka
    if show == "All":
      plt.figure(figsize=(12, 6))
      plt.pcolormesh(t, f, mag_log, shading="gouraud", cmap="turbo")
      plt.ylabel("Frekvencia [Hz]")
      plt.xlabel("Čas [s]")
      plt.title(f"Spektrogram nahrávky bez úpravy (v dB)")
      plt.colorbar(label="Intenzita [dB]")
      plt.show()

    # 3. Orezanie frekvencii (pasmovy propust)
    # idx_v_spektre = pozadovana_freq / f_rozlisenie
    from_idx = int(round(from_freq / (Fs / 512)))
    to_idx = int(round(to_freq / (Fs / 512)))
    mag_cut = mag_log[from_idx:to_idx, :]
    f_cut = f[from_idx:to_idx]
    # Ukazka
    if show == "All":
        # print(f"Cut: From index: {from_idx}  To index: {to_idx}")
        plt.figure(figsize=(12, 6))
        plt.pcolormesh(t, f_cut, mag_cut, shading="gouraud", cmap="turbo")
        plt.ylabel("Frekvencia [Hz]")
        plt.xlabel("Čas [s]")
        plt.title(
            f"Spektrogram nahrávky orezaný ({f_cut[0]:.0f}Hz - {f_cut[-1]:.0f}Hz)  (v dB)"
        )
        plt.colorbar(label="Intenzita [dB]")
        plt.show()

    # 4. Normalizacia
    # **Priemernym spektrom
    if type == "P":
        # priemerna hodnota kazdej danej freq v case
        # (axis=1 znamena priemer cez casovu os)
        mean_spectrum = np.mean(mag_cut, axis=1, keepdims=True)
        std_dev_spectrum = np.std(mag_cut, axis=1, keepdims=True)
        # Odcitanie premeru a delenie smerodajnou odchylkou
        mag_norm = (mag_cut - mean_spectrum) / (std_dev_spectrum + 1e-10)
    # **Z-score
    elif type == "Z":
        mean = np.mean(mag_cut)
        std_dev = np.std(mag_cut)
        # Odcitanie premeru a delenie smerodajnou odchylkou
        mag_norm = (mag_cut - mean) / (std_dev + 1e-10)
    else:
        mag_norm = 0  # kvoli pythonu
    # Ukazka
    if show == "All" or show == "Norm":
        plt.figure(figsize=(12, 6))
        plt.pcolormesh(t, f_cut, mag_norm, shading="gouraud",cmap="turbo")
        plt.ylabel("Frekvencia [Hz]")
        plt.xlabel("Čas [s]")
        if type == "Z":
            plt.title(
                f"Spektrogram nahrávky orezaný ({f_cut[0]:.0f}Hz - {f_cut[-1]:.0f}Hz) a normalizovaný (Z-score)"
            )
        elif type == "P":
            plt.title(
                f"Spektrogram nahrávky orezaný ({f_cut[0]:.0f}Hz - {f_cut[-1]:.0f}Hz) a normalizovaný (Priemerným spektrom)"
            )
        else:
            print("ZLY ARGUMENT")
        plt.colorbar(label="Intenzita")
        plt.show()
    return mag_norm
# Ukazkove spustenie na Booty booty nahravach:

Nahrávka known_signals[144]

_ = get_spectrogram(known_signals[144], 250, 4000, "All", "Z")
_ = get_spectrogram(known_signals[144], 250, 4000, "Norm", "P")

Nahrávka eval_signals[0]:

_ = get_spectrogram(eval_signals[0], 250, 4000, "All", "Z")
_ = get_spectrogram(eval_signals[0], 250, 4000, "Norm", "P")

2.4 Hlavná funkcia na porovnanie

Pre určovanie podobnosti nahrávok funkcia compute_similarity_matrix() volá funkciu na ziskanie apltitúdovej zložky pre káždú známu a porovnávanávanú nahrávku (rovnakým spôsobom ako som popísal vyššie) a uloží ich do dvoch polí. Ako parameter som pridal aj typ normalizácie aby som jednoducho mohol demonštrovať pri porovnávaní oboch typov.

Ako prvá myšlienka mi napadlo posúvať testovaný signál cez známy, keďže je kratší (akoby som posúval obrázky cez seba a postupne pre každý časový okmažik vypočítal ich koreĺaciu/podobnosť). Tento prístup síce fungoval ale bol veľmi algoritmicky náročný využitím python cyklov.

Optimalizáciu som teda rozhodol riešit vypočítavaním korelácie pomocou konvolucie, keďže konvolúcia je "otočená" korelácia. Výsledok je rovnaký ako pri posúvaní pomocou cyklov - 2D matica korelácie.

Po získaní amplitúdových zložiek sa vykoná konvolúcia známeho signálu a otočeného neznámeho signálu pomocou signal.convolve(). Parameter mode="valid" značí, že sa má počítať bez zero-paddingu teda len na nenulovom intervale nahrávky a parameter method="fft" akceleruje výpočet pomocou Rýchlej Fourierovej trasnformácie. Výsledkom je 2D matica korelácie jednotlivých časových a frekvenčných posunov, ktorej výpočet je možné matematicky zapísať ako:

$$R[f, t] = \sum_{i=0}^{H_f-1} \sum_{j=0}^{W_t-1} S_{known}[f+i, t+j] \cdot S_{cmp}[i, j] = S_{known} \star S_{cmp} = S_{known} * S^{-}_{cmp}$$

$$\text{Podobnosť} = \frac{\max(R)}{H_f \cdot W_t} $$

• $R[f, t]$ je vysledná korelačná matica pre posun o frekvenciu $f$ a čas $t$
• $t$ je čas (index posunu)
• $f$ je frekvencia (index posunu)
• $H_f$ je výška matice porovnávanej nahrávky
• $W_t$ je šírka matice porovnávanej nahrávky

Pre výpočet podobnosti sa delí veľkosťou matice aby som dostal číslo na intervale $[0-1]$. Toto číslo sa zapíše do výslednej matice podobností pre kombináciu dvoch nahrávok. Takto sa postupne prejde každá nahrávka z neznámych s každou nahrávkou zo známych a zapíše do výslednej matice podobnosti.

# Hlavna funkcia na porovnavanie
# Rozhodol som sa vyuzit posuvanie
def compute_similarity_matrix(to_cmp_signals, known_signals, type, verbose=True):
    N_to_cmp = len(to_cmp_signals)
    N_known = len(known_signals)

    # Vypocet amplitud signalov
    if verbose:
        print("Spracovanie nahrávok...")
    known_mag_arr = []
    for i in range(N_known):
        known_mag_arr.append(get_spectrogram(known_signals[i], 250, 4000, "No", type))
    to_cmp_mag_arr = []
    for i in range(N_to_cmp):
        to_cmp_mag_arr.append(get_spectrogram(to_cmp_signals[i], 250, 4000, "No", type))

    # Vytvorenie pozadovanej matice podobnosti
    if verbose:
        print(f"Vytváranie matice {N_to_cmp}x{N_known} ...")
    similarities = np.zeros((N_to_cmp, N_known))

    if verbose:
      print("Počítanie podobnosti...")
    # Kazdy neznamy signal
    for i in range(N_to_cmp):
        to_cmp_sig = to_cmp_mag_arr[i]
        # Porovna s kazdym znamym
        for j in range(N_known):
            known_sig = known_mag_arr[j]
            # Vypocet korelacie pre kazdy casovy okamzik
            # mode='valid' znamena, ze nepocita s zero paddingom teda, 
            # ze prejde iba cez maticu známeho signálu a nie mimo
            # method='fft' pouziva rychlu fourierovu transformaciu aby bol vypocet rychlejsi
            corr_array = signal.convolve(known_sig, to_cmp_sig[::-1, ::-1], mode="valid", method="fft")

            # Vysledkom je najvacsia hodnota correlacie vydelena velkostou matice (normalizuje na rozsah 0-1)
            similarities[i, j] = np.max(corr_array) / to_cmp_sig.size
        if verbose:
            print(f"Spracovávaný porovnávaný signal {i+1}/{N_to_cmp} s {N_known} známymi")
    if verbose:
        print("Spracovanie dokončené")
    return similarities

2.5 Funkcie na porovnanie metód normalizácií

Vytvoril som si dve jednoduché funkcie na demonštráciu a porovnávanie úspešnosti metód normalizácie pre known_signals[114] a eval_signals[0]:

• compare_two_recordings() - ktorá zavolá funkciu na výpočet podobnosti vyššie, ale len pre dve nahrávky pre konkrétny typ normalizácie zadaný parametrom.
• visualize_convolution_3d() - ktorá pekne vizualizuje maticu 2D maticu korelácie pre dve nahrávky

# Test porovnavania
def compare_two_recordings(cmp_idx, known_idx, to_cmp_source, known_source, type):
    # 1. Vyberie signaly
    to_cmp_sig = to_cmp_source[cmp_idx]
    known_sig = known_source[known_idx]
    # 2. Zabali do arr aby brala funkcia compute_similarity_matrix a spusti ju
    to_cmp_signals = np.array([to_cmp_sig])
    known_signals = np.array([known_sig])
    result = compute_similarity_matrix(to_cmp_signals, known_signals, type, verbose=False)
    # Vyprintuje vysledok
    print(f"Výsledná podobnosť: {result[0, 0] * 100:.2f} %")
    return

# Vizualizacia vypocitanie korelacie/konvolucie
def visualize_convolution_3d(known_sig, to_cmp_sig):
    corr = signal.convolve(known_sig, to_cmp_sig[::-1, ::-1], mode='same', method='fft')
    # X os = t
    # Y os = f
    ny, nx = corr.shape
    x = np.arange(nx)
    y = np.arange(ny)
    X, Y = np.meshgrid(x, y)

    # Vykreslenie
    fig = plt.figure(figsize=(14, 8))
    ax = fig.add_subplot(111, projection='3d')
    surf = ax.plot_surface(X, Y, corr, cmap='jet', linewidth=0, antialiased=False)

    # Vrchol - finalna zhoda
    max_idx = np.unravel_index(np.argmax(corr), corr.shape)
    ax.scatter(max_idx[1], max_idx[0], corr.max(), color='red', s=100, label='Maximálna zhoda')

    ax.set_title("3D Vizualizácia Korelácie")
    ax.set_xlabel('Časový posun')
    ax.set_ylabel('Frekvenčný posun')
    ax.set_zlabel('Miera zhody (korelácia)')
    fig.colorbar(surf, shrink=0.5, aspect=5)
    plt.legend()
    plt.show()

# Booty booty test:

Ukážka pre normalizáciu Z-score:

visualize_convolution_3d(get_spectrogram(known_signals[114], 250, 4000, "No", "Z"), 
                         get_spectrogram(eval_signals[0], 250, 4000, "No", "Z"))
compare_two_recordings(0, 114, eval_signals, known_signals, "Z")

Výsledná podobnosť: 80.97 %

Ukážka pre normalizáciu priemerným spektrom:

visualize_convolution_3d(get_spectrogram(known_signals[114], 250, 4000, "No", "P"),
                          get_spectrogram(eval_signals[0], 250, 4000, "No", "P"))
compare_two_recordings(0, 114, eval_signals, known_signals, "P")

Výsledná podobnosť: 68.30 %

3 Finálna evaluácia metód

3.1 Porovnanie metód normalizácie

Síce sa môže zdať, že vyššie percento podobnosti znamená, že Z-score normalizácia je lepšia, opak je ale pravdou. Z-normalizácia síce priniesla vyššie percento podobnosti ale pri porovnaní obrázku s normalizáciou priemerným spektrom je vidno, že konzistencia výsledokov je ďaleko lepšia pri priemernom spektre. Naznačuje to aj vyššie spomenutý test evaluácie ktorý bol úspešnejší pri normalizácii priemerným spektrom. Rozhodol som sa preto použiť túto normalizáciu pri finálnom odovzdávaní výstupu porovnania mojich login nahrávok eval.txt.

3.2 Spúštanie kontrolného vyhodnotenia dát

# vyhodnotenie - funkcia vracia presnost Top-1 a Top-5 na validacnych datach
def eval(scores, key):
  indices = np.flip(np.argsort(scores), axis=-1) # chceme od najvyssieho po najnizsie ...
  #print(scores[0,key[0]], key[0], indices)
  top1acc = np.sum(key == indices[:,0]) / indices.shape[0]
  top5acc = 0
  for ii in range(5):
    top5acc += np.sum(key == indices[:,ii])
  top5acc /=  indices.shape[0]
  return top1acc, top5acc

key = np.loadtxt("valid/key.txt", delimiter = ',', usecols=(1), dtype ='int')

# tu prebieha vypocet a vyhodnotenie na validacnych datach ...
scores_matrix = compute_similarity_matrix(valid_signals, known_signals, "P",verbose=False)
#print (key)
top1, top5 = eval(scores_matrix, key)
print("Top 1 accuracy ", top1 * 100, "%, Top 5 accuracy ", top5 * 100, "%" )

Top 1 accuracy  100.0 %, Top 5 accuracy  100.0 %

3.3 Spúštanie finálneho vyhodnotenia nahrávok loginu

# Tento proces zaberie vela casu, ak je potrebny zmente na True.
if False:
  scores_matrix_final = compute_similarity_matrix(eval_signals, known_signals, "P")
  # na exportovanie
  np.savetxt("eval.txt", scores_matrix_final)