{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib notebook"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Adding random noise to frequency-domain CSEM data\n\nAdding random noise to frequency-domain CSEM data is not a trivial task, and\nthere are many different ways how it can be done. The problem comes from the\nfact that we live in the time domain, we do our measurements in the time\ndomain, our noise is therefore time-domain noise, but we want to add\ncomplex-valued noise in the frequency domain.\n\nHere we are going to look at some possibilities. However, keep in mind that\nthere are more possibilities than the ones shown here.\n\n## 1. Theory\n\nLet's assume we have complex-valued data $d=x+\\text{i}y$. We can add\nrandom noise to the data in the following way,\n\n\\begin{align}:label: generalnoise\n\n    \\tilde{d} = d + \\sigma \\left[(1 + \\text{i})\\,\\mu + \\mathcal{R} \\right] \\, ,\\end{align}\n\nwhere $\\tilde{d}$ is the data with added noise, $\\sigma$ is the\nstandard deviation, $\\mu$ is the mean value of the randomly distributed\nnoise, and $\\mathcal{R}$ is the random noise. We define the standard\ndeviation as\n\n\\begin{align}:label: stdev\n\n    \\sigma = \\sqrt{\\epsilon_\\text{n}^2 + \\left(\\epsilon_\\text{r}|d|\\right)^2 }\n    \\, ,\\end{align}\n\nwhere $\\epsilon_\\text{n}$ is the noise floor and\n$\\epsilon_\\text{r}$ is the relative error.\n\nWe compare here three ways of computing the random noise $\\mathcal{R}$.\nOf course there are other possibilities, e.g., one could make the non-zero mean\na random realization itself.\n\n\n1. Adding random uniform phases but constant amplitude\n\n   .. math::\n       :label: uniform\n\n       \\mathcal{R}_\\text{wn} = \\exp[\\text{i}\\,\\mathcal{U}(0, 2\\pi)] \\, ,\n\n   where $\\mathcal{U}(0, 2\\pi)$ is the uniform distribution and its\n   range. This adds white noise with a flat amplitude spectrum and random\n   phases.\n\n\n2. Adding Gaussian noise to real and imaginary parts\n\n   - Adding correlated random Gaussian noise\n\n      .. math::\n          :label: cgaussian\n\n          \\mathcal{R}_\\text{gc} = (1+\\text{i})\\,\\mathcal{N}(0, 1) \\, ,\n\n      where $\\mathcal{N}(0, 1)$ is the standard normal distribution of\n      zero mean and unit standard deviation.\n\n   - Adding uncorrelated random Gaussian noise\n\n      Above is the correlated version. Noise could also be added completely\n      uncorrelated,\n\n      .. math::\n          :label: ugaussian\n\n          \\mathcal{R}_\\text{gu} =\n          \\mathcal{N}(0, 1) + \\text{i}\\,\\mathcal{N}(0, 1) \\, .\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import empymod\nimport numpy as np\nimport matplotlib.pyplot as plt\nplt.style.use('bmh')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "### Noise computation\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Initiate random number generator.\nrng = np.random.default_rng()\n\n\ndef add_noise(data, ntype, rel_error, noise_floor, mu):\n    \"\"\"Add random noise to complex-valued data.\n\n    If `ntype='white_noise'`, complex noise is generated from uniform randomly\n    distributed phases.\n\n    If `ntype='gaussian_correlated'`, correlated Gaussian random noise is added\n    to real and imaginary part.\n\n    If `ntype='gaussian_uncorrelated'`, uncorrelated Gaussian random noise is\n    added to real and imaginary part.\n\n    \"\"\"\n\n    # Standard deviation\n    std = np.sqrt(noise_floor**2 + (rel_error*abs(data))**2)\n\n    # Random noise\n    if ntype == 'gaussian_correlated':\n        noise = rng.standard_normal(data.size)*(1+1j)\n    elif ntype == 'gaussian_uncorrelated':\n        noise = 1j*rng.standard_normal(data.size)\n        noise += rng.standard_normal(data.size)\n    else:\n        noise = np.exp(1j * rng.uniform(0, 2*np.pi, data.size))\n\n    # Scale and move noise; add to data and return\n    return data + std*((1+1j)*mu + noise)\n\n\ndef stack(n, data, ntype, **kwargs):\n    \"\"\"Stack n-times the noise, return normalized.\"\"\"\n    out = add_noise(data, ntype, **kwargs)/n\n    for _ in range(n-1):\n        out += add_noise(data, ntype, **kwargs)/n\n    return out"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## 2. Graphical illustration\n\nThe following is a graphical illustration. Please note that the relative\nerror is **very** high (20%)! This is only for illustration purposes.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Inputs\nd = np.array([6+2j])         # observed data point\nmean = 0.3                   # Non-zero mean\nrelative_error = 0.2         # Very high relative error\nstd = relative_error*abs(d)  # std (without noise floor)\n\n# Create figure\nfig, axs = plt.subplots(2, 1, figsize=(8, 10), sharex=True, sharey=True)\nax1, ax2 = axs\n\n# Titles\nfig.suptitle(r\"Random noise with $\\epsilon_n = 0, \"\n             f\"\\\\epsilon_r={relative_error}, \\\\mu={mean}$\", y=1, fontsize=20)\nax1.set_title('Theoretical distributions')\nax2.set_title('Random realizations')\n\n# Plot data point\nfor ax in axs:\n    ax.plot(np.r_[0., d.real], np.r_[0., d.imag], '--', c='.5')\n    ax.plot(d.real, d.imag, 'ko', ms=10, label='$d^{obs}$', zorder=10)\n\n\n# Mean and standard deviation\nax1.plot(d.real+np.r_[0, std*mean], d.imag+np.r_[0, std*mean],\n         'C8', label=r'Scaled mean $\\sigma (1+i)\\mu$', zorder=9)\nax1.plot(d.real+np.r_[std*mean, std*(1+mean)],\n         d.imag+np.r_[std*mean, std*mean],\n         'C1', label=r'Standard deviation $\\sigma$')\n\n\n# Random uniform phase\nuniform_mean = std * ((1+1j)*mean + np.exp(1j*np.linspace(0, 2*np.pi, 301)))\nax1.plot(d.real+uniform_mean.real, d.imag+uniform_mean.imag,\n         'C0', label='Random uniform phase')\n\n\n# Gaussian\nfor i in range(1, 3):\n    # Correlated\n    ax1.plot(d.real + np.r_[-std, +std]*i + std*mean,\n             d.imag + np.r_[-std, +std]*i + std*mean,\n             'C3-', lw=6-2*i,\n             label=f'Gaussian $\\\\pm {i} \\\\sigma$, correlated')\n\n    # Uncorrelated\n    ax1.plot(d.real + np.r_[-std, -std, +std, +std, -std]*i + std*mean,\n             d.imag + np.r_[-std, +std, +std, -std, -std]*i + std*mean,\n             'C2:', lw=6-2*i,\n             label=f'Gaussian $\\\\pm {i} \\\\sigma$, uncorrelated')\n\n\n# Plot random realizations\ndata = np.ones(300, dtype=complex)*d\nshape = data.shape\nrng = np.random.default_rng()\nls = ['C0x', 'C3+', 'C2x']\nntypes = ['white_noise', 'gaussian_correlated', 'gaussian_uncorrelated']\nfor i, ntype in enumerate(ntypes):\n\n    # Add random noise of ntype.\n    ndata = add_noise(data, ntype, relative_error, 0.0, mean)\n    ax2.plot(ndata.real, ndata.imag, ls[i], label=ntype)\n\n\n# Axis etc\nfor ax in axs:\n    ax.axhline(c='k')\n    ax.axvline(c='k')\n    ax.legend(framealpha=1, loc='upper left')\n    ax.set_aspect('equal')\n    ax.set_ylabel('Imaginary part')\n    ax.set_xlim([-4, 10])\nax2.set_xlabel('Real part')\n\n# fig.tight_layout()\nfig.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Intuitively one might think that the Gaussian uncorrelated noise is the\n\"best\" one, as it looks truly random. However, it is arguably the least\n\"physical\" one, as real and imaginary part of the electromagnetic field are\nnot independent - if one changes, the other changes too. The uniformly\ndistributed phases (blue circle) is the most physical noise corresponding to\nwhite noise adding random phases with a constant amplitude.\n\nTo get a better understanding we look at some numerical examples where we\nplot amplitude-vs-offset for a fixed frequency, and amplitude-vs-frequency\nfor a fixed offset; for single realizations and when we stack it many times\nin order to reduce the noise.\n\n## 3. Numerical examples\n\n### Model\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Model parameters\nmodel = {\n    'src': (0, 0, 0),  # Source at origin\n    'depth': [],       # Homogenous space\n    'res': 3,          # 3 Ohm.m\n    'ab': 11,          # Ex-source, Ex-receiver}\n}\n\n# Single offset and offsets\noffs = np.linspace(1000, 15000, 201)\noff = 5000\n\n# Single frequency and frequencies\nfreqs = np.logspace(-3, 2, 201)\nfreq = 1\n\n# Responses\noresp = empymod.dipole(\n    rec=(offs, offs*0, 0),  # Inline receivers\n    freqtime=freq,\n    **model\n)\nfresp = empymod.dipole(\n    rec=(5000, 0, 0),     # Inline receiver\n    freqtime=freqs,\n    **model,\n)\n\n# Relative error, noise floor, mean of noise\nrel_error = 0.05\nnoise_floor = 1e-15\nn_stack = 1000\n\n# Phase settings: wrapped, radians, lag-defined (+iw)\nphase = {'unwrap': False, 'deg': False, 'lag': True}"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "### Plotting function\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "def error(resp, noise):\n    \"\"\"Return relative error (%) of noise with respect to resp.\"\"\"\n    return 100*abs((noise-resp)/resp)\n\n\ndef figure(x, data, reim, comp):\n    fig, axs = plt.subplots(2, 4, constrained_layout=True,\n                            figsize=(14, 8), sharex=True)\n\n    axs[0, 0].set_title('|Real| (V/m)')\n    axs[0, 0].plot(x, abs(data.real), 'k')\n    axs[0, 0].plot(x, abs(reim.real), 'C0')\n    axs[0, 0].plot(x, abs(comp.real), 'C1--')\n    axs[0, 0].set_yscale('log')\n\n    axs[1, 0].plot(x, error(data.real, reim.real), 'C0')\n    axs[1, 0].plot(x, error(data.real, comp.real), 'C1--')\n    axs[1, 0].set_ylabel('Rel. Error (%)')\n\n    axs[0, 1].set_title('|Imaginary| (V/m)')\n    axs[0, 1].plot(x, abs(data.imag), 'k', label='Data')\n    axs[0, 1].plot(x, abs(reim.imag), 'C0', label='Noise to Re; Im')\n    axs[0, 1].plot(x, abs(comp.imag), 'C1--', label='Noise to Complex')\n    axs[0, 1].set_yscale('log')\n    axs[0, 1].legend(fontsize=12, framealpha=1)\n\n    axs[1, 1].plot(x, error(data.imag, reim.imag), 'C0')\n    axs[1, 1].plot(x, error(data.imag, comp.imag), 'C1--')\n\n    axs[0, 2].set_title('Amplitude (V/m)')\n    axs[0, 2].plot(x, data.amp(), 'k')\n    axs[0, 2].plot(x, reim.amp(), 'C0')\n    axs[0, 2].plot(x, comp.amp(), 'C1--')\n    axs[0, 2].set_yscale('log')\n\n    axs[1, 2].plot(x, error(data.amp(), reim.amp()), 'C0')\n    axs[1, 2].plot(x, error(data.amp(), comp.amp()), 'C1--')\n\n    axs[0, 3].set_title('Phase (rad)')\n    axs[0, 3].plot(x, data.pha(**phase), 'k')\n    axs[0, 3].plot(x, reim.pha(**phase), 'C0')\n    axs[0, 3].plot(x, comp.pha(**phase), 'C1--')\n\n    axs[1, 3].plot(x, error(data.pha(**phase), reim.pha(**phase)), 'C0')\n    axs[1, 3].plot(x, error(data.pha(**phase), comp.pha(**phase)), 'C1--')\n\n    return fig, axs"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "### 3.1 Offset-range for single frequency\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "def offset_single(mu):\n    \"\"\"Single frequency, many offsets, one realization.\"\"\"\n    # Add noise\n    inp = {'rel_error': rel_error, 'noise_floor': noise_floor, 'mu': mu}\n    onoise_reim = add_noise(oresp, 'gaussian_correlated', **inp)\n    onoise_comp = add_noise(oresp, 'white_noise', **inp)\n\n    fig, axs = figure(offs, oresp, onoise_reim, onoise_comp)\n    fig.suptitle(f\"Inline $E_{{xx}}$; $s_z=r_z=0$; $f=${freq} Hz; \"\n                 f\"fullspace of {model['res']} $\\\\Omega\\\\,$m; $\\\\mu=${mu}\",\n                 fontsize=20)\n\n    for i in range(3):\n        axs[0, i].set_ylim([1e-19, 3e-10])\n\n    for i in range(4):\n        axs[1, i].set_xlabel('Offset (m)')\n        axs[1, i].set_yscale('log')\n        axs[1, i].set_ylim([1e-2, 1e6])"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "offset_single(mu=0.0)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "offset_single(mu=0.5)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "### 3.2 Offset-range for single frequency - STACKED\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "def offset_stacked(mu):\n    \"\"\"Single frequency, many offsets, stacked.\"\"\"\n    # Stack noise\n    inp = {'rel_error': rel_error, 'noise_floor': noise_floor, 'mu': mu}\n    sonoise_reim = stack(n_stack, oresp, 'gaussian_correlated', **inp)\n    sonoise_comp = stack(n_stack, oresp, 'white_noise', **inp)\n\n    fig, axs = figure(offs, oresp, sonoise_reim, sonoise_comp)\n    fig.suptitle(f\"STACKED {n_stack} times; $\\\\mu=${mu}\", fontsize=20)\n\n    for i in range(3):\n        axs[0, i].set_ylim([1e-19, 3e-10])\n\n    for i in range(4):\n        axs[1, i].set_xlabel('Offset (m)')\n        axs[1, i].set_ylim([-5, 40])"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "offset_stacked(mu=0.0)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "offset_stacked(mu=0.5)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "### 3.3 Frequency-range for single offset\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "def frequency_single(mu):\n    \"\"\"Single offset, many frequencies, one realization.\"\"\"\n    # Add noise\n    inp = {'rel_error': rel_error, 'noise_floor': noise_floor, 'mu': mu}\n    fnoise_reim = add_noise(fresp, 'gaussian_correlated', **inp)\n    fnoise_comp = add_noise(fresp, 'white_noise', **inp)\n\n    fig, axs = figure(freqs, fresp, fnoise_reim, fnoise_comp)\n    fig.suptitle(f\"Inline $E_{{xx}}$; $s_z=r_z=0$; offset$=${off/1e3} km; \"\n                 f\"fullspace of {model['res']} $\\\\Omega\\\\,$m; $\\\\mu=${mu}\",\n                 fontsize=20)\n\n    for i in range(3):\n        axs[0, i].set_ylim([1e-18, 1e-11])\n\n    for i in range(4):\n        axs[0, i].set_xscale('log')\n        axs[1, i].set_xlabel('Frequency (Hz)')\n        axs[1, i].set_yscale('log')\n        axs[1, i].set_ylim([1e-1, 1e5])"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "frequency_single(mu=0.0)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "frequency_single(mu=0.5)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "### 3.4 Frequency-range for single offset - STACKED\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "def frequency_stacked(mu):\n    \"\"\"Single offset, many frequencies, stacked.\"\"\"\n    # Stack noise\n    inp = {'rel_error': rel_error, 'noise_floor': noise_floor, 'mu': mu}\n    sfnoise_reim = stack(n_stack, fresp, 'gaussian_correlated', **inp)\n    sfnoise_comp = stack(n_stack, fresp, 'white_noise', **inp)\n\n    fig, axs = figure(freqs, fresp, sfnoise_reim, sfnoise_comp)\n    fig.suptitle(f\"STACKED {n_stack} times; $\\\\mu=${mu}\", fontsize=20)\n\n    for i in range(3):\n        axs[0, i].set_ylim([1e-18, 3e-11])\n\n    for i in range(4):\n        axs[0, i].set_xscale('log')\n        axs[1, i].set_xlabel('Frequency (Hz)')\n        axs[1, i].set_ylim([-5, 40])"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "frequency_stacked(mu=0.0)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "frequency_stacked(mu=0.5)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "empymod.Report()"
      ]
    }
  ],
  "metadata": {
    "kernelspec": {
      "display_name": "Python 3",
      "language": "python",
      "name": "python3"
    },
    "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.8.6"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}