{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# BirdCLEF 2023 - Audio Processing - First Steps\n\nThis notebook collects my first steps regarding to understand audio processing. \n\nCore reference is the starter notebook from the 2021 competition [BirdCLEF2021: Processing audio data](https://www.kaggle.com/code/stefankahl/birdclef2021-processing-audio-data). The goal was to understand the details of the respective transformation and techniques there.\n\n**All comments welcome!**\n\n## Table of Contents\n- [Preparation](#Preparation)\n- [Pick example](#Pick-example)\n- [Wavefront](#Wavefront)\n- [Spectrogram](#Spectrogram)\n    - [Without centering](#Without-centering)\n    - [With centering](#With-centering)\n- [Recap: Discrete Fourier Transform](#Recap:-Discrete-Fourier-Transform)\n    - [Real-valued DFT](#Real-valued-DFT)\n    - [Frequency](#Frequency)\n- [Amplitude to dB](#Amplitude-to-dB)\n- [Mel spectrogram](#Mel-spectrogram)\n    - [Mel basis](#Mel-basis)","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19"}},{"cell_type":"code","source":"import os\nimport numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\nimport librosa\nprint(\"librosa:\", librosa.__version__)","metadata":{"execution":{"iopub.status.busy":"2023-03-20T20:37:01.276254Z","iopub.execute_input":"2023-03-20T20:37:01.276740Z","iopub.status.idle":"2023-03-20T20:37:01.326030Z","shell.execute_reply.started":"2023-03-20T20:37:01.276702Z","shell.execute_reply":"2023-03-20T20:37:01.324455Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Preparation ","metadata":{}},{"cell_type":"code","source":"base_dir = \"/kaggle/input/birdclef-2023/\"\npath_train = base_dir + \"train_metadata.csv\" \ntrain_sound_dir = \"/kaggle/input/birdclef-2023/train_audio/\"\nsampling_rate = 32_000","metadata":{"execution":{"iopub.status.busy":"2023-03-20T20:37:02.278701Z","iopub.execute_input":"2023-03-20T20:37:02.279145Z","iopub.status.idle":"2023-03-20T20:37:02.287404Z","shell.execute_reply.started":"2023-03-20T20:37:02.279102Z","shell.execute_reply":"2023-03-20T20:37:02.284967Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train = pd.read_csv(path_train)\ntrain[\"path_ogg\"] = train_sound_dir + train[\"filename\"]\ntrain.head()","metadata":{"execution":{"iopub.status.busy":"2023-03-20T20:37:02.960811Z","iopub.execute_input":"2023-03-20T20:37:02.961293Z","iopub.status.idle":"2023-03-20T20:37:03.158619Z","shell.execute_reply.started":"2023-03-20T20:37:02.961214Z","shell.execute_reply":"2023-03-20T20:37:03.157503Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Pick example\nLet us pick a single example and analyze.","metadata":{}},{"cell_type":"code","source":"from IPython.display import Audio","metadata":{"execution":{"iopub.status.busy":"2023-03-20T20:37:05.478963Z","iopub.execute_input":"2023-03-20T20:37:05.480399Z","iopub.status.idle":"2023-03-20T20:37:05.485932Z","shell.execute_reply.started":"2023-03-20T20:37:05.480346Z","shell.execute_reply":"2023-03-20T20:37:05.484352Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# lets make it reproducible, otherwise use rec = train.sample(1).iloc[0]\nrec = train.loc[11470]\npath_ogg = rec[\"path_ogg\"]\nrec","metadata":{"execution":{"iopub.status.busy":"2023-03-20T20:37:06.044892Z","iopub.execute_input":"2023-03-20T20:37:06.045339Z","iopub.status.idle":"2023-03-20T20:37:06.058686Z","shell.execute_reply.started":"2023-03-20T20:37:06.045294Z","shell.execute_reply":"2023-03-20T20:37:06.057661Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"librosa.get_samplerate(path_ogg)","metadata":{"execution":{"iopub.status.busy":"2023-03-20T20:37:06.659114Z","iopub.execute_input":"2023-03-20T20:37:06.660436Z","iopub.status.idle":"2023-03-20T20:37:13.352802Z","shell.execute_reply.started":"2023-03-20T20:37:06.660386Z","shell.execute_reply":"2023-03-20T20:37:13.351154Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"Audio(path_ogg, rate=sampling_rate)","metadata":{"execution":{"iopub.status.busy":"2023-03-20T20:37:13.355718Z","iopub.execute_input":"2023-03-20T20:37:13.356424Z","iopub.status.idle":"2023-03-20T20:37:13.381786Z","shell.execute_reply.started":"2023-03-20T20:37:13.356378Z","shell.execute_reply":"2023-03-20T20:37:13.380432Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sig, sr = librosa.load(path_ogg, sr=sampling_rate, duration=15)\nsig, sr, len(sig) / sampling_rate","metadata":{"execution":{"iopub.status.busy":"2023-03-20T20:37:13.383729Z","iopub.execute_input":"2023-03-20T20:37:13.384481Z","iopub.status.idle":"2023-03-20T20:37:20.410176Z","shell.execute_reply.started":"2023-03-20T20:37:13.384344Z","shell.execute_reply":"2023-03-20T20:37:20.408511Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Wavefront\n- The librosa utility function seems to support higher resolution than the matplotlib one.","metadata":{}},{"cell_type":"code","source":"fig, (ax0, ax1) = plt.subplots(2, 1, figsize=(12, 6))\nax0.plot(sig)\nax0.set_title(\"plot\")\nlibrosa.display.waveshow(sig, sr=sr, ax=ax1)\nax1.set_title(\"librosa.display.waveshow\")\nplt.tight_layout()\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-03-20T20:37:20.412762Z","iopub.execute_input":"2023-03-20T20:37:20.413153Z","iopub.status.idle":"2023-03-20T20:37:21.741201Z","shell.execute_reply.started":"2023-03-20T20:37:20.413116Z","shell.execute_reply":"2023-03-20T20:37:21.739542Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Spectrogram\n- Note, pick the same FFT and window function as used by librosa!\n- Using the inital ``scipy.signal.windows.hann`` was always giving something slightly different (though filter looks very similar).\n- Spectrogram is an image of the Fourier transform for different time frames of the signals. Those frames usually partially overlap.\n- There is an additional filter in place to avoid border effects, so one computes the Discrete Fourier Transform (DFT) of ``window * frame``.","metadata":{}},{"cell_type":"code","source":"from librosa import get_fftlib\nfrom librosa import util\nfrom librosa.filters import get_window\nfrom scipy.signal.windows import hann as _hann","metadata":{"execution":{"iopub.status.busy":"2023-03-20T20:37:23.411028Z","iopub.execute_input":"2023-03-20T20:37:23.411520Z","iopub.status.idle":"2023-03-20T20:37:23.422596Z","shell.execute_reply.started":"2023-03-20T20:37:23.411473Z","shell.execute_reply":"2023-03-20T20:37:23.421000Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fft = get_fftlib()\nfft","metadata":{"execution":{"iopub.status.busy":"2023-03-20T20:37:25.527466Z","iopub.execute_input":"2023-03-20T20:37:25.527901Z","iopub.status.idle":"2023-03-20T20:37:25.536623Z","shell.execute_reply.started":"2023-03-20T20:37:25.527861Z","shell.execute_reply":"2023-03-20T20:37:25.535116Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"n_fft = 2048\nhop_length = int(n_fft // 4)\nwin_length = n_fft","metadata":{"execution":{"iopub.status.busy":"2023-03-20T20:37:27.018046Z","iopub.execute_input":"2023-03-20T20:37:27.019280Z","iopub.status.idle":"2023-03-20T20:37:27.025980Z","shell.execute_reply.started":"2023-03-20T20:37:27.019187Z","shell.execute_reply":"2023-03-20T20:37:27.024471Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"w = get_window(\"hann\", Nx=win_length)\n_w = _hann(win_length)\nfig, axs = plt.subplots(ncols=3, figsize=(14, 3))\naxs[0].plot(w)\naxs[0].set_title(\"librosa\")\naxs[1].plot(_w)\naxs[1].set_title(\"scipy\")\naxs[2].plot(w - _w)\naxs[2].set_title(\"diff\")\nplt.tight_layout()\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-03-20T20:37:27.728514Z","iopub.execute_input":"2023-03-20T20:37:27.728916Z","iopub.status.idle":"2023-03-20T20:37:28.217554Z","shell.execute_reply.started":"2023-03-20T20:37:27.728881Z","shell.execute_reply":"2023-03-20T20:37:28.216384Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Without centering\n- This is an easier slicing of the individual frames on which the FFT is being computed.\n- More precisely, the frames start with ``frame[0] = sig[0:n_fft], frame[1] = sig[hop_length:hop_length + n_fft], ...`` and don't pad beyond.","metadata":{}},{"cell_type":"code","source":"stop = None\nspec = librosa.stft(sig[:stop], n_fft=n_fft, hop_length=hop_length, win_length=win_length, center=False)\nprint(spec.shape)","metadata":{"execution":{"iopub.status.busy":"2023-03-20T20:37:33.688201Z","iopub.execute_input":"2023-03-20T20:37:33.688721Z","iopub.status.idle":"2023-03-20T20:37:33.727003Z","shell.execute_reply.started":"2023-03-20T20:37:33.688664Z","shell.execute_reply":"2023-03-20T20:37:33.725515Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# manual size calculation\nnrows = 1 + n_fft // 2\nncols = 1 + (len(sig[:stop]) - n_fft) // hop_length \nnrows, ncols","metadata":{"execution":{"iopub.status.busy":"2023-03-20T20:37:34.177020Z","iopub.execute_input":"2023-03-20T20:37:34.177455Z","iopub.status.idle":"2023-03-20T20:37:34.185242Z","shell.execute_reply.started":"2023-03-20T20:37:34.177414Z","shell.execute_reply":"2023-03-20T20:37:34.184314Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# frame check with [0, 1, ..., 9]\nutil.frame(np.arange(10), frame_length=3, hop_length=2)","metadata":{"execution":{"iopub.status.busy":"2023-03-20T20:37:34.832939Z","iopub.execute_input":"2023-03-20T20:37:34.834202Z","iopub.status.idle":"2023-03-20T20:37:34.845185Z","shell.execute_reply.started":"2023-03-20T20:37:34.834035Z","shell.execute_reply":"2023-03-20T20:37:34.843660Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# pick a column and compute manual\nidx = 15\nstart = idx * hop_length\nstop = start + n_fft\ns0 = spec[:, idx]\ns1 = fft.rfft(w * sig[start:stop]).astype(s0.dtype)\nprint(\"librosa:\", (s0.shape, s0.dtype))\nprint(\"manual:\", (s1.shape, s1.dtype))\nprint(\"allclose:\", np.allclose(s0, s1))","metadata":{"execution":{"iopub.status.busy":"2023-03-20T20:37:35.519057Z","iopub.execute_input":"2023-03-20T20:37:35.519528Z","iopub.status.idle":"2023-03-20T20:37:35.528366Z","shell.execute_reply.started":"2023-03-20T20:37:35.519483Z","shell.execute_reply":"2023-03-20T20:37:35.527324Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig, axs = plt.subplots(nrows=3, figsize=(14, 5), sharex=\"col\", gridspec_kw={\"hspace\": 0})\naxs[0].plot(np.abs(s0), label=\"librosa\")\naxs[0].legend(loc=\"upper right\")\naxs[1].plot(np.abs(s1), label=\"manual\")\naxs[1].legend(loc=\"upper right\")\naxs[2].plot(np.abs(s0 - s1), label=\"diff\")\naxs[2].legend(loc=\"upper right\")\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-03-20T20:37:36.011668Z","iopub.execute_input":"2023-03-20T20:37:36.013206Z","iopub.status.idle":"2023-03-20T20:37:36.653511Z","shell.execute_reply.started":"2023-03-20T20:37:36.013155Z","shell.execute_reply":"2023-03-20T20:37:36.652114Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### With centering\n- Centering pads the signal with zeros beyond border so that column ``i`` is centered at ``i * hop_length``.","metadata":{}},{"cell_type":"code","source":"stop = None\nspec = librosa.stft(sig[:stop], n_fft=n_fft, hop_length=hop_length, win_length=win_length, center=True)\nprint(spec.shape)","metadata":{"execution":{"iopub.status.busy":"2023-03-20T20:37:37.030623Z","iopub.execute_input":"2023-03-20T20:37:37.031060Z","iopub.status.idle":"2023-03-20T20:37:37.066846Z","shell.execute_reply.started":"2023-03-20T20:37:37.031019Z","shell.execute_reply":"2023-03-20T20:37:37.065330Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# manual size calculation\nnrows = 1 + n_fft // 2\nncols = 1 + (len(sig[:stop]) - 1) // hop_length \nnrows, ncols","metadata":{"execution":{"iopub.status.busy":"2023-03-20T20:37:37.714194Z","iopub.execute_input":"2023-03-20T20:37:37.715664Z","iopub.status.idle":"2023-03-20T20:37:37.723826Z","shell.execute_reply.started":"2023-03-20T20:37:37.715614Z","shell.execute_reply":"2023-03-20T20:37:37.722156Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# pick a column and compute manual\nidx = 1\nmid = idx * hop_length\nwidth = n_fft // 2\nstart = mid - width\nstop = mid + width\npad = (max(0, -start), 0)\nstart = max(0, start)\ntot_len = len(sig[start:stop]) + pad[0]\npad = pad[0], n_fft - tot_len\nh = np.pad(sig[start:stop], pad_width=pad)\ns0 = spec[:, idx]\ns1 = fft.rfft(w * h).astype(s0.dtype)\nprint(\"start, stop, pad:\", (start, stop, pad))\nprint(\"librosa:\", (s0.shape, s0.dtype))\nprint(\"manual:\", (s1.shape, s1.dtype))\nprint(\"allclose:\", np.allclose(s0, s1))","metadata":{"execution":{"iopub.status.busy":"2023-03-20T20:37:38.350077Z","iopub.execute_input":"2023-03-20T20:37:38.350508Z","iopub.status.idle":"2023-03-20T20:37:38.363773Z","shell.execute_reply.started":"2023-03-20T20:37:38.350466Z","shell.execute_reply":"2023-03-20T20:37:38.362379Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig, axs = plt.subplots(nrows=3, figsize=(14, 5), sharex=\"col\", gridspec_kw={\"hspace\": 0})\naxs[0].plot(np.abs(s0), label=\"librosa\")\naxs[0].legend(loc=\"upper right\")\naxs[1].plot(np.abs(s1), label=\"manual\")\naxs[1].legend(loc=\"upper right\")\naxs[2].plot(np.abs(s0 - s1), label=\"diff\")\naxs[2].legend(loc=\"upper right\")\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-03-20T20:37:39.043498Z","iopub.execute_input":"2023-03-20T20:37:39.043941Z","iopub.status.idle":"2023-03-20T20:37:39.492841Z","shell.execute_reply.started":"2023-03-20T20:37:39.043897Z","shell.execute_reply":"2023-03-20T20:37:39.491352Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Recap: Discrete Fourier Transform\n- Notation: A column vector is denoted as $\\left| v\\right>$, its complex conjugated row vector as $\\left< v \\right|$ and $\\left<u \\mid v\\right>$ denotes the dot product. The vectors $\\left| n \\right>$ with $n=0, \\dots, N-1$ are used to denote standard basis vectors, so any arbitrary signal is expressed as $\\left| s \\right> = \\sum_{n=0}^{N-1} s_n \\left| n \\right>$ with $s_n$ being the coordinates.\n- This precise Fourier transform used here is the Discrete Fourier Transform (DFT).\n- In the DFT the signal $\\left| s \\right>$ is projected onto (or expanded in) the following basis: $\\left|b_k\\right> = \\sum_n e^{i \\frac{2 \\pi}{N} kn }\\left|n \\right>$ with $k=0,\\dots,N-1$.\n- This is an un-normalized basis $\\left< b_k \\mid b_l \\right> = N \\delta_{kl}$.\n- The DFT has values $x_k = \\left< b_k \\mid s\\right> = \\sum_n e^{-i \\frac{2 \\pi}{N} kn} s_n$ for all $k$.\n- The $N$ factor is being used in the back trafo, i.e., $\\left| s\\right> = \\frac{1}{N}\\sum_k \\left|b_k\\right>\\left< b_k \\mid s\\right> $.\n\n### Real-valued DFT\n- If the signal is real-valued, $s_n = s_n^*$ then we have $\\left< b_k \\mid s\\right> = \\left< b_{N-k} \\mid s\\right>^*$ or $x_k = x_{N-k}^*$.\n- With this property the first set of coefficients uniquely determine all coefficients:\n    - $x_0 = \\sum_n s_n$\n    - $x_k = \\Re(x_k) + i \\Im(x_k) = \\frac{1}{2}(x_k + x_{N-k}) + i \\frac{1}{2i}(x_k - x_{N-k}^*)$ with $k \\leq \\frac{N}{2}$.\n- For even $N=2K$, this means for the coefficient $x_k$ that:\n    - $x_0$ real\n    - $x_k$ for $k=1, \\dots, K-1$ complex\n    - $x_K$ real\n    - First $1 + K - 1 + 1 = 1 + K$ coefficients contain all $1 + 2 (K -1) + 1 = 2K = N$ real parameters in the original vector.\n- For odd $N=2K + 1$, this means for the coefficient $x_k$ that:\n    - $x_0$ real\n    - $x_k$ for $k=1, \\dots, K$ complex\n    - First $1 + K$ coefficients contain all $1 + 2 K = N$ real parameters in the original vector.\n    \n### Frequency\n- Basis vector $\\left| b_k\\right>$ with $k=0$ corresponds to the maximum period of $T=\\infty$.\n- Basis vector with $k=1$ corresponds to a period of $N$ steps or $T = \\frac{N}{N_0}$ with $N_0$ being steps per second (sampling rate). This is central frequency unit then $f_1 = \\frac{N_0}{N}$.\n- More general, basis vector $k$ corresponds to $f_k = k f_1$.\n- Note that this goes till $k \\leq \\frac{N}{2}$, then one switches to negative frequencies, turning clockwise if seeing the complex coefficients $\\left< n \\mid b_k\\right> = e^{i \\frac{2 \\pi}{N} kn}$, see [Negative frequency](https://en.wikipedia.org/wiki/Negative_frequency).","metadata":{}},{"cell_type":"code","source":"N = 10\ns = np.random.randn(N)\nwith np.printoptions(precision=3, suppress=True):\n    display(fft.fft(s))\n    display(fft.fftfreq(N))\n    display(fft.rfft(s))\n    display(fft.rfftfreq(N))","metadata":{"execution":{"iopub.status.busy":"2023-03-20T20:18:30.880159Z","iopub.execute_input":"2023-03-20T20:18:30.880843Z","iopub.status.idle":"2023-03-20T20:18:30.897948Z","shell.execute_reply.started":"2023-03-20T20:18:30.880804Z","shell.execute_reply":"2023-03-20T20:18:30.896691Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"N = 11\ns = np.random.randn(N)\nwith np.printoptions(precision=3, suppress=True):\n    display(fft.fft(s))\n    display(fft.fftfreq(N))\n    display(fft.rfft(s))\n    display(fft.rfftfreq(N))","metadata":{"execution":{"iopub.status.busy":"2023-03-20T20:18:30.900228Z","iopub.execute_input":"2023-03-20T20:18:30.901165Z","iopub.status.idle":"2023-03-20T20:18:30.922537Z","shell.execute_reply.started":"2023-03-20T20:18:30.901119Z","shell.execute_reply":"2023-03-20T20:18:30.921012Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Amplitude to dB\n- Spectrum gets squared and converted to log scale, $S \\to 10 \\left( \\log_{10}{|S|^2} - \\log_{10}{S_{ref}^2}\\right)$\n- $S_{ref}$ is a reference value and below it is computed as the maximal value on full spectrum.\n- Note, this also has an effect on the relative top dB value.","metadata":{}},{"cell_type":"code","source":"spec = librosa.stft(sig, n_fft=n_fft, hop_length=hop_length, win_length=win_length, center=False)\nspec_db = librosa.amplitude_to_db(np.abs(spec), ref=np.max)\nspec.shape, spec_db.shape","metadata":{"execution":{"iopub.status.busy":"2023-03-20T20:18:30.923955Z","iopub.execute_input":"2023-03-20T20:18:30.924392Z","iopub.status.idle":"2023-03-20T20:18:30.992820Z","shell.execute_reply.started":"2023-03-20T20:18:30.924358Z","shell.execute_reply":"2023-03-20T20:18:30.991699Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.figure(figsize=(12, 6))\nlibrosa.display.specshow(\n    spec_db, sr=sampling_rate, x_axis=\"time\", y_axis=\"hz\", cmap=plt.get_cmap(\"viridis\")\n)\nplt.tight_layout()\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-03-20T20:18:30.994353Z","iopub.execute_input":"2023-03-20T20:18:30.994727Z","iopub.status.idle":"2023-03-20T20:18:32.165420Z","shell.execute_reply.started":"2023-03-20T20:18:30.994692Z","shell.execute_reply":"2023-03-20T20:18:32.163992Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# pick a column and compute manual\nidx = 30\nv0 = spec_db[:, idx]\nref = np.max(np.abs(spec))\nv1 = np.abs(spec[:, idx])\nv1 = np.maximum(10 * np.log10(v1**2) - 10 * np.log10(ref**2), -80)\nprint(\"librosa:\", (v0.shape, v0.dtype))\nprint(\"ref value:\", ref)\nprint(\"manual:\", (v1.shape, v1.dtype))\nprint(\"allclose:\", np.allclose(v0, v1))","metadata":{"execution":{"iopub.status.busy":"2023-03-20T20:18:32.167387Z","iopub.execute_input":"2023-03-20T20:18:32.168086Z","iopub.status.idle":"2023-03-20T20:18:32.185706Z","shell.execute_reply.started":"2023-03-20T20:18:32.168045Z","shell.execute_reply":"2023-03-20T20:18:32.184371Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig, axs = plt.subplots(nrows=3, figsize=(14, 5), sharex=\"col\", gridspec_kw={\"hspace\": 0})\naxs[0].plot(v0, label=\"librosa\")\naxs[0].legend(loc=\"upper right\")\naxs[1].plot(v1, label=\"manual\")\naxs[1].legend(loc=\"upper right\")\naxs[2].plot(np.abs(v0 - v1), label=\"diff\")\naxs[2].legend(loc=\"upper right\")\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-03-20T20:18:32.186955Z","iopub.execute_input":"2023-03-20T20:18:32.187310Z","iopub.status.idle":"2023-03-20T20:18:32.618480Z","shell.execute_reply.started":"2023-03-20T20:18:32.187276Z","shell.execute_reply":"2023-03-20T20:18:32.617237Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Mel spectrogram\n- Similar built as spectrogram but with other frequency scaling.\n- This mel scale is a **perceptual scale** of pitches, see [Mel scale](https://en.wikipedia.org/wiki/Mel_scale).\n    - Based how we perceive a tone.\n    - _A double a high perceived tone should have a double as high value_. This is not the case on the frequency scale. \n- Implemented as basis change on Fourier transform, i.e., columns of the spectrogram.\n- Parameters ``n_fft, hop_length`` have the same meaning as before.\n- Only frequency output is transformed to given number of mels and min/max frequency range.","metadata":{}},{"cell_type":"code","source":"from librosa.filters import mel","metadata":{"execution":{"iopub.status.busy":"2023-03-20T20:18:32.620005Z","iopub.execute_input":"2023-03-20T20:18:32.620359Z","iopub.status.idle":"2023-03-20T20:18:32.625763Z","shell.execute_reply.started":"2023-03-20T20:18:32.620324Z","shell.execute_reply":"2023-03-20T20:18:32.624534Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"spec_height = 128\nspec_width = 256\nseconds = 5\npower = 2\nn_fft = 1024\nfmin = 500\nfmax = 12_500","metadata":{"execution":{"iopub.status.busy":"2023-03-20T20:18:32.627611Z","iopub.execute_input":"2023-03-20T20:18:32.628420Z","iopub.status.idle":"2023-03-20T20:18:32.636076Z","shell.execute_reply.started":"2023-03-20T20:18:32.628324Z","shell.execute_reply":"2023-03-20T20:18:32.634834Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"n_mels = spec_height\nhop_length = (seconds * sampling_rate - n_fft) // (spec_width - 1) \nprint(\"hop_length:\", hop_length)\nprint(\"manual width:\", 1 + (seconds * sampling_rate - n_fft) // hop_length)\nprint(\"spec_height, spec_width:\", (spec_height, spec_width))","metadata":{"execution":{"iopub.status.busy":"2023-03-20T20:18:32.637990Z","iopub.execute_input":"2023-03-20T20:18:32.638531Z","iopub.status.idle":"2023-03-20T20:18:32.647471Z","shell.execute_reply.started":"2023-03-20T20:18:32.638480Z","shell.execute_reply":"2023-03-20T20:18:32.646332Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"start = 0\nstop = start + seconds * sampling_rate\nmel_spec = librosa.feature.melspectrogram(\n    y=sig[start:stop], \n    hop_length=hop_length,\n    sr=sampling_rate, \n    n_fft=n_fft, \n    n_mels=n_mels,\n    center=False,\n    fmin=fmin,\n    fmax=fmax,\n    power=power,\n)\nmel_spec.shape","metadata":{"execution":{"iopub.status.busy":"2023-03-20T20:18:32.649645Z","iopub.execute_input":"2023-03-20T20:18:32.650113Z","iopub.status.idle":"2023-03-20T20:18:33.906827Z","shell.execute_reply.started":"2023-03-20T20:18:32.650064Z","shell.execute_reply":"2023-03-20T20:18:33.905481Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Mel basis","metadata":{}},{"cell_type":"code","source":"mel_basis = mel(sr=sampling_rate, n_fft=n_fft, n_mels=n_mels, fmin=fmin, fmax=fmax)\nmel_basis.shape","metadata":{"execution":{"iopub.status.busy":"2023-03-20T20:18:33.916645Z","iopub.execute_input":"2023-03-20T20:18:33.918373Z","iopub.status.idle":"2023-03-20T20:18:33.933582Z","shell.execute_reply.started":"2023-03-20T20:18:33.918304Z","shell.execute_reply":"2023-03-20T20:18:33.931862Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# compute frequency bounds\nfreqs = fft.rfftfreq(n_fft, 1 / sampling_rate)\nidxmin = np.where(fmin <= freqs)[0].min()\nidxmax = np.where(freqs <= fmax)[0].max()\nprint(f\"min valid freq: {freqs[idxmin]} at {idxmin}\")\nprint(f\"max valid freq: {freqs[idxmax]} at {idxmax}\")","metadata":{"execution":{"iopub.status.busy":"2023-03-20T20:18:33.935879Z","iopub.execute_input":"2023-03-20T20:18:33.936724Z","iopub.status.idle":"2023-03-20T20:18:33.949894Z","shell.execute_reply.started":"2023-03-20T20:18:33.936666Z","shell.execute_reply":"2023-03-20T20:18:33.948310Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"S = librosa.stft(y=sig[start:stop], n_fft=n_fft, hop_length=hop_length, center=False)\nS = np.abs(S)**power\nS.shape","metadata":{"execution":{"iopub.status.busy":"2023-03-20T20:18:33.952070Z","iopub.execute_input":"2023-03-20T20:18:33.952760Z","iopub.status.idle":"2023-03-20T20:18:33.985616Z","shell.execute_reply.started":"2023-03-20T20:18:33.952696Z","shell.execute_reply":"2023-03-20T20:18:33.983838Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# manual size calculation\nnrows = 1 + n_fft // 2\nncols = 1 + (len(sig[start:stop]) - n_fft) // hop_length \nnrows, ncols","metadata":{"execution":{"iopub.status.busy":"2023-03-20T20:18:33.988708Z","iopub.execute_input":"2023-03-20T20:18:33.990075Z","iopub.status.idle":"2023-03-20T20:18:34.003445Z","shell.execute_reply.started":"2023-03-20T20:18:33.990001Z","shell.execute_reply":"2023-03-20T20:18:34.001646Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# manual basis transform\nmal_spec_manual = mel_basis @ S\nprint(\"allclose:\", np.allclose(mel_spec, mal_spec_manual))","metadata":{"execution":{"iopub.status.busy":"2023-03-20T20:18:34.006765Z","iopub.execute_input":"2023-03-20T20:18:34.008360Z","iopub.status.idle":"2023-03-20T20:18:34.020085Z","shell.execute_reply.started":"2023-03-20T20:18:34.008278Z","shell.execute_reply":"2023-03-20T20:18:34.018046Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.matshow(mel_basis)\nplt.axvline(x=idxmin)\nplt.axvline(x=idxmax)\nplt.xlabel(\"freq index\")\nplt.ylabel(\"mel index\")\nplt.title(\"mel basis transform\")\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-03-20T20:18:34.023660Z","iopub.execute_input":"2023-03-20T20:18:34.025051Z","iopub.status.idle":"2023-03-20T20:18:34.431314Z","shell.execute_reply.started":"2023-03-20T20:18:34.024980Z","shell.execute_reply":"2023-03-20T20:18:34.429732Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.matshow(mel_basis @ mel_basis.transpose())\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-03-20T20:18:34.432854Z","iopub.execute_input":"2023-03-20T20:18:34.433217Z","iopub.status.idle":"2023-03-20T20:18:34.752237Z","shell.execute_reply.started":"2023-03-20T20:18:34.433180Z","shell.execute_reply":"2023-03-20T20:18:34.750383Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# to dB\nmel_spec_db = librosa.power_to_db(mel_spec, ref=np.max)\nmel_spec_db.shape","metadata":{"execution":{"iopub.status.busy":"2023-03-20T20:18:34.754896Z","iopub.execute_input":"2023-03-20T20:18:34.756079Z","iopub.status.idle":"2023-03-20T20:18:34.766559Z","shell.execute_reply.started":"2023-03-20T20:18:34.755910Z","shell.execute_reply":"2023-03-20T20:18:34.765710Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"librosa.display.specshow(\n    mel_spec_db, \n    sr=sampling_rate, \n    hop_length=hop_length, \n    x_axis='time', \n    y_axis='mel',\n    fmin=fmin, \n    fmax=fmax,\n    cmap=plt.get_cmap('viridis')\n)\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-03-20T20:18:34.767941Z","iopub.execute_input":"2023-03-20T20:18:34.768531Z","iopub.status.idle":"2023-03-20T20:18:35.049854Z","shell.execute_reply.started":"2023-03-20T20:18:34.768496Z","shell.execute_reply":"2023-03-20T20:18:35.048699Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}