{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Excitation Signals for Room Impulse Response Measurement\n",
    "\n",
    "### Criteria\n",
    "- Sufficient signal energy over the entire frequency range of interest\n",
    "- Dynamic range\n",
    "- Crest factor (peak-to-RMS value)\n",
    "- Noise rejection (repetition and average, longer duration)\n",
    "- Measurement duration\n",
    "- Time variance\n",
    "- Nonlinear distortion\n",
    "\n",
    "#### _References_\n",
    "* Müller, Swen, and Paulo Massarani. \"Transfer-function measurement with sweeps.\" Journal of the Audio Engineering Society 49.6 (2001): 443-471.\n",
    "[link](http://www.aes.org/e-lib/browse.cfm?elib=10189)\n",
    "\n",
    "* Farina, Angelo. \"Simultaneous measurement of impulse response and distortion with a swept-sine technique.\" Audio Engineering Society Convention 108. Audio Engineering Society, 2000.\n",
    "[link](http://www.aes.org/e-lib/browse.cfm?elib=10211)\n",
    "\n",
    "* Farina, Angelo. \"Advancements in impulse response measurements by sine sweeps.\" Audio Engineering Society Convention 122. Audio Engineering Society, 2007.\n",
    "[link](http://www.aes.org/e-lib/browse.cfm?elib=14106)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from scipy.signal import chirp, max_len_seq, freqz, fftconvolve, resample\n",
    "import sounddevice as sd\n",
    "import tools"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def crest_factor(x):\n",
    "    \"\"\"Peak-to-RMS value (crest factor) of the signal x\n",
    "\n",
    "    Parameter\n",
    "    ---------\n",
    "    x : array_like\n",
    "        signal\n",
    "    \"\"\"\n",
    "\n",
    "    return np.max(np.abs(x)) / np.sqrt(np.mean(x**2))\n",
    "\n",
    "\n",
    "def circular_convolve(x, y, outlen):\n",
    "    \"\"\"Circular convolution of x and y\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    x : array_like\n",
    "        Real-valued signal\n",
    "    y : array_like\n",
    "        Real-valued signal\n",
    "    outlen : int\n",
    "        Length of the output\n",
    "    \"\"\"\n",
    "\n",
    "    return np.fft.irfft(np.fft.rfft(x, n=outlen) * np.fft.rfft(y, n=outlen), n=outlen)\n",
    "\n",
    "\n",
    "def plot_time_domain(x, fs=44100, ms=False):\n",
    "\n",
    "    time = np.arange(len(x)) / fs\n",
    "    timeunit = 's'\n",
    "    if ms:\n",
    "        time *= 1000\n",
    "        timeunit = 'ms'\n",
    "\n",
    "    fig = plt.figure()\n",
    "    plt.plot(time, x)\n",
    "    plt.xlabel('Time / {}'.format(timeunit))\n",
    "    return\n",
    "\n",
    "\n",
    "def plot_freq_domain(x, fs=44100, khz=False):\n",
    "\n",
    "    Nf = len(x) // 2 + 1\n",
    "    freq = np.arange(Nf) / Nf * fs / 2\n",
    "    frequnit = 'Hz'\n",
    "    if khz:\n",
    "        freq /= 1000\n",
    "        frequnit = 'kHz'\n",
    "\n",
    "    fig = plt.figure()\n",
    "    plt.plot(freq, db(np.fft.rfft(x)))\n",
    "    plt.xscale('log')\n",
    "    plt.xlabel('Frequency / {}'.format(frequnit))\n",
    "    plt.ylabel('Magnitude / dB')\n",
    "    return\n",
    "\n",
    "\n",
    "def compare_irs(h1, h2, ms=False):\n",
    "    t1 = np.arange(len(h1)) / fs\n",
    "    t2 = np.arange(len(h2)) / fs\n",
    "    timeunit = 's'\n",
    "    if ms:\n",
    "        t1 *= 1000\n",
    "        t2 *= 1000\n",
    "        timeunit = 'ms'\n",
    "    fig = plt.figure()\n",
    "    plt.plot(t1, h1, t2, h2)\n",
    "    plt.xlabel('Time / {}'.format(timeunit))\n",
    "    return\n",
    "\n",
    "\n",
    "def compare_tfs(h1, h2, khz=False):\n",
    "    n1 = len(h1) // 2 + 1\n",
    "    n2 = len(h2) // 2 + 1\n",
    "    f1 = np.arange(n1) / n1 * fs / 2\n",
    "    f2 = np.arange(n2) / n2 * fs / 2\n",
    "    frequnit = 'Hz'\n",
    "    if khz:\n",
    "        freq /= 1000\n",
    "        frequnit = 'khz'\n",
    "    fig = plt.figure()\n",
    "    plt.plot(f1, db(np.fft.rfft(h1)), f2, db(np.fft.rfft(h2)))\n",
    "    plt.xscale('log')\n",
    "    plt.xlabel('Frequency / {}'.format(frequnit))\n",
    "    plt.ylabel('Magnitude / dB')\n",
    "    return\n",
    "\n",
    "\n",
    "def pad_zeros(x, nzeros):\n",
    "    \"\"\"Append zeros at the end of the input sequence\n",
    "    \"\"\"\n",
    "    return np.pad(x, (0, nzeros), mode='constant', constant_values=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fs = 44100\n",
    "dur = 1\n",
    "L = int(np.ceil(dur * fs))\n",
    "time = np.arange(L) / fs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## White Noise"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Generate a random signal with normal (Gaussian) amplitude distribution. Use `numpy.random.randn` and normalize the amplitude with `tools.normalize`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's listen to it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot the signal in the time domain and in the frequency domain."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Is the signal really white?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What is the crest factor of a white noise?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now feed the white noise to an unkown system `tools.blackbox` and save the output signal."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "How do you think we can extract the impulse response of the system?\n",
    "Try to compute the impulse response from the output signal.\n",
    "Compare it with the actual impulse response which can be obtained by feeding an ideal impulse to `tools.blackbox`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Maximum Length Sequence\n",
    "\n",
    "> Maximum-length sequences (MLSs) are binary sequences that can be generated very easily with an N-staged shift register and an XOR gate (with up to four inputs) connected with the shift register in such a way that all possible 2N states, minus the case \"all 0,\" are run through. This can be accomplished by hardware with very few simple TTL ICs or by software with less than 20 lines of assembly code.\n",
    "\n",
    "(Müller 2001)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nbit = int(np.ceil(np.log2(L)))\n",
    "mls, _ = max_len_seq(nbit)  # sequence of 0 and 1\n",
    "mls = 2*mls - 1  # sequence of -1 and 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Take a look at the signal in the time domain."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Examine the properties of the MLS\n",
    "* frequency response\n",
    "* crest factor\n",
    "* simulate the impulse response measurement of `tools.blackbox`\n",
    "* evaluate the obtained impulse response"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In practice, the (digital) signal has to be converted into an analog signal by an audio interface?\n",
    "Here, the process is simulated by oversampling the signal by a factor of 10.\n",
    "Pay attention to the crest factor before and after upsampling."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "upsample = 10\n",
    "mls_up = resample(mls, num=len(mls) * upsample)\n",
    "time = np.arange(len(mls)) / fs\n",
    "time_up = np.arange(len(mls_up)) / fs / upsample\n",
    "\n",
    "plt.figure(figsize=(10, 4))\n",
    "plt.plot(time_up, mls_up, '-', label='Analog')\n",
    "plt.plot(time, mls, '-', label='Digital')\n",
    "plt.legend(loc='best')\n",
    "plt.xlabel('Time / s')\n",
    "plt.title('Crest factor {:.1f} -> {:.1f} dB'.format(\n",
    "    tools.db(crest_factor(mls)), tools.db(crest_factor(mls_up))))\n",
    "\n",
    "plt.figure(figsize=(10, 4))\n",
    "plt.plot(time_up, mls_up, '-', label='Analog')\n",
    "plt.plot(time, mls, 'o', label='Ditigal')\n",
    "plt.xlim(0, 0.0025)\n",
    "plt.legend(loc='best')\n",
    "plt.xlabel('Time / s')\n",
    "plt.title('Crest factor {:.1f} -> {:.1f} dB'.format(\n",
    "    tools.db(crest_factor(mls)), tools.db(crest_factor(mls_up))))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Linear Sweep"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Generate a linear sweep with `lin_sweep`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def lin_sweep(fstart, fstop, duration, fs):\n",
    "    \"\"\"Generation of a linear sweep signal.\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    fstart : int\n",
    "           Start frequency in Hz\n",
    "    fstop  : int\n",
    "           Stop frequency in Hz\n",
    "    duration : float\n",
    "           Total length of signal in s\n",
    "    fs : int\n",
    "           Sampling frequency in Hz\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    array_like\n",
    "         generated signal vector\n",
    "\n",
    "    Note that the stop frequency must not be greater than half the\n",
    "    sampling frequency (Nyquist-Shannon sampling theorem).\n",
    "\n",
    "    \"\"\"\n",
    "    if fstop > fs / 2:\n",
    "        raise ValueError(\"fstop must not be greater than fs/2\")\n",
    "    t = np.arange(0, duration, 1 / fs)\n",
    "    excitation = np.sin(\n",
    "        2 * np.pi * ((fstop - fstart) /\n",
    "                     (2 * duration) * t ** 2 + fstart * t))\n",
    "    # excitation = excitation - np.mean(excitation)  # remove direct component\n",
    "    return excitation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fs = 44100\n",
    "fstart =\n",
    "fstop =\n",
    "duration =\n",
    "\n",
    "lsweep ="
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Examine the properties of linear sweeps\n",
    "* spectrogram (Use `pyplot.specgram` with `NFFT=512` and `Fs=44100`)\n",
    "* frequency response\n",
    "* crest factor\n",
    "* simulate the impulse response measurement of `tools.blackbox`\n",
    "* evaluate the obtained impulse response"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exponential Sweep"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Generate a exponential sweep with `exp_sweep`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def exp_sweep(fstart, fstop, duration, fs):\n",
    "    \"\"\"Generation of a exponential sweep signal.\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    fstart : int\n",
    "           Start frequency in Hz\n",
    "    fstop  : int\n",
    "           Stop frequency\n",
    "    duration : float\n",
    "           Total length of signal in s\n",
    "    fs : int\n",
    "           Sampling frequency in Hz\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    array_like\n",
    "           Generated signal vector\n",
    "\n",
    "    Note that the stop frequency must not be greater than half the\n",
    "    sampling frequency (Nyquist-Shannon sampling theorem).\n",
    "\n",
    "    \"\"\"\n",
    "    if fstop > fs / 2:\n",
    "        raise ValueError(\"fstop must not be greater than fs/2\")\n",
    "    t = np.arange(0, duration, 1 / fs)\n",
    "    excitation = np.sin(2 * np.pi * duration *\n",
    "                        fstart / np.log(fstop / fstart) *\n",
    "                        (np.exp(t / duration * np.log(fstop / fstart)) - 1))\n",
    "    # excitation = excitation - np.mean(excitation)  # remove direct component\n",
    "    return excitation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fs = 44100\n",
    "fstart = \n",
    "fstop = \n",
    "duration = \n",
    "\n",
    "esweep = "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Examine the properties of linear sweeps\n",
    "* spectrogram (Use `pyplot.specgram` with `NFFT=512` and `Fs=44100`)\n",
    "* frequency response\n",
    "* crest factor\n",
    "* simulate the impulse response measurement of `tools.blackbox`\n",
    "* evaluate the obtained impulse response"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "mycomac",
   "language": "python",
   "name": "mycomac"
  },
  "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.9.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}