{ "cells": [ { "cell_type": "markdown", "metadata": { "id": "EPtqtPIUvdc2" }, "source": [ "# FMCW Radar 101 - Intro\n", "[![](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/matt-chv/mmWrt/blob/main/docs/FMCW-Radar-101_Intro.ipynb)\n", "\n", "Goal:\n", "\n", "> Explicit in one single place code and maths needed to understand FMCW radar.\n", "\n", "Status:\n", "\n", "* this work book runs in any browser and correctly computes distance and speed for a single target.\n", "\n", "Howto:\n", "\n", "> you can start with the maths or with the code, you should finish the otherway round and check they both match.\n", "\n", "1. if needed check the introduction to Google Colab to run notebooks in your web browser [here](https://colab.research.google.com/)\n", "\n", "2. Navigate to the cell you want to execute and press execute\n", "3. make changes to chirp ramp-time, distances, speeds and observe the changes (keep in mind that all measures are related to the FFT bin sizes).\n", "\n", "Next:\n", "\n", "* click to access other workbooks:\n", " * [FMCW 102 - CFAR](https://colab.research.google.com/gist/matt-chv/33e98a23d4b9d90dd27c1bf7f0a54781/fmcw-radar-102-cfar.ipynb) : CFAR or how to detect objects of interest from range FFT.\n", " * [FMCW 103 - AoA](https://colab.research.google.com/gist/matt-chv/d81f7e2166009a623a36781a0773ae47/fmcw-radar-103-aoa.ipynb) : angle of arrival (CAPON vs Bartlett)\n", " * [FMCW 104 - increased resolution vs FFT bin](https://colab.research.google.com/gist/matt-chv/0b25dbc4673f2d7d63804cc6241643b9/fmcw-radar-104-1-fft-freq-estimation.ipynb)increase accuracy option compared to standard FFT.\n", " * Also available on github as gist for forking:\n", " * [fmcw 101](https://gist.github.com/matt-chv/bdd8b835c5cb7e739bb8b68d00257690)\n", " * [fmcw 102](https://gist.github.com/matt-chv/33e98a23d4b9d90dd27c1bf7f0a54781)\n", " * ...\n", "\n", "History:\n", "\n", "* 2022-Dec-28: Added 2D FFT - 3 targets - 2 targets same range, 2 same speed\n", "* 2022-Dec-23: Added 2D FFT example - 1 targets\n", "* 2022-Dec-15: Clean-up to have MRE for range and speed as standalone cells\n", "* 2022-Apr-28: creation\n", "\n", "Related ressources:\n", "\n", "* How to run a google Colab notebook (requires a Google account to login): [Colab intro](https://colab.research.google.com/)\n", "* OpenRadar notebook on range: [Range example](https://github.com/PreSenseRadar/OpenRadar/blob/master/Presense%20Applied%20Radar/basics/Range%20-%20COMPLETED.ipynb)\n", "\n", "\n", "\n" ] }, { "cell_type": "markdown", "metadata": { "id": "oV8JqYMFNYd-" }, "source": [ "## FMCW Maths 101" ] }, { "cell_type": "markdown", "metadata": { "id": "kxWYcEQHig7p" }, "source": [ "\n", "### linear chirps and IF frequencies\n", "\n", "Signal model\n", "According to (Barriok 1973; Stove 1992; Komarov and Smolskiy 2003; Winkler 2007) the\n", "transmitted signal of an FMCW radar system can be modeled as\n", "where\n", "\n", "$$ y_T (t) = A_T \\cdot cos\\left( 2 \\pi \\cdot ( f_{0min} \\cdot t +\\int_0^tf_T({\\tau}) d\\tau)\\right)$$\n", "\n", "Given for a linear chirp that\n", "$$ f_T(\\tau) = \\frac{B}{T} \\cdot \\tau = s \\cdot \\tau $$\n", "We derive the phase, given\n", "$$ \\int_0^t f_T({\\tau}) d\\tau = \\frac{s}{2} \\cdot \\tau^2 - K$$\n", "\n", "which can be written as:\n", "$$ y_T (t) = A_T \\cdot cos \\left( 2 \\pi \\cdot ( f_{0min} \\cdot t + \\frac{s}{2} \\cdot t^2 )+\\Phi_0 \\right) $$\n", "\n", "Where:\n", "* $f_{0min}$ is the start frequency at the begining of the raising frequency of the chirp.\n", "* s is the slope at which the frequency is ramped ( $S = \\frac{B}{T}$)\n", "* B is the total bandwdith of the chirp\n", "* T is the total time of the chirp\n", "\n", "Considering a reflected signal with a time delay $\\delta = 2 · \\frac{R0+ v\\cdot t}{c}$ and Doppler shift $f_D = −2 · \\frac{f_c \\cdot v}{c}$\n", "\n", "Where:\n", "\n", "* c is the speed of light\n", "* $f_D$ is the doppler shift\n", "* R0 is the nominal distance to the target\n", "* $\\Delta$ is the time of flight (to and from the target)\n", "* v is the velocity of the target\n", "\n", "The receive signal $y_R(t)$ can be written as :\n", "\n", "$$y_R(t) = A_R \\cdot cos ( 2\\pi \\cdot (f_{0min} + \\frac{s}{2} \\cdot (t-\\delta)) \\cdot (t-\\delta))$$\n", "\n", "\n", "$y_{IF}(t)$ is the IF signal (after mixer) which is obtained by multiplication in the time domain, and passed to a low-pass filter (LPF)\n", "\n", "This can be done easily when remembering the trigonometric relation:\n", "\n", "$$cos(\\alpha) \\cdot cos(\\beta) = \\frac{cos(\\alpha + \\beta) + cos(\\alpha - \\beta)}{2}$$\n", "\n", "$$y_{MIX}(t) = y_R(t) \\cdot y_T(t)$$\n", "\n", "Noticing that the element which sums the elements will be higher frequency and will be filtered by the LPF, it remains that:\n", "\n", "$$y_{IF}(t) = \\frac{A_t \\cdot A_r}{2} \\cdot cos(2 \\pi \\cdot [f_{0min} \\cdot \\delta + s \\cdot \\delta \\cdot t - \\frac{s}{2} \\cdot \\delta^2])$$\n", "\n", "Where:\n", "\n", "* $f_{0min}$ the starting frequency of the chirp\n", "* s is the slope of the chirp\n", "* $\\delta$ is the total time of flight between antennas and target\n", "* At, Ar: Amplitude of the RX and TX waves\n", "\n", "## phase information\n", "\n", "$$ y_{IF}(t) = \\frac{A_t \\cdot A_r}{2} \\cdot cos(2 \\pi \\cdot [f_{0min} \\cdot \\delta + s \\cdot \\delta \\cdot t - \\frac{s}{2} \\cdot \\delta^2])$$\n", "\n", "could also be written as\n", "\n", "$$ y_{IF}(t) = A_{IF} \\cdot e^{j \\cdot \\Phi(t)} $$\n", "\n", "with\n", "$$ \\Phi(t) = 2 \\cdot \\pi \\cdot [f_{0min} \\cdot \\delta + s \\cdot \\delta \\cdot t - \\frac{s}{2} \\cdot \\delta^2] $$\n", "\n", "let's say $$ \\Phi_0 = \\Phi(\\tau) $$ and\n", "$\\Delta t$ later\n", "$$ \\Phi_1 = \\Phi(\\tau+\\Delta t) $$\n", "\n", "$$ \\Delta \\Phi = \\Phi_1 - \\Phi_0$$\n", "$$ \\iff $$\n", "$$ \\Delta \\Phi = \\Phi_1 - \\Phi_0 = f_0⋅\\Delta t−K \\cdot \\Delta t−K\\tau \\cdot \\Delta t-\\frac{K}{2}\\Delta t^2 $$\n", "\n", "which can be simplified as\n", "$$ \\Delta \\Phi = 2 \\cdot \\pi \\cdot f \\cdot \\Delta t $$\n", "with $c= \\lambda \\cdot f$\n", "\n", "$$ \\Delta \\Phi = 4 \\cdot \\pi \\cdot \\frac{v \\cdot T_c}{\\lambda} $$\n", "\n", "More ressources:\n", "* Eq 8 [FMCW Radar system](https://uwaterloo.ca/centre-for-intelligent-antenna-and-radio-systems/sites/ca.centre-for-intelligent-antenna-and-radio-systems/files/uploads/files/fmcwradarsystem.pdf)\n", "* Eq 7, 10 [spyy005 on ti.com](https://www.ti.com/lit/wp/spyy005a/spyy005a.pdf)\n", "* Eq 1, 2, 3, 4 [Design of an FMCW radar baseband\n", "signal processing system for automotive\n", "application](https://link.springer.com/content/pdf/10.1186/s40064-015-1583-5.pdf)\n", "* Eq 2.7 in [Object Detection with\n", "Automotive Radar\n", "Sensors using CFARAlgorithms](https://www.jku.at/fileadmin/gruppen/183/Docs/Finished_Theses/Bachelor_Thesis_Katzlberger_final.pdf)" ] }, { "cell_type": "markdown", "metadata": { "id": "BbRdF5nfR_Kh" }, "source": [ "### Distance from IF FFT" ] }, { "cell_type": "markdown", "metadata": { "id": "6quZUR4fSCJx" }, "source": [ "### Speed from IF FFT" ] }, { "cell_type": "markdown", "metadata": { "id": "5o6JdjegAhC1" }, "source": [ "## Minimum Reproductible Examples\n", "\n", "Below is code generating IF for single target and using the IF to comupte distance and speed.\n", "\n" ] }, { "cell_type": "markdown", "metadata": { "id": "WiJNA2I2Amt_" }, "source": [ "## Distance" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "colab": { "base_uri": "https://localhost:8080/" }, "id": "6NJQ_925A4DN", "outputId": "4378e4c9-d0ae-4d61-c746-6dff79f0fa30" }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "target at:12 is computed to be at 12\n", "target at:21 is computed to be at 21\n", "target at:40 is computed to be at 40\n" ] } ], "source": [ "from numpy import abs ,angle, arange, arcsin, cos, linspace, pi, sqrt, tan\n", "from scipy.fft import fft\n", "\n", "def y_IF(f0_min, slope, T, antenna_tx, antenna_rx, target, v=3e8):\n", " \"\"\" This function implements the mathematical IF defined in latex as\n", " y_{IF} = cos(2 \\pi [f_0\\delta + s * \\delta * t - s/2* \\delta^2])\n", " into following python code\n", " y_IF = cos (2*pi*(f_0 * delta + slope * delta * T - slope/2 * delta**2))\n", " Parameters:\n", " -----------\n", " f0_min: float\n", " the frequency at the begining of the chirp\n", " slope: float\n", " the slope with which the chirp frequency inceases over time\n", " T: ndarray\n", " the 1D vector containing time values\n", " antenna_tx: tuple of floats\n", " x, y, z coordinates\n", " antenna_rx: tuple of floats\n", " x, y, z coordinates\n", " target: tuple of floats\n", " x, y, z coordinates\n", " v: float\n", " speed of light in considered medium\n", " Returns:\n", " --------\n", " YIF: ndarray\n", " vector containing the IF values\n", " \"\"\"\n", " tx_x, tx_y, tx_z = antenna_tx\n", " rx_x, rx_y, rx_z = antenna_rx\n", " t_x, t_y, t_z = target\n", " # distance tx antenna to target\n", " distance = sqrt((tx_x-t_x)**2 + (tx_y-t_y)**2 + (tx_z-t_z)**2)\n", " # distance target to rx antenna\n", " distance += sqrt((rx_x-t_x)**2 + (rx_y-t_y)**2 + (rx_z-t_z)**2)\n", " # usually delta_t = 2*d/c, but\n", " # distance is already 2*D (TX to target + distance target to RX)\n", " # so delta = distance/v\n", " delta = distance/v\n", " YIF = cos(2 *pi *(f0_min * delta + slope * delta * T - slope/2 * delta**2))\n", " return YIF\n", "\n", "f0_min = 60e9\n", "c = 3e8\n", "# lambda ~5mm at 60GHz\n", "lambda0_max = 3e8/f0_min\n", "n_rx = 2\n", "Distance = 10\n", "k = 10e12\n", "n_samples = 512\n", "f_if = 2*k*Distance/c\n", "fs = 50e6\n", "ts = 1/fs\n", "\n", "# antenna_tx = (-lambda0_max/2,0,0)\n", "antenna_tx = (0,0,0)\n", "antenna_rx = (0,0,0)\n", "T = arange(0, n_samples*ts+ts, ts)\n", "\n", "for d in [12, 21, 40]:\n", " target = (0, d, 0)\n", "\n", " # sanity check\n", " f_if = 2*k*d/c\n", " assert f_if < fs/2\n", "\n", " YIF = y_IF(f0_min, k, T, antenna_tx, antenna_rx, target)\n", " FT = fft(YIF)\n", " MAG = abs(FT)[0: n_samples//2]\n", " ANG = angle(FT)[0: n_samples//2]\n", "\n", " # now find the peak\n", " amplitude_peak = sorted(MAG, reverse = True)[0]\n", " i_peak = list(MAG).index(amplitude_peak)\n", " # max un-ambiguous speed is fs*c/2/k\n", " distances = linspace(0, fs*c/2/k, n_samples)\n", " d_calc = distances[i_peak]\n", " print(f\"target at:{d:.2g} is computed to be at {d_calc:.2g}\")" ] }, { "cell_type": "markdown", "metadata": { "id": "WiTVHL7RAoBq" }, "source": [ "## Speed w/ phase (no FFT)" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "colab": { "base_uri": "https://localhost:8080/" }, "id": "yEyC0CprDZPd", "outputId": "9dcc61fc-1e74-478d-d76c-eee0a59a988d" }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "if we computed speed with FFT:\n", "speed resolution: 2.1e+03\n", "v_max: 1e+03\n", "speed is:1e+02 is computed to be at 1.1e+02\n", "speed is:2e+02 is computed to be at 2.2e+02\n", "speed is:3e+02 is computed to be at 3.4e+02\n" ] } ], "source": [ "from numpy import abs ,angle, arange, arcsin, cos, pi, sqrt, tan\n", "from scipy.fft import fft\n", "\n", "def y_IF(f0_min, slope, T, antenna_tx, antenna_rx, target, v=3e8):\n", " \"\"\" This function implements the mathematical IF defined in latex as\n", " y_{IF} = cos(2 \\pi [f_0\\delta + s * \\delta * t - s/2* \\delta^2])\n", " into following python code\n", " y_IF = cos (2*pi*(f_0 * delta + slope * delta * T - slope/2 * delta**2))\n", " Parameters:\n", " -----------\n", " f0_min: float\n", " the frequency at the begining of the chirp\n", " slope: float\n", " the slope with which the chirp frequency inceases over time\n", " T: ndarray\n", " the 1D vector containing time values\n", " antenna_tx: tuple of floats\n", " x, y, z coordinates\n", " antenna_rx: tuple of floats\n", " x, y, z coordinates\n", " target: tuple of floats\n", " x, y, z coordinates\n", " v: float\n", " speed of light in considered medium\n", " Returns:\n", " --------\n", " YIF: ndarray\n", " vector containing the IF values\n", " \"\"\"\n", " tx_x, tx_y, tx_z = antenna_tx\n", " rx_x, rx_y, rx_z = antenna_rx\n", " t_x, t_y, t_z = target\n", " # distance tx antenna to target\n", " distance = sqrt((tx_x-t_x)**2 + (tx_y-t_y)**2 + (tx_z-t_z)**2)\n", " # distance target to rx antenna\n", " distance += sqrt((rx_x-t_x)**2 + (rx_y-t_y)**2 + (rx_z-t_z)**2)\n", " # usually delta_t = 2*d/c, but\n", " # distance is already 2*D (TX to target + distance target to RX)\n", " # so delta = distance/v\n", " delta = distance/v\n", " YIF = cos(2 *pi *(f0_min * delta + slope * delta * T - slope/2 * delta**2))\n", " return YIF\n", "\n", "f0_min = 60e9\n", "c = 3e8\n", "# lambda ~5mm at 60GHz\n", "lambda0_max = 3e8/f0_min\n", "n_rx = 2\n", "Distance = 10\n", "k = 10e12\n", "n_samples = 512\n", "f_if = 2*k*Distance/c\n", "fs = 50e6\n", "ts = 1/fs\n", "\n", "antenna_tx = (0,0,0)\n", "antenna_rx = (0,0,0)\n", "T = arange(0, n_samples*ts+ts, ts)\n", "t_chirp_to_chirp = 1.2e-6\n", "print(\"if we computed speed with FFT:\")\n", "print(f\"speed resolution: {lambda0_max/(2*t_chirp_to_chirp):.2g}\")\n", "print(f\"v_max: {lambda0_max/(4*t_chirp_to_chirp):.2g}\")\n", "\n", "for d in [10]:\n", " for v in [100, 200, 300]:\n", " target_t0 = (d, 0, 0)\n", " f_if = 2*k*d/c\n", " # sanity check\n", " assert f_if < 1/ts/2\n", "\n", " YIF0 = y_IF(f0_min, k, T, antenna_tx, antenna_rx, target_t0)\n", "\n", " target_t1 = (d+v*t_chirp_to_chirp, 0, 0)\n", " YIF1 = y_IF(f0_min, k, T, antenna_tx, antenna_rx, target_t1)\n", "\n", " # since we have a real FFt, only want to see the first peak\n", " FT0 = fft(YIF0[:n_samples//2])\n", " FT1 = fft(YIF1[:n_samples//2])\n", " MAG0 = abs(FT0)\n", " amplitude_peak0 = sorted(MAG0, reverse = True)[0]\n", " i_peak0 = list(MAG0).index(amplitude_peak0)\n", "\n", " MAG1 = abs(FT1)\n", " amplitude_peak1 = sorted(MAG1, reverse = True)[0]\n", " i_peak1 = list(MAG1).index(amplitude_peak1)\n", " assert i_peak0 == i_peak1\n", "\n", " ANG0 = angle(FT0)\n", " ANG1 = angle(FT1)\n", "\n", " ph0 = ANG0[i_peak0]\n", " ph1 = ANG1[i_peak1]\n", " # d, v, ph0 10 100.0 -1.5646604036459237 -1.2627114022166515\n", " v_est = lambda0_max*(ph1-ph0)/(4*pi*t_chirp_to_chirp)\n", " print(f\"speed is:{v:.2g} is computed to be at {v_est:.2g}\")" ] }, { "cell_type": "markdown", "metadata": { "id": "C0xtKfC_FbkJ" }, "source": [ "## Range-Doppler FFT (2D FFT)" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "colab": { "base_uri": "https://localhost:8080/", "height": 261 }, "id": "Y8lC-Q4GFeQA", "outputId": "96a4f77e-02ba-4e10-e2ad-a04f3da3b60c" }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "512\n" ] }, { "data": { "text/plain": [ "Text(0.5, 1.0, 'Velocity-Range 2D FFT')" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkYAAADRCAYAAAApIX2+AAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjMsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvZiW1igAAAAlwSFlzAAAPYQAAD2EBqD+naQAAQvlJREFUeJztnQWYFWX7xp8NYgmlkW4p6ZIQKQmRMkBKUNK/tCIgEoKCoIgISH0S0iEg4gcISLc00giS0iq1sDH/637YOd+cs2eXPVsn9v5d18CemTlz3vfEzD1P+hmGYQghhBBCCBF/dw+AEEIIIcRToDAihBBCCImAwogQQgghJAIKI0IIIYSQCCiMCCGEEEIioDAihBBCCImAwogQQgghJAIKI0IIIYSQCCiMCCGEEEIioDAixEfx8/OToUOHJsixz507p8efOXNmghyfEELcBYURIW6mcePGkipVKrlz506U+7Ru3VqSJ08uN2/eFE/lv//9b4IIMRwTIsxckiVLJnnz5pUePXrI33//Lb7C0qVLpUWLFpI/f379PhQuXFjef/99p3O0vh+BgYGSIUMGKVeunPTs2VOOHj0a49fE+2g9lnUJDg7WfSB+o9qnf//+UqNGjSi3W5eEEumExDeB8X5EQohLQPT89NNPsmzZMnnrrbcibb9//778+OOPUr9+fcmYMaN4Anny5JEHDx6oSLEKo4kTJybYBXDSpEmSJk0auXfvnqxfv17Gjx8v+/btk61bt4ov0LlzZ8mePbu0adNGcufOLYcPH5YJEybo+4p5BgUF2e3/0ksv6fcF7S7/+ecfOXjwoMyaNUu+/fZbGTVqlPTp0ydGr1u6dGkVYI5AiFsZNmyY5MuXz27dc889J7Vr15aOHTva1u3Zs0e++eYb+eijj6Ro0aK29SVLlozxe0GIW0ETWUKI+7h//76RNm1ao169ek63z5s3D42ejQULFrh0XDxnyJAhRmLx3nvv6WvGN5gDjnv9+nW79S1atND1u3btMnyBDRs2RFo3a9YsneO0adPs1mMd3m9Hbty4YVSuXFm3//zzz098zTx58hgNGzaMdp8ZM2bo8fbs2ROjeSxevFj3dzYfQrwButIIcTOwBLz66qtqBbl27Vqk7fPmzZO0adOqyw3AtdKrVy/JlSuXpEiRQgoWLKgWgvDw8Ce+1v79+6VBgwby1FNPqfUFd/s7d+6MtB9eo3fv3upqwWvkzJlTrRM3btxwGmPUvn17tRYBq/sE13Aco0mTJpFeA66ap59+Wrp06RKLd03khRde0P/PnDljW3fr1i354IMPpESJEjo/zBPzhTXFysaNG3V8ixYtks8++0znlzJlSn0/Tp8+Hem1MDe4uPBZVaxYUbZs2aIuJCxWHj58KEOGDNHPBO8bPqMPP/xQ1z8Jx2OBZs2a6f/Hjh2L0XsCi+KCBQvUvYZ5EUJch640QjzEnQY3CC7U3bp1s7vQr1mzRlq2bKkXZbjVXnzxRbl06ZIKCrhctm/fLgMGDJArV67I119/HeVr/P777yomIBZwsYYbbMqUKXpB3rRpk1SqVEn3u3v3ru6Hi/E777wjZcuWVUG0YsUKuXjxomTKlCnSsTGWy5cvy9q1a2X27Nm29RAfcA2NHj1a54JYGBO4D//991/dHhsgzkD69Olt6/744w9Zvny5vPHGG+r2uXr1qs4R7xlib+CqsvL555+Lv7+/iim4ozBOfBa7du2yc+HhM8F7ArGI123atKm+LgSVCYQpxCtce3CLwY0Ed9jYsWPl5MmTOi5X+euvv/R/Z+95VOA7gflu2LBB31983tEREhJiE7wmiHHCYgXvj+N+royLEK/B3SYrQohhhIaGGtmyZVM3iJXJkyerW2LNmjX6ePjw4Ubq1KmNkydP2u3Xv39/IyAgwDh//nyUrrSmTZsayZMnN86cOWNbd/nyZXXjVa9e3bZu8ODB+tylS5dGGmd4eLj+f/bsWd0HbpYnudJOnDih6ydNmmS3vnHjxkbevHltx3ySKw3HgTvt3LlzxvTp042goCAjc+bMxr1792z7BgcHG2FhYXbPx1hTpEhhDBs2zLYObh4cs2jRosbDhw9t68eNG6frDx8+rI+xLWPGjEaFChWMkJAQ234zZ87U/V588UXbutmzZxv+/v7Gli1bnH6G27ZtM1ylQ4cO+rk6ft5RudJMevbsqfscPHjwia407Oe4WL83pivN2eIMutKIt0NXGiEeQEBAgLz55puyY8cOmyXEdKNlzZpVXTxg8eLFarmAtQJ37+ZSp04dCQsLk82bNzs9Prb98ssvaumAS8gkW7Zs0qpVK7VywLoAfvjhBylVqpTNjWMFFiBXefbZZ9UaNXfuXNs6WI9WrVql1pmYHhNZWpkzZ1bXHCxZcFfhGFbLBtxXsACZc0YWH1xqeC4CmB15++237YKMTfccLE/gt99+02N06tRJ3VMmGLfVUmV+NrASFSlSxO6zqVWrlm6HBccV8Nl/9913GhhdqFAhl56LOYPoMh1N8NnA0mddnCUBwJ3ouB8hvghdaYR4CLjYwu2CCyIyeuC2QiwL0tIhnMCpU6fk0KFDKhCc4SxGCVy/fl3dcBAIjuBiDjfQhQsXpHjx4hqz89prr8Xr3HChhTvqzz//1Iw2iAi4cNq2bavbHz16pGLJCuZoztsUbHALYS7Iejp79mykTC3MY9y4cZqZhe0QRybOMvrgdrJiip3bt2/r/xgvgAizApEEgWYFnw3cj65+Ns7A596hQwepV69erGKF4A4FiE17EnCHQVg/CcRWlS9f3uWxEOJtUBgR4iGgDg2sDfPnz1dhhP/hNYFgsl74kaaNGKGorDOeCKxhiM+B1QhzmzNnjl5kTaGGOKmaNWvaPQfCxio+qlevbotpadSokQZY473Zu3evzUo0YsQIGTRokFqUhg8frjFN2IZgdWfB6VbhZeWxt8o1cHyM6auvvnK6HYHYMQGB4ohVQir8kiVL7CxVMeXIkSM6N8f0ekLIk6EwIsSDwIUeF3ZYhWA5ggulQoUKtu0FChRQa0BM7vCtwIoBl9OJEycibTt+/LiKB/PCjdfAhdVVonOJQaA0bNhQhRHmuG3bNrtAcbjuHF0zzzzzTLSuImR/wRWGgHUILwAhAYEFF5Rjll1sAoVh3QLIVLMKt9DQUHV5Wmvz4H2DqIHbMzYuRwBrHepVZcmSResXmS4xVzh//rwG01euXDlGFiNCiD2MMSLEgzCtQ4MHD5YDBw7YWYtA8+bNNQ4JmWqO4OKPC7YzYD2oW7euFoq0xjAhawsCrFq1arbsJbjRcIFHwUlXLCmpU6e2jcMZcJshM6xv3762mCqrCwtiz7ogfT468N4gKwylCqzzdBwj3HbI4osNsGrBBTdt2jS79xYCz3S3WT8bvA72dQTFMFGY8kkZaPiMIFLx+UblkosOuCORwQgX4sCBA11+PiGEFiNCPAq4PqpUqaICBjgKI4gKpM2/8sorWjsI7jdccJEWDmsJRE9UlpFPP/1UrTIQQf/3f/+nLhqksqPGDtLUra+BYyHlHS4pvAYuuHjdyZMnq3XHGdgPICYKsTGO4gcWI4gMCBXUFoJVJC6g3ABaYGC8q1evVksL3hdUaIYlCe8j3heIGGvAuSsgMBuVvLt3765B1BA/eI9RvwkWIqtlCMIP1quuXbtqoHXVqlVVoMAih/UQO9HF6GD8CPqGmxTB8NaK3gjAhwvVCkoAwCUJIYjAeYhZvLewKMKdh+MRQmKBu9PiCCH2TJw4UdOdK1as6HT7nTt3jAEDBhgFCxbU9PtMmTIZVapUMb788kvj0aNH0Va+3rdvn1bYTpMmjZEqVSqjZs2axvbt2yO9xs2bN41u3boZOXLk0NfImTOn0a5dO62sHFW6PkoOdO/eXVPo/fz8nKZz/9///Z+uRzXvuFa+Bv/884/x9NNP29Lmka7//vvva+kDpPNXrVrV2LFjh263ptab6fpILbfibF7gm2++0dR2pP3jc0Hqfbly5Yz69evb7Yf3f9SoUUbx4sV13/Tp0+t+n3zyiY41OqJKiXcsC+C4L0oEpEuXzihTpoym6f/+++8xfm9Z+ZqQyPjhn9gIKkIIcRUEYCP+B24jxwKC3gQCreHqQsVyZ64zQoj3whgjQkiigBYgcP0ghsmbRBHG7Xj/+P3336t70VkbD0KId8MYI0JIgoL6PevWrdO4JRRLRFyQN4FecrB0IeYKMVIoFAmrF9LpsY4Q4ltQGBFCEhRkoiGIHMHWKMxYunRp8SZQSwmlDDB2s98bClaiz5q1ajYhxDfw6RgjlLD/4osvNJ4BmTTjx4/X6q2EEEIIIUkqxmjhwoXSp08fLQIH0zeEEVKIXSnLTwghhJCkhc9ajNAYERWDJ0yYYMsigTkc9Uj69+/v7uERQgghxAPxyRgjNKRE/6QBAwbY1qGaLKrpomqwM1DkDosJhBTiCRBsGdvy/oQQQghJXGDvuXPnjmTPnt3WR1GSujC6ceOGVpxFtVgreIwqtM4YOXKkfPLJJ4k0QkIIIYQkJBcuXNC2Qa7ik8IoNsC6hJgkk3/++Udy584t1eRlCZRkbh0beQK4I3DSOZ0QQhIEnnM8mlAJka3y31g3UfZJYYReUejThAaZVvA4qo7dKVKk0MURiKJAPwojT8bPP1AMCYP91N1DIYQkBfz8RfwojDyWiEtBbMNgfDIrDbVF0NBy/fr1djFDeFy5cmW3jo0QQgghnotPWowA3GLt2rXTbtaoXfT1119rF3J03SaEEEIISVLCqEWLFnL9+nUZPHiwFnhEtd3Vq1dHCsgmhBBCCPH5OkZx5d9//5Wnn35aakgTxhh5OH4BgWKEM8aIEJJIMPjaowk1QmSj/KhJVE899ZTLz/fJGCNCCCGEkNhAYUQIIYQQEgGFEfF+WJmcEEKILwojVJ9GfzMUZcqSJYs0bdpUTpw4YbdPly5dpECBAhIUFCSZM2eWJk2aRKpm3aNHD03XR10iBF0TQgghhHidMNq0aZO89957snPnTlm7dq2EhIRI3bp1Nc3eBIJnxowZcuzYMVmzZo32RME+aAFi5Z133tHMNEIIIYQQn8hKQ7o9LEcQTNWrV3e6z6FDh6RUqVJy+vRptSRZGTp0qCxfvlwOHDjg8mszK81L8PMTv8BAMUJDmZVGCEkcmJXm01lpHl3HCJMCGTJkcLodliRYj/Llyye5cuWK02s9fPhQF6swIl4SX8QYI0IIIb7oSrOCFh69evWSqlWrynPPPWe37dtvv5U0adLosmrVKnW7oQ1IXOObYCEyl7gKLUIIIYR4Hx4rjBBrdOTIEVmwYEGkba1bt5b9+/eri+3ZZ5+V5s2bS3BwcJxeb8CAAWqhMpcLFy7E6XgkMaHFiBCSiNBt79O47EqDu2nXrl3y559/yv379zUzrEyZMurOii+6desmK1eulM2bN0vOnDkjbTetOoUKFZLnn39e0qdPL8uWLZOWLVvG+jWRwYaFeB9+/n5iwJ3GkxUhhJDEEkbbtm2TcePGyU8//aTZYhAmSJm/deuWiqX8+fNL586dpWvXrppuHxsQB969e3cVORs3boyR2MJzsFjjg0gSw4wxojgihBCSGK60xo0ba+p73rx55ZdffpE7d+7IzZs35eLFi2o1OnXqlHz88ceyfv16dW0h5ie27rM5c+bIvHnzVFyh+SuWBw8e6PY//vhDY4H27t0r58+fl+3bt8sbb7yhAu3ll1+2HQcZashEM5+Lv7E8evQoVuMing5daYQQQhIxXX/KlClaFyhZsienrR89elSuXLkitWvXdn0wUWQXIfOsffv2cvnyZenYsaMKo9u3b0vWrFk1jX/w4MFSuHBh2/41atTQ+CNHzp49q+IuJjBd30vw9xf/oCAJR4wZ0mdpMSKEJDS0Tvt0ur5H1zFyJxRGXnJygjBKlUqM+w/ECA/jyYoQkvBQGPm0MHI5Kw3ZWnChmezevVvT6qdOneryixMSV2Bl9Avw/583jTWNCCGExAGXhVGrVq1kw4YN+jdieF566SUVRwMHDpRhw4bFZSyExC7gOiAQqWkURYQQQhJfGKG2UMWKFfXvRYsWafFFBEHPnTtXZs6cKfHJ559/rhYBWKSs7NixQ2rVqiWpU6dWMxnijMwA7XPnzkmHDh00ow1B2WgTMmTIEAZe+wp+DqIIpfmTJxNRq5FlPSGEEJIYdYyQqm/W+1m3bp1mrIEiRYpo0HV8sWfPHg36LlmyZCRRVL9+fS3IOH78eAkMDJSDBw+KPy6QInL8+HGtmo3nFixYUIVcp06dtH3Il19+GW/jI4mMo9iBC83PX/wCAkRSpRS/fwNEQkMf1zNyhLEAhBBCEkoYFS9eXCZPniwNGzbUtPzhw4fremSMZcyYUeKDu3fvanXradOmyaeffmq3rXfv3tKjRw/p37+/bZ01Iw2iCYsJ6iudOHFCJk2aRGHkDaiuicLiY7EIQRQJRFHKFPIoa1pJ/s+dx1lpEEfW5o7QRFEZkCiYCCGExNWVNmrUKLXGICUelabR2R6sWLHC5mKLK6hnBOFVp04du/XXrl3TqttZsmSRKlWqaLr+iy++KFu3bo32eIhMj6oRrQkKRCITzbp4ZTNVdy6w2sVkgaBxWPwCAsUvMPB//2NJlkz8kiUXv+QpxD9lSk3L90+dWgKeSiv+GdKJX7bMkqp0BunXd4f8WymHSI6s4p8xgwQ8/ZT4p0n9eP+gFOKXIoX4JU/++HiBkV/H2XhiPBcs7n7fCSGEJL7FCIUcU6VKpYLoxo0bKhzQisMEVa+xPa6gN9q+ffvUleYICjyCoUOHqvWndOnS8v3332vNJLjM0CLEERR7hMvtSdYiFI785JNPxGvxBOuHszFEunDjceT9DKwzrPWscNG3rRAjHMLLEPELEwnDEq5HCr4ZIDN+LyVB10LE794DkeCHYoSEqPVILUcYk22xvtYT3i/H7Z7w/hJCCPEci1GmTJnklVde0bR8CCOrKAIonAhLTlxAKYCePXtqIHfKlCkjbUfsEOjSpYu8/fbb2qNt7Nix6kqbPn16pP0vXbqkbjVUx0acUXSwiWwCYSdMsIQ/dnlFsxhhYWKEhUb8HyZGKJZQMUIeifHokRgPH4nx4IGE/3tHwi/9K1d/SSXJzt8Q4587En7vnhjBDyUc+4WEiBES+vj5OE744+OqsHrCGCKNmxBCSJIgxsIIQc316tXTTLQ8efJIpUqV5LPPPpPDhw/H22BQ0RrusrJly2pQNRZUsP7mm2/0b7jOQLFixeyeV7RoUW0RYgUxTzVr1lSXW0xqLCGgHBlu1oW4GauY0uXx48cCJxyZACIPHkqKy3fEuHv/sWiKEED/Ezf2zyWEEELiRRjlzp1bG7wiE+3q1auaQg9R9MILL2iAMx7/+uuvEoa78VgClxiOafY2w1K+fHkNxMbfeJ3s2bNrMLWVkydPqlizWorg8itXrpy2EzEz1oiXYxU3RvjjhyGhIveCIwSRQ0sQCiFCCCEJnZUG0CoDgddYkL6Pgo8//fSTurfQYBYxPRAzroLGsaiLZAW1ipDtZq7v27ev1iVC0DdijGbNmqXWrCVLltiJIgglxBVdv37ddqxnnnkmNtMlngYEj5bkjxBCETFFkfYhhBBCEkMYWUFj2bp16+oCQbR//34JDQ2VhAKWqeDgYE3bv3XrlgoklA1AIUeAvxFwjSVnzpx2z2VbOB8VR6FwnVmtSfycCSGExI5YNZGFMDl06JDGA5kB0XowPz9p1KiR+AJsIusFoJ6Rf4D4pQpiE1lCSOLBJrI+3UTWZYvR6tWr5a233tLMNEcgjOISY0SIqxhGONS9u4dBCCHER3A5KhkB2Eh/R/sPWIusC0URSVRMQYQUf7jUCCGEkMQWRshI69Onjy11nhC3Q4sRIYQQdwmj119/XTZu3CgJBbLK2rRpo5loQUFBUqJECfntt9/shFn79u01bR+VtlHA8dSpU06PhfCpBg0aqItv+fLlCTZm4l4MBF7rHxRIhBBC4obLMUYTJkxQV9qWLVtUtCArzQoavMaW27dvS9WqVbUw46pVqyRz5swqeswq2xA6TZs21df88ccfNajqq6++0p5qR48e1dR+K19//XVEewlCCCGEkAQQRvPnz5dffvlFW3bAcmQVHvg7LsIIDWpz5cqlRRlN8uXLZ/sbImnnzp3aF6148eK6btKkSVqfCOPq2LGjbV8UhBwzZoxam7JlyxbrMREvwKxnRAghhCS2K23gwIHabBVpcOfOnZOzZ8/aFrPJa2xZsWKFVrqGRQp919ALbdq0abbtDx8+1P+tfdRQ1RrtPLZu3WrX8LZVq1YyceLEGBd1xLGRom9dCCGEEJK0cFkYPXr0SFq0aJEgbTYgrGABKlSokKxZs0beffddtUChujUoUqSItiZBw1e43TAWWJkuXryoWXImKP6IHmlNmjSJ8WuPHDlS6xaZCyxXxAvQpq+0FhFCCIkfXFY37dq1k4ULF0pCgJR/NJAdMWKEWos6d+4snTp1ksmTJ+t2xBYtXbpUe6NlyJBBg6/RjgQB1qZQg9UJPdsQX+QKEFuwgpnLhQsXEmSOhBBCCPGhGCPUKho9erRadEqWLBkp+BrB0LEFsUDFihWzW1e0aFH54YcfbI/RGBbxQxAvsBghQLtSpUrqggMQRWfOnJF06dLZHee1117ThrdRZdTBHYeFEEIIIUkXl4XR4cOH1ZoDEARtJa4ZYMhIO3HihN06WIfQENYRuLvMgGwEWA8fPlwf9+/f3y4IGyB7buzYsT7TroQQQgghHiKM4LpKKMzYILjSmjdvLrt375apU6fqYrJ48WK1EiHWCCKtZ8+emsKPJrYAwdbOAq6xvzXDjfgQzEgjhBDiLmGUkFSoUEGWLVum8T7Dhg1TIYNYodatW9v2QZA1Km+j0CNcb+jbNmjQILeOmxBCCCG+gZ+BqolPoGvXrvLxxx9Lzpw5n3hABGaHhobaiRlvBOn6cNfVkCYS6GcfR0U8C7+AQDHCw2g5IoQkDggb4fnGYwk1QmSj/KixyCgEnSAWI7iuUFARMUCI00GgM1pyoJ4Q0uZRdRp1hBYsWKDrra4vQgghhBCfshgBuK7+85//qPiBELKSNm1abcuBoGf0LvMFaDHyHmgxIoQkKrQY+bTFKMbCyAqsROfPn5cHDx5IpkyZpECBAj7Xk4zCyHugMCKEJCoURj4tjGJVvhpNXUuVKiXPP/+8FCxYMF5F0Z07d6RXr16aoh8UFKRZanv27NFtISEh0q9fP02/R8NYuO0QfH358mW7Y+TNm1fHZF0+//zzeBsj8TR4giKEEOKDWWkA7jjUR5o9e7YKnzlz5qibDu67NGnSyL59+zQLDcIMliuk6zdu3FhrGVlBVhuqZlvdfcQ3iYXRkxBCCPF8YQTXHKpc//jjj1K9enVdN3ToUPnpp5+0h9qnn34qa9eutXvOhAkTpGLFiuraQ60iqxCKaQNZs4ms2aQWsImsl0FxRAghJB6I/06wcQBp/mg5gmw3K3CpIevNGfAhwlXm2AIErrOMGTNqle4vvvhCjx0dbCJLCCEkRvhYTC2JQ/A1dkVz1SxZskQSL/EFYoqSJ08u8+bNk6xZs8r8+fO1cS1imRzbhQQHB2sJgSJFisjcuXPt+rWhGS0azW7fvl0LRr799tvR9nFzZjGCOGLwtReABsLh4e4eBSEkqcBzjkeTqFlp4eHhKoh+//13KVSokCQEaAD7zjvvyObNmyUgIEAFzrPPPit79+6VY8eO2fZDIDYaw168eFEbw0Y3+enTp0uXLl3k7t27MW4Uy6w0L4InKUJIYsJzjkeTqFlp/v7+Kohu3rwpCQVS/zdt2qQiBtYp9EuDCMqfP79tHzxGL7U///xTY46eNPFKlSqpK+3cuXMJNm5CCCGEJMEYI8Tu9O3bVzPHEhKk46MXGjLP1qxZI02aNLETRadOnZJ169ZpHNGTOHDggIo6uAAJIYQQQuItKw11g+7fv6/p8ogFQmC0lVu3bklcgAiCd69w4cJy+vRpFWGIIUKMEETR66+/rin7K1eu1EDtv/76S5+HeCKMZ8eOHbJr1y6pWbOmZqbhce/evaVNmzZaf4kQQgghJN6EEbrdJyTwCSJYGrFDEDuII/rss88kWbJk6gpbsWKF7le6dGm7523YsEFq1KihMURoW4I0fwRT58uXT4VRnz59EnTchBBCCPF+YtUSJCnA4GsvgoGQhJDEhOccj8YtLUGQOfbxxx9Ly5Yt5dq1a7pu1apVmq1GCCGEEOKtuCyMkDGGXmWI41m6dKlmj4GDBw/KkCFDEmKMhBBCCCGeKYz69+9va82BYGeTWrVqyc6dO+M0mKtXr0r79u21R1qqVKmkfv36mn3mDHgAGzRooFWvly9fblsPgQZLFoozIjC8aNGiMm7cuDiNixBCCCFJA5eDrw8fPqxVqR1BKvyNGzdiPRAInaZNm2qQNXqlwS+IStVmA1mk7zsGgUMUOYJCkBgLms9CHKHydefOnbVYZLdu3WI9PkIIIYT4Pi4LI/Qku3LlimZ7Wdm/f7/kyJEj1gOBZQgWJ9RHKl68uK5D41g0gkVbkI4dO9rVJRozZoz89ttvWuvICqpmW0FhSKTsw+1HYUQIIYSQeHWlvfnmm9KvXz+tHwSLDdqEbNu2TT744AOtcRRbzD5l1h5sKMqI9HtrA1nUUGrVqpVMnDhRRVNMQGQ6Uv+f9PrIRLMuhBBCCElauCyMRowYoQUX4aZC4HWxYsWkevXq2vwVmWqxBcfMnTu31jBCtetHjx7JqFGjtJ4RLFQmqEmE1zIrYT8JuNIWLlyo7rToGDlypKbnmwvmRwghhJCkhcvCCAHX06ZN05R9VJ9GLM/x48dl9uzZGscTU+bOnStp0qSxLXCjwd118uRJte4g+BpFGxFgDcsRQHHHX3/9NcZFJuGWg4BCtlzdunWj3ReCDJYlc0GfNkIIIYQkLVyOMTKBdQdLbGncuLE2dzVBfBKyyBA/BGECi1HmzJl1n/Lly+s+EEUQZIhzsoLq2C+88IJs3LjRtg4B27Vr11ZLUUwsWXDZYSFeCGuUEkIISczK166000AmWXyBgGy42FA8EhYfxDU5Zr6hphLS8Rs1amQLCEehSZQPaNeunYwePTpWr83K114EshMpjgghiXW+wcLK1z5b+TpGFiNknMUEZ+nzrrB48WK1EsEShbIAPXv21BR+0w2GYGtnAdfY3xRFcJ9BFNWrV08FndlkFm4+HJsQQgiJC35+/mIIhZGvEiNhhFifxABB1hAzKPSINHxkuQ0aNMilYyxZskSuX7+usU9YTPLkyaNNaAkhhBBCEqSJrBmg7IsZXHSleRF0pRFCEgs/P/HzDxAjLNTdIyGe0kQ2NDRUrTgQDXnz5tUFfyPAOSQkxOUBEEIIIYR4bVZa9+7dNa0eQc2VK1fWdagsPXToULl586ZWqyaEEEIISRLCCH3SFixYoPWFTEqWLKnuNDRvpTAihBBCiLfisisNtX7gPnMEWWEo/hgXkNXmbPniiy9s+3z22Wda+RoFIB3rGZmcP39eGjZsqPugoWzfvn3VBUgIIYQQEq/CCI1Yhw8fbuttBvA3BEtcm7QiK826TJ8+XYURCjiaoPDjG2+8Ie+++67TY4SFhakown5oBzJr1iyZOXOmDB48OE5jI4QQQhT/uJWmIT6WldasWTNZv369Wo5KlSql6w4ePKhCBJWmrSAWKS6ghtGdO3f09RyB2OnVq5f8/fffdutRDPKVV16Ry5cvS9asWXXd5MmTtfEt0vhjatViVpoXwaw0QkhiAU9GYKAYTDZK2gUercB9ZbXgJFS6PmoZ/fzzz2rxcQUEgqMatimKAIo9wsKEithlypRx+jxYvaxWMAgjQgghhCQtXBZGM2bMkMQAgiht2rTy6quvuvQ8VLq2iiJgPjarYDtj5MiR8sknn8RytIQQQghJkjFGJteuXZMtW7bogr9dZe7cuZImTRrbguNYQXxR69atJWXKlJIYDBgwQM1u5mIWrySEEEKc9ksjPonLFiO4mN577z1N2Uegs9mHrEWLFjJx4kSNy4kJjRs3lkqVKtke58iRw/Y3RNKJEydk4cKFrg5Pe6nt3r07klvO3BYViJnCQgghhJCki8sWo06dOsmuXbtk5cqVGviMBX//9ttv0qVLlxgfB26yggUL2pagoCDbtu+++07KlStnC+52BRSdRANaqxVr7dq1GoBVrFgxl49HCCGE2IClyB+XTlqMfBWXLUYQQWvWrJFq1arZBTdPmzZN6tevH+cBwSK1ePFiGTNmTJQ1im7duqX/w2J14MABXQ9xBZdc3bp1VQC1bdtWq3MjrgjtSmDlokWIEEIIIfEqjDJmzOjUXYZ16dOnl7gCFx0qCKCKtjNQj8iaqWZmmW3YsEFq1Kihbj2IN2ShwXqUOnVqadeunQwbNizOYyOEEJKEiYgrQn09A3+ySohP4nIdo6lTp6pFZ/bs2baYHVhlID6QQeaKO82TYR0jL4J1jAghiRh07R8UJOEPHoiEh7t7RMQT6hihF9rp06cld+7cugC4teCmQgHFKVOm2Pbdt2+fywMihBBCPBU/P///ZaTxpswnCYxNNWpCCCEkSbcEoSjyWVwWRkOGDJHEoGvXrmp9Gjt2rLb+ABs3bpSaNWs63R8p+hUqVNC/Fy1aJCNGjJCTJ09K5syZtYcbGskSQgghscK0Evn7iV9AgPiJH0OMfBSXhZHJ3r175dixY/p38eLFo2y1ERuWLVsmO3fulOzZs9utr1KlijaXtTJo0CDtpVa+fHlbrzQUhhw/frxmqGGMKDGAcgBxbXJLCCEkCWJ1ncGVFhj42GqEECNajnwOl4UR6gO9+eabar1B3zSAWkaw5CCjDBaauHDp0iXp3r27lgRo2LCh3TY0gLUWaQwJCZEff/xR90eWAEBQONx9sDiB/Pnza1XrUaNGacq+uR8hhBDyROyuGX7iB0GUPNn/RJIRTnGU1As8QoSg4z0asqKeEJYjR45oFlePHj3iNJjw8HCtPwS3F6xQT2LFihVy8+ZNefvtt23r0AjWsY0IrEUXL16UP//8M8pj4XmYg3UhhBCShLGKIj9/8fOPsBalChK/wMDHm7XYI9uEJGlhtHr1avn222+laNGitnUoqIh2IHBjxQVYdQIDA2MssFAhG8Ulc+bMaVuHx0uXLlX3GoQW4ozMYpGObjjHJrJIzzeXXLlyxWkuhBBCvAxT3MASBMETUeVaY4oCA8QveTLxS5NagnOlEQlKIX7Jkomff4Autuf6W55LoZQ0hBHERrJkkev6YB22xbaJ7KZNm2TcuHEyc+bMGLm7YAGCu61Dhw526xFPhFiiV155RV1vzz//vLr+gL+p7J3AJrKEEOJrAscicqJbIHwCAh8LIBVBgY+X5MnFP3kK8UuZUvxTpxa/dE/Lo3wZpcDLf4t/jnTi9/RT4gfrUcoUj/dLnvzx86zH0r8RkxSDcdiNm6LKa4RRrVq1pGfPnnL58mW7uKDevXtL7dq1Y3wcNJFFOw9z2b59u8YvoTYSrEZY4Pp6//33JW/evJGeP2PGDK3CjeNYgaiC5enu3bv6fBSfrFixoi3eKCpQhwmFoKwLIYSQKHC8iCf64h/1on3MIhbVF34O6/y0HhFcY48tPhAmEccMCBAJiLAS4W9YipIFalyREZRcgtMFyKt5jknAUwEiKZI9thrBWJAsYn/zWDbBg78Rm+T/+DVV8DhbnI3T2fzc/b77ia/jcuVrWFIgRhBjZLqbsO65557TmB+rW8sVECvk6OqCWwwxR4ghKly4sG09hlygQAGttP3ll18+8dhvvfWWFqWE+IoprHztRTDwkRASW6K80Jui6vE+KmxgDUqdSkKezS6Bx86LcT9YJCxcDARgm+egqM5FPEf5buVriCFUtF63bp0cP35c1yHeqE6dOhIXYP3B4uieQxaaVRSBX3/9Vc6ePSsdO3aMdJwbN27IkiVLtG9acHCwWpbQwgSuOkIIISRmgsV43AstQjgZYWGPHwc/lGRX74gR/EgkNCxmooj4fh0jmAJfeuklXdwBgq5R06hIkSJOt6PJ7AcffKCWJTSSRWkB051GCCGExBiL2IEI8gsNFbkfrELJJoooiJKuKw3B1QiORtbXuXPnVCDly5dPXn/9dXV5+VKNILrSvAi60gghCU1EfA0Cqf2fSivh//wrRljo4208//iUKy3GwdfQT4gtgvsKwdYlSpTQWkMIcG7fvr00a9bM5RcnhBBCvAa1DoWLqFuN7jNJ6q40WIo2b96s9YEc+5Uh5gfVpr///nsNdCaEEEJ8CgggtU4/jjdyMW+JeBExthjNnz9fPvroI6dNXJHC379/f61NRAghhPgqBpSR1uxjbJEkdWF06NAhqV+/fpTbGzRoIAcPHoz1QND3rF+/fuqiS506tTaQhfXJWi/JsYVH6dKlNa4JdZCsoPAjCjumTZtWe7e99tprGhNFCCGExAkGW/s8MRZG6ImWNWvWKLdj2+3bt2M9kPv372sZgEGDBun/CPA+ceJEpAKOJh9++KGKJ0eQxt+kSRO1YkEwQSQhhR81jwghhJA4Q2Hk08Q4xigsLEyrUUdFQECAhCKNMZYgA2zt2rV26yZMmKBp9ufPn9eK2CboyfbLL7/IDz/8EKk/2969e3Wsn376qa0FCFL3IZZglXLWzoQQQgiJeZxRRI0jkrSFEQLNkH2G1hlRubbiG6TawVWWLl0627qrV69qP7Tly5dLqlSpIj2nXLlyKohQ2BHjRWuQ2bNnawHK6EQRxm+dA9L1CSGEEDsMQ4xwqCIqI0nqrrR27dpJlixZ7DrQWxdsi8+MNFStRsxRy5YtbXUITHHWtWtXKV++vNPnoa4SrEkIFIeIg6hCw9lFixZF+3ojR460m4/Z7oR4AT5UP4sQ4g1QFPkyLvdKiy+QwdalSxfbY7jEXnjhBf0bLi8ETEPQoGq1KYy++eYbFTho7wHXHQKqIYT279+vgdgATWOrV6+u5QMgqu7cuSODBw9WNyBcdVEVoXRmMYI4YoFHLwAuU80SIYSQBAZFHpMlEyMkhLFGHkqi90qLLxBUXalSJdvjHDly2ERR8+bNtXAk6iNZJ4XHO3bsiOTOg/WodevW2gpk4sSJavEZPXq0bfucOXNU5OzatUuz1ZyBY0blJiSEEEJsMDPNp3GbMEIqPRYrpig6deqUbNiwIVJTWViMEFRtglT+evXqycKFC20iC9ltZtC1CaxLZksTQgghhBCPE0aOQBSh5xpS9VeuXKmZZXCLgQwZMkjy5MntMtNAmjRp9P8CBQpIzpw59e+GDRvK2LFjZdiwYTZXGuKN8uTJI2XKlHHDzAghhBDic8HXCQ36r61YsULjihAvlC1bNtuyffv2GB8H9YvmzZunWWsQQihKCRfZ6tWrJSgoKEHnQAghhBDvxm3B154Ogq8Rq8Tgay+AwdeEkMQMvvYPECMs9nX7iGcHX3uMxYgQQgghxN1QGBFCCCGEeKswQg+1unXrasaaswayJkjrR7wRGtLClIbaRg8ePEj08RJCCCHEe/A6YXTv3j2pVq2ajBo1Ksp9IIoQdA0BtXv3btmzZ49069YtUho/IYQQQohHpuvHlLZt2+r/qHodFb1795YePXpI//79besKFy6cKOMjhBBCiPficyaUa9euaYVr9G6rUqWKZM2aVV588UXZunVrtM9DOxBkolkXQgghhCQtfE4Y/fHHH/r/0KFDpVOnTlq/qGzZslK7dm2tqB0VbCJLCCHkiRiGGAbLg/gyHu1Ki67RbFSYbT/wvLffflv/RqHH9evXy/Tp01UAOWPAgAHSp08f22PUP0Cl7VBBo8B4mhBJGAx/EZ6oCCGJBc85Ho1et1XDGr4njKJqNBsdqJQNihUrZre+aNGicv78+Sif59hE9saNG/r/VvlvrMZOEhGenwghiQnPOV4BWoLBA+RTwshZo9knkTdvXsmePbucOHHCbv3JkyelQYMGMT4O+rMBiKnYvLGeCmKn4Ca8cOFCrCqCejK+OjfOy7vw1Xn58tw4L9+al2EYKoqgBWKDRwsjZ9y6dUvFyuXLl/WxKYCeeeYZXVDbqG/fvjJkyBApVaqU9l2bNWuWHD9+XJYsWRLj1zFT+yGKfOkLZYI5+eK8fHlunJd34avz8uW5cV7eRXTziotBw+uEERrNmrFD4M0339T/IYQQcA169eolwcHBmrYPIQWBtHbtWilQoIDbxk0IIYQQz8frhFH79u11eRKoYWStY0QIIYQQkuTS9eMLBGLDCmUNyPYFfHVevjw3zsu78NV5+fLcOC/vIkUCz8vPiG0+GyGEEEKIj0GLESGEEEJIBBRGhBBCCCERUBgRQgghhERAYUQIIYQQEgGFURRMnDhRq2inTJlS25Ls3r1bvIWwsDAZNGiQ5MuXT4KCgrR+0/Dhw219Y0JCQqRfv35SokQJSZ06tVYHfeutt2xFMz2JzZs3S6NGjXSMKN65fPnySPscO3ZM28egoBfmU6FCBaftXzB/VD+P6jiJCXr2YZyo7J4lSxZp2rRppGrtMR33zJkzpWTJkvpdxbHee+89cReTJk3SsZiF1ypXrqw9DgFqinXv3l0KFy6s30v0IuzRo4f2JbSyZ88ebfqcLl06SZ8+vdSrV08OHjwonsTnn3+unwdqpplMnTpVatSoofPGtr///jvS8/AetG7dWvfB/Dp06CB3794VT57XmTNnpFmzZpI5c2Ydd/PmzeXq1at2z8PvD58nvoNoy9S2bVu3n09Q1w5zsS5FihTRbefOnYu0zVwWL15sO4az7QsWLBB3c+nSJWnTpo1kzJhRf0s4l//222+27Shp4zju+vXr27Zv3Lgxyvnv2bPHLXPC9dbZeMzzGfqf4lqG+eK72KRJEy3cbAU9UatUqaLnVRR8xnUuNDTU5bFQGDlh4cKF2lAW6YD79u3TApE4OV+7dk28gVGjRukFasKECSoa8Hj06NEyfvx43X7//n2dF8QT/l+6dKlelHFy8zTu3bun7z+EqjNw0q5WrZqe8PBjP3TokM4LJ2hHvv76a/2heQKbNm3SH/zOnTu1+CjEat26dXW+roz7q6++koEDB2rNrt9//13WrVun31V3kTNnTr247t27V0/UtWrV0hMYxoYLJZYvv/xSjhw5ooJu9erVKg5MIBJwAsdFdteuXbJ161Y9yWFOeI88AVw4pkyZogLQCn5XGPtHH30U5XMhivBe4DNfuXKlCv/OnTuLp84L30d8L/H9+/XXX2Xbtm3y6NEjvVkxG3aDmjVryqJFi/Q88sMPP+jv8vXXX3fTTP5H8eLF5cqVK7YF3yeAdhLW9Vg++eQTSZMmTaTWUTNmzLDbDzcx7uT27dtStWpVSZYsmd50HD16VMaMGaM3EVbwXbSOe/78+bZtEA+O8+/YsaPeTJcvX95t3z/rePAbAW+88Yb+X65cOf0scE1bs2aN3jDiuwlDAMDN08svv6zz3r9/v17HURA6VvUMka5P7KlYsaLx3nvv2R6HhYUZ2bNnN0aOHGl4Aw0bNjTeeecdu3Wvvvqq0bp16yifs3v3bpiTjD///NPwVDC+ZcuW2a1r0aKF0aZNmyc+d//+/UaOHDmMK1euOD2Ou7l27ZqOa9OmTTEe961bt4ygoCBj3bp1hieTPn164z//+Y/TbYsWLTKSJ09uhISE6OM9e/boPM+fP2/b59ChQ7ru1KlThru5c+eOUahQIWPt2rXGiy++aPTs2TPSPhs2bNDx3r5922790aNHdT3maLJq1SrDz8/PuHTpkuGJ81qzZo3h7+9v/PPPP7Z9//77bx0z9o2KH3/8Ufd59OiR4S6GDBlilCpVKsb7ly5dOtJ50xPPFf369TOqVasW7T7t2rUzmjRpEuNj4nPKnDmzMWzYMMNTwHewQIECRnh4uNPtBw8e1M/n9OnT+njAgAFG+fLl7fZZsWKFkTJlSuPff/916bVpMXIAd0O4261Tp45d3zQ83rFjh3gDuBuASRGNc00ljTul6Jrowp2Bu0KY970F3LH+/PPP8uyzz6pFAW4kuD0d3U24k2/VqpVanWBe9URMd5LZvDgm48YdFd4DmNWLFi2q1hq4OdBY0RPAnRzcDrA6wKUW1bzhngkMfFyEH242uAe+++47/S0+ePBA/8b8YGp3N7DyNWzY0O78EFNw/sDvy3pHjuPg/ALrmCfO6+HDh3pesBbSgzUWYzatL87chXPnztXzEKwa7uTUqVPqhs+fP79a65y52AHO+QcOHLCzXlrfm0yZMknFihVl+vTptpAEdwErCL5DsKTgnFemTBmZNm1apP1gQcd2/KbeffdduXnzZrTHxHZruy13gt/+nDlz5J133nFqLcc5BdYjWLhg/TO/q46eArjd0B4Mn69LxEHQ+SS4c8Pbsn37drv1ffv2VUuSNwALF+4qcMcWGBio/48YMSLK/R88eGCULVvWaNWqleHJON69mVaUVKlSGV999ZVaV2DVw3w3btxo269z585Ghw4dojyOJ3xesPJVrVrVbv2Txo25JkuWzChcuLCxevVqY8eOHUbt2rX18cOHDw13AQtP6tSpjYCAAOPpp582fv75Z6f7Xb9+3cidO7fx0Ucf2a0/fPiw3inCUoEF8zl37pzhbubPn28899xz+nsBrlqMPvvsM+PZZ5+NtD/u1L/99lvDE+cFS+ZTTz2lj+/du2fcvXvX6Natm84P308rH374of4Wse355583bty4YbiT//73v2qRhGUBv4/KlSvr982Z9eDdd981ihYtGmk9LChbt2419u3bZ3z++edGihQpjHHjxhnuBGPAAgsJxjVlyhS1isycOdPuM4XVDr9FnDMwtwoVKhihoaFOj9mgQQNdPIWFCxfq+cPRkjpx4kQ9t+A7hvOCaS2yWjfnzZun87x48aLxwgsv6L5Y5woURj4ojPCjyJkzp/6PH8b3339vZMiQwe6HYzWhNmrUyChTpoydudwTcRQG5mfVsmVLu/0wnzfffFP/xsmhYMGC6iqI6jjupmvXrkaePHmMCxcu2NbFZNy40GIdTggmuJDh5IALgbuAKIPb67fffjP69+9vZMqUyfj999/t9sF3Db+n+vXr27lb7t+/r+vfeustde9C7L322mtG8eLFdZu7gGsvS5YsepE18QVhFJN54fuVP39+veHAxQqua9xI4XvrKHRPnDhh/PLLLyryX3755SjdIO4AnwdEnqNbF98rCPgvv/zyiccYNGiQnlvdCW6GIPKsdO/eXcVoVJw5c0a/k87c7jjv4JyxZMkSw1OoW7eu8corr0RaDzfuyZMnNeQA53l8D01BD8aMGaOfMb6nEOm4ecS8FyxY4NLrUxg5OanjTXW8cOJE3bhxY8MbwA93woQJduuGDx+uCtsKLkhNmzY1SpYs6fa7u5jgKAzwWcEihrk53rlWqVJF/8YJ3jyhmwuOgxMBLgDuBrFs+Lz++OMPu/UxGff06dN1nVVQAVzopk6dangKsGJZrQu4Y8eJHeutJzWAixbGDyua9XPGSQ5C313ge4f32vHzMD8j6514VMLou+++M9KlS2e3DrFVeP7SpUsNT58XhI85p6xZsxqjR4+O8rj4Tjq7wXQ3iEGBWLeCG0eIDdxUPImVK1fqvIKDgw13AauX1ZIMIKwRBxsduEGZPHmyU6sYxLk748GswDqM89zy5cuN6DDPC47WIIhx3DRD8JpxfbjJcoXHjn1iI3ny5Br9jhgdM/sAcRx43K1bN/EGEJuCGAArAQEBdlkkyPBBPAp88Bs2bNC4Dm/8rJDy7pjmjtiqPHny6N/ISEC2hRWkto4dO1Yza9wFdB5S15ctW6axAPCVW4nJuJGZAjB/xBeZ8R03btywzd8TwPcO/n/w77//ajwYYlYQ1+AYE2B+d61xBeZj6/c3sUH5gMOHD9utQzwGsiGREozf15NAnBVS+BHvgHMMQKYX5oXYOE+fF+JszDEjQze6LFbzszI/d08AGY/IlkMpASuIYcNckAL+JBCHhOwvdzZlxe8+unOeMy5evKgxRCil4HgeQqwOyrW4Ox7MBONBbBRi3qIjwrAT6TuGcwXiygAy8RCDVLZs2WiP5ezgxAGY3eDDhesJihN3u7jT++uvvwxvABkJyGTC3c3Zs2f1bhR3C7CkANwZwPoFS8WBAwc0Vsdc3Bmb4gy4khA7hAVfVzOWyMyew9xwtwcLCdw348eP1zvdLVu2RHlMT3ClIaYB5nvEQlnf/+jcRc7GjcwTuJm2bdumsTkwPxcrVsxtd3+4G4eZG987uHHxGNYHuFfgPqtUqZJRokQJjQ2wztu0TBw7dkx/e3h/8Ns7cuSIum7wXl2+fNnwJBxdTpgHvpvTpk3Tz2rz5s36+ObNm7Z94DqE23rXrl0au4JMMEdXsKfNC5ZJuDTxmc2ePVvd8n369LFt37lzp/7uMFfc7a9fv14ttogTc6dl5f3339ffF76L+H3UqVNHz4NWyxDOGfh+IjvQEWQ04bPE7wr7wSoDC8XgwYMNdwLrByzlcM1iXHPnztVxzZkzx3bO/OCDD/Qzw9zhPoPLCd81x88D2/Bdxe/OEwgLC1OLGGJkHV2BiJOFex7nfnyecKXhu3j16lXbfrBi4ryD8wYsYbg2xOZcT2EUBfih4wNCKjFiHvDj9xbgqsCJDeNHUB7iAwYOHGgTPfix4MfgbIEbwJMw3RKOC8Sf1UWBeBzMFem5TzLBeoIwiur9nzFjhkvjhthAijGEO04SzZo1s0t1T2wwFsRL4XcD8zzcZRBF0X2WWPCdNDFjVCCGkOpfq1YtPcl7Go4CAunhT/pMIZIghNKkSaOxEG+//bZdHJknzgsXKbjOcJHBxRVxHNbYIVyIatasqd8/iNq8efNq/BGCX90JSnlky5ZNv4u4UcRja7AuQABzrly57Fy3JhBLSOHHZ4WAX5xb4Ipytm9i89NPP2nAPN7vIkWK2LnOcXOFGB38/vCZ4ffYqVMnpzf2+C6aYQeewJo1a/Q3g1g1K3CNITgcbnbMCTf1SBY6fvy43X74HuK8gWsBbsIQgB8b/PCPazYmQgghhBDfhHWMCCGEEEIioDAihBBCCImAwogQQgghJAIKI0IIIYSQCCiMCCGEEEIioDAihBBCCImAwogQQgghJAIKI0IIIYSQCCiMCCHEBQYNGiSdO3eO0zGOHj2q/e3u3bsXb+MihMQPFEaEkEShffv22uARCxpWonHuhx9+KMHBweIt/PXXXzJu3DgZOHBgnI5TrFgxef755+Wrr76Kt7ERQuIHCiNCSKJRv359uXLlivzxxx8yduxYmTJligwZMkS8hf/85z9SpUqVaDuZxxR0sZ80aZKEhobGy9gIIfEDhREhJNFIkSKFPPPMM5IrVy5p2rSp1KlTR9auXWvbfvPmTWnZsqXkyJFDUqVKJSVKlJD58+fbHaNGjRrSo0cPtTZlyJBBjzd06FC7fY4fPy7VqlWTlClTqnVm3bp1aqlavny5bZ8LFy5I8+bNJV26dHqcJk2ayLlz56Id/4IFC6RRo0aRxtO9e3fp1auXpE+fXrJmzSrTpk1TNxnET9q0aaVgwYKyatUqu+e99NJLcuvWLdm0aVOs3ktCSMJAYUQIcQtHjhyR7du3S/LkyW3r4FYrV66c/Pzzz7odsTxt27aV3bt32z131qxZkjp1atm1a5eMHj1ahg0bZhNYYWFhKrogrLB96tSpkVxfISEhUq9ePRUtW7ZskW3btkmaNGnUovXo0SOn44WIQWxQ+fLlI23DeDJlyqTjhEh699135Y033lDr0r59+6Ru3bo6j/v379ueg3mXLl1aX58Q4kEYhBCSCLRr184ICAgwUqdObaRIkcLA6cff399YsmRJtM9r2LCh8f7779sev/jii0a1atXs9qlQoYLRr18//XvVqlVGYGCgceXKFdv2tWvX6ustW7ZMH8+ePdsoXLiwER4ebtvn4cOHRlBQkLFmzRqn49i/f78e4/z583brHccTGhqqc2zbtq1tHcaC5+7YscPuuc2aNTPat28f7fwJIYlLoLuFGSEk6VCzZk2Nq4GbCTFGgYGB8tprr9m2w9ozYsQIWbRokVy6dEmtNw8fPlTrj5WSJUvaPc6WLZtcu3ZN/z5x4oS66uBiM6lYsaLd/gcPHpTTp0+rxcgKLFZnzpxxOvYHDx7o/3DPOWIdT0BAgGTMmFHdgCZwrwFzjCZBQUF2ViRCiPuhMCKEJBpwfyHeBkyfPl1KlSol3333nXTo0EHXffHFF5r19fXXX6uwwP6I3XF0byGrzQrih8LDw2M8jrt376rLbu7cuZG2Zc6c2elz4CoDt2/fjrSPs/FY1+ExcBwj3HMFChSI8bgJIQkPY4wIIW7B399fPvroI/n4449t1hjE+iAIuk2bNiqa8ufPLydPnnTpuIULF9bA6qtXr9rW7dmzx26fsmXLyqlTpyRLliwq1KzL008/7fS4EDBPPfWUxhnFF4ijKlOmTLwdjxASdyiMCCFuAwHKcD1NnDhRHxcqVEiDqBGUfezYMenSpYudwIkJyPaCiGnXrp0cOnRIxRbEl9Vy07p1a7UAQYQh+Pns2bOyceNGzXa7ePFilEIOWXRbt26V+AAZcHAX4piEEM+BwogQ4jYQY9StWzfNLEPcEQQMrDnIGEMaPOKEkGHmChBaSMuHu6xChQrSsWNHW1aaGR+EmKXNmzdL7ty55dVXX5WiRYuqOw8xRrAKRQWOhZR9V9x2UYEyBMhWi4+aSISQ+MMPEdjuHgQhhCQksBqhrhECruMS04PTZaVKlaR3795abym2IGYK1rF58+ZJ1apVY30cQkj8w+BrQojPsWzZMq1LBPEBMdSzZ08VIHENdIYrDnWRDh8+HKfjnD9/XuOrKIoI8TxoMSKE+Bzff/+9fPrppypAEEuEOJ4xY8ZoGj0hhEQHhREhhBBCSAQMviaEEEIIiYDCiBBCCCEkAgojQgghhJAIKIwIIYQQQiKgMCKEEEIIiYDCiBBCCCEkAgojQgghhJAIKIwIIYQQQuQx/w/+3QjdkaeK5gAAAABJRU5ErkJggg==", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "from numpy import abs ,angle, arange, arcsin, cos, concatenate, pi, sqrt, tan, zeros, linspace\n", "from scipy.fft import fft, fft2\n", "import matplotlib.pyplot as plt\n", "\n", "\n", "def y_IF(f0_min, slope, T, antenna_tx, antenna_rx, target, v=3e8):\n", " \"\"\" This function implements the mathematical IF defined in latex as\n", " y_{IF} = cos(2 \\pi [f_0\\delta + s * \\delta * t - s/2* \\delta^2])\n", " into following python code\n", " y_IF = cos (2*pi*(f_0 * delta + slope * delta * T - slope/2 * delta**2))\n", " Parameters:\n", " -----------\n", " f0_min: float\n", " the frequency at the begining of the chirp\n", " slope: float\n", " the slope with which the chirp frequency inceases over time\n", " T: ndarray\n", " the 1D vector containing time values\n", " antenna_tx: tuple of floats\n", " x, y, z coordinates\n", " antenna_rx: tuple of floats\n", " x, y, z coordinates\n", " target: tuple of floats\n", " x, y, z coordinates\n", " v: float\n", " speed of light in considered medium\n", " Returns:\n", " --------\n", " YIF: ndarray\n", " vector containing the IF values\n", " \"\"\"\n", " tx_x, tx_y, tx_z = antenna_tx\n", " rx_x, rx_y, rx_z = antenna_rx\n", " t_x, t_y, t_z = target\n", " # distance tx antenna to target\n", " distance = sqrt((tx_x-t_x)**2 + (tx_y-t_y)**2 + (tx_z-t_z)**2)\n", " # distance target to rx antenna\n", " distance += sqrt((rx_x-t_x)**2 + (rx_y-t_y)**2 + (rx_z-t_z)**2)\n", " # usually delta_t = 2*d/c, but\n", " # distance is already 2*D (TX to target + distance target to RX)\n", " # so delta = distance/v\n", " delta = distance/v\n", " YIF = cos(2 *pi *(f0_min * delta + slope * delta * T - slope/2 * delta**2))\n", " return YIF\n", "\n", "f0_min = 60e9\n", "c = 3e8\n", "# lambda ~5mm at 60GHz\n", "lambda0_max = 3e8/f0_min\n", "n_rx = 2\n", "Distance = 10\n", "k = 10e12\n", "n_samples = 512\n", "f_if = 2*k*Distance/c\n", "fs = 50e6\n", "ts = 1/fs\n", "\n", "antenna_tx = (-lambda0_max/2,0,0)\n", "antenna_rx = (0,0,0)\n", "T = arange(0, n_samples*ts, ts)\n", "t_chirp_to_chirp = 1.2e-6\n", "n_chirps = 128\n", "print(len(T))\n", "\n", "cube2D = zeros((n_chirps, n_samples))\n", "for d in [160]:\n", " for v in [460]:\n", " for chirp_i in range(n_chirps):\n", " d_i = d + v*t_chirp_to_chirp*chirp_i\n", " f_if = 2*k*d_i/c\n", " target_t_i = (0, d_i, 0)\n", "\n", " # sanity check\n", " assert f_if < 1/ts/2\n", " assert v < lambda0_max/4/t_chirp_to_chirp\n", "\n", " YIFi = y_IF(f0_min, k, T, antenna_tx, antenna_rx, target_t_i)\n", " cube2D[chirp_i, :] = YIFi\n", "\n", "Z_fft2 = abs(fft2(cube2D))\n", "\n", "Data_fft2 = Z_fft2 # [0:n_chirps//2,0:n_samples//2]\n", "\n", "# change scale for 2D plot\n", "# range formula\n", "# ranges from 0 to max distance with n_samples values\n", "# max un-ambigous range is fs*c/2/k\n", "ranges = linspace(0, fs*c/2/k, n_samples)\n", "no_labels = 10 # how many labels to see on axis x\n", "step_x = int(n_samples / (no_labels - 1)) # step between consecutive labels\n", "x_positions = arange(0, n_samples, step_x) # pixel count at label position\n", "x_labels = ranges[::step_x] # labels you want to see\n", "# rounding up the values for easier display\n", "x_labels = [int(d) for d in x_labels]\n", "# allocate labels to ticks\n", "plt.xticks(x_positions, x_labels)\n", "\n", "# Max un-ambigous speed is lambda0_max/4/t_chirp_to_chirp\n", "# if we `assume` no negative speed then\n", "# Max non-un-ambiguous is 2x lambda0_max/4/t_chirp_to_chirp\n", "# i.e. from 0 to 2x lambda0_max/4/t_chirp_to_chirp\n", "# instead of [- lambda0_max/4/t_chirp_to_chirp ; lambda0_max/4/t_chirp_to_chirp]\n", "no_labels_y = 10 # how many labels to see on axis x\n", "step_y = int(n_chirps / (no_labels_y - 1)) # step between consecutive labels\n", "no_negative_speeds = False\n", "if no_negative_speeds:\n", " speeds = linspace(0, lambda0_max/4/t_chirp_to_chirp*2, n_chirps)\n", "else:\n", " # labels as we haven't done freqshift\n", " speeds_pos = linspace(0, lambda0_max/4/t_chirp_to_chirp, n_chirps//2)\n", " speeds_neg = linspace(-lambda0_max/4/t_chirp_to_chirp, 0, n_chirps//2)\n", " speeds = concatenate((speeds_pos, speeds_neg))\n", "y_positions = arange(0, n_chirps, step_y) # pixel count at label position\n", "y_labels = speeds[::step_y] # labels you want to see\n", "y_labels = [int(dop) for dop in y_labels]\n", "plt.yticks(y_positions, y_labels)\n", "\n", "# change scale for 2D plot\n", "\n", "plt.imshow(Data_fft2)\n", "plt.xlabel(\"Range (m)\")\n", "plt.ylabel(\"Doppler (m/s)\")\n", "plt.title('Velocity-Range 2D FFT')" ] }, { "cell_type": "markdown", "metadata": { "id": "S66vx0NVPbsz" }, "source": [ "## 2D FFT - 3 targets\n", "\n", "> 3 targets 2 same range, 2 same speed" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "colab": { "base_uri": "https://localhost:8080/", "height": 367 }, "id": "FHn7s0BiPglc", "outputId": "214e40f1-259f-4982-dfa8-37bc0b2ed1b2" }, "outputs": [ { "data": { "text/plain": [ "Text(0.5, 1.0, 'Velocity-Range 2D FFT')" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkQAAAFNCAYAAADhB3APAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjMsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvZiW1igAAAAlwSFlzAAAPYQAAD2EBqD+naQAAPItJREFUeJzt3QmcTfX/+PH3nX3s+5adCmUnaxFKkohvxVcliTZKlNJCSbT8opQoeyRFkVL2LCFrSpQtImt2BrOe/+P96X/v994xxsw9M3PvzHk9H49r3LN+zj33nvu+n/fn8zkuy7IsAQAAcLCQQBcAAAAg0AiIAACA4xEQAQAAxyMgAgAAjkdABAAAHI+ACAAAOB4BEQAAcDwCIgAA4HgERAAAwPEIiACHcLlc8sorr2TKtvfu3Wu2P3ny5EzZPgBkNgIiIMjceeedkitXLjl79uxll+natatERETI8ePHJVh99913mRKA6TY1+HI/wsPDpXz58vLkk0/KqVOnJKf46quv5N5775WKFSua98O1114r/fv3T/EYvV+PsLAwKVSokNStW1eeeuop2bZtW5r3qa+j97a8HxcvXjTLaNB7uWWef/55ad68+WXnez8yKzgH/BXm95oAMoUGO998843Mnj1bHnjggUvmnz9/Xr7++mu57bbbpHDhwhIMypUrJxcuXDDBiXdANHr06Ez74hszZozkyZNHYmJiZMmSJfL+++/Lpk2b5Mcff5ScoFevXlKqVCm57777pGzZsrJlyxb54IMPzOuqxxkdHe2z/C233GLeL3p7ytOnT8svv/wiU6ZMkQ8//FDefPNN6devX5r2W6tWLRN4JacBuLchQ4ZIhQoVfKZdf/310rJlS3n44Yc909avXy+jRo2SF154QapWreqZXqNGjTS/FkCW0Ju7Agge58+ft/LmzWu1bt06xfnTp0/XGzJbM2bMSNd2dZ3BgwdbWeWJJ54w+8xoegy63X/++cdn+r333mumr1271soJfvjhh0umTZkyxRzjuHHjfKbrNH29kzt27JjVqFEjM3/evHlX3Ge5cuWstm3bprrMpEmTzPbWr1+fpuOYOXOmWT6l4wGCCSkzIMjoL/+OHTuaWo+jR49eMn/69OmSN29ek1pTmkLp27evlClTRiIjI6Vy5cqmRiApKemK+/r555+lTZs2ki9fPlPbor/uf/rpp0uW0308/fTTJqWi+yhdurSpjTh27FiKbYgefPBBUzukvNMk+t2t22jfvv0l+9CUTP78+eWRRx7x41UTufHGG83f3bt3e6adOHFCnnnmGalevbo5Pj1OPV6tPfG2bNkyU74vvvhCXn/9dXN8UVFR5vXYtWvXJfvSY9NUlp6rG264QVauXGlSRfrwFhsbK4MHDzbnRF83PUcDBgww068k+bbUXXfdZf7+/vvvaXpNtAZxxowZJo2mxwXg8kiZAUGaNtN0h35B9+7d2+cLfsGCBdKlSxfzZazps2bNmsmBAwdMIKGpldWrV8vAgQPl0KFD8u677152H1u3bjVBhAYJ+iWt6a6PPvrIfBEvX75cGjRoYJY7d+6cWU6/hB966CGpU6eOCYTmzp0rf//9txQpUuSSbWtZDh48KIsWLZKpU6d6pmvQoSmgt956yxyLtnVx0zThmTNnzHx/aFCmChYs6Jn2559/ypw5c+Tuu+826Z0jR46YY9TXTNvWaErK2xtvvCEhISEmiNK0k5ZTz8XatWt9UnV6TvQ10SBR99uhQwezXw2k3DQg1aBVU3ia/tJ0kaa9Ro4cKTt27DDlSq/Dhw+bvym95pej7wk93h9++MG8vnq+UxMfH+8JdN20DZM+vOnrk3y59JQLCDqBrqICcKmEhASrZMmSJt3hbezYsSb9sGDBAvP8tddes3Lnzm3t2LHDZ7nnn3/eCg0Ntfbt23fZlFmHDh2siIgIa/fu3Z5pBw8eNOm6m266yTNt0KBBZt2vvvrqknImJSWZv3v27DHLaDrlSimz7du3m+ljxozxmX7nnXda5cuX92zzSikz3Y6mzfbu3WtNnDjRio6OtooWLWrFxMR4lr148aKVmJjos76WNTIy0hoyZIhnmqZzdJtVq1a1YmNjPdPfe+89M33Lli3muc4rXLiwVb9+fSs+Pt6z3OTJk81yzZo180ybOnWqFRISYq1cuTLFc7hq1SorvXr06GHOa/LzfbmUmdtTTz1llvnll1+umDLT5ZI/vN837pRZSo+UkDJDdkHKDAhCoaGh0rlzZ1mzZo2n5sOdLitevLhJ5aiZM2eamgqtndBf6+5Hq1atJDExUVasWJHi9nXewoULTc2Gpn7cSpYsKf/9739NrYbWJqgvv/xSatas6UnXeNMan/S65pprTO3Tp59+6pmmtUXff/+9qY1J6za111XRokVNCk5rrjQtpdvwrsnQNJXW+LiPWXvlaepM19WGycl1797dp/GwOw2nNU1qw4YNZhs9e/Y0aSg3Lbd3zZT73GitUJUqVXzOTYsWLcx8rbFJDz33EyZMMA2er7766nStq8esUuu56KbnRmv2vB8pNe7XtGHy5YDsjJQZEKT0S1bTK/pFqD10ND2lbVW0e7kGTGrnzp3y66+/msAgJSm1QVL//POPSbdpYJCcfolrumf//v1y3XXXmTY5nTp1ytBj0y9YTTv99ddfpoeaBg+aqrn//vvN/Li4OBMkedNjdB+3O1DT9I8ei/Zi2rNnzyU9r/Q43nvvPdPTSudrUOSWUg89TS95cwc5J0+eNH+1vEqDL28aHGlg5k3PjaYZ03tuUqLnvUePHtK6dWu/2gJp2lNp27Mr0bSXBtRXom2n6tWrl+6yAMGKgAgIUjqOjNYufPbZZyYg0r+aHdFAyfsLX7tbaxugy9XGBCOt/dL2N1pLpMc2bdo08+XqDtC0HdTNN9/ss44GNN5Bx0033eRps9KuXTvTcFpfm40bN3pqhYYNGyYvv/yyqUF67bXXTJslnaeN0FNqdO4dcHn7NyuVPrp9LdOIESNSnK8NrNNCG4BrWyTt0j5r1iyfmqm0+u2338yxJe8mD+B/CIiAIKZf8PqFrrVAWlOkqZL69et75leqVMn8+k/LL3pvWmuhqaXt27dfMu+PP/4wQYP7C1v3oV+o6ZVa6ksDk7Zt25qASI9x1apVPg3ANUWXPAVTokSJVFNC2ptLU17aEF0DLqUBhAZWmmpK3mvOnwbAWpultOeZd8CWkJBgUpveY+vo66bBjKY3/UktKq2d0/GmihUrZsYfcqe+0mPfvn2mkXyjRo3SVEMEOBVtiIAg5q4NGjRokGzevNmndkjdc889pp2R9jxLTr/09Ys6JVpbcOutt5oBHr3bKGkvLA28mjZt6umNpOky/WLXgSLTU3OSO3duTzlSoukx7en17LPPetpMeaeqNMjzfmg3+NToa6O9vHTIAe/jTF5GTc9przx/aC2WptrGjRvn89pqYOdOq3mfG92PLpucDmKpA0peqUeZniMNTvX8Xi71lhpNO2qPRE0Vvvjii+leH3ASaoiAIKYpjsaNG5vARSUPiDSY0O7vd9xxhxn7R9Ns+kWr3bu1dkSDncvVhAwdOtTUwmjw8/jjj5tUjHZJ1zFytLu59z50W9p1XVNPug/9otX9jh071tTmpESXU9rmSdu+JA96tIZIgwsNUHRsIK0FsUOHDdBbVWh558+fb2pW9HXREZW15khfR31dNHjxbkieHtrgWkfe7tOnj2kcrUGPvsY6/pLWCHnXBGnAp7VVjz76qGlA3aRJExOYaA2cTtcgJ7U2OFp+bcyt6VBt5O49Arc2rNdUqTftyq+pRw0AtUG8BrH62moNoqbtdHsAUhHobm4AUjd69GjTbfmGG25Icf7Zs2etgQMHWpUrVzbd6IsUKWI1btzY+r//+z8rLi4u1ZGqN23aZEbEzpMnj5UrVy7r5ptvtlavXn3JPo4fP2717t3buuqqq8w+SpcubXXr1s2MhHy5bvc6dECfPn1MV3iXy5Vit+zHH3/cTNfRt+2OVK1Onz5t5c+f39P9Xbvd9+/f3wxhoN3ymzRpYq1Zs8bM9+4i7+52r13EvaV0XGrUqFGmi7p239fzol3o69ata912220+y+nr/+abb1rXXXedWbZgwYJmuVdffdWUNTWX69qevHt/8mW1q3+BAgWs2rVrm+72W7duTfNry0jVcDKX/pNawAQAmUUbVmv7Hk0PJR/4LzvRBtSa0tIRxlNKkQEIfrQhAhAQeqsOTfFoG6XsFAxpuZP/jvzkk09MGjGl220AyB5oQwQgS+n4O4sXLzbtknSQQ233k53ovd60ZkvbVGkbKB3gUWu5tFu8TgOQPREQAchS2rNMG4drI2odULFWrVqSnehYSDokgZbdfT82HWhS74PmPco1gOyFNkQAAMDxaEMEAAAcj4AIAAA4Hm2I0tGt9uDBg2boe3+H4QcAAFlLWwadPXtWSpUq5bnPYUoIiNJIg6G03owRAAAEl/3795vb+1wOAVEauW+K2FRulzAJD3RxkNMFshbSZS+THhLh/2XFFRVpa99J5877va6VmCjZllP7xgTwc+IKDbW3ARvrW3Fx9vbtsPdLgsTLj/LdFW9uTECURu40mQZDYS4CImSy7BwQ2fh8uFz2uq0nueL9XteyedyB5awvuKAIiFw2AyIb61vmTjh2OOz9Yv3750rNXbLzFQAAACBDOCogGj16tBlULSoqSho0aCDr1q0LdJEAAEAQcExA9Pnnn0u/fv1k8ODBZqj9mjVrSuvWrc1tBAAAgLM5JiAaMWKE9OzZU7p37y7VqlWTsWPHmhtKTpw4McXlY2Nj5cyZMz4PAACQMzkiIIqLi5ONGzdKq1atPNN0LAJ9vmbNmhTXGT58uOTPn9/zoMs9AAA5lyMComPHjkliYqIUL17cZ7o+P3z4cIrrDBw4UE6fPu156PgFAAAgZ6Lb/WVERkaaBwAAyPkcUUNUpEgRCQ0NlSNHjvhM1+clSpQIWLkAAEBwcERAFBERIXXr1pUlS5b43JtMnzdq1CigZQMAAIHnmJSZdrnv1q2b1KtXT2644QZ59913JSYmxvQ6AwAAzuaYgOjee++Vf/75RwYNGmQaUteqVUvmz59/SUNrAADgPC7Lcthd3vyk4xBp9/vm0p57mSHzZed7mUWEB/DmrjF+r8vNXbMhbu7q5wac9X5JsOJlmXxteozny5fvsss5poYIQBYJsRFQhdu7uatITMACQbGSxJEBdHZm55zbDIjsBFTOCmeyjiMaVQMAAKSGgAgAADgeAREAAHA8AiIAAOB4BEQAAMDxCIgAAIDjERABAADHIyACAACOR0AEAAAcj4AIAAA4HgERAABwPAIiAADgeAREAADA8QiIAACA44UFugBA0HK5JNty+f9bxxUaam/XEeH+r5s72t6+z571f+WEBFv7tpL4fZnd2Hmvh0RG2tu5nc9JXJytXVuJibbWt8WyJFjxCQYAAI5HQAQAAByPgAgAADgeAREAAHA8AiIAAOB4BEQAAMDxCIgAAIDjERABAADHIyACAACOR0AEAAAcj4AIAAA4HgERAABwPAIiAADgeAREAADA8QiIAACA44UFugDZjSssTFyuALxsrgDGriEuW6u7XDbWD7F53Db27QoNtbfvcP/fJ66ICFu7tnJH+71uQvH8tva9664o/9f971hb+27c71G/182z/6KtfYeei7W1vis+0f+VE5Mku14frHD/P2eJuSNt7ftCcf/XP9LF3vtlQI2Ffq87tV87W/vO9ddpv9d1nT1va9/WRRufk/g4//ZpxYmcuvJy1BABAADHIyACAACOR0AEAAAcL6gDouHDh0v9+vUlb968UqxYMenQoYNs377dZ5mLFy/KE088IYULF5Y8efJIp06d5MiRIz7LPPnkk1K3bl2JjIyUWrVqZfFRAACAYBfUAdHy5ctNsPPTTz/JokWLJD4+Xm699VaJiYnxLPP000/LN998IzNnzjTLHzx4UDp27HjJth566CG59957s/gIAABAdhDUvczmz5/v83zy5Mmmpmjjxo1y0003yenTp2XChAkyffp0adGihVlm0qRJUrVqVRNENWzY0EwbNWqU+fvPP//Ir7/+GoAjAQAAwSyoa4iS0wBIFSpUyPzVwEhrjVq1auVZpkqVKlK2bFlZs2aNrX3FxsbKmTNnfB4AACBnyjYBUVJSkvTt21eaNGki119/vZl2+PBhiYiIkAIFCvgsW7x4cTPPbvul/Pnzex5lypSxtT0AABC8sk1ApG2JfvvtN5kxY0aW7G/gwIGmRsr92L9/f5bsFwAAZL2gbkPk1rt3b/n2229lxYoVUrp0ac/0EiVKSFxcnJw6dcqnlkh7mek8O7RHmj4AAEDOF9Q1RJZlmWBo9uzZsnTpUqlQoYLPfO1KHx4eLkuWLPFM0275+/btk0aNGgWgxAAAIDsKC/Y0mfYg+/rrr81YRO52QdqmJzo62vzt0aOH9OvXzzS0zpcvn/Tp08cEQ+4eZmrXrl1y7tw5s/6FCxdk8+bNZnq1atVMGyQAAOBsQR0QjRkzxvxt3ry5z3TtWv/ggw+a/48cOVJCQkLMgIzaM6x169by4Ycf+iz/8MMPmzGK3GrXrm3+7tmzR8qXL58FRwIAAIJZWLCnzK4kKipKRo8ebR6Xs2zZsgwuGQAAyElcVlqiDphxiDRF11zaS5grPNDFATKPy2Vv9TD/Px+uKHsdGawLF/xfN8nmpdBKsre+U7n8b8rqCnEFbN8huaPt7dtGp52kE6ds7dpKTMye73PLv89oghUvy+Rr02Ncm9Zky0bVAAAAWYGACAAAOB4BEQAAcDwCIgAA4HgERAAAwPEIiAAAgOMREAEAAMcjIAIAAI5HQAQAAByPgAgAADgeAREAAHA8AiIAAOB4BEQAAMDxCIgAAIDjERABAADHCwt0AQAEGcuyt3pCvP/rxiTa2rck2VwfWc/y/5xZlsvevl3+v9eTLly0t+v4BL/XtRJtvs+tpIBdH4IZNUQAAMDx0lVDdOrUKZk9e7asXLlS/vrrLzl//rwULVpUateuLa1bt5bGjRtnXkkBAAACWUN08OBBefjhh6VkyZIydOhQuXDhgtSqVUtatmwppUuXlh9++EFuueUWqVatmnz++eeZVVYAAIDA1RBpDVC3bt1k48aNJuhJiQZJc+bMkXfffVf2798vzzzzTEaXFQAAIFO4LOvKLaSOHz8uhQsXTvNG07t8dnDmzBnJnz+/NJf2EuYKD3RxgODlstHQ1WWzWSONqp3FznvN5vvNFW6vT5IrLCxgDbqd1qg6wYqXZfK1nD59WvLly3fZ5dL0bkhvcJPTgiEAAJCzpTs8njJlisybN8/zfMCAAVKgQAHToFobWgMAAOT4gGjYsGESHR1t/r9mzRoZPXq0vPXWW1KkSBF5+umnM6OMAAAAmSrdSUxtMF25cmXzf21E3alTJ+nVq5c0adJEmjdvnhllBAAACK4aojx58phG02rhwoWmu72KiooyPc0AAAByfA2RBkA6JpF2xd+xY4fcfvvtZvrWrVulfPnymVFGAACA4Koh0jZDjRo1kn/++Ue+/PJLT48yHaOoS5cumVFGAACAwI9DpCZOnCh33nmnaTztRIxDBKQR4xAhqzAOkX8YhyhFaX43TJs2zdymQ7vXv/nmm/LHH3+kdVUAAICgluYQdenSpXLy5EkzBtHcuXPl9ddfl+LFi5tao/bt20vTpk0lJMTmrzsAGfPL19a+7X2O7fxydkVE2Np3Usz5wPxqzqa/nLM926+5jXOeaK820lbJea9minRd+QoWLCj33XeffPHFF3Ls2DF5//33Tc+yrl27SrFixeSBBx6QWbNmSUxMTOaUFgAAIBP4/VMwIiJCbrvtNvnwww/N2ETz5883vcxee+01GTFiRMaWEgAAIBgaVadHfHy8hIfnrIbHNKpGliJl5hdSZsiqz5krNNTevm2sb8XF2du3w96rCWlsVJ3uK5fGT5oW++GHH+To0aOSlPS/i4jL5TJd8XNaMAQAAHK2dP8U7Nu3r9x///2yZ88eM2q11pq4H6lFXhnhjTfeMEGXlsHt4sWL8sQTT5jxkLQ8eiuRI0eOeObrqNqa2itVqpRERkZKmTJlpHfv3qbGBwAAwK8aoqlTp8pXX33lGaE6q6xfv14++ugjqVGjhs90vaGs9nybOXOmCco02OnYsaOsWrXKzNeeb9oLbujQoVK0aFHZtWuXCaBOnDgh06dPz9JjAAAAOSQg0qCjYsWKkpXOnTtnerKNGzfOBDZumg+cMGGCCWxatGhhpk2aNEmqVq0qP/30kzRs2ND0jHvsscc865QrV04ef/xxefvtt1PdZ2xsrHm4UaMEAEDOle6U2SuvvCKvvvpqlt7IVWt02rZtK61atfKZrrcL0Qbc3tOrVKkiZcuWlTVr1qS4rYMHD5oarmbNmqW6z+HDh/ukAzXVBgAAcqZ0B0T33HOPGaBRxx2qXr261KlTx+eR0WbMmCGbNm0yAUpyhw8fNt3/CxQo4DNdB4zUed70Pmu5cuWSq666yrR1Gj9+fKr7HThwoKmBcj90aAEAAJAzpTtl1q1bN1MzowM0auChjZwziwYhTz31lCxatEiioqJsbWvkyJEyePBg2bFjhwl2+vXrZ8ZQuhxtgK0PAACQ86U7INIGzAsWLDC36shsGnhp137vmqfExERZsWKFfPDBB6YccXFxcurUKZ9aIu1lVqJECZ9t6XN9aEqtUKFCcuONN8rLL78sJUuWzPTjAAAAOSwg0rY0md293q1ly5ayZcsWn2ndu3c3Qc1zzz1nyqJjHi1ZssR0t1fbt2+Xffv2SaNGjS67XffYSd6NpgEAgHOlOyB65513ZMCAATJ27Fhzq47MlDdvXrn++ut9puXOnduMOeSe3qNHD5P+0lofDdT69OljgiHtYaa+++47U2NUv359M07R1q1b5dlnn5UmTZpkevkBAEAODYi07dD58+elUqVKppFy8lGpdXyfrKRtg3SsIa0h0hqf1q1b+7QNio6ONt31dbwina+1SjpO0fPPP5+l5QQAADnoXmZTpky5YqPrnIh7mSFLcS8zv3AvM6QL9zJzhITMupdZTg14kEMFMrAIJBtBjSvE3mtmJ6gJyZ3L1r6ti/63C7QSbe1awzEJGId9wQXD62Yl2XvNXS7OWbBJ01UzJiYmXRtN7/IAAABBHxBVrlzZ3Fj10KFDl11GM286XlCbNm1k1KhRGVlGAACATJWmlNmyZcvkhRdeMLftqFmzptSrV8/cPV4HS9RRq7dt22ZulREWFmYGPXzkkUcyt9QAAABZHRBde+218uWXX5rxffSu8itXrpTVq1eb+5kVKVJEateubXpyae1QqN2GZgAAAMHey8yp6GWWTdGoOusbVUdHB6xRdeLxk36vayXabFVtt5earX1zGc9yIfZ+/NvppWYlxNvat9PeLwlp7GVmr38tAABADkBABAAAHI+ACAAAOB4BEQAAcDwCIgAA4HjpDoj0DvFDhgwxXfABAAAcGRD17dtXvvrqK6lYsaLccsstMmPGDHMXeQAAAEcFRJs3b5Z169ZJ1apVpU+fPlKyZEnp3bu3bNq0KXNKCQAAEIxtiOrUqWPuWXbw4EEZPHiwjB8/XurXry+1atWSiRMnmnubAQAA5Jhbd6QkPj5eZs+eLZMmTTI3dW3YsKH06NFD/v77b3Pfs8WLF8v06dMlR4587NTRjx02YnMg2R0t2tZI1eF+XxaMkHx5/V43qWgBW/t2xZz3f+W4OFv7FpsjXVtJNn5E2r0kBXKUbTsC+cPb9msWwNtcBfI7zAreypJ0X/k0LaZB0GeffSYhISHywAMPyMiRI6VKlSqeZe666y5TWwQAAJAdpDsg0kBHG1OPGTNGOnToIOHhl97Xq0KFCtK5c+eMKiMAAEBwBUR//vmnlCtXLtVlcufObWqRAAAAsoN0NzS4+eab5fjx45dMP3XqlOmKDwAAkOMDor1790piCo0HdSyiAwcOZFS5AAAAgi9lNnfuXM//FyxYIPnz5/c81wBpyZIlZhRrAACAHBsQaQNq5XK5pFu3bj7ztGG1BkPvvPNOxpcQAAAgWAKipKQkTw+y9evXS5EiRTKzXAAAAMHby2zPnj2ZUxIAAIBgDoj0Fh29evWSqKgo8//UPPnkkxlVNgAAgCzhstJw0zFNk23YsEEKFy5s/n/ZjblcZpyinOjMmTOmIXlzVwcJc106GCWCFLfuyPpbdxQsELBbd8ju/X6vamXnW3fYxa07svz2F64w/79HrIR4ybasrD9nCVa8LJOv5fTp05IvX77LLheW3jQZKTMAAJDTZM+fzwAAAIEMiDp16iRvvvnmJdPfeustufvuuzOqXAAAAMEbEK1YsUJuv/32S6a3adPGzAMAAMhu0t168ty5cxIREXHJdB2cURse53QhkRES4ooIwI5Dsm3jQTtld9net431Q0Pt7TvM/8bJrhQ+Y+lh5Yrye934onls7Xtfm2i/1331PzNs7fujp/7j97rR++1dv1znzttaX2L9b9RtJSTY27ed9W02Bk9Dv55Ma8huh2Vz3y4714fQwH0fWIFsyO4nlxUiEnvl5dL9qlavXl0+//zzS6bPmDFDqlWrlt7NAQAABFy6Q9SXX35ZOnbsKLt375YWLVqYaXofs88++0xmzpyZGWUEAAAIroCoXbt2MmfOHBk2bJjMmjVLoqOjpUaNGrJ48WJp1qxZ5pQSAAAgE/mViGzbtq2sWrVKYmJi5NixY7J06dJMC4YOHDgg9913nxkUUoMvTdnpIJHe+cxBgwZJyZIlzfxWrVrJzp07U9xWbGys1KpVy7RL2bx5c6aUFwAAZD9+t8zauHGjTJs2zTx+/vlnyQwnT56UJk2amAbb33//vWzbtk3eeecdKViwoE93f72dyNixY2Xt2rWSO3duad26tVy8ePGS7Q0YMEBKlSqVKWUFAAAOSpkdPXpUOnfuLMuWLZMCBf4dZv/UqVNy8803m4bVRYsWzbDC6XhHZcqUkUmTJnmmed86RGuH3n33XXnppZekffv2Ztonn3wixYsXN2k9LaebBlQLFy6UL7/80vwfAADA7xqiPn36yNmzZ2Xr1q1y4sQJ8/jtt99Ml/uMvrHr3LlzpV69embAx2LFiknt2rVl3LhxPrcROXz4sEmTuen9xho0aCBr1qzxTDty5Ij07NlTpk6dKrly5UrTvjW9psfk/QAAADlTugOi+fPny4cffihVq1b1TNPu9qNHj87wmhe9UeyYMWPk6quvlgULFshjjz1mgq4pU6aY+RoMKa0R8qbP3fO0FunBBx+URx991ARXaTV8+HATXLkfWlMFAABypnQHRElJSaZNT3I6TedlJN1enTp1TI82rR3q1auXqenR9kJp9f7775sarYEDB6Zr37q83hnX/di/3/+7aAMAgBwWEOnYQ0899ZQcPHjQpyfY008/LS1btszQwmnPseSDPWrN1L59+8z/S5Qo4UmJedPn7nnaA07TZ5GRkRIWFiaVK1c207W2qFu3bpfdty6fL18+nwcAAMiZ0h0QffDBB6Y9Tfny5aVSpUrmoQ2ddZrWxmQk7WG2fft2n2k7duyQcuXKmf/rfjXw0YEh3bQc2tusUaNG5rn2QPvll19MN3t9fPfdd2a6jrb9+uuvZ2h5AQCAQ3qZaVuaTZs2mYEY//jjD0+tjXfD5oyitU6NGzc2KbN77rlH1q1bJx9//LF5KB1PqG/fvjJ06FDTzkgDJB1JW7vWd+jQwSxTtmxZn23myfPvfZo0kCtdunSGlxkAAGQ/ft1dTgORW265xTwyU/369WX27NmmPc+QIUNMwKPd7Lt27eoztpAOEKnti7T7f9OmTU3D76go/29uCQAAnCVNAZGmndIqo7ve33HHHeaRWnCmwZI+0kJTfdnxbr0AACDAAdHIkSPTtDENTjI6IAIAAAiKgEgHQMS/kmLjJMnlsBoml993eAk4V4grkDv3f12b5XaF+ZUNN8KPRdrad8VTxfxed/SGe2ztO+/P/l+rrPMXbO3bio+3t76dmuskm9ckK2OHTEnfri1HXluS4uIyrCxInWWl7bPp9xU7Li7O9ABLSEjwdxMAAABBId0B0fnz56VHjx7mFhjXXXedZ0wgvaXHG2+8kRllBAAACK6ASHt86bg+enNX755c2u1ex/YBAADIbtLd0EDvIq+BT8OGDU0jajetLdq9e3dGlw8AACD4aoj++ecfc+f55HQsIO8ACQAAIMcGRHoPsHnz5nmeu4Og8ePHe26XAQAAkKNTZnobjTZt2si2bdtMD7P33nvP/H/16tWyfPnyzCklAABAMNQQ/fbbb+av3hpDb5KqwVD16tVl4cKFJoWmd5SvW7duZpYVAAAgsDVENWrUMPcWe/jhh6Vz584ybty4zCkRAABAsNYQaTpMe5L1799fSpYsKQ8++KCsXLkyc0sHAAAQTAHRjTfeKBMnTpRDhw7J+++/b27n0axZM7nmmmvkzTfflMOHD2duSQEAAIKll1nu3Lmle/fupsZox44dcvfdd8vo0aOlbNmycuedd2ZOKQEAADKRrbt2Vq5cWV544QV56aWXJG/evD7d8QEAALILv2+JvWLFCpNC+/LLLyUkJETuuecec48zAACAHB0QHTx4UCZPnmweu3btksaNG8uoUaNMMKSpNEewLP1HHMVKDNy+bY5+Hsiii8vG+8RuuRP934CVkGBr1yFHjvu9bt54e/u2zp7zf12b+xYryebqAbyu2Cy7E1lJIQHcucO+g7Lo9UpzQKSDMS5evFiKFCkiDzzwgDz00ENy7bXX2ikiAABAUEhzQBQeHi6zZs2SO+64Q0JDQzO3VAAAAMEYEM2dOzdzSwIAABAgAUyCAgAABAcCIgAA4HgERAAAwPEIiAAAgOMREAEAAMcjIAIAAI5HQAQAAByPgAgAADgeAREAAHA8AiIAAOB4BEQAAMDxCIgAAIDjpfnmrkBAWFbg9u1y2VvfSsqokvixa/9/67gk0d7OY2P93/e587Z2nZSYGLDzZSVZ2fb9Aj9wvnIcaogAAIDjERABAADHIyACAACOF9QBUWJiorz88stSoUIFiY6OlkqVKslrr70mlle7Ev3/oEGDpGTJkmaZVq1ayc6dOy/Z1rx586RBgwZmmYIFC0qHDh2y+GgAAECwCupG1W+++aaMGTNGpkyZItddd51s2LBBunfvLvnz55cnn3zSLPPWW2/JqFGjzDIaOGkA1bp1a9m2bZtERUWZZb788kvp2bOnDBs2TFq0aCEJCQny22+/BfjoAABAsHBZ3tUtQeaOO+6Q4sWLy4QJEzzTOnXqZGp5pk2bZmqHSpUqJf3795dnnnnGzD99+rRZZ/LkydK5c2cT/JQvX15effVV6dGjR5r3HRsbax5uZ86ckTJlykhzaS9hrvAMPlIEJbu9zALJZaOXWYi94w7Jlcv/fefNY2vficeO21jZXu86epkhywTv13ZQSrDiZZl8beKDfPnyZc+UWePGjWXJkiWyY8cO8/yXX36RH3/8Udq0aWOe79mzRw4fPmzSZG5ae6SpsTVr1pjnmzZtkgMHDkhISIjUrl3bpNZ0/SvVEA0fPtxsy/3QYAgAAORMQR0QPf/886aWp0qVKhIeHm4Cmr59+0rXrl3NfA2GlNYIedPn7nl//vmn+fvKK6/ISy+9JN9++61pQ9S8eXM5ceLEZfc9cOBAE026H/v378/EIwUAAIEU1G2IvvjiC/n0009l+vTppg3R5s2bTUCkabJu3bqlaRtJSf9WQ7/44osm3aYmTZokpUuXlpkzZ8ojjzyS4nqRkZHmAQAAcr6gDoieffZZTy2Rql69uvz1118mnaUBUYkSJcz0I0eOmFSYmz6vVauW+b97erVq1TzzNdCpWLGi7Nu3L4uPCAAABKOgTpmdP3/etP3xFhoa6qn10V5lGhRpOyPvxs9r166VRo0amed169Y1AdD27ds9y8THx8vevXulXLlyWXYsAAAgeAV1DVG7du3k9ddfl7Jly5qU2c8//ywjRoyQhx56yMx3uVwmhTZ06FC5+uqrPd3uNaXmHmdIW5Q/+uijMnjwYNMwWoOgt99+28y7++67A3p8AAAgOAR1QPT++++bAOfxxx+Xo0ePmkBH2/zoQIxuAwYMkJiYGOnVq5ecOnVKmjZtKvPnz/eMQaQ0AAoLC5P7779fLly4YHqhLV261DSuBgAACOpxiIKJpuK0+z3jEDkI4xD5hXGI/N6AvfXhHHxtZ8o4REFdQwQ49qJjNxiz8eVqJdlrWmglJPi/cly8rX3bCWocHdDwBeusH0zIfo2qAQAAsgIBEQAAcDwCIgAA4HgERAAAwPEIiAAAgOMREAEAAMcjIAIAAI5HQAQAAByPgAgAADgeAREAAHA8AiIAAOB4BEQAAMDxCIgAAIDjERABAADHIyACAACOFxboAgBIgWXZW9/lkoBJSvJ/3YQECRjLRrmz+/sFADVEAAAABEQAAMDxCIgAAIDjERABAADHIyACAACOR0AEAAAcj4AIAAA4HgERAABwPAIiAADgeAREAADA8QiIAACA4xEQAQAAxyMgAgAAjkdABAAAHC8s0AUAEGSsJHurW5bf67oSE+3tO8n/fQecjdcNgH3UEAEAAMcjIAIAAI5HQAQAABwvoAHRihUrpF27dlKqVClxuVwyZ86cS9oiDBo0SEqWLCnR0dHSqlUr2blzp88yJ06ckK5du0q+fPmkQIEC0qNHDzl37pzPMgsWLJCGDRtK3rx5pWjRotKpUyfZu3dvlhwjAAAIfgENiGJiYqRmzZoyevToFOe/9dZbMmrUKBk7dqysXbtWcufOLa1bt5aLFy96ltFgaOvWrbJo0SL59ttvTZDVq1cvz/w9e/ZI+/btpUWLFrJ582YTHB07dkw6duyYJccIAACCn8uy0yUkA2kN0ezZs6VDhw7muRZLa4769+8vzzzzjJl2+vRpKV68uEyePFk6d+4sv//+u1SrVk3Wr18v9erVM8vMnz9fbr/9dvn777/N+rNmzZIuXbpIbGyshIT8G/998803JkjSaeHh4Wkq35kzZyR//vzSXNpLmCtt6wAB43IFbtcREX6vGxIZaWvfiediAta7zrbguBQjG3zGeK+kT4IVL8vkaxNDaDYp27Uh0pqdw4cPmzSZmwYkDRo0kDVr1pjn+lfTZO5gSOnyGvhojZKqW7eueT5p0iRJTEw0L8jUqVPNcqkFQxosaRDk/QAAADlT0AZEGgwprRHyps/d8/RvsWLFfOaHhYVJoUKFPMtUqFBBFi5cKC+88IJERkaaAEprj7744otU9z98+HATgLkfZcqUyeAjBAAAwSJoA6KMooFRz549pVu3bia1tnz5comIiJD//Oc/qQ4gN3DgQFOb5H7s378/S8sNAACyTtCOVF2iRAnz98iRI6aXmZs+r1WrlmeZo0eP+qyXkJBgep6519cG21rDow203aZNm2ZqfDStpr3PUqK1SfoAAAA5X9DWEGmqS4OaJUuWeKZpOx4NYho1amSe699Tp07Jxo0bPcssXbpUkpKSTFsjdf78eU9jarfQ0FDzV5cDAAAIaECk4wVpV3h9uBtS6//37dtnep317dtXhg4dKnPnzpUtW7bIAw88YHqOuXuiVa1aVW677TaTElu3bp2sWrVKevfubXqg6XKqbdu2JlU2ZMgQM4bRpk2bpHv37lKuXDmpXbt2IA8fAAAEiYAGRBs2bDBBiTsw6devn/m/DsaoBgwYIH369DHjCtWvX98EUNqtPioqyrONTz/9VKpUqSItW7Y03e2bNm0qH3/8sWe+jj80ffp0M+ijblsDKE2F6XZ0sEcAAICgGYco2DEOEbIVxiFKP8YhQnowDlGOG4coaBtVBxt33Jgg8SK8FxH0AhgQWf7vO8TGuirRivd/ZQIipAsBUXZhvre9vscvh4Aojc6ePWv+/ijfBboowJUF8noZG6B1gaxETJItv8c103M5pMzSSHukHTx40NwgVht8J0+naTd+Hasoteo4BA/OWfbC+cp+OGfZy5kcfL40zNFgSDtbJe917o0aojTSF7F06dKpLqNvopz2RsrpOGfZC+cr++GcZS/5cuj5Sq1mKOjHIQIAAMgqBEQAAMDxCIgygI5rNHjwYG71kY1wzrIXzlf2wznLXiI5XzSqBgAAoIYIAAA4HgERAABwPAIiAADgeAREAADA8QiIMsDo0aOlfPnyEhUVJQ0aNJB169YFukiOtGLFCmnXrp0ZjVRHE58zZ84ly/z+++9y5513mkG6cufOLfXr15d9+/b5LLNmzRpp0aKFma8DlN10001y4cKFLDwSZxgzZozUqFHDMxBco0aN5Pvvv/fMf+SRR6RSpUoSHR0tRYsWlfbt28sff/zhmX/8+HG57bbbzPnWnjE6ym7v3r3NiLvIHK+88or5bHk/qlSp4pl/8eJFeeKJJ6Rw4cKSJ08e6dSpkxw5csRnG/p5a9u2reTKlUuKFSsmzz77rCQkJATgaJzhwIEDct9995lzop+l6tWry4YNG3yW4br4LwIimz7//HPp16+f6a64adMmqVmzprRu3VqOHj0a6KI5TkxMjHn9NUBNye7du6Vp06bmAr5s2TL59ddf5eWXXzaBrPeHXr9kb731VhPYrl+/3nzJpjbcO/yjI7+/8cYbsnHjRnOB1outBj1bt2418+vWrSuTJk0yF+sFCxaY4ff1vCQmJpr5ek50+blz58qOHTtk8uTJsnjxYnn00UcDfGQ523XXXSeHDh3yPH788UfPvKefflq++eYbmTlzpixfvtzc7qhjx46e+XruNBiKi4uT1atXy5QpU8x5GzRoUICOJmc7efKkNGnSRMLDw82PjW3btsk777wjBQsW9CzDddGLdruH/2644QbriSee8DxPTEy0SpUqZQ0fPjyg5XI6fWvPnj3bZ9q9995r3Xfffamu16BBA+ull17K5NLhcgoWLGiNHz8+xXm//PKLOa+7du267PrvvfeeVbp06UwsobMNHjzYqlmzZorzTp06ZYWHh1szZ870TPv999/NOVuzZo15/t1331khISHW4cOHPcuMGTPGypcvnxUbG5sFR+Aszz33nNW0adNUl+G6+D85LLzLWvorR3/dtmrVyjNNI2Z9rhE1guvmvPPmzZNrrrnG1OBpVb2mN73Talqrt3btWjOvcePGUrx4cWnWrJnPL2BkDq05mDFjhqnl09RZcjpda4sqVKhgUmMp0dqIr776ypwzZJ6dO3eaNGXFihWla9euntSKXgvj4+N9roda61C2bFnP9VD/aspGP1tu+nnUNKe7ZhAZR2tP69WrJ3fffbe5rtWuXVvGjRvnmc910RcBkQ3Hjh0zF3LvD7fS54cPHw5YuXAp/VCfO3fOpGi06nfhwoVy1113mep8rdpXf/75p6edRM+ePWX+/PlSp04dadmypfkSQMbbsmWLaWuibYA01TV79mypVq2aZ/6HH35o5utDq/wXLVokERERPtvo0qWLaY9y1VVXmbYN48ePD8CROIN+WWqKSz8b2gZsz549cuONN5o7ies1T89NgQIFLns91L8pXS/d85Cx9Jqm5+nqq682aefHHntMnnzySZOqVFwXk/GqLUI6HThwwFQHr1692mf6s88+a1JpCJ6UmftcdenSxWe5du3aWZ07dzb/X7VqlVlm4MCBPstUr17dev7557Oo5M6iaZKdO3daGzZsMK9xkSJFrK1bt/qkYXbs2GEtX77cnKs6depYFy5c8NnGoUOHTGrm66+/tqpVq2Y99thjATgSZzp58qRJd2ma89NPP7UiIiIuWaZ+/frWgAEDzP979uxp3XrrrT7zY2JizOdO02nIWJrCbNSokc+0Pn36WA0bNjT/57roixoiG4oUKSKhoaGX9KLQ5yVKlAhYuZDyuQoLC/OpfVBVq1b1VPmXLFnS/E1tGWQsrVGoXLmyaUA9fPhw0yj+vffe88zXXi/661Z7tMyaNcv0MtNaJG/6WdPUjPaS+eijj8wvYm3si8yntUGabtm1a5c5D9qM4NSpU5e9HurflK6X7nnIWHpNS+16xnXRFwGRzYu5XsiXLFnik5PV5ym1g0Bgz5V2Jd2+fbvPdO2dVK5cOfN/HTpB20aktgwyl35+YmNjU5ynFX/6uNx89/oqtWWQcTTdor2U9EtTr4Xam8n7eqifJf3SdF8P9a+mSb174WoaVFOdyb9wYZ/2MEvtesZ1MZlkNUZIpxkzZliRkZHW5MmTrW3btlm9evWyChQo4NOLAlnj7Nmz1s8//2we+tYeMWKE+f9ff/1l5n/11VemCvnjjz82aZr333/fCg0NtVauXOnZxsiRI00KQHvK6DLasyIqKirVnk3wj1a3aypsz5491q+//mqeu1wua+HChdbu3butYcOGmVSanj+tttdq/EKFCllHjhwx68+bN8+aOHGitWXLFrONb7/91qpatarVpEmTQB9ajtW/f39r2bJl5vXWc9KqVSuT5jx69KiZ/+ijj1ply5a1li5das6dpmu8UzYJCQnW9ddfb9JmmzdvtubPn28VLVr0knQMMsa6deussLAw6/XXXzfXM01r5sqVy5o2bZpnGa6L/0NAlAH0DaQXAc2fa9uhn376KdBFcqQffvjBBELJH926dfMsM2HCBKty5crmw6zdh+fMmXPJdnTIBO26rRcOvZh7XxiQcR566CGrXLly5nOjX4otW7Y0wZC7bUObNm2sYsWKmYu1no///ve/1h9//OFZX7909fzkz5/fnM+rr77adDPWdi3IHNpFu2TJkuacXXXVVea595eitu96/PHHzfAJ+vm56667TBsvb3v37jXnNjo62gRTGmTFx8cH4Gic4ZtvvjFBqP5wr1Kligl8kuO6+C+X/pO81ggAAMBJaEMEAAAcj4AIAAA4HgERAABwPAIiAADgeAREAADA8QiIAACA4xEQAQAAxyMgAgAAjkdABAA2vPzyy9KrVy9b29i2bZuULl1aYmJiMqxcANKHgAhAQDz44IPicrnMQ28KWqFCBRkwYIBcvHhRsovDhw/Le++9Jy+++KKt7eiNTRs2bCgjRozIsLIBSB8CIgABc9ttt8mhQ4fkzz//lJEjR8pHH30kgwcPluxi/Pjx0rhx4wy563f37t1lzJgxkpCQkCFlA5A+BEQAAiYyMlJKlCghZcqUkQ4dOkirVq1k0aJFnvnHjx+XLl26yFVXXSW5cuWS6tWry2effeazjebNm8uTTz5papcKFSpktvfKK6/4LPPHH39I06ZNJSoqytTGLF682NRMzZkzx7PM/v375Z577pECBQqY7bRv31727t2bavlnzJgh7dq1u6Q8ffr0kb59+0rBggWlePHiMm7cOJMO06Anb968UrlyZfn+++991rvlllvkxIkTsnz5cr9eSwD2EBABCAq//fabrF69WiIiIjzTNH1Wt25dmTdvnpmvbXXuv/9+Wbdunc+6U6ZMkdy5c8vatWvlrbfekiFDhngCq8TERBNsaUCl8z/++ONLUlzx8fHSunVrE6ysXLlSVq1aJXny5DE1WHFxcSmWV4MXbftTr169S+ZpeYoUKWLKqcHRY489JnfffbepTdq0aZPceuut5jjOnz/vWUePu1atWmb/AALg/9/1HgCyVLdu3azQ0FArd+7cVmRkpKWXo5CQEGvWrFmprte2bVurf//+nufNmjWzmjZt6rNM/fr1reeee878//vvv7fCwsKsQ4cOeeYvWrTI7G/27Nnm+dSpU61rr73WSkpK8iwTGxtrRUdHWwsWLEixHD///LPZxr59+3ymJy9PQkKCOcb777/fM03LouuuWbPGZ9277rrLevDBB1M9fgCZIywQQRgAqJtvvtm0m9F0krYhCgsLk06dOnnma+3OsGHD5IsvvpADBw6Y2prY2FhT2+OtRo0aPs9LliwpR48eNf/fvn27SclpKs3thhtu8Fn+l19+kV27dpkaIm9aQ7V79+4Uy37hwgXzV9NwyXmXJzQ0VAoXLmzSfW6aRlPuMrpFR0f71BoByDoERAACRtNc2p5GTZw4UWrWrCkTJkyQHj16mGlvv/226cX17rvvmoBCl9e2OcnTWNpLzZu2D0pKSkpzOc6dO2dSc59++ukl84oWLZriOpoSUydPnrxkmZTK4z1Nn6vkZdQ0XKVKldJcbgAZhzZEAIJCSEiIvPDCC/LSSy95al+0LY82br7vvvtMsFSxYkXZsWNHurZ77bXXmgbTR44c8Uxbv369zzJ16tSRnTt3SrFixUyA5v3Inz9/itvVwCVfvnymHVFG0XZStWvXzrDtAUg7AiIAQUMbHmuKafTo0eb51VdfbRpHa2Pr33//XR555BGfwCYttPeWBi/dunWTX3/91QRZGnR519R07drV1Pho8KWNmvfs2SPLli0zvdf+/vvvywZw2ivuxx9/lIygPdo0LajbBJD1CIgABA1tQ9S7d2/TU0zbFWngorU32gNMu7NrOyDtMZYeGmBp93pNi9WvX18efvhhTy8zd/sfbZO0YsUKKVu2rHTs2FGqVq1q0nbahkhrgS5Ht6Vd79OTnrscHU5Ae59lxJhGANLPpS2r/VgPALItrSXScYm0IbWdNjt6+WzQoIE8/fTTZrwkf2mbKK0Nmz59ujRp0sTv7QDwH42qAeR4s2fPNuMKadChQdBTTz1lAg+7DZg15abjGm3ZssXWdvbt22faTxEMAYFDDRGAHO+TTz6RoUOHmsBD2wppO5133nnHdIcHAEVABAAAHI9G1QAAwPEIiAAAgOMREAEAAMcjIAIAAI5HQAQAAByPgAgAADgeAREAAHA8AiIAACBO9/8ADSvtRwXIh2MAAAAASUVORK5CYII=", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "from numpy import abs ,angle, arange, arcsin, cos, pi, sqrt, tan, zeros\n", "from scipy.fft import fft, fft2\n", "import matplotlib.pyplot as plt\n", "\n", "\n", "def y_IF(f0_min, slope, T, antenna_tx, antenna_rx, target, v=3e8):\n", " \"\"\" This function implements the mathematical IF defined in latex as\n", " y_{IF} = cos(2 \\pi [f_0\\delta + s * \\delta * t - s/2* \\delta^2])\n", " into following python code\n", " y_IF = cos (2*pi*(f_0 * delta + slope * delta * T - slope/2 * delta**2))\n", " Parameters:\n", " -----------\n", " f0_min: float\n", " the frequency at the begining of the chirp\n", " slope: float\n", " the slope with which the chirp frequency inceases over time\n", " T: ndarray\n", " the 1D vector containing time values\n", " antenna_tx: tuple of floats\n", " x, y, z coordinates\n", " antenna_rx: tuple of floats\n", " x, y, z coordinates\n", " target: tuple of floats\n", " x, y, z coordinates\n", " v: float\n", " speed of light in considered medium\n", " Returns:\n", " --------\n", " YIF: ndarray\n", " vector containing the IF values\n", " \"\"\"\n", " tx_x, tx_y, tx_z = antenna_tx\n", " rx_x, rx_y, rx_z = antenna_rx\n", " t_x, t_y, t_z = target\n", " # distance tx antenna to target\n", " distance = sqrt((tx_x-t_x)**2 + (tx_y-t_y)**2 + (tx_z-t_z)**2)\n", " # distance target to rx antenna\n", " distance += sqrt((rx_x-t_x)**2 + (rx_y-t_y)**2 + (rx_z-t_z)**2)\n", " # usually delta_t = 2*d/c, but\n", " # distance is already 2*D (TX to target + distance target to RX)\n", " # so delta = distance/v\n", " delta = distance/v\n", " YIF = cos(2 *pi *(f0_min * delta + slope * delta * T - slope/2 * delta**2))\n", " return YIF\n", "\n", "f0_min = 60e9\n", "c = 3e8\n", "# lambda ~5mm at 60GHz\n", "lambda0_max = 3e8/f0_min\n", "n_rx = 2\n", "Distance = 10\n", "k = 10e12\n", "n_samples = 64\n", "f_if = 2*k*Distance/c\n", "fs = 50e6\n", "ts = 1/fs\n", "\n", "antenna_tx = (-lambda0_max/2,0,0)\n", "antenna_rx = (0,0,0)\n", "T = arange(0, n_samples*ts, ts)\n", "t_chirp_to_chirp = 1.2e-6\n", "n_chirps = 32\n", "\n", "# define empty 2D array to be filled by 2D FFT\n", "cube2D = zeros((n_chirps, n_samples))\n", "\n", "# define 3 targets\n", "# one target is a 2-uple (distance, speed)\n", "targets = [(160, 600), (160, 200), (300, 200)]\n", "\n", "for chirp_i in range(n_chirps):\n", " d0, v0 = targets[0]\n", " d0_t = d0 + v0*t_chirp_to_chirp*chirp_i\n", "\n", " d1, v1 = targets[1]\n", " d1_t = d1 + v1*t_chirp_to_chirp*chirp_i\n", "\n", " d2, v2 = targets[2]\n", " d2_t = d2 + v2*t_chirp_to_chirp*chirp_i\n", "\n", " f_if = 2*k*d0_t/c\n", "\n", " target_0_t = (0, d0_t, 0)\n", " target_1_t = (0, d1_t, 0)\n", " target_2_t = (0, d2_t, 0)\n", "\n", " # sanity check\n", " assert f_if < 1/ts/2\n", " assert v0 < lambda0_max/4/t_chirp_to_chirp\n", "\n", " YIFi0 = y_IF(f0_min, k, T, antenna_tx, antenna_rx, target_0_t)\n", " YIFi1 = y_IF(f0_min, k, T, antenna_tx, antenna_rx, target_1_t)\n", " YIFi2 = y_IF(f0_min, k, T, antenna_tx, antenna_rx, target_2_t)\n", "\n", " YIFi = YIFi0 + YIFi1 + YIFi2\n", "\n", " cube2D[chirp_i, :] = YIFi\n", "\n", "Z_fft2 = abs(fft2(cube2D))\n", "\n", "Data_fft2 = Z_fft2[0:n_chirps//2,0:n_samples//2]\n", "\n", "# change scale for 2D plot to display range and velocity values\n", "# https://stackoverflow.com/a/53746824\n", "# range formula\n", "ranges = linspace(0, fs*c/2/k *2, n_samples)\n", "no_labels = 10 # how many labels to see on axis x\n", "step_x = int(n_samples / (no_labels - 1)) # step between consecutive labels\n", "x_positions = arange(0, n_samples, step_x) # pixel count at label position\n", "x_labels = ranges[::step_x] # labels you want to see\n", "x_labels = [int(x) for x in x_labels]\n", "plt.xticks(x_positions, x_labels)\n", "\n", "# we plot up to 2x max un-ambigous speed as we have no negative speeds\n", "# to display so the positive speeds wrap-up and are displayed above\n", "# the max un-ambigous speed\n", "# hence speeds from 0 to max-speed *2\n", "speeds = linspace(0, lambda0_max/4/t_chirp_to_chirp *2,\n", " n_chirps)\n", "no_labels_y = 10 # how many labels to see on axis x\n", "step_y = int(n_chirps / (no_labels_y - 1)) # step between consecutive labels\n", "y_positions = arange(0, n_chirps, step_y) # pixel count at label position\n", "y_labels = speeds[::step_y] # labels you want to see\n", "y_labels = [int(y) for y in y_labels] # remove decimlas for cleaner plotting\n", "plt.yticks(y_positions, y_labels)\n", "\n", "# change scale for 2D plot\n", "\n", "plt.imshow(Data_fft2)\n", "plt.xlabel(\"Range (m)\")\n", "plt.ylabel(\"Velocity (m/s)\")\n", "plt.title('Velocity-Range 2D FFT')" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [], "source": [ "# non regression hook, check the index and max value haven't changed\n", "from numpy import argmax, unravel_index\n", "flat_index = argmax(Data_fft2)\n", "\n", "# convert to 2D index\n", "i, j = unravel_index(flat_index, Data_fft2.shape)\n", "assert (i,j) == (3,14)\n", "assert Data_fft2[i, j]==879.8465949086351" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "colab": { "name": "FMCW-Radar-101_Intro.ipynb", "provenance": [], "toc_visible": true }, "kernelspec": { "display_name": "Python 3 (ipykernel)", "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.10.9" } }, "nbformat": 4, "nbformat_minor": 4 }