exercise05.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "489c7bcd-d5b9-4a10-beb9-e58120f4758b",
   "metadata": {},
   "source": [
    "Sascha Spors,\n",
    "Professorship Signal Theory and Digital Signal Processing,\n",
    "Institute of Communications Engineering (INT),\n",
    "Faculty of Computer Science and Electrical Engineering (IEF),\n",
    "University of Rostock,\n",
    "Germany\n",
    "\n",
    "# Data Driven Audio Signal Processing - A Tutorial with Computational Examples\n",
    "\n",
    "Winter Semester 2023/24 (Master Course #24512)\n",
    "\n",
    "- lecture: https://github.com/spatialaudio/data-driven-audio-signal-processing-lecture\n",
    "- tutorial: https://github.com/spatialaudio/data-driven-audio-signal-processing-exercise\n",
    "\n",
    "Feel free to contact lecturer frank.schultz@uni-rostock.de"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "201797c2-3cc5-4cce-b1b9-684286f3b759",
   "metadata": {},
   "source": [
    "# Exercise 5: SVD Matrix on Multitrack Audio \n",
    "\n",
    "## Objectives\n",
    "\n",
    "TBD"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "66703710",
   "metadata": {},
   "outputs": [],
   "source": [
    "# note: we will delete (own) and save files on the computer,\n",
    "# so we should check what\n",
    "# wavfile.write, wavfile.read, rmtree, mkdir and plt.savefig ist doing\n",
    "# in the cells below. I tested the notebook on two machines and it did what it\n",
    "# is intended to do (i.e. deleting or creating left_singular_vectors folder,\n",
    "# deleting or creating reduced_rank_mixdown folder, write *.wav files into there,\n",
    "# write *.png graphic files of output plots into the chosen folder 'path').\n",
    "# So, harmless stuff actually. Anyway, I am not responsible for data loss on\n",
    "# other machines, if something went wrong. So, please, check what you are doing\n",
    "# on your machine\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import os\n",
    "import pyloudnorm as pyln\n",
    "import shutil\n",
    "from scipy.io import wavfile\n",
    "from numpy.linalg import cond, matrix_rank\n",
    "from scipy.linalg import svd, norm, diagsvd, inv, pinv, null_space"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8c7c76a4",
   "metadata": {},
   "outputs": [],
   "source": [
    "# some 'global' variables that should appear here at top for convenience\n",
    "target_lufs = -32  # target loudness\n",
    "flag_channel = \"Mono\"  # 'Left', 'Right' or 'Mono'"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "72e22fb2",
   "metadata": {},
   "source": [
    "## Data Preparation\n",
    "\n",
    "This notebook gains most fun, if we have multi-track audio to work with, ideally on our own computer. Maybe we even have some own recorded multi-tracked music ;-)\n",
    "\n",
    "We can get some nice multi-track audio material from [https://www.cambridge-mt.com/ms/mtk/](https://www.cambridge-mt.com/ms/mtk/) .\n",
    "\n",
    "This data base supports Mike Senior's 2nd edition of the book \"[The Mixing Secrets For The Small Studio](https://cambridge-mt.com/ms/main/)\".\n",
    "\n",
    "The data base is provided for educational purpose and we will use it precisely this way...maybe in slightly different mixing approach than initially intended, but we will doing a **mixing job** here...**using the SVD**.\n",
    "\n",
    "So, if we need some multi-track material, we can download one of the zip files. Below we use 4 examples that worked out very illustratively.\n",
    "\n",
    "To create a matrix where the columns represent the multi-tracked files (either the left channel, or the right channel or summed to a mono channel) we need equal length *.wav files to work with the code below (I was too lazy to program this nicely).\n",
    "\n",
    "A reasonable work flow is to use a digital audio workstation (DAW) software, such as [Reaper](https://www.reaper.fm/purchase.php) (please support these guys, non-commercial license is very cheap).\n",
    "\n",
    "In Reaper we would do:\n",
    "- grabbing all multi-track files from a certain zip into an empty project\n",
    "- making them as separate tracks\n",
    "- choose an appropriate time selection\n",
    "- do some simple level mixing according to actual taste\n",
    "- selecting all tracks\n",
    "- open the `render to file` dialog\n",
    "    - choose Source: `Selected tracks (stems)`\n",
    "    - choose Bounds: `Time selection`\n",
    "    - Directory: create a `multi-track` folder as subfolder where the raw wav files are stored\n",
    "    - File name: `multi-track`\n",
    "    - make sure that original sampling rate is used\n",
    "    - Format: `Wav`\n",
    "    - WAV bit depth: `32 Bit FP` (uses more storage, but then we don't need to pay attention to numbers larger/smaller than 1 if some mixing was performed)\n",
    "    - click `Dry Run (no output`)\n",
    "    - check if this simulated rendering comes with the expected channels into the expected number of files with naming convention `multi_track-xxx.wav`\n",
    "    - we are ready to render by clicking `Render xxx files``\n",
    "- maybe it's a good idea to save this Reaper project for further reference, such as checking what different mixes do to the SVD data..."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "66d3ab71",
   "metadata": {},
   "source": [
    "## Data Generation\n",
    "\n",
    "Instead of setting up multi-track audio described above, we can create some data to play with, in order to make this notebook self-contained and running out of the box."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "62816cdf",
   "metadata": {},
   "outputs": [],
   "source": [
    "if True:\n",
    "    path_gen = \"svd_multitrack_audio/test_data/\"\n",
    "    # clear/del path and re-create new to start fresh\n",
    "    try:  # for left_singular_vectors/\n",
    "        shutil.rmtree(path_gen)\n",
    "    except OSError as e:\n",
    "        print(\"Error: %s : %s\" % (path_gen, e.strerror))\n",
    "    os.makedirs(path_gen)\n",
    "    path_multi = \"multi_track\"\n",
    "    os.mkdir(path_gen + path_multi)\n",
    "\n",
    "    fs = 48000\n",
    "    N = 48\n",
    "    k = np.arange(N * 1000)\n",
    "    dOm = 2 * np.pi / N\n",
    "\n",
    "    for track in range(16):\n",
    "        x = np.sin(dOm * (track + 1) * k)\n",
    "        x = np.tile(x, (2, 1)).T\n",
    "        tmp = \"%0*d\" % (3, track + 1)  # multitrack number with leading zeros\n",
    "        wavfile.write(\n",
    "            path_gen + path_multi + \"/\" + path_multi + \"-\" + tmp + \".wav\",\n",
    "            fs,\n",
    "            x.astype(np.float32),\n",
    "        )"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b1d460c8",
   "metadata": {},
   "source": [
    "## Choose Multitrack Example"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "781e056b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# choose one of these illustrative multitrack projects:\n",
    "\n",
    "path = \"svd_multitrack_audio/cfx_Mathematician_Full/\"\n",
    "# c:fx, 'Mathematician', https://soundcloud.com/c-fx\n",
    "# multitrack for educational use only:\n",
    "# https://mtkdata.cambridgemusictechnology.co.uk/MTK005/cfx_Mathematician.zip\n",
    "\n",
    "path = \"svd_multitrack_audio/MaurizioPagnuttiSextet_AllTheGinIsGone/\"\n",
    "# Maurizio Pagnutti Sextet, 'All The Gin Is Gone', https://www.artesuono.it/album.aspx?id=76&p=0\n",
    "# for educational use only:\n",
    "# https://multitracks.cambridge-mt.com/MaurizioPagnuttiSextet_AllTheGinIsGone.zip\n",
    "\n",
    "path = \"svd_multitrack_audio/cryonicPAX_Excessive/\"\n",
    "# cryonicPAX, 'Excessive'\n",
    "# multitrack for educational use only:\n",
    "# https://mtkdata.cambridgemusictechnology.co.uk/MTK015/cryonicPAX_Excessive.zip\n",
    "\n",
    "path = \"svd_multitrack_audio/Fin_Echoes/\"\n",
    "# FIN, 'Echoes'\n",
    "# multitrack for educational use only:\n",
    "# https://multitracks.cambridge-mt.com/Fin_Echoes.zip\n",
    "\n",
    "path = \"svd_multitrack_audio/test_data/\"\n",
    "# we use the audio that is generated above in order to make this notebook\n",
    "# self-contained"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ba965fe2",
   "metadata": {},
   "outputs": [],
   "source": [
    "# set up path names\n",
    "pathr = path + \"multi_track/\"\n",
    "pathu = path + \"left_singular_vectors/\"  # i.e. the column space of matrix\n",
    "patha = path + \"reduced_rank_mixdown/\"\n",
    "\n",
    "# clear/del path and re-create new to start fresh\n",
    "try:  # for left_singular_vectors/\n",
    "    shutil.rmtree(pathu)\n",
    "except OSError as e:\n",
    "    print(\"Error: %s : %s\" % (pathu, e.strerror))\n",
    "os.mkdir(pathu)\n",
    "\n",
    "try:  # for reduced_rank_mixdown/\n",
    "    shutil.rmtree(patha)\n",
    "except OSError as e:\n",
    "    print(\"Error: %s : %s\" % (pathu, e.strerror))\n",
    "os.mkdir(patha)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c89d1f31",
   "metadata": {},
   "source": [
    "## Read Audio Files into Matrix X"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "93ae1832",
   "metadata": {},
   "outputs": [],
   "source": [
    "# read in multitrack wavs and store in a matrix X\n",
    "\n",
    "files = sorted(os.listdir(pathr))\n",
    "print(files)\n",
    "flag_conc = True\n",
    "\n",
    "for i in files:\n",
    "    # consider only .wav (this excludes .wav.reapeaks files from Reaper)\n",
    "    if i[-4:] == \".wav\":\n",
    "        # consider only multi_track-\n",
    "        if i[0:12] == \"multi_track-\":\n",
    "            # so we read in all stuff with name 'multi_track-xxxxxxxx.wav'\n",
    "            # if we set up the folder multi_track properly, we only find\n",
    "            # multi_track-001.wav, multi_track-002.wav, multi_track-003.wav, ...\n",
    "            fs, tmp = wavfile.read(pathr + i)  # get raw data from file\n",
    "            # we assume that the files are stereo\n",
    "            if flag_channel == \"Mono\":\n",
    "                # make mono and make it (xxx, 1) dimension\n",
    "                x = np.expand_dims((tmp[:, 0] + tmp[:, 1]) / 2, 1)\n",
    "            elif flag_channel == \"Left\":\n",
    "                # take left channel and make it (xxx, 1) dimension\n",
    "                x = np.expand_dims(tmp[:, 0], 1)\n",
    "            elif flag_channel == \"Right\":\n",
    "                # take right channel and make it (xxx, 1) dimension\n",
    "                x = np.expand_dims(tmp[:, 1], 1)\n",
    "            else:\n",
    "                print(\"!!! unknown flag_channel !!!\")\n",
    "                break\n",
    "            # if we got here, we can print some log data\n",
    "            print(i, fs, \"Hz\", x.shape, x.dtype)\n",
    "            # non-elegant way to stack all stems (single channel tracks) into matrix\n",
    "            if flag_conc:  # very first x into non-exisiting X\n",
    "                X, flag_conc = x, False\n",
    "            else:  # X is already existing and axis=0 is assumed to be const\n",
    "                # for all concatenated x,\n",
    "                # so multi_track-001.wav, multi_track-002.wav, multi_track-003.wav, ...\n",
    "                # must have exactly same length\n",
    "                X = np.concatenate((X, x), axis=1)\n",
    "# we could ignore the warning  Chunk (non-data) not understood, skipping it.\n",
    "\n",
    "# since we know sampling frequency (assuming that all multitracks have same),\n",
    "# we can instantiate a lufs meter (to measure loudness), that is later to be used\n",
    "lufs_meter = pyln.Meter(fs)\n",
    "\n",
    "print(\"\\nmulti-track matrix X\", X.shape, X.dtype)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "26f8f8e2",
   "metadata": {},
   "source": [
    "## SVD of X\n",
    "\n",
    "Let us assume tall, thin, full column rank matrix $\\mathbf{X}$ with rank $r=N$ here, which is a fair assumption for our multi-track data. It is probably hard to set up a rank deficient tall, thin matrix with different multi-track audio channels (except two guitars play exactly the same stuff ;-) ).\n",
    "\n",
    "Full SVD \n",
    "\n",
    "$$\\mathbf{X}_{M \\times N} = \\mathbf{U}_{M \\times M} \\mathbf{S}_{M \\times N} (\\mathbf{V}_{N \\times N})^\\mathrm{H}$$\n",
    "\n",
    "Economy SVD\n",
    "\n",
    "$$\\mathbf{X}_{M \\times N} = \\mathbf{U}_{M \\times r} \\mathbf{S}_{r \\times r} (\\mathbf{V}_{N \\times r})^\\mathrm{H}$$\n",
    "\n",
    "with diagonal matrix\n",
    "\n",
    "$$\\mathbf{S}_{r \\times r}=\n",
    "\\begin{bmatrix}\n",
    "\\sigma_1  & 0  & 0 & 0\\\\\n",
    "0 & \\sigma_2 & 0 & 0\\\\\n",
    "0 & 0 & : & 0\\\\\n",
    "0 & 0 & 0 & \\sigma_r\n",
    "\\end{bmatrix}\n",
    "$$\n",
    "\n",
    "containing singular values $\\sigma_1 > \\sigma_2 > ... > \\sigma_r$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ce98c398",
   "metadata": {},
   "outputs": [],
   "source": [
    "# SVD of X and some checks\n",
    "[U, s, Vh] = svd(X, full_matrices=False)\n",
    "print(\"U shape: \", U.shape)\n",
    "print(\"Vh shape: \", Vh.shape)\n",
    "print(\"sing vals: \", s)\n",
    "print(\n",
    "    s.size,\n",
    "    \"well behaved singular values ==\",\n",
    "    X.shape[1],\n",
    "    \"tracks ?\",\n",
    ")\n",
    "print(\"condition number: \", cond(X), s[0] / s[-1])\n",
    "\n",
    "# check SVD synthesis of X vs. stored X\n",
    "print(norm(U @ np.diag(s) @ Vh - X, \"fro\"))\n",
    "print(norm(U @ np.diag(s) @ Vh - X, 2))\n",
    "print(norm(U @ np.diag(s) @ Vh - X, \"nuc\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f73630d8",
   "metadata": {},
   "source": [
    "## Weighted Left Singular Vectors of the Column Space\n",
    "\n",
    "We set up\n",
    "\n",
    "$$\\mathbf{H}_{M \\times r} =  \\mathbf{U}_{M \\times r} \\mathbf{S}_{r \\times r}$$\n",
    "\n",
    "as a matrix, that contains the left singular vectors weighted by their corresponding non-zero singular values. We consider the **economy SVD version** from above and are still under the assumption of full column rank, tall/thin matrix $\\mathbf{X}$, so $r=N$.\n",
    "\n",
    "In literature which discusses principal component analysis (PCA), this matrix $\\mathbf{H}$ is often termed as \n",
    "- principal component signals\n",
    "- principal component **scores**\n",
    "\n",
    "The right singular value matrix $\\mathbf{V}_{N \\times r}$ is then referred to as\n",
    "- principal component coefficients\n",
    "- principal component **loadings**\n",
    "\n",
    "So, there is a strong link of SVD and PCA, in fact one way to derive a prove of the PCA is via the SVD. Note, however, that PCA typically assumes matrices $\\mathbf{X}$, where the **columns** are assumed or enforced (by `scipy.stats.zscore()` ) to be **mean free** and very often also having **unit variance**.\n",
    "\n",
    "Besides that fundamental difference, the main motivation of the PCA and our little toy example is about the same: what can the weighted left singular vectors $\\mathbf{H}$ that span the column space tell us about our data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "81f169e6",
   "metadata": {},
   "outputs": [],
   "source": [
    "# weight the left singular (unit) vectors by their corresponding\n",
    "# singular values\n",
    "H = U @ np.diag(s)\n",
    "H.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b0389586",
   "metadata": {},
   "source": [
    "## Plot Column Space Vectors and Render Audio of it"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fea31549",
   "metadata": {},
   "outputs": [],
   "source": [
    "# here we\n",
    "# 1. we plot the weighted left singular vectors\n",
    "# 2. we create wav files of these signals such that we can listen to them\n",
    "# (loudness normalization using BS1770 for convenience)\n",
    "#\n",
    "# do not harm your ears!!!\n",
    "#\n",
    "# this handling indicates:\n",
    "# 1. the relative weighting of the signals in terms of their singular values\n",
    "# can be seen in the plots\n",
    "# 2. the weighted left singular vector signals in wave files are encoded such\n",
    "# that we perceive about same loudness when listening to all of them\n",
    "\n",
    "N = H.shape[0]\n",
    "t = np.arange(N) / fs\n",
    "\n",
    "nr = 3  # if we have more than 18 tracks we should\n",
    "nc = 6  # choose other nr, nc to fit into a subplot matrix\n",
    "# plt.figure\n",
    "fig, axs = plt.subplots(nr, nc, figsize=(16, 10))\n",
    "cnt = 0\n",
    "for r in range(nr):\n",
    "    for c in range(nc):\n",
    "        # we start plotting/wav writing with index 1\n",
    "        # in order to match rank 1...R wavfiles\n",
    "        axs[r, c].plot(t, H[:, cnt])\n",
    "        axs[r, c].set_ylim(-1.5, 1.5)\n",
    "        axs[r, c].set_ylabel(\n",
    "            r\"$\\sigma_{{{0:d}}} \\cdot U_{{{0:d}}}$\".format(cnt + 1)\n",
    "        )\n",
    "        axs[r, c].set_xlabel(\"t / s\")\n",
    "        axs[r, c].set_title(\n",
    "            r\"$\\sigma_{{{0:d}}}$={1:3.2f}\".format(cnt + 1, s[cnt])\n",
    "        )\n",
    "        axs[r, c].grid(True)\n",
    "        # calc 1770 loudness\n",
    "        lufs = lufs_meter.integrated_loudness(H[:, cnt])\n",
    "        # adapt to target_lufs\n",
    "        tmp = H[:, cnt] * 10 ** ((target_lufs - lufs) / 20)\n",
    "        # get info if clipping occurs\n",
    "        if np.max(np.abs(tmp)) > 1.0:\n",
    "            print(\"!!! wav file clipping !!!, decrease target_lufs to:\")\n",
    "            print(target_lufs - 20 * np.log10(np.max(np.abs(tmp))))\n",
    "        # check if we got target_lufs\n",
    "        lufs = lufs_meter.integrated_loudness(tmp)\n",
    "        print(\"track\", cnt + 1, lufs, \"dB LUFS\")\n",
    "        # write wav of column space signals\n",
    "        # for convenient listening all files exhibit about equal loudness\n",
    "        wavfile.write(\n",
    "            pathu + \"left_singular_vector_\" + str(cnt + 1) + \".wav\", fs, tmp\n",
    "        )\n",
    "        cnt += 1  # prep for next column space signal\n",
    "        # in case we have more subplots than col space signals we stop:\n",
    "        if cnt == H.shape[1]:\n",
    "            break\n",
    "plt.tight_layout()\n",
    "plt.savefig(path + \"plot_left_singular_vectors.png\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b652d667",
   "metadata": {},
   "source": [
    "## Apply Equal Mixing Gain "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e3cf5c9c",
   "metadata": {},
   "source": [
    "We could assume that the multi-track data is somehow pre-processed, such that we obtain a reasonable mixdown to a mono signal $\\mathbf{y}$ when applying equal mixing gains stored in the vector $\\mathbf{g} = [g, g, g, ..., g]^\\mathrm{T}$.\n",
    "\n",
    "This linear model is written as\n",
    "\n",
    "$$\\mathbf{y} = \\mathbf{X} \\mathbf{g}$$\n",
    "\n",
    "and we should assume unity gain $g=1$, so we have $\\mathbf{g} = [1, 1, 1, ..., 1]^\\mathrm{T}$\n",
    "\n",
    "With our nice SVD factorization above, we can write\n",
    "\n",
    "$$\\mathbf{y} = \\mathbf{X}\\mathbf{g} = \\mathbf{U} \\mathbf{S} \\mathbf{V}^\\mathrm{H} \\mathbf{g}$$\n",
    "\n",
    "and with the weighted column space matrix $\\mathbf{H} = \\mathbf{U} \\mathbf{S}$ we get\n",
    "\n",
    "$$\\mathbf{y} = \\mathbf{X}\\mathbf{g} = \\mathbf{H} \\cdot \\mathbf{V}^\\mathrm{H} \\mathbf{g}$$\n",
    "\n",
    "Therefore, the vector $\\mathbf{V}^\\mathrm{H} \\mathbf{g}$ defines the mixing weights for the linear combination of the column space signals.\n",
    "\n",
    "We should check the polarity and the level of these mixing weights."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "eb1a505f-55f1-4ff6-8dc7-cd8e09b3ebe0",
   "metadata": {},
   "outputs": [],
   "source": [
    "# check the polarity and level of the mixing weights\n",
    "# unity gain as weights in terms of the V space vectors (inner products)\n",
    "mixing_weight = Vh @ np.ones((s.size, 1))\n",
    "ui = np.arange(mixing_weight.size) + 1\n",
    "level = 10 * np.log10(np.abs(mixing_weight) ** 2)  # dB\n",
    "level_pos = np.copy(level)\n",
    "level_pos[mixing_weight <= 0] = 0\n",
    "level_neg = np.copy(level)\n",
    "level_neg[mixing_weight > 0] = 0\n",
    "\n",
    "plt.figure(figsize=(6, 3))\n",
    "plt.plot(\n",
    "    ui[np.squeeze(mixing_weight <= 0)],\n",
    "    level_neg[np.squeeze(mixing_weight <= 0)],\n",
    "    \"C0o-\",\n",
    "    label=\"polarity -\",\n",
    ")\n",
    "plt.plot(\n",
    "    ui[np.squeeze(mixing_weight > 0)],\n",
    "    level_pos[np.squeeze(mixing_weight > 0)],\n",
    "    \"C3o-\",\n",
    "    label=\"polarity +\",\n",
    ")\n",
    "\n",
    "plt.xticks(ui)\n",
    "plt.legend()\n",
    "plt.xlabel(\"U column index\")\n",
    "plt.ylabel(\"mixing weight in dB\")\n",
    "plt.grid(True)\n",
    "plt.tight_layout()\n",
    "plt.savefig(path + \"mixing_weights.png\")\n",
    "# print(mixing_weight)\n",
    "# print(20*np.log10(np.abs(mixing_weight)))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "63af16ff",
   "metadata": {},
   "source": [
    "## Low-Rank Approximation of X and Mixdown\n",
    "\n",
    "The best rank $q<=r$ approximation for $\\mathbf{X}$ is given as\n",
    "$$\\mathbf{X}_q = \\sum_{i=1}^{q} \\sigma_i \\mathbf{u}_i \\mathbf{v}_i^\\mathrm{H},$$\n",
    "well known as Eckart–Young theorem.\n",
    "\n",
    "It is interesting to listen to a mixdown using the equal mixing gains just as above, here now with different best rank approximations.\n",
    "\n",
    "Instead of creating the best rank approximation $\\mathbf{X}_q$ straightforwardly with a for loop for the sum above, we could also treat the singular value matrix as\n",
    "$$\\mathbf{S}^*_{r \\times r}=\n",
    "\\begin{bmatrix}\n",
    "\\sigma_1  & 0  & 0 & 0 & 0\\\\\n",
    "0 & \\sigma_2 & 0 & 0 & 0\\\\\n",
    "0 & 0 & \\sigma_q & 0 & 0\\\\\n",
    "0 & 0 & 0 & \\sigma_{r-1}=0 & 0\\\\\n",
    "0 & 0 & 0 & 0 & \\sigma_r=0\n",
    "\\end{bmatrix},\n",
    "$$\n",
    "thus setting all singular values $\\sigma_{>q}=0$, and then re-calculating\n",
    "$$\\mathbf{H}_{M \\times r} =  \\mathbf{U}_{M \\times r} \\mathbf{S}^*_{r \\times r}$$\n",
    "and\n",
    "$$\\mathbf{y} = \\mathbf{H} \\cdot \\mathbf{V}^\\mathrm{H} \\mathbf{g}$$\n",
    "using the economy SVD as above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bd3d24ed-5a43-4bbd-ba47-fc34510b3ec6",
   "metadata": {},
   "outputs": [],
   "source": [
    "# reconstruct X from SVD with rank 1...r\n",
    "# and mixdown to a mono signal\n",
    "# loudness normalization of mixdown\n",
    "for r in range(s.size):\n",
    "    s_reduced = np.copy(s)  # get original singular values\n",
    "    s_reduced[r + 1 :] = 0  # throw away lower right diagonal parts\n",
    "    # print(s_reduced)\n",
    "    mixing_weight = Vh @ np.ones((s.size, 1))\n",
    "    H = U @ np.diag(s_reduced)\n",
    "    tmp = H @ mixing_weight  # this is the mixdown\n",
    "    lufs = lufs_meter.integrated_loudness(tmp)  # calc 1770 loudness\n",
    "    tmp *= 10 ** ((target_lufs - lufs) / 20)  # adapt to target_lufs\n",
    "    if np.max(np.abs(tmp)) > 1.0:  # clipping might occur\n",
    "        print(\"!!! wav file clipping !!!, decrease target_lufs to:\")\n",
    "        print(target_lufs - 20 * np.log10(np.max(np.abs(tmp))))\n",
    "    lufs = lufs_meter.integrated_loudness(tmp)  # check if we got target_lufs\n",
    "    print(\"rank\", r + 1, \"approx\", lufs, \"dB LUFS\")\n",
    "    # write wav of rank reduced mixdown, for convenient listening files have equal lufs\n",
    "    wavfile.write(patha + \"rank_\" + str(r + 1) + \".wav\", fs, tmp)\n",
    "# by that we can listen how the reduced SVD-factorization of matrix X\n",
    "# acting on equal mixing input vector (i.e. all X matrix columns get unity gain)\n",
    "# sounds like"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b96bff94-6ef3-40e3-8273-2aa218d2edb7",
   "metadata": {},
   "source": [
    "## Copyright\n",
    "\n",
    "- the notebooks are provided as [Open Educational Resources](https://en.wikipedia.org/wiki/Open_educational_resources)\n",
    "- the text is licensed under [Creative Commons Attribution 4.0](https://creativecommons.org/licenses/by/4.0/)\n",
    "- the code of the IPython examples is licensed under the [MIT license](https://opensource.org/licenses/MIT)\n",
    "- feel free to use the notebooks for your own purposes\n",
    "- please attribute the work as follows: *Frank Schultz, Data Driven Audio Signal Processing - A Tutorial Featuring Computational Examples, University of Rostock* ideally with relevant file(s), github URL https://github.com/spatialaudio/data-driven-audio-signal-processing-exercise, commit number and/or version tag, year."
   ]
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "1743232f157dd0954c61aae30535e75a2972519a625c7e796bafe0cd9a07bf7e"
  },
  "kernelspec": {
   "display_name": "myddasp",
   "language": "python",
   "name": "myddasp"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}