DFT.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# DFT\n",
    "\n",
    "DFT(Discrete Fourier Transform)，离散傅里叶变化，可以将离散信号变换到频域，它的公式非常简单:\n",
    "\n",
    "$$\n",
    "X[k] = \\sum_{n=0}^{N-1} x[n] e^{-j2\\pi kn/N}\n",
    "$$\n",
    "\n",
    "$X[k]$：离散频率下标为k时的频率大小\n",
    "\n",
    "$x[n]$: 离散时域信号序列\n",
    "\n",
    "$N$: 信号序列的长度，也就是采样的个数\n",
    "\n",
    "如果你刚接触DFT，并且之前没有信号处理的相关经验，那么第一次看到这个公式，你可能有一些疑惑，为什么这个公式就能进行时域与频域之间的转换呢？\n",
    "这里，我不打算去解释它，因为我水平有限，说的不清楚。相反，在这里我想介绍，作为一个程序员，如何如实现DFT"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 从矩阵的角度看DFT\n",
    "\n",
    "DFT的公式，虽然简单，但是理解起来比较麻烦，我发现如果用矩阵相乘的角度来理解上面的公式，就会非常简单，直接上矩阵：\n",
    "$$\n",
    "\\begin{bmatrix}\n",
    "s_0^0 & s_0^1 & \\cdots & s_0^{N-1} \\\\\n",
    "\\vdots & \\vdots  & \\vdots &  \\vdots\\\\\n",
    "s_k^0 & s_k^1 & \\cdots & s_k^{N-1} \\\\\n",
    "\\vdots & \\vdots  & \\ddots &  \\vdots\\\\\n",
    "s_{N-1}^0 & s_{N-1}^1 & \\cdots & s_{N-1}^{N-1} \\\\\n",
    "\\end{bmatrix}\n",
    "\\begin{bmatrix}\n",
    "x[0] \\\\\n",
    "x[1] \\\\\n",
    "\\vdots\\\\\n",
    "x[n] \\\\\n",
    "\\vdots \\\\\n",
    "x[N-1]\n",
    "\\end{bmatrix} = \\begin{bmatrix}\n",
    "X[0] \\\\\n",
    "X[1] \\\\\n",
    "\\vdots\\\\\n",
    "X[k] \\\\\n",
    "\\vdots \\\\\n",
    "X[N-1]\n",
    "\\end{bmatrix}\n",
    "$$\n",
    "\n",
    "$S$矩阵中的每一行都是一个$S_k$向量，$S_k = e^{-j2\\pi kn/N}, n=0,1,\\cdots,N-1$，进一步简化上面的表示，得到：\n",
    "$$\n",
    "\\begin{bmatrix}\n",
    "\\cdots & S_0 & \\cdots \\\\\n",
    "       & \\vdots &     \\\\\n",
    "\\cdots & S_k & \\cdots \\\\\n",
    "       & \\vdots &     \\\\\n",
    "\\cdots & S_{N-1} & \\cdots \\\\\n",
    "\\end{bmatrix}\n",
    "\\begin{bmatrix}\n",
    "x[0] \\\\\n",
    "x[1] \\\\\n",
    "\\vdots\\\\\n",
    "x[n] \\\\\n",
    "\\vdots \\\\\n",
    "x[N-1]\n",
    "\\end{bmatrix} = \\begin{bmatrix}\n",
    "X[0] \\\\\n",
    "X[1] \\\\\n",
    "\\vdots\\\\\n",
    "X[k] \\\\\n",
    "\\vdots \\\\\n",
    "X[N-1]\n",
    "\\end{bmatrix}\n",
    "$$\n",
    "\n",
    "OK，通过上面的表示，我们很容易将DFT理解成为一种矩阵相乘的操作，这对于我们编码是很容易的。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# IDFT\n",
    "IDFT(Inverse Discrete Fourier Transform), 傅里叶逆变换，可以将频域信号转换到时域中"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Talk is cheap, show me the code\n",
    "\n",
    "根据上面的理解，我们只需要构建出$S$矩阵，然后做矩阵相乘，就等得到DFT的结果\n",
    "\n",
    "在这之前，我们先介绍如何生成正弦信号，以及如何用scipy中的fft模块进行DFT操作，以验证我们的结果是否正确"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 正弦信号\n",
    "$$\n",
    "x[n] = A\\cos(2\\pi fnT + \\phi)\n",
    "$$\n",
    "\n",
    "A: 幅度\n",
    "\n",
    "f: 信号频率\n",
    "\n",
    "n: 时间下标\n",
    "\n",
    "T: 采样间隔, 等于 1/fs，fs为采样频率\n",
    "\n",
    "$\\phi$: 相位\n",
    "\n",
    "下面介绍如何生成正弦信号"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX8AAAD8CAYAAACfF6SlAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJztvXl0I9d95/v5ASAALgB3kGyyu9mrpNZut2Ut3i1nZCexnDd2Yk9mosxzrHmT58zMy/LGec7xmXGSc5ztJZM3nkwcZ1EySRzbWazITmxZ3jdJbWtttXpnd7NJAtwBLtjv+6OqQIpNNhfUcguszzl9mgCKqPqxbn3v7/7u7/6uKKUICAgICNhdhLy+gICAgIAA9wnEPyAgIGAXEoh/QEBAwC4kEP+AgICAXUgg/gEBAQG7kED8AwICAnYhgfgHBAQE7EIC8Q8ICAjYhQTiHxAQELALiXh9ARvR09OjhoeHvb6MgICAAF/x/e9/f0op1bvZcdqK//DwMCdOnPD6MgICAgJ8hYhc2spxQdgnICAgYBcSiH9AQEDALiQQ/4CAgIBdSCD+AQEBAbuQQPwDAgICdiG2iL+IPCAip0XknIh8aJ3P94nIV0XkGRF5XkTeYcd5AwICAgJ2Rt3iLyJh4OPA24FjwPtE5Niaw34F+LRS6k7gvcD/qPe8AQEBAQE7xw7P/y7gnFLqglKqCHwKeHDNMQpImj+3A2M2nHddlFL8+udf4qmLM06dQitOT+T4o29c4Gw65/WluMJ3z0/zx9+6yGSu4PWlOI5Sis8/P85fPnmJpWLZ68txnHypwl8/dZnPPXuVSrXxt5edWijwp9++yDfPTnpyfjsWeQ0CV1a9HgVeu+aY/wJ8SUR+DmgF7l/vi0TkYeBhgH379u3oYkZnl/mrJy/zR9+8yG+/53be/eqhHX2PH3jywjT/5o+folip8juPn+avP3A3d+7r9PqyHOMvn7zEh//+RQD+8OvneeznXkcqGff4qpzjI587yV98z1iv85kTo3zq4buJN4U9vipnKFeq/NSfPFVz2r55dorffs/tHl+Vc0wtFHjXx7/N6OwyAB/5kWP876874Oo12OH5yzrvre223wf8mVJqCHgH8Bcics25lVKfUEodV0od7+3ddHXyuuztauHpX7mfuw928V8ePdmwHmKxXOU//+3zDHY28/n/8Dq6W2P80mefb1iPKZ3N89F/fIk3HO3l0//uHrL5Eh997CWvL8sxvnN+ir/43iX+7X3D/Lf33sGzV+b402+PeH1ZjvHIdy/x1MUZPva/3cq/f9MhPvv9Ub76csbry3KM3/inl0ln8/zVB17L24718bF/epnL00uuXoMd4j8K7F31eohrwzrvBz4NoJT6LhAHemw497q0RCP82rtuZbFYrnlOjcYXXhhnZHqJX/nhm7h5Tzv/zztu4lxmgS+enPD60hzhk9+8QKWq+LUHb+GuA128/3UHeOz5cS5MLnh9aY7wP756nlQixn9+4EYevGOQt9yY4g+/cZ58qeL1pdlOsVzlE984z72HuvmJ1+zl/7r/KHu7mvn4V895fWmOcHVumb975io/dc8w9x7q4dfedQsAn/jmeVevww7xfxo4IiIHRCSKMaH76JpjLgNvBRCRmzDE39FA1+FUG2++IcWnnrrckN7wXz11mQM9rbz5hhQAD9zSz2BHM3/z9JVNftN/FMtVPvv9UX7o5j72dbcA8NC9w4RDwt+caDx7M9k83z4/xb967b5amOf9rzvA3FKJf36x8Tr3r53OkM4W+JnXH0BEiEZCPHTPMCcuzXKmAeeyPnPiClWl+Ol7hwHoS8b5kdsH+Idnxlzt3OsWf6VUGfgg8EXgFEZWz0kR+aiIvNM87BeAD4jIc8BfAz+tlHJckX/szkEyuQJPXpx2+lSuMrNY5MTIDO+8fQ+hkBF1C4eEd925h2+enWR6obFCXU9enGZ2qcSP3bkyf5NKxHnDkR6+8MI4LjQlV/mnFydQCn7ktoHae/cc7GagPc5jz497eGXO8OhzY3S3RnnDkZVQ7zvv2ENI4LHnHMsN8YzHX0rz6n2d7O1qqb33L181xEKh7Gqoy5Y8f6XUF5RSR5VSh5RSv26+9xGl1KPmzy8ppe5TSt2ulLpDKfUlO867GW+9KUU0HOJrp72ZTXeKr76coaoM+1bztmP9VBV8+3xjdXZfOz1JNBLivsPdr3j/LTemuDKzzMWpRY+uzBk+/8I4R/vaOJxK1N4LhYR/cXM/3zg7yXKxcUI/lariG2cmeetNKSLhFTlKJeIcH+7iy6caK+4/Pr/MybEs9x/re8X7dx/sJhmP8PUz7mlVQ6/wbYlGuHNfB98+N+X1pdjKEy+nSSVi3LKn/RXv3zrYTntzE9/yKHXMKb52OsPdB7tpib4yOe2NR43Oz80Hxmnmloo8PTLDA7cMXPPZG2/opViu8oPLsx5cmTO8cHWebL7M645cm+Dx+sM9nJrIMrtY9ODKnMHqzO6/6ZXiHw4J9xzq5ptnp1wbyTa0+APcd7iHk2NZZhqkAZUrVb5xZoq33JiqhXwswiHhvsPuNiCnuTKzxPnJRd58w7XisK+7hYM9rQ01sjsxMotScN+h7ms+O76/k3BI+N6FxhnZWY7ZvevYe+/hbpSioez9+ukM+7tbONTbes1nrzvSy9W5ZUZcyvrZBeJvNKrvNkgo5HQ6x0KhzD3rPCwArz/Sy/h8nvMNkgVj3bfXH1k/OewNR3v53oVpSpWqm5flGE+PzNAUFm7f23HNZ4l4E7cMtvPkhcZZwPits1PcNJCkpy12zWe3DXXQEg3znQZ5dpVS/ODyHHcNdyFybYb86w8bbdytkXvDi/9tQx20RsN890JjhH6euTwHwJ1711/MdfdBo1M4MdIYoYFnrsyRjEc42NO27uev3t9JoVzl9ERjZIU8NTLDbUMdGy7muvtgF89emWuIuH+xXOX7l2fX9foBmsIh7jrQxXfON8azOzK9xMxikVftX//Z3d/dwkB7nCddqk7Q8OLfFA5x61A7L4zOe30ptvDM5Tl62qLs7Wpe9/Ph7haS8QjPX20Me5+9MsftezuuCXFZ3GF6yM+Nzrl5WY6wXKzwwug8rxnu2vCY1x7oolipNoS9Z9I5iuVq7R6ux/H9nZyfXCSbL7l4Zc7wg0uGQ/aqDVbhiwi3D3XwokvPbsOLPxje/6lxo6H5nWcuz3LH3s51h41gNKDbhjp4vgHEYalY5kw6d11xGOpspqs1ynNX/G/vC1fnKVcVrxneuETHrYPG38ItgXCSF0wbbhtq3/CYWwaNz05ezbpyTU7yg8uzJGIRjqTWH8UC3DrUzsj0EvNLznd2u0L8bxlsp1ip+n7ByNxSkQtTi9y5b2MxBKMBvTye8/1q0BevZqlU1XXF3/CW2nnuiv/F8OSYYYMleOvRm4jRl4xxcsz/Yvj86DzJeIR9q/Ld12L9LRqhs3vm8hx37Nt4FAsrHeELLti7K8T/1gZpQC+a3s/tQ9cX/9uH2ilXFafG/S0Qlje/3uTnam7f28GZjDER7mdOjWfpbo2SSlw7+bmaW/a0+74tA7xwdY7bhjo2HMUC9LTFGGiPuyKGTlIsVzmbyV23Ywe4zRzZPX/V+ZHsrhD//V0tJOIR3zeglycMMb9xIHHd4241Owe/23tqPEt/Mr5uJshqbh1sRyk4PeHvzu6l8Sw3DSSvK4YAN+9Jcn5ywdeTvvlShdMTOW69TsjH4pbBdl4c83dbvjC1QKmiuLH/+s9ue0sT+7tbXJmj3BXiHwpJQ3hLpydy9LRFNxXDPe1x2pubeNnnGTCn0zmObvKwABztM445PeHf9NZSpcqZiQWO7UlueuzNg+1UFZzycWd3LmOI4c1bsPfWwXYuTi36emRnZaPd2L+5vcf3d+HGMp1dIf4AN/QnOJtZ8PXip9PpHDdsQQxFhCOpNl9v8FKpKs5mFrihb+PJMYvBjmZao2Ffz+lcmFykWKlybGAL4m8Kpp/j/mczxr26oW/z9nxjfwKl8HV7PjWeoyksHFxncddafufHb+d//ptXO35Nu0b8D6faWCpWGJvPe30pO6JSVZxJ57ihb3NxADjSl+BM2r+d3aXpRYrlas2rvx6hkJj2+lccXho3RqU3bUH8rc7unI/tPZteIBIS9ndvLoZHzDZwNuPfkd3LE1kO9bbRFNZHcvW5EoexRMSv3sPlmSXypeqmMUOLo31tzC+XfLuZjTVM3spIBwwP0s/ib4nhVjxDEeFQqo1zPl7FfTazwIGeVqKRzSVob2cz0UiI8z4W/9MTuS117G6ya8Tfyq0959MGtF0xPOpzb+l0OoeIMWLbCkf7E0wtFJnyaTnrC5OL7Otu2bJneLi3zbdtGQwn7MgWQnoAkXCIgz2tvm3L88slxufzWxrFusmuEf/O1ig9bVHOpv3ZgKxaPVsVQ+vB8qs3fCadY19XyzWVPDfiqGWvTye5z08ucKh3a/cW4FCqjXS2QM6HK1/zpQqXZ5Y4ktq6GB7pS9TmCfyGVXJ8q8+uW9gi/iLygIicFpFzIvKhDY75cRF5SUROishf2XHe7XI41ebbBnRhcpG+ZIzW2NbEsLctRkdLE2d82tldmFzclhhaQnLeh7X9y5UqI9Pbs9cSkvOT/rP3wuQiVcWWPX8wRjqjs8ssFf2X8WNtNXqgZ/OQnpvULf4iEgY+DrwdOAa8T0SOrTnmCPDLwH1KqZuB/1TveXfCkZR/M34uTi1sq/GICId623y5x221qhiZXtyWvX3JGM1NYS76UAxHZ5cpVdSW4v0Wh30cxrQcsO15/m0oZXQcfuPi1CLhkFx3JbMX2OH53wWcU0pdUEoVgU8BD6455gPAx5VSswBKKU+25znY20ouX/Zlbf+LU4sc2KCy5UYMd7dyyaXa4HaSzuXJl6oMb7Oz29/dwsi0/8TBCultx/Pf39VCU1h8Kf6XppcQMapYbhXrb+PHUuUXphZrk9Y6YcfVDAKrd9EeNd9bzVHgqIh8W0S+JyIP2HDebTNsppW5tVmCXcwuFpldKnFwm8PGAz0tTGTzvlsJannv27X3YG8rIz4M+6yI/9btjYRD7O9u5eKU/8RwZGqRgWR8w7LV62F1FH50Zi5Mbm8U6xZ2iP96a9HXxlUiwBHgTcD7gE+KyDUFW0TkYRE5ISInJift39BgpQH5SyAumte7nbAAUPOc/eYNW/Zux/MHo3O/PLNE2Wcbu1yYXKS7NUpHS3Rbv7e/q8WXYjgyvbil/P7VxJvC9CfjvmvL1apiZGqRg9sY1bmFHeI/Cuxd9XoIGFvnmM8ppUpKqYvAaYzO4BUopT6hlDqulDre23vttn31MtTZQkjwnXdoxTm36z1YIx3fdXaTi8QiIQaS8W393nBPK+WqYnR22aErc4bLM0vs20YIxGJfdwuXZ5Z8N4d1aXqJ4Z7t27u/u4XLPuvs0rk8y6VKw3r+TwNHROSAiESB9wKPrjnmH4A3A4hID0YY6IIN594W0UiIwc5m34V9Lk4tEA4Je7c5YWR5zhen/GXvyPQiw92t1y19ux5WmOiizzq7yzNLO5oMHO5uZalYYdJHaxuy+RLTi8Vte/5g2Ou7Z3eHIUw3qFv8lVJl4IPAF4FTwKeVUidF5KMi8k7zsC8C0yLyEvBV4JeUUp5szGlMgvpNHJYZ7Gje9tLwtliEnraY/0Y6U4s78gxrYS4f2VuqVBmfz+9I/K3Rgp+84UumIzK8A/Hf39PC1ELBVwXerswa9m7XcXMDW6aflVJfUEodVUodUkr9uvneR5RSj5o/K6XUzyuljimlblVKfcqO8+4EIyPEPw8LwJWZJYY619+2cTOGfZYBUzXDNjsRh+7WKG2xiK/i4ONzeSpVxd7OHYRBTEHxU3seqc3n7GykA/4KY16ZWSYcEgbatxfCdAO9co9cYLi7lfnlErM+SvccnV3ekTgA7PdZuufUQoFiubqjzk5EGOps9lXMvx7P0JrDuuwjMbSEeycjHT9m/FyZXWKgPU5Eo4JuFvpdkcNYscZLM/5oQMvFClMLhQ03bN+Moc5m0rm8b/YvtsRwaIed3VBnC6Oz/ri3YMT7gR3d32gkxJ4Of81hXZlZpjcR23LZjtXs7/Zf9lo9jpvT7DrxH+wwHrKrPvEOr87VK4bNKAUTPillbXntOw1zDXU2c3V22TcZMJdnloiEhIH2ndm7v7vFN44MwOjcUu0Z3C5tsQidLU3+GtnNLO3YcXOa3Sf+pqhYoqo7V2aMhr7TBmTZ6xdv2HqwB+sQ/1yhTHbZH5OC1nxOeJuZTRZ7O1t848iA4XTttGMHo134xd58qUImV9ix4+Y0u07825ubSMQjvmlAo3WGQawhp1+8pdHZJbpbozsKC8DK3+mKTzo7wzPcuTgMdjQztVAgX9J/FXe1qhiby++4YwfD3qtzfmnL9TluTrPrxB+MBuQXMbwyu0w0EqJ3k317N6K/PU5IYNRHD8xQHWI4VBvp+MPeK7PLdYn/HjOEMuaD+zu5UKBYqTK0w7APwGBHi2/CepbjFsT8NWKo00/egxEW2O6CJ4umcIj+ZNw3YZ960lph9UhHf3sXCkaRwXrEYSWMqX97rncUa/xuM8ulCrNL+u9jcKXm+Qfirw1DPoqTXplZrjtmaGTA6G9vtaq4OldfTDjZHCERi/jC3ivmRG09pX79lMBQ73zO6t/1hb0zS3WN2p1mV4r/YIcxKTi/rL/3MDq7xN46HhZYyYDRnUyuQKmi6ursRITBzmZfeP6XbRB/K6znD8/fFP+6wj7+SWC4MrvEUMfOR+1OszvF3ycZMAuFMrNLpbo9/8HOZiayee2rXa6EBert7Pwx0qk3rRVWwnp+6Nyvzi3T2dK05d3o1mPIV2Gu+uavnGZ3ir9Phsq1CaM6swWGOpupVBXjmuf617Ij6g5zGRP6uk8Kjs8tE28K0dHSVNf37PFJBszobP0hzPbmJlqjYV907ldm6h+1O8muFH+/eA9Wjr8dMX/wg712ef7NLPggrDc2v8ye9mZE6gsLDPokgeHq7M4XeFlYYT3d7V00R+31zG84za4U/67WKPGmkPae//i8cX17OuorCrUSJ9Xb3tHZZXraYtva4Wk9hnyytmFsLs9AnfcWjPs7MW8UiNMVpeqfzLfwQ8KGNcqut7Nzkl0p/iLii1z/8fk8TWGhp7W+bIGBjjgi+s9xjM0vM2iDGA75ZE5nfH55x2UdVjPY2Uy5qkhn9Q3rTS8WyZeqtnjCfljoZZVT6d/mhkRusivFH0zvwQcNqC8ZrztbIBYJ05eI+6Kzs0MMrYVPOs9xlCpVMrlC7VrroTaHpXF7vmpDpo/FYGcz88slrev6j5mjdjvas1PsWvH3Q9zQ8Azt8RwGOuLaF3ebmM/Tb4O9nS1NRCMhre1NZ/MoBXtssHfIB7nv1rNml+cPettrtb2+dj1z/MEm8ReRB0TktIicE5EPXee4d4uIEpHjdpy3HgY7mplZLGpdE8UQQ3s8h4H2eG0OQUeyecOTq3d+A4ywnmGvvuJvXduADZ7wHh94/pa9e2wKc4HexRnH5/P0tEWJReqbv3KSusVfRMLAx4G3A8eA94nIsXWOSwD/AXiy3nPagRWL09U7VEqZYRB7PP/+ZDPj83lt0x9rMVKbOrv+pN6dnVWLxw7PvyVqlDrWWfwn5peJRepPawVqz8TEvL57F0/ML9syinUSOzz/u4BzSqkLSqki8CngwXWO+1XgNwEt1Na6MROaTpLNLpUolKu2TRgNtMdZKlbIaRonrXnCNj0wezqad43nD0anqasjA9QcmXrTWgF622KExBBYXRmfz9Of1DfeD/aI/yBwZdXrUfO9GiJyJ7BXKfWYDeezhZr4a/rAjNcmjGzy/DW313qQ7ers+tvjpLN5qpqmP47PLZOIR2irY7XravqTMW3vLRhzHHZ5wpFwiFRC77DeRNa+UbtT2CH+63XltSdORELA7wK/sOkXiTwsIidE5MTk5KQNl7Yxlsjo2oAmbPYMrYaoq71jc3lEoM/GkU6popjWdK/msfm8LfFvC6uz0xW7Mrks+trj2o7al4sV5pZKuyLsMwrsXfV6CBhb9ToB3AJ8TURGgLuBR9eb9FVKfUIpdVwpdby3t9eGS9uY1liEZDyi7dDR7jDIiuevp70T83l62mJEI/YkoOk+pzM2t2zLAi+LvmSc6cUihbJ+CQxVcw2CnWI4kNQ3e83uUbtT2PGkPQ0cEZEDIhIF3gs8an2olJpXSvUopYaVUsPA94B3KqVO2HDuuujX2HuYmM8TDgk9NpWDTSWMhV66ev7jNg+TLS9zTNPOzm5P2PrbZbL6TYJOLxYpVZStC5762/UV/5XkhQYXf6VUGfgg8EXgFPBppdRJEfmoiLyz3u93Ep0nycbn8/QlYjve23Ut0UiInjZ948IT88u2i4PxvfrZmy9VmFks2pLpY2GFy3QM/Tghhv3tcXKFspYLvVZG7XpP+Noy26SU+gLwhTXvfWSDY99kxzntYCAZ5+XxrNeXsS4TWftTxXTOfR+fz3PvoR7bvq+7NUo0HNLSXrszfUDv7DXrmuwd2a107odTbbZ9rx04Ya8T7NoVvmA8MJMLBUoa1rm3OywARhxcR094oVAmly/b2tmFQkJfe0zLOY4Jm+dzQO85jlom1y6xd3ze2Leg3gKFTrPrxV8pYwcpnVBK2VbqYDX9mq7ynXBogmwgqWeufyZnXFMqYd/S//bmJmKRkJZhn/H5PJFQ/QUKV9Nfy17TsT3btzLfSXa9+IN+3kM2X2apWLFdDPvb42TzZRY1i5NaAm13BcR+TcNck6azkbLRXhExExj0cmTAvgKFq9F5jsPOlflOsqvFf0BT8XcqW2BA07hwre6LzbXPB8yMEN1KWqSzeWKREMm4PQu8LPqScdKatWVwRgzjTWG6WqNadu7jDozanWBXi//KQi+9ho5O5Qlby8117exSSXsrIPa3xylWqtot9MrkCqSSMVtKHaymPxlnPKtXWwZ7V/eupk/DOSwrk2tA4zr+Frta/Nubm4g36Vf61+4iZxa6rvLN5PJ0tjTZXgHR6tx1y33PZAv0JewXB2OVb0GrkY7dBQpXM6DhOh0rDBV4/ppjlP5t1q4Bjc8bpQ7snBCElQapW5w0nS2QckAMrZh6OqeZvbm87aMcMDzhYrnK7JI+exdnl8sslyq2le1YjY4LvfyS4w+7XPwB+jQsiDUxn6e3LUZT2N7bE28Kk4xHyGgm/lYYxG76zO/Uzd5Jhzo7HeewrI7XEfE3S1rotCfHSk2uwPPXnr5kXLtUz3GHYqRgeMNpzcIgk9k8vTaPcoDad+pk71KxTK5QdszzB71Gdta1OCX+sJI9pQMTDtprN7te/FOJmLmlnj5x0rSZGucEfclYLc9cB6pVxeRCwRF7Y5EwnS1NWtlrzT844fnruMp3xV4HOndrZKfR/U1n87RGw7aV6naSXS/+fck4hXKVbF6f3PdMLu/IwwKG6OjkCc8uGUW/nLK3T7ORjjXKdMLeVCKGiF5hn5q9Dox0rL+hThP6RghTf68fAvGvhQZ0iQtbE3ZOeIZgPISTOX0yQlbE0Lkwly73Fla8VCdGOk3hEN2tMa3CPplcnrZYhJao/Z6w9TfUKWw7mS04EsJ0gl0v/ro1oMkF5zwlMES2WKkyp0lGiJOeIVhhPT3uLazMPzg10ulvj+kV9nFoMh+gqyVKJCTadXZO3Vu72fXiXxs6ahI3tLxU58Iglr16CKJlrxN572DYO7lQ0GY7x0wuTzRsz0bm66Fb8b5M1jkxDJn7XejSlsHs7Bxqy3YTiH8tQ0KPBuR4GCShV0aI055/XzJOparPdo5WWMDu1b0WxhyHHvcWnBdDI4FBj2d3oWDU5OpzqC3bza4X/7ZYhNZoWJtJI+fFUD/PPxGPOFb+VrfOzqkFXhb9yTizSyUtct+VUmSyBUfDIL0JfeZ0aqP23ST+IvKAiJwWkXMi8qF1Pv95EXlJRJ4XkSdEZL8d57WLVDKuzSrQyWyekBibkTiBbmJoeIbOPSwpzdIBnRZDy96pBe8794WCsbrXSTHUyfNPO5jG6wR1i7+IhIGPA28HjgHvE5Fjaw57BjiulLoN+Czwm/We105SiRiTGnn+3W0xIjav7rVojoZJxCPaLIxJZ51b0wCrJvQ1ur9O2msJjw6CaF2D0/bOLBYplr3fkMmJfRqcxA6FuQs4p5S6oJQqAp8CHlx9gFLqq0qpJfPl94AhG85rG6lkXBvPMO3gBJlFKqHPQi+nPf/eNn1W+eZLFeaXSw6HQfTJfbdGl06mPlqjikkNRjqTDs/X2Y0d4j8IXFn1etR8byPeD/zTeh+IyMMickJETkxOTtpwaVvDSgfUIffdaTEEfRY+KaUcXxQTjYTobo1qEdZzQxxqYrhL7NWpflMmVyAaCZFs1n91L9gj/uulLayroiLyr4HjwG+t97lS6hNKqeNKqeO9vb02XNrW6EvGWC5VWNBghys3UsV08fyzy2WK5arjnV1vIqaJODg/IdjdGiMkmoR9ss4mL4BmYS5z1O5UJpfd2CH+o8DeVa+HgLG1B4nI/cCHgXcqpby/U6tYmQT19rIqVcX0gnOLYiwsz9/rkU66JobOdna6FO9zY0IwHBK622JahH0yuTzxphAJB+vcpDRaoe/0fI7d2CH+TwNHROSAiESB9wKPrj5ARO4E/hBD+DM2nNNWdMkImV4oUFXOTxj1JmIUy1Wyy96OdJws+rWavqQeJQ/cSgXUZWSXzhpi6KQn3N2m0UjHhZCtndQt/kqpMvBB4IvAKeDTSqmTIvJREXmnedhvAW3AZ0TkWRF5dIOv8wTLE/M6A8ZqwL1Oh3002eTEreyIvmScyVyBiserfDO5ApGQ0NXiTBqvhSH+Ooih88kLYXOVrw6duxvJGnZiy3hMKfUF4Atr3vvIqp/vt+M8TmF5Yl43IDdiwgB9qzJCjvYlHD3X9VhZ0Ob8HEdVwfSit0vv0+bq3lDI2ZhwKhHnxbGso+fYCplcgZv6k46fJ6VBrn++VCGXL/umoicEK3wBSMQiNDd5v8rXOr/TccOVkhbednZu1T5PaZLr71bRr1QyxvSC9yMdtypc9iXb+sKeAAAgAElEQVTi3t/brDVq94/nH4g/xl6+OngPtbBPm/Mx4dXn8wq3ap/rssPVpEv21kY6Hua+O7lj2VpSGmxQ5LcFXhCIfw1rRy8vSWfzdLY0EY04e1taYxHaYhHP7XXNM0zqsdDLrZhwrwbpj07uWLaWVMLYy7dc8W6Vr9MFGZ0gEH+TlDkp6CVuloO1NnXxkrRLYZCeNmOHKy87O6c36VmNDtlrK6Ud3PH8lYKpBe8qt9ZKk/ukqBsE4l9DB8/fyY0v1uK1vVbFRzfyopvCIbpaop56wlb5AVfEMGGt8vXQ86+FQdzx/MHbzt3K5Op0OJPLTgLxN+lLxlksVlj0cJXvZDbv2oSR1wufahUfXbLXGNl5Jw5pF8v96lDfx601HKBHmfJMzp1MLjsJxN/E60lQpRSTCy6GfUzP36tVvk7vW7AWr3Pf3YyBxyJhOlqaPBdDJ3csW81KiQdvO3c/TfZCIP41vB46zi6VKFWUaw2oLxmnUK6SzXsz0ql5wi52dl56wpMureGw8HqVb8YcxbpR56anLWrO6Xh5fwuOL860m0D8TbweOloPqlu1QVZCA94IxKSLE4JgTnB7mPuezhbMTXrcsbfX65GOi/NXkXCI7taYp2E9N+21i0D8TWpDR4/E0I0KiKupbXLiVWdXWxTjludv7OU749Fevplcnp62GGGXYsIpjxc+ubWgzcLLkV2xXGVmsRiEffxKsjlCLBLy0PN3b4Js9Xm8CnNlcnlikRDJuDu1z1fmdLyy192Kj6mEkcrr5ZyOmznvqWTMs1pVU7VMriDs40tqq3w9FENwMQbuseefzhrDZLdqn2thr4ueYW8iRrFSZX655No5LfKlCnNLzu5YthYvSzy47bjZRSD+q0glvNvhKpMtGDWGomFXztcWi9AS9a6eUSaXp89Nz9DKfffI3slc3tWYsJed3aTLmVzWuaY8mtPJuJy8YBeB+K/CywyJTC5Pr8sTRn0e7l3s9gRZr4dhrlKlyvRi0d0wiIe5/m5Va12Nl/WM0h50dnYQiP8qvFz4lHE5LADeZoRMZt2NCcebwrQ3e5P7PrVQQCmXPWEP5zgmPShy5mU9o8lsHhHobvXP6l6wSfxF5AEROS0i50TkQ+t8HhORvzE/f1JEhu04r930JmLk8mXypYrr53Z7ggxWJgXdxqr46Hb5W69Gdm4u8LLwMuzjRZGzPg/rGWVyBbpbY0TC/vKl675aEQkDHwfeDhwD3icix9Yc9n5gVil1GPhd4DfqPa8TeDVUVkq5nhoHVjqgd2LodnaEVyM7N4ucWXg5p5PO5o29hF30hL3cs8Fv2zda2NFV3QWcU0pdUEoVgU8BD6455kHgEfPnzwJvFQ23uPdqe8NcoUy+VHVdDFPJGIvFCgsu1zPyKjvCq1xwt1czW6QSsVpBOTfJZAv0tEVdrXNj7YHhRcJGJpf3VTVPCzvEfxC4sur1qPneuseYe/7OA902nNtWvPL83V7gZbFir7udnZtFzlbTm/Qm9z2TKyBilCFwE89Gdi6vaQCIRkJ0tjR5FtbzW6YP2CP+63Xva5+urRyDiDwsIidE5MTk5KQNl7Y9vJoks87nfgzcm7iwVxtfpBJxipUqc0vu5r5P5vKexIStzs5tvAqDpBLuh/UqVcXUgv9KO4A94j8K7F31eggY2+gYEYkA7cDM2i9SSn1CKXVcKXW8t7fXhkvbHp0tUSIhcb0BTXolhh7VM8rk8jSFhU4XKj6uxqvKrW4v8LLwqpLpZC7vSZEzL7ZinV4sUFX+W+AF9oj/08ARETkgIlHgvcCja455FHjI/PndwFeUV+vOr0MoJJ7Ehb0Kg3gV9rHSPN2e9vFyZOeFZ9ibiLFQKLNUdG9Op1SpMrXgTZ2bVCLOpMtt2e0aVXZSt/ibMfwPAl8ETgGfVkqdFJGPisg7zcP+GOgWkXPAzwPXpIPqQq8HC58y2QLxphCJmDt1bizam439gt0ODaRz7m1as5o+jzJCMtmCq6uZLVaKFbpnr1XnxovOzqrcWnVxlW/G5VLddmKL2iilvgB8Yc17H1n1cx54jx3ncppUIsaVmSVXz2nl+LvtCYsIvW3uD5Uz2QIHelpdPSd4E+byMia8Osw17NLfu5bG60lnF6NUUcwuFeluc+fv7eaOZXbjr1UJLuDF3rZe5Phb9CXdX/jkRTYIQEs0Qlss4ur9nV7wLibsxUbubu/QthovEhisc3kxkq2XQPzXkErEmV0qUSxXXTunlxtBuF33PV+qML/sbsXH1bi9qtmLOjcWXoR9vFrTAN6M7DK5PB0tTcQi7hRktJNA/NdgNSA3F8e4XedmNW5nSHhR8XE1vS6XeFgRQ/ft7WxpoinsbvaaV2sawJsEBq/mc+wgEP81uN2ArDo33nn+MeaXS67VM3J734K1pFwu8eCl578yp+OeGHq1pgG8Cfukfbh9o0Ug/mtwe3tDLyfIYOWBcSsU4tVqZos+M5XXrUzjWiqgSxOQa+lNxl0Nc3m1pgGgORomEY+4au9kNu/L1b0QiP81uL0QyMsJMqC2h4Bb3qEVBvFqy7tUMsZyqULOpXpG6VyertYo0Yg3j5rb61a8WtNg4WbCRrWqfLlxu0Ug/mvobosRElxbLOK5GLpczyiTKxAJCV0t3tQ+d3sS1It9GlbjdnE3r2PgbpZ4mF0qUq4q+nyY6QOB+F9DOCR0t8Vcqw7o9f6fbsdJ09kCvYmYqxUfV+P2Kl9j+0ZvxXBmsehK9poOdW5SLqYuezmfYweB+K+Dm5t+ZLJ5opEQ7c3u1rmx6G6NEg6Jiw+Md2saYFU2l4udnQ72Trng/Xu5psEi5eKcjpeZXHYQiP86uFkQy6qA6NX2BqGQ0NMWdTcM4qUn7GKJh6rlCXsshuDOyE4HT7gvGadQrpLNOz+ns7JJT+D5Nwxuxg3TWW89YXDXXq89/0QsQrwp5MpIZ8aKCXsshuDOxvUrabze3V9rpe2kC/fXSgf34+peCMR/XVLJGNMLBSouFIjyqtTBatwa6RTKFWaXSp7aKyKkEnFX5nR0CAvUVr26IP7prPeevzWH5cb9zeQKtDc3EW/y3+peCMR/XVLJOFVlxDCdJqOD55+Mu+IpTXo8uW3h1pyODmGQ7tYY4ZC4I4Yer2kAd+sZeZ3JVS+B+K+DW3HSfKlCNl/2PFsglYgxvVikXHE2I8QSIM9HOi6VtMjU0ni9E4hwyFjl61bYx8s1DeBu6nLa4zUN9RKI/zpYDcjpB0aXcrCpZAylYGqh6Oh5rNGF1zFSY9MPN8I+elR87EvGSLvQ2Xmd2QTQFovQ3BR2qXP3b10fCMR/XVIulXhI1zaC8Nrzd2dSUCfPP+fCDlfprOEJe13xMZV0ZyN3r9c0gDmn48LITinFZK5QWyHvR+oSfxHpEpHHReSs+X/nOsfcISLfFZGTIvK8iPxEPed0Aytm6fTQsVbXx+MG5FaYK5PLG4voWr1Z3Wvh1ipfrzYyX0tf0q2wjyb2Jpzv7OaWShQr1V3t+X8IeEIpdQR4gvW3Z1wCfkopdTPwAPB7ItJR53kdJRoJ0dnS5PikkZe1z1fj1iRZOlugt8271b0WrnV2We89YTDEcHapRKHsXOXWatXwhHUQ/14XPH+va3LZQb3i/yDwiPnzI8C71h6glDqjlDpr/jwGZIDeOs/rOG7kvmdyBZrCQmeLN6t7LXraYoi45Alr8LC42dnpUPdlJd3Tufurw5oGC2OV7+5w3OqhXvHvU0qNA5j/p653sIjcBUSB83We13HciBsaC57c37t3LU3hEF0tUXc8YQ0eFjfCPpWqYnLB+zUc4M4cli7JC8Y1xFksVlh0sHLryupe7+3dKZtu4C4iXwb61/now9s5kYgMAH8BPKSUWjenUEQeBh4G2Ldv33a+3nZSiTjnM1OOniNjFjnTgd5EzPFc/0yuwKv2XzMt5Dpu7HA1vWgsEtRBHPpqnZ1z93clecF7e1eH9Q7ENpW4HdEInv+mfxml1P0bfSYiaREZUEqNm+Ke2eC4JPB54FeUUt+7zrk+AXwC4Pjx4+7strEBludfrSrHYtSZXJ4DPa2OfPd2cXqHq2K5ysxiUYsJMmuVr5Nhn4wGq10trA7IyUnfyZrn7729q1c1O/V8TeYKJOIRmqP+XN0L9Yd9HgUeMn9+CPjc2gNEJAr8PfDnSqnP1Hk+10glYpSritkl53Lf0x7u3bsWpzf9sGrK6+AZgrmXr4P2Wh2LDmGfzpYoTWFxNNc/o8kaDnBnNz6va1TZQb3i/zHgbSJyFnib+RoROS4inzSP+XHgDcBPi8iz5r876jyv4zhd5z5fqjC/XNIiLACG+E8tGCMdJ0hrsNp1NU6XeEhrFAMPhax6Rs7aq0udGzcWaerkuO2UugJiSqlp4K3rvH8C+Bnz5/8F/K96zuMFKxkhBW4asP/7V+rc6NGA+pJxylXFzFKRHgdqs2Q0CguAcX+fGplx7PvTmlV8TCWdH+no0NEBtDc3EY2EHN2zIZPL8+p93s9f1UOwwncDVmqEOOM9ZDSaIAPna6LoZ2+cOQdz39PZAj1tUZrCejxifQ57/jpUp7UQMeoZOTVqV0oZnr8m9u4UPVqmhjgd9klr6AnDStaG3WSyBUJiVJnUASv85JR3OJnTI63VwulVvrpVuHRyO8fscpliuaqVvTshEP8NaI6GScQjjomDNaLQyRMGHCt4ls7m6U0Y5YV1wOm67+lsQZv5DTCyjrL5MstF+0c6Ota5cTKBIaNJTa56CcT/OqQSznlLmVyBSEjoavG2zo1Fby032jl7dfKEnd7xKa3JgjYLJzeu17HOTZ+Dqcs6TebXQyD+18HJEg9pc4GX13VuLOJNYZLxiIP25jXzhJ2r71OuVJla0MvzX9nO0X57dVrgZZFKxJhfLpEv2T/S0SmNtx4C8b8OTsYNMxqUv11Lf3uciXln7J3MFejVyDPsbo0Rcqie0fRikarSKyywkvtu//3VLZMLVoUxHejcA89/F9CXNPZ6Vcr+3HfdJsgA+tubmXAgzFWqVJleLGrlCYdDYiz0clAMdfIMV1b52i+GVpsZaNfH3l4Hi/dlcnlao2FaHSod4RaB+F+H/mScYrnK7FLJ9u9Oa5QXbTGQjDPugOev25oGC6fCerotaIOV3HcnUpet0aJOYR+rI5qYd+b+9mnU0e2UQPyvw0oDsveByZcqzC2V2NPRbOv31kt/e5yphQLFsr17+a4UwdJHHMCa0HcuBq6T5y8ijqV7js/n6dZgx7LV9Jt/+/H5Zdu/e3w+r9UoZ6cE4n8drN59ImtvA7I6k36NxAGMzk4p+4fKluAMdOhlbyrpzMKndLaACJ7vWLYWY6GXQ56wZm3ZKDURcmQOa2I+T39SL8dtJwTifx2s3t3uUIj1fbp5D/0OjXRW7NXrgRlojzOzWLQ9IyQ9n6e3LUZEk9W9Fn3JuCOL+HT0hEWEAQfmsCpVRSZX0M7enaBX69SM3jYjIyRtsxhaI4l+zRqQJc52PzAT8/na1pg6Yf397c74GZtfZkCzkB44V99H1xh4f9L+7LWpBWOfBt2e3Z0QiP91iIRD9CZijnn+ujUgJz3/gXbvdyxby8rIzv6w3oBmYRAwPP+FQpkFG3e4ypcqzCwWtbS3v93+BIZxTUO2OyEQ/01wIv1xYj5PMh6hJapXqphxTWHbHxgjRqrfw+LkSEe3jh2c2dTFGknoaG9/uzGnY2eZ8ol5PUftOyEQ/00YcGDoODGf1y7+DUac1ImFXuPZZS1jpP0OzOnk8iVyhbKe9pqTlHbe33GNxXCg3ShTPr1o34ZME5rO1+2EQPw3wQkxnMjq6RmCMZy1MwxSrSrS8wX6Nezs2mIREvEI43P22VsTBw1j/nvMbKsxO+3VcIGXhTXatLWzy+aJhkN0aZbJtRPqEn8R6RKRx0XkrPn/hrsbiEhSRK6KyH+v55xu098eJ2dznFTH7AgLuzu76cUixUpVW3sHbI4L65rJBc6MdKy2oluqJ6yE9ex0ZqyQnm7zVzuhXs//Q8ATSqkjwBPm6434VeDrdZ7Pdexe6FUsG0W/dPX8B9rjpHNGRoMdTGg6uW1h95xOLQyioRjGImF62qK2iuH4fN4cQemVyQXQ127Mcdh7f/Wcv9oJ9Yr/g8Aj5s+PAO9a7yAReTXQB3ypzvO5jt1Dx0wuj1J6eoZgiGGlqphasCcl0BIaXe21u6TF+HweET09YTC84bE5++zVrVrranpaY0RCYutINq1xyHa71Cv+fUqpcQDz/9TaA0QkBPwO8EubfZmIPCwiJ0TkxOTkZJ2XZg+19EebvIcVT1i/mDBQS9mzSxBXYsJ62mt3SYuJ+Tw9bTGiET2n04wwl72ev673NhQS+mxM2FBKaR2y3S6btlAR+bKIvLjOvwe3eI6fBb6glLqy2YFKqU8opY4rpY739vZu8eudpa/m+dvzwFhiqOvQcSXX3x57x+fzNIVFu1IHFnaXtBjTXBz2dDQzbrvnr6+9ds7pzC6VKJarDeP5b5porpS6f6PPRCQtIgNKqXERGQAy6xx2D/B6EflZoA2IisiCUup68wPaEG8K09Uatc8T1jwGbndJi4l5Qxx02bRmLasXtg11ttT9fRPzywx3t9b9PU4xYCYw5PKluuP05UpV+1IHfe1xTo1lbfku3UOY26XesemjwEPmzw8Bn1t7gFLqJ5VS+5RSw8AvAn/uF+G3sHPoOD6fpyVq7JqlI12tUaJh+wpijc/rmeNvYVVWHbPx/upsr5WCakfnPumDUgfWnI4de3KslOrW197tUK/4fwx4m4icBd5mvkZEjovIJ+u9OF0Y7Ihz1abcaGu1q66pYtZCLzs9f13nN8DeMNdCoUwuX9Yyx99iT7t9uf5XZ43vGOzU196BjmaWzRLq9WJNlOs6x7Fd6hJ/pdS0UuqtSqkj5v8z5vsnlFI/s87xf6aU+mA95/SCwY7mWkOvl9G5Za0fFrAv198PE2SJWIRWm0paTPggLGBnrr/lEA1p3NkNmtdmh/M2OrtMU1i025dip+iZkqAZg53N5Apl5pfr9x6uzi7XGqSuDLTHGbdhD4OZxSKFsr4LvGDVSMeGSVA/FP3qS8YRwZZVzaM+8PyHzGsbtcF5uzq3zJ6OZm3nr7ZLIP5bwJoIrNf7z5cqTC0Uag1SV/rb46TnC3UXxLIeODsmUp1kqLPFFs/QCqXotkPbaprCIVKJmC1zHFfnlulsadKuQOFqVsR/qe7vujq7pL3jth0C8d8Cdg0drd/X2VMCYxhfrFSZrHOh14r4a25vZ7Mt4jA6u0w4JFqPdMCIWduR6391Vv8QZntzE63RsG1hn0D8dxlWA79ap0D4yROG+r0l6/d1F4ihzhZml0p112+6MrNEfzKu3Q5ea9nTYU+Y6+qc/mIoIgx21j9nVyhXyOQK2j+720HvVqoJ3a1R4k2h+j1/K0aq+QOzt8u4visz9dk7OrtMe3MTSQ3rvqxmqNa512+v7qMcgD3tzYzNL9eV/qiUMuev9BfDwY7mup9dq7PU3ZHZDoH4bwERYY8NDWh0domIueRcZ6wH2g7P3w9iaFdc2BB//cVwb1cL+VJ9Yb3ZpRLLpYovxHCos6XuCd9Rnzhu2yEQ/y0y2NFcdwO6OrfMQEecsObZAs1Ro/qjHfb6Q/ytzm7n9hbKFdK5fG3UpDN2jOz8MooFw1ufX64vrHd1znAM/NCet0og/ltkyIa44ejsMkM+GCaDIYhX6vCElVK+8YR72qLEIqG6PP/xOaNaqx/s3ddV/8jOT2JYS9io4/m9OrtMSPQty7ITAvHfIoMdzUwvFlkuVnb8HX7IjrAwMmB2/rDMLpVYKlZ8IQ4iUre9fslsgpUO6vL0zsXfT2GQ2pzOXH329ifjNGk+mb8dGscSh6ll/Oww7l8sV0nn8r4QBzAEYmxuecebutQyfXwgDlB/XNiy1w/3N94UJpWI1TWyuzq3TEs0TEeL3pP5sPLs1nV/fbAyf7sE4r9FrEnQnYr/+PwySvlHDPd2NVOqqFoxq+3il7RWi3pz/a/MLhEOidare1ezt6uFyzN1iL+Z865rjarV9LQa+yvUG/bxS1veKoH4b5HBOtMB/SeG9U2C+iXH36LeXP/R2WX2dOif42+xt7O5vglfH3nCoZAYCRs7dNzKlSoT2bxvHLet4o+WqgH9yTiRkOx4qHzVRzFhMMQBjIVLO+Hq7DLJeIT2Zv3DAlB/rr+fJvPBmPQdn1+mVNn+DmZKKS7PLLHXJ44M1DeHNZHNU6kq33R2WyUQ/y0SDokxVN7hJNmV2SVfZQtY9Wl27vn7a5hcb66/X9Y0WAx1tVBVOyvtPLtUIpcvs7/bP/d3b1fLjh0ZKzzmp85uKwTivw32d7cwMr24o98dmV5iqLPFN9kC1qTgTsXw8ozPxNB8sHciEIVyhXTWX0v/99bs3b74XzKfAZ13LFvLcHcLM4vFHVXmvWQ6fMM9/rm/W8EfSqQJw92tXJpe2tGy+EvTiwz3+OdhAdNb2oH4V6uKSzNLvrK3py1KSzTMyA5GdtboyA8LvCz2mV77Tu6vH8Vwv9lRXdqB8zYytUg0HGqYTVws6hJ/EekSkcdF5Kz5f+cGx+0TkS+JyCkReUlEhus5r1fs725hoVBmerG4rd9TSnFxapFhHw2TwbD30g7EcGx+mWK56ivPUEQY7m7d0chuZMr0hH3U2Rk567KjjJ9L00uI+Cd5AVZGKTvp3EemF9nb1az9yvztUq/n/yHgCaXUEeAJ8/V6/DnwW0qpm4C7WH+jd+0ZrnkP22tAKzFS/4gDwIHuVsbn8ywVt5cBMzLlP88Q4EBPa03It8NF83cO+kj8wyGjXtVOwlyXphcZSMaJN4UduDJnsFY1X9rB/b00vcQBH93brVKv+D8IPGL+/AjwrrUHiMgxIKKUehxAKbWglKq/eLoHWBNc2x06WuJwwG9i2Gt6S1Pbu10Xpy17/fXADPe0cGV2+xkwF6cW6WhpoqMl6tCVOcP+nY50phdrYSO/0BwNM9Ae37bnr5RiZHrRd47bVqhX/PuUUuMA5v+pdY45CsyJyN+JyDMi8lsisq7LICIPi8gJETkxOTlZ56XZz1BnCyHZ/tDR6iz81oAs8b64TW/p0tQi8aYQfQl/ZDZZDHe3UqmqbXvDF6cWfdfRgTFSuTi5uO05rMszS74K6VkYYcztteVMrkC+VPVdyHYrbCr+IvJlEXlxnX8PbvEcEeD1wC8CrwEOAj+93oFKqU8opY4rpY739vZu8evdIxoJMdjZvO3QwMi0kebpt1Qx6wG/OLWwrd8bmV5kuLvVd3udWgK+XW/44tQiB3wohod6W1ksGplKWyWXLzG1UPSd5w+Yczrb69itZ91vjttW2HTzTaXU/Rt9JiJpERlQSo2LyADrx/JHgWeUUhfM3/kH4G7gj3d4zZ5yoKeNC9sVw6lF9nQ0E434K7mqNRahPxnnwjY7u4tTixxJJRy6KudYGelsXSCWixXG5/O+9PwP9bYBcH5yYcvrT85ljLZ/2PxdPzHc08rUQoFsvrTlDYasUa8fRzqbUa8aPQo8ZP78EPC5dY55GugUEcuVfwvwUp3n9YwjqTbOZRa2tbn5ucwCh1P+e1jAiINvZ6RTrlS5MrPsq8wXi67WKIl4ZFsjHWuU4Ed7D5oCfmFy6/Za4n+kz3+du9VhWTZshbOZBeJNoYZb3Qv1i//HgLeJyFngbeZrROS4iHwSQClVwQj5PCEiLwAC/FGd5/WMw6k28qXqlgu8VaqK85MLHPGp+Bsjna3HhS/NLFGsVH3Z2YkIh3rbOJ/Zeme3MpnvP/HvS8ZojYY5P7l1e89lFoiGQ7XyH37iSJ8p/untif+h3raGS/OELYR9rodSahp46zrvnwB+ZtXrx4Hb6jmXLlgifi6zwN6uzeOeV2aWKJSrvvSUwOjs5paMOG9vIrbp8WcmcgDc4FN7b+hL8OVT6S0ff3oiR0jwbWd3sLeN89v0/A/0tPqmgN1qhjpbiEVCnM3ktvw759I5Xnuw28Gr8g7/3UGPsR7yrTagM2njOL96/jf2GyJ+emKr9i4gPhVDMLzD6cUiU1vc3/b0RI7h7lZf5byv5mBvKxe24/lPLnC4z5/3NhwyRnZntxj2yeVLjM3nfduWNyMQ/23S0RKlpy3G2S0OHc/6OEYKcIMp/i9PZLd0/Jl0jr2dLTRH/SmGlr1Wp70ZZ9I5jvr03gIc7UtwdW6ZbH7zmjf5UoXLM0u+nOy1ONLXtuVntza/EYh/gMWRVBvntjhUPpvOsac9TlusrgibZ/S0xehuje4aMbTCVWe2MNLJlyqMTC9ytN+/9t40YHbu45vbe2FyEaVWYud+5EiqjatzyyxuYd8GvztumxGI/w442tfGmYncljJ+zqQXOOzzxnNDf2JLYZ9iucrFqUWO+lgcehMxOlqaOLOF0MC5zAJVtRIa8yPHBtoBODW++cjuJfMYv87nABxObX1kdzadIxrx5+T2VgjEfwfcPNjOYrGyaf57oVzhbCbHTT4WBzDE/0x68/TWC1MLlKvK156/iHA0tbXOzjrGz/b2JWN0tjRtSfxfvDpPc1O4liLqR27ekwTgxbGt2Jvlpv6ELye3t0JjWuUwtw4a3tLJsfnrHndmYoFSRXHbUIcbl+UYN/QlWDbjvdfjhVHj73HLYNKNy3KMY3uSnBrPbrp5/anxLNFIyNdL/0WEmwaSWxL/k2PzHNuT9HXa41BnM50tTbw4ev1nt1pVvDg2zy3ms96IBOK/A46k2ohFQjWx24jnr84BcNuQvxvQreb1Pzc6d93jnhudoy0W4WCPfz1DgNv3trNUrGy6GOjZK3Pcsifpe8/wpoEkp9O563Z21ari5Fi25vj4FRHhlsF2Xrh6/Wf30swSuXzZ98/u9fB3q/WISDjEjQPJTRvQC6PzdLQ0+WpHq/W4oS9Bc1OYZy5fX/yfHyZsqKgAAAt6SURBVJ3n1sF239X0Wcvt5kjtuSsb21uqVHlxbJ479q67hYWvuHlPknypet305YvTiywVK7WwiZ+5dbCdM+kc+VJlw2OsZzvw/AOu4dbBJC+NZa8bB7fEUMTfYhgJh7h1qJ1nriOGhXKFU+NZbt/r7xAXGHVckvHIde09PZEjX6pyxz7/2/vq/UYHdmJkdsNjXmwgMbxtqJ1yVfHydeZ1XhidIxoJ+Xo+ZzMC8d8htw11kCuUN1wwslyscCada5hh4537Ojg1lqVQXt9bOjWeo1RR3N4A9oZCwm1DHdf1/J81P7uzATq7fV0t9CZifP/SxuL/9MgMLdFwQyx4utUc2T17eWN7nxud56aBpG/23N4JjWuZw9xjLvn+7vmpdT8/cWmGclVxfLjLzctyjDv3dlKsVDm5QZaEJRyN4PmDEfc/nc5tmA/+7JU5ulujvg/pgREHf81wJ0+PzGx4zHfPT/Oa4a6GEMPBjmaGOpv5zvnpdT9fKpZ59vIcdx9ojGd3I/x/Jz1ib1cLe7s2bkDfOjdFU1h4bYM0ICs08N0N7P3GmUkO9rSyp8P/Yghw76EeKlW1rr1KKb59borXDHf5PqRncXx/F6Ozy4zPX1uwMJPNc35ykXsONU6Nm/sO9fC9C9PrTnI/PTJLsVLl3sM9HlyZewTiXwf3Hty4AX3r7BSv2tdJS9SfK3vX0puIcctgkq++fO2WDflShScvTvOGo/ptwLNTjg930hIN8/Uz1+4odzqdY3w+z5tvbBx77zKdlO+cu7az++4F4717GqjA2b2Hu8nmy+uma3/n3BTRcIjXDPt/Mv96BOJfBxs1oJnFIifHsryuwTyHt9yQ4geXZ5ldLL7i/adHZsiXqrzhaOPYG4uEufdQN187k7mmnPXXThsdwptuWG/XUn9ybCDJQHucfz45cc1n3zk3TSIWaYhMHwtrFPOtc9eGbb91boo793U0jOO2EYH418Hrj/QSDgn//OIrH5gnzJLAr28gTxjgzTemqCr4xtlXesNffXmSaDjE3Q3kGQK88WgvV2aWr1nJ/ZWXMxwbSNKX9NcexdcjFBL+xc39fOPM5CvmOcqVKo+fSvOmG1O+X8+wmlQizs17knzp5CvLd4/OLnFyLNtQo9iNqOtuikiXiDwuImfN/9cdJ4nIb4rISRE5JSK/Lw0SKO1qjXLvoW4efW7sFaGfzz07xr6ulobIfFnN7UMdpBIx/vG5sdp75UqVR58b40039Dacp/S2Y/2EBP7hmau1967OLfP0yAz3H+vz8Mqc4e239FMoV/nKqtDeN89NMbNY5IdvHfDwypzhR2/fw7NX5mob8oDx7AL86G17vLos16i3K/8Q8IRS6gjwhPn6FYjIvcB9GJu53IKxifsb6zyvNrz3NfsYnV2ueftn0jm+dW6K97x6qGEmAy1CIeE9x4f4yssZLpnbF37+hXGmFgq85/hej6/Ofvrb47zhaC9//dSV2oKgv/juJQB+/PiQl5fmCMeHu9jTHq/ZCPBn3x6hpy3GW25snBCXxY/dOUhTWHjkOyOAsXDvL793ibsPdvlyg/rtUq/4Pwg8Yv78CPCudY5RQByIAjGgCdj6Vkma80M397G/u4WP/fPLzC+X+NXHXqI1GuZf373f60tzhIfuGSYWCfPRf3yJ+aUSv/2l0xzta+OtDSgOAP/+jYeYWijw/33lLOcyCzzynRF+9LY9DHU2njiEQ8IH3nCQp0ZmeOz5Mb78Upqvn5nk/a87QDTSOCEfi75knB+7c5C/fPISL09k+Z9fO8/YfJ7/442HvL40V5Ct7s267i+LzCmlOla9nlVKXRP6EZHfxtjWUYD/rpT68Gbfffz4cXXixIkdX5ubfOvsFA/96VNUlUIp+PUfu4WffG1jij/An3zrIh997CVCYuSI/83DdzfMeob1+MXPPMdnvz9KSIzNfB77udc1TErrWkqVKu/+g+/wnFm36qaBJH//s/f6dqeyzZjMFXjH73+T6YUCVQU/fNsAH/9Xr/L6supCRL6vlDq+6XGbib+IfBnoX+ejDwOPbCb+InIY+G/AT5hvPQ78Z6XUN9Y518PAwwD79u179aVLl9Yeoi3fOTfFo8+Ncd/hHn709saOFyql+PtnrvLUxRn+5auHeE0DCz8Ygvin377IlZll/u19w74uabwVsvkSf/C181SV4uHXH6S7bfO9m/3MpelF/vTbI/QmYnzg9Qd9P8qxTfw3Oclp4E1KqXERGQC+ppS6Yc0xvwTElVK/ar7+CJBXSv3m9b7bT55/QEBAgC5sVfzr7eIeBR4yf34I+Nw6x1wG3igiERFpwpjsPVXneQMCAgIC6qBe8f8Y8DYROQu8zXyNiBwXkU+ax3wWOA+8ADwHPKeU+sc6zxsQEBAQUAd1JWYrpaaBt67z/gmMCV6UUhXg39VznoCAgIAAe/H3zEZAQEBAwI4IxD8gICBgFxKIf0BAQMAuJBD/gICAgF1IIP4BAQEBu5C6Fnk5iYhMAvUs8e0B1t9jsfHYTbZCYG8js5tsBWfs3a+U2rQmtbbiXy8icmIrq9wagd1kKwT2NjK7yVbw1t4g7BMQEBCwCwnEPyAgIGAX0sji/wmvL8BFdpOtENjbyOwmW8FDexs25h8QEBAQsDGN7PkHBAQEBGxAw4m/iDwgIqdF5JyIXLOnsB8RkT8RkYyIvLjqvS4ReVxEzpr/d5rvi4j8vmn/8yLiq22JRGSviHxVRE6JyEkR+Y/m+41qb1xEnhKR50x7/6v5/gERedK0929EJGq+HzNfnzM/H/by+neCiIRF5BkRecx83ci2jojICyLyrIicMN/Toi03lPiLSBj4OPB24BjwPhE55u1V2cKfAQ+see9DwBNKqSPAE+ZrMGw/Yv57GPgDl67RLsrALyilbgLuBv5P8x42qr0F4C1KqduBO4AHRORu4DeA3zXtnQXebx7/fmBWKXUY+F3zOL/xH3nlnh6NbCvAm5VSd6xK6dSjLSulGuYfcA/wxVWvfxn4Za+vyybbhoEXV70+DQyYPw8Ap82f/xB433rH+fEfxgZBb9sN9gItwA+A12Is/ImY79faNfBF4B7z54h5nHh97duwcQhD8N4CPIaxr3dD2mpe9wjQs+Y9LdpyQ3n+wCBwZdXrUfO9RqRPKTUOYP6fMt9vmL+BOcy/E3iSBrbXDIM8C2Qw9rg+D8wppcrmIattqtlrfj4PdLt7xXXxe8D/DVTN1900rq0ACviSiHzf3KMcNGnLdW3moiGyznu7LZ2pIf4GItIG/C3wn5RSWZH1zDIOXec9X9mrjA2P7hCRDuDvgZvWO8z837f2isiPABml1PdF5E3W2+sc6ntbV3GfUmpMRFLA4yLy8nWOddXeRvP8R4G9q14PAWMeXYvTpEVkAMD8P2O+7/u/gbnX898Cf6mU+jvz7Ya110IpNQd8DWOuo0NELOdstU01e83P24EZd690x9wHvFNERoBPYYR+fo/GtBUApdSY+X8Go2O/C03acqOJ/9PAETN7IAq8F2OT+UbkUeAh8+eHMGLj1vs/ZWYO3A3MW0NMPyCGi//HwCml1P+76qNGtbfX9PgRkWbgfozJ0K8C7zYPW2uv9Xd4N/AVZQaIdUcp9ctKqSGl1DDGs/kVpdRP0oC2AohIq4gkrJ+BHwJeRJe27PWEiAMTLO8AzmDETT/s9fXYZNNfA+NACcM7eD9G7PMJ4Kz5f5d5rGBkPJ0HXgCOe33927T1dRhD3eeBZ81/72hge28DnjHtfRH4iPn+QeAp4BzwGSBmvh83X58zPz/otQ07tPtNwGONbKtp13Pmv5OWHunSloMVvgEBAQG7kEYL+wQEBAQEbIFA/AMCAgJ2IYH4BwQEBOxCAvEPCAgI2IUE4h8QEBCwCwnEPyAgIGAXEoh/QEBAwC4kEP+AgICAXcj/D8Gd0Dh/a6YyAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "def generate_sinusoid(N, A, f0, fs, phi):\n",
    "    '''\n",
    "    N(int) : number of samples\n",
    "    A(float) : amplitude\n",
    "    f0(float): frequency in Hz\n",
    "    fs(float): sample rate\n",
    "    phi(float): initial phase\n",
    "    \n",
    "    return \n",
    "    x (numpy array): sinusoid signal which lenght is M\n",
    "    '''\n",
    "    \n",
    "    T = 1/fs\n",
    "    n = np.arange(N)    # [0,1,..., N-1]\n",
    "    x = A * np.cos( 2*f0*np.pi*n*T + phi )\n",
    "    \n",
    "    return x\n",
    "\n",
    "N = 511\n",
    "A = 0.8\n",
    "f0 = 440\n",
    "fs = 44100\n",
    "phi = 0\n",
    "\n",
    "x = generate_sinusoid(N, A, f0, fs, phi)\n",
    "\n",
    "plt.plot(x)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7ffa5a71ae48>]"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYYAAAD8CAYAAABzTgP2AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzsvXmUZNdd5/n5xZ77XntJVSWXF3mT5ELsDHgBe2baMmC6bTYxA6NhwDPMMDDYh27gsJwxzZkDhxkaEGAwS3vBDVh9sNvYxgYabypZshbLkkqlKlXWmntmZGbsd/5470ZGZsbylntvRFHve06eyox48eJXd/l9f9v9PVFKkSBBggQJEmik+i1AggQJEiQYLCTEkCBBggQJdiEhhgQJEiRIsAsJMSRIkCBBgl1IiCFBggQJEuxCQgwJEiRIkGAXEmJIkCBBggS7kBBDggQJEiTYhYQYEiRIkCDBLmT6LUAUzM7OqhMnTvRbjAQJEiS4qfDII48sKqXmel13UxLDiRMnOHv2bL/FSJAgQYKbCiJyMch1SSgpQYIECRLsQkIMCRIkSJBgFxJiSJAgQYIEu5AQQ4IECRIk2IWEGBIkSJAgwS4YIQYReZ+I3BCRJzu8LyLy2yJyTkQeF5F7Wt67X0Se83/uNyFPggQJEiSIDlMew58Ab+7y/luA0/7PA8DvAojINPCLwNcD9wK/KCJThmRKkCBBggQRYIQYlFL/CCx3ueQ+4E+Vhy8AkyJyGPgu4JNKqWWl1ArwSboTTCw89JUrfPSxy7ZuHxilap0PfOlFzi8U+y0Kzy8U+cCXXqRcq/dVDqUUH3viKl88v9RXOQBWtyr8xRcvcn291G9ReOzSKn/z6GUajf4+grfeUPzVl+d5Yn6tr3IAXFsr8WdfuMh6qdpvUfin5xb41Fev91sM43B1wO0ocKnl73n/tU6v74OIPIDnbXDbbbdFEuKvvjzPZ59ZYLyQ5TtefiDSPUzg//7Y07z/8xc5PFHgsz/77eQz6b7IsVWp8W9+//MsFitcWNrkPW95RV/kAPi7r17nJ/7iywB8/Ke+lVccHu+bLP/nh7/Cp792gw8/fIm/+clvRkT6IsfVtW3e/rufo9ZQrJeq/PA3nuiLHAB//M8v8Kt/+zT5TIp/+rnv4MBYoS9yKKX40fc/zFNX1vn884v8hx94XV/kAHh8fpUf+qMvAfBH95/hDa842DdZTMNV8rndzlJdXt//olIPKqXOKKXOzM31PNHdFg/+0BkOjRf48y8EOvxnBduVOn/5yDxzY3murpX4+6dv9E2W//LkNRaLFWZHc3zgiy9S76NV+pFH5slnUmTTwge+9GLf5Li8us2nv3aDubE8X5lf44nL/bOQP3J2nlpDMTua7+uaVUrxZ1+4yOxonnKtwUcfvdI3WZ6+usFTV9YZzWf4+JPXWNgo902W//jFF8llUowVMn2dHxtwRQzzwPGWv48BV7q8bgW5TIrvfOVBPn9+iVq9YetruuKRiytsVeq893tezUguzef7GDr57DMLHBjL82//uztZL9V46kp/lOBmucY/PLvAD3z97fw3L53jH55d6IscAJ9/3puP//eddwPwX88t9k2Wf3pukdcem+CBbzvJs9eLfVOCF5e2uLi0xU+98TSn5kb43PP9G5OPP3mVdEr4vR98HUrBP/ZxrXzu+SW+/aVzfO89x/jc80tU+6RTbMAVMTwE/LBfnfQNwJpS6irwCeA7RWTKTzp/p/+aNbzu9im2KnWeub5h82s64rFLKwCcOTHNPbdP8aUXuqVm7OLx+VVed/sU33THDABf6BNJPXVlnUqtwbeenuWe26e4uLTF6lalL7I8cnGZ8UKGe09Mc2p2hC9fXO2LHNV6g8cvr3LP7VO87vbppmz9wFfmvTE4c/sUX39yhrMXVvrmXT764ip3Hh7nG++YYTiX5vH5/szP9fUSLy5vce9Jbx+Xaw2eudYfnWIDpspVPwB8HniZiMyLyI+KyI+LyI/7l3wMOA+cA/4A+AkApdQy8CvAw/7PL/uvWcM9t3lFT1++uGLzazrisUtrnJobYWIoy9edmOaZ6xts9CGJtrZV5cLSFq8+NsGB8QK3TQ/z6Iv92WRP+uGaVx4Z565jkwA83qck5yMXV7jn9ilSKeGu2yZ57NIKSrlXgs9c26BUbXD3bVO86ug4+UyKsxf6s2a/cmmNQjbF6QOj3Htyio1yjeduuFeCSimeuLzGq46Ok04Jrzo6wVf6uE7AM/DuPu6t2Ucv9Wf/2ICR5LNS6p093lfAT3Z4733A+0zIEQTHpoaYHM7ydJ/Y/cnLa3yjb6G//NAYSsG5G0Xuvs1tle6TftjoNUe9Rf3Sg2Ocu9GfKqmnrqwzO5rnwHiBQs5LxD8+v8q3vTRaLikqKrUGzy9s8qY7vSTia49N8ldfvsz19TKHJtwmW59qzs8E+Uya0wdHebZP8/PE5VVeeWSCTDrFyw56RQHPXS/y8kNuCwTmV7ZZ267yyiMTgDc2f/qFi1TrDbJpt2d1n72+gYi3h/OZFNMjOZ4cgIotU7jlTj6LCKdmR/pSKrpZrnFtvcRLDowCNP99fmHTuSzP+aG0lx7yZDh9cJQXFjf7Eif92rV17jziKZnxQpbDEwXO92FMXlzepN5QzXm5Y8779/yi+7VyfnGTXDrF8enhpizP94kYnl/Y5KUHxwA4NTeCCH0xInSoRlesvezQGJVag8sr285leX5hk2NTQxSyaUSEl8x5++dfCm45YgA4NTfaF8WjF87J2REAjk8Pk00Lz/eBpC4sbTGSSzM3mgfg9IFRag3FxSW346KU4uLSFqf8MQFvfF5wLAfsEPSpWY8QTs55MvVjw7+wsMntM8OkU17h3h1zo1xe3Wa74va8ydp2leXNCidnPYIqZNMcnxruCzFcXN4CdvbPCf/fvqyVG8Wm4eDJMsz5hBhubpyaG+HGRtl5bP+Cv4BPzHgLOptOcfvMSH822dImt8+MNGv0tZV87obbxb28WaFYrnGbbxmDt+Ev9GGTaWPhlE8Ih8cLFLKpvhkRJ1vIUsvk2nvR83D7zI4sLzkw2hdj5uLSJmOFDFPDWWBnH110vFYaDcX5xd3EcHJ2lMVieSAO3ZnArUkMvkV4YXHL6ffqTXZitkUJzoxwadmtHOCVILYqnuNTnkzzK25l0Vbg7TM7Y3JyZoSVrarzyqQLi5vMjeUZK3iKJ5USTsyMOPcY6g3Pi9IeC+xYyc7X7NJuLxfgtulhLi1vOU/KX1ja4kSLMTM7mmMkl+bCktsxWSiWKVUbnGhds9p76YMRYQO3JDEcmxoCvMNMLvHi8hYHxvIM53Zy/kcnC85jpPWG4tLKFre1LOzJ4SzDubTzMdGhq1Zi0L+73vBX1rY5Mjm067Xj08POyXJho0yl3miSNcCxSe/3q2uu58f7v7d6dMemhtis1Fnbdmsdv7i0uWvNiojnXToOJV3x98jhiZ21otfsJcdrxRZuSWI46m/+K46V4NW1Eof3VLccnRpio1xz6oIubJSp1lVzHMDbZEcnh5h3TFKXlr3vO9aiBLVyvuZYCV5bK3F4fM/8TA5xddVtzySt/FvXyvhQhpE+EPfVtRIzIzkK2Z22LXrduFwrjYbi8ur2LrLUsvRjHwMcntyZnyM+SbheK7ZwSxLD5HCWQjblfEFdXy/tK3vUStCl13DNbw53aI8SPDY15Nx7uba+X/EcaRK32012bW3//ByeKDgnbt28r1UWEeFIH5Tg9fUSB/eS5ZR7w2p5q0K1rvYZVocnCk1F7Qr6+460eAzjQxmGc2muODZmbOGWJIbmJnM8iVfXSvuUcT+8l2tr+xUPeBvetUV6Y73EgT1jMjWcJZ9JNQnMBTZKVTbKNY5M7lE8k+4twaZFOrE7rOURg3uy3OflTroPxWqyPDie3/X64ckhNko1iuWaM1murm5TyKaY9JPg4OmUwxOFxGO42XFkwu0m2yzX2CjVOLRns/d3k+21voZY266yVXG3ya6vl/dtdr3J+kOWe+fHGyOXRsS1tRK5TKpZfdOUZapPHsMeYpgeyTn3uDuvWe9vl2FHLyQ8tK/r7pHJIec5IFu4dYlh0rHiaYYHdivBmdE86ZRwY91dg7Rr6yWyaWFmJLfr9QNjnmyLG+6qga6v7/eiwCMplyGCK00rfT9ZgnuP4dB4YZ/iOTxeYGmz4uzZGeVanaXNyr75ERHmxvJOm/pd9/fHXmLQsrlcK1fXtvetE/DWzhXHYS1buGWJYW4sz9JmxdkDUK5ri3R8t0WaTnkK2u0mK3FgrEAqtVvxzPnEsFB0s7hr9QaLxfK+UBJ4m+yaw012rU3CF3bG5MaGQ1na5KIAZn1ZlopuiFsbK+2Ie240z0LRoTGzVkJkZz40jvQp1Lc3zAeet7lYLPetc7NJ3LrEMJqn3lCsOKqVv9ohrg8wO5pn0eEm8xKK+X2vN4nBEUktFis01P64sZZloVh2Vit/ZdVTPHsfQJNNeyEdl8TdLq4PNE+pu5JFe7l7Q0mAc4/hxkaJmZH8vp5IO8aMG1lq9QbX10v7clEAc6M5lPIS5Tc7blli0NbXoiPrq1MlEOwoQVe4tra/0kTLAXDD0YZvxo3bPA1sdjRPpdZgw1FS8dqap3hymf1bwiVxK6XaVkeBe+Ju5l06rVnHZLk3DAtei46xfKYPxkz7NQtuQ7G2cOsSg2Pr6/p6iYmhLEO5/Y/x7Ee8tt3CnhnJkxL3Fmn7sImX/1h0tuHLzRzLXsyN5Z0ZECtbVSr1RnuybBozbom7fSip4MlacxM2ub5ebjsmADOjOZY23cyPHnutP1ox69h7sYlblhjmHG+yhY3yvvhoqyyLxbKTfMdm2Svta6eM0ylhZtQdSd3wFc+BNqGkpvXlSCEvbVaYGc21fc+lx7CkFU+btTLry+dsfjbK5DMpxof2d+fXa3lp0x1JtQtpgT8/jsZEE9Bsm7Wy4zEkxACAiLxZRJ4RkXMi8u427/+miDzm/zwrIqst79Vb3nvIhDxB4NpjWNqsMD3SXvHMjeap1pWTFgPdLB4ti7tQUtlPvu+XRb/mTCFvljuOyaxDsmwqnjZrJZ9JM17IOLNIF4vemOytjgK3Ya1avcHSZqWZY9mLfhD3TDuPwScLl/lCW4j9oB4RSQO/A7wJ7xnOD4vIQ0qpr+prlFL/R8v1/ytwd8sttpVSd8WVIyzGCxlymZSzSVzerPDSg6Nt32tNoE11IA9T0Ipnb6mqxoFxt0pwajjXbC3dCh1KWnK24TsT9+xYjq1Kna1KbVefK1tyAEx38F60d+kCy12MmQMOiUEnc9tZ6eCFkr50wZFn6c9PO+9yNJ8h71Cn2IQJj+Fe4JxS6rxSqgJ8ELivy/XvBD5g4HtjQUScltwtFcudFY9DF3RZK54u3osrYljeLDM9km373vRwDhFYcBBK2qrU2KrUu4aSwE1ScdkPzbTzosBtPqobMbj0GFY2PU96usOYzI7mWdmqOCkTXdwsk0unGMvvNxBExPdekuQzwFHgUsvf8/5r+yAitwMngb9vebkgImdF5Asi8jYD8gTG7Kib8wP1hmJ1u9pxYbssuVve7EEMDvMd3RRPJp1iejjnxPrSVuDsAMyPVip7Tz1rOA1rFTvnXWYc5jt0HmOqgxExq8tEHSSg9Zi0C6+BlxtKPAYP7Uaok1Z5B/ARpVTr0c3blFJngO8HfktE7mj7JSIP+ARydmFhIZ7EPlxVm6xsVVCqsyvcLBN1cPq5GUrquOHz1BrKSdO45c1KR8vYkyXnxovqMSZzo+7yHcubFaaGs2Q6PMPYZYXU0ma5Y8gxn0kzMZR1asx0WisuCxWWiuWO6wS8swwuKwxtwQQxzAPHW/4+BlzpcO072BNGUkpd8f89D3yW3fmH1useVEqdUUqdmZsz85B4V9bXUo/wzXghQzYtTg7GLG+WKWRTHWPlOrTjwvpa3qx0tALBXVJRW6TtEopaDnBnHXdaJ1qWYrlGqWq3LcZWpUap2ujo5YK3nl2tE/197eCyjHephzGThJJ28DBwWkROikgOT/nvqy4SkZcBU8DnW16bEpG8//ss8M3AV/d+1hamR3KsblWsn67ViqfTwhYRJodzrLhwhXss7KlhT0bbJ8J7hdfA3SbT39HJOtYWootWFEvFwZifpR5jAl77+tUt+57lUo/wmpbRVdixm8cwO5pnebNM3VGrHVuITQxKqRrwLuATwNPAh5VST4nIL4vIW1sufSfwQbVbC78COCsiXwE+A7y3tZrJNqaGc9Qayvrp2mYMu4NF6smSddKeo1tc35PDVzybdjf8qh9e66Z4ZkbdWKTdKk3Aa4sxls84mZ9u5ylgRznaHpdeVjp4BQKuPIbJLuE1Vx6DUqprWTPA1EiOhsL58+RNw0jtnVLqY8DH9rz2C3v+/qU2n/sc8GoTMkSB7qe+ulllvNA5pBEXQTbZ5HCOFQfWVy9i0O/ZDmvpMelWnjs1nKNYrlGpNdq2qjCFpWKZoWy6aynq5EjWyTOoexK3/55tS725ZruQ1ORwjqevrluVA7y1OD3cWY6xfIZMSqyPyValTqna6GrMaOJe2aoy2UXmQccte/IZHLrlmxVEdr6vHaZdhZKKlZ7hAcC6Eux1ngJ2Ntnqtn1ZuiljTxb7xK2bOnbKdWg5wM2ahe7zMz2SdZMX63LGBHQoNmt9fnrlCsHd/NjGrU0MI5rdbcdryx0PcrXKMggew2jeT4RbDiWtBPCiJvxNtmZ5XFa2ehPD5HDOOlmuBAivtVqkNqEPFvbyckvVBtsVu4nwXmtWy2J7frSB0s3Ac2VY2cYtTQyTjth9ZavSMXHWKovtRPh2pc52td41POAqEb4UgBhcKcG17equxzR2ksW2HFqZdAuv6TW76iDHkEunGG1zkEtDz50L76Vb3gW8+bEdStL377ZWXOXobOOWJgZXk+gpnu4Le9pBIly7/d3itfp925u9mWPoIosrt3xtq8rEUC9isD8mWvF0kyWX8ZS17RDO6laVqZFsx4Nc4GZ+lFKBPLqJIQfzsx2CGBKP4ebFxFAWEftu3+pWlckeiqc1EW5Pjsqu7+oEL6xlnxjG8pmuSWVXbvlqAI9hcjjLRqlmte2CbqIYZK1Yt463KwHI0vfoLK7Z9VKNekN1NSC0LLbHZM1fh+NdxmWskCEl9osDbOOWJoZ0SpgYchEiCGaRgl1LQyuebgtby2K7DHF9u8pEgPAa2N1kjYZidavC5FCPGPaQToTbJO7eFim48V7Wtqs9x8RFKGk94JqdHM5aL1LQ+6fbXk75OsW2LLZxSxMDuNlkQZTglMNN1mvDT43kHFikvclyJJcmmxarxF2s1Gio7psdWstE7c3PagDFAzipwFndqgZQxu6Mmd5j4iXCbZ4IX92qMpxLk8/sf9hWK1xUsNnGLU8Mtt3yWt17PGVgt9ym4tEx7ECJ1orVRnprAYhBJ8JtKuO1gGPiwntZ264iAmM9ztRMOajAWQ8YXgO7oaSgxDDlYH5Wt3uHhEHrlMRjuKlh22NYL3nJ5F4LykUiPMwmayisNtILQgyeLHaJO2hc30WF1NpWhfFCtmtZM7jpURTEo8umU4wV7J4IXwuQ8G1937ZhNdEj1wG+Tkmqkm5u2PYYtOXQyyKdGMqSEvtueSYljLR57vReWQDWt+1VSAUlhknL1SY7cf3eB9zA/vwEGhPLifBqvcFWpR5YFqvhtQCVWloOsD0/lYAeg32PzjZueWKw7TGsBYzrp1LCWCHbzAPYgLYCu5Ugws4mtPWoUaVUKCVoNzwQrFLLRYVUkOoo2PFubM1PUCsdvLWivWKbsgQNJdk8DLm6FWx+XJx5sY1bnhgmh7JsVepULVlfqwGrKrxrMlaf+xxUGWtZbYWSStUGlVqjpxcF9om76TH0GJdRvx+PzQ0fpHoNdrxPWwo5qJWur7G9ZrNpYSjb3cuddBDqC0rcUyM5tqt1663RbeKWJ4amErS0uNcHyPoKUh2l5QD7Fmkgj2Eka7VENGgJr4hXhmjTo1sPStwFu2s2zPyMW/ZyPWOm8xPTNGyH+pRS/kHI3jkG24aVCyTEMOQd+bdtfQWJTY4X7FpfQS3S8QEihvFClkrNXhni2naVQjZFoYdFCt642CTuoBapfeKu7PqeXrLYXLMeWfZuAl3IpslnUna93Hoj4Jr1dYrFHJ1t3PLE4Mo6DhJKsm2RBg0lTVj2okIRg2XrK8jhtqYshYy1MQmTd7E9JmHnx6oxE+AEdqsstpRx0FwU2NcpLnDLE4Ntt3x1q+p3LO091PY9hmBVFSO5NOmUDIjHYNf6CppQBO0x2BmTYtlr/RCEpGwrnqCVWlqWsmWPLjAxWCTuUJ5/EkryICJvFpFnROSciLy7zfs/IiILIvKY//NjLe/dLyLP+T/3m5AnDFx4DEEX9sSwPcXT8Bv0BZFFRLxNZtFKh+ChCrDoMYRRPBat46CHD6HVmLFDlk0vt9A7hNMkboveS9AH3tgk7mjzc/MSQ+wnuIlIGvgd4E3APPCwiDzU5hGdH1JKvWvPZ6eBXwTOAAp4xP/sSly5gsK+Wx7CFS5kKFUblGv1nsfuw2KjVEOpYCEt0LFju4oniHVsO9+xvl3ltunhQNd6iVa7YxJkrRSyKXLplFWSGstnOj5KsxU7xRs1DoyZlyVI59umLAV7Zyp03iXYmrWbt3QBEx7DvcA5pdR5pVQF+CBwX8DPfhfwSaXUsk8GnwTebECmwLBtfQVN+ILdg2U7MdJg1pfNfMd6s/VDEIvUfqgveCjJnhcV9AQ2+B6dRVnWt3v3SdKwSdz1hmK9VAsli/UiklvEYzBBDEeBSy1/z/uv7cX3isjjIvIRETke8rPWYNv6CvIQGA2bmyyMRaplsTkmY/kMqR6tHzw5LFeNbVdCkaWtCqkwoQqwm48KFf606HFvlEKuWZs5hhBl580KqVucGNrt7r3d1/4zcEIp9RrgU8D7Q3zWu1DkARE5KyJnFxYWIgvb5r5Wra+wMWyws8miEIPNqqQwChDsWF+lqvdw9zChCluyhAmvgd35CVo2C67GJDhJrZeqVp6CGPSgnYbNfIcLmCCGeeB4y9/HgCutFyillpRSZf/PPwBeF/SzLfd4UCl1Ril1Zm5uzoDYO7B1SKdZgjgAmyyMKww7m8wGwpClzfp0fc8gSVawS9xhyiG1LLa8qCgew6B4udW6olQ138Vgo1RlrNC7nUxTlkLmlj/H8DBwWkROikgOeAfwUOsFInK45c+3Ak/7v38C+E4RmRKRKeA7/decwlbYpFzzWz/cjJvMD1XYsr6CygH2rOMNX7H2anOtYXN+1rdr5NLBDtppWQYl7wJ2PYawhpWN+dko1QLlxDRsH/yzjdhVSUqpmoi8C0+hp4H3KaWeEpFfBs4qpR4C/jcReStQA5aBH/E/uywiv4JHLgC/rJRajitTWNgiBr1Zgioem/H0sMQw0WJ9DfXoxhpFliMTQ4Gvt2V97RBDQI/B4pkKzyINvh1txtPXS9Wmku2FfCZNIWsnRxemZxO07p8qhyYKRmVZ3w45P0NZ663RbSI2MQAopT4GfGzPa7/Q8vt7gPd0+Oz7gPeZkCMqxgsZ5pe3jN9XK/jAoQqLoaT1UjWURao32dp21TgxhKl68WSxE9bSyc3gxG0z0RrNIlVKBQ5vBEG5VqdSa4SWxQZZrodOPtvbPxulGmP5EGu2kOXC4qZxOVzhlj/5DPbcvh3FE2yT2axm2CjVGA252cGOElzfrjWJJwhs5YDCegw2Q0k6hh0U40NZag3FtuEKqbDhNbBXIRXao7NM3GHWrO2GmLaREAM7FqnpeHqUTWaLpIoRLFIwrwRL1TqVeiNwqELLYmOThSXusYK9eHoUjwHMz09YZaxlseXRpVMhKoGsh/rCEHfGWo7OBRJiYHc83SSK5fCbzGbYJFwM21c8hvvb6zEZzYcZEzvx9LDErePpdkgqHDHYOpgZNrwG9nJ0ekyChsrsnqkIPz/1hmKrcnM+kyEhBuxVM0TZZPbCWuFipLY2WRSL1FaFlFbwYUjKVjVQFIsUBmN+rK7ZEHKMWcox1P0+Y2HJEm7eRnoJMTBYm8xmBc5ghCqiWaR24uneCex0gBPYTVksxtMjzY9hjy5seA3sVUhtlKqMhjBmcpkUQ9m08fnRXm7QIhJw89x0m0iIAXtKcL1UQwRGc4NgfYWzSHWiesNw2CSqxwA2wibhlDHYCfU1GopiJaRFWrBjke5U0oUj7o1yjUbDfI4u/PyYN6yikeXN/UyGhBiwV+a2UaoymgvWE6gpi60cQzncJsumPetrw5IXFTZ8AzY8unBkqWUxrXiKFa/z7ViEMbGVfA4zP2OFDErBZsU8cYex0sGvYLPm+UcI9SXEcPPCVjwwisUzms9QLNWMxtMbDUUxJDGAt+HNewy6DUX/N1kki7SQGYhKIHseXXXX/YNAK0zjspTDE7cNw2ojihdlyaNzhYQY2NmQNjZZ2IU9VvDi6SYrpDa1RRqWpAoZNsqDE0qyoZDDk2XWghcVPu+STacoZFPN+Lc5WWoMZdOBnjioYW//RCPugQgl3eSP90yIgR23eRAW9s4mM7egorjC+npbOYYwFqk9jy4KcWf8hx6Z8+iikKV3vXmSCnveBXb2T7FsThallHcoM0RIC+zk6NYjEIO+tniTHnJLiAHvxHEunRosYjBoCUZVPOOFjHnFU66GtkibimcA5me0kKHWUJRr5jy6KBapvt5O+CY8QYHZHl+laoN6Q0UyZmx4UfreQaE9OtMetyskxOBjzIISjGqRep81t6C0JRfW+tL5DpMI25oDWk4cWyGG8IrHk2UAPLq8jRxQ+DEZt7Bmo5LlqL+PB8GjG82b97hdISEGH6OFjBVLI6wS1HXbJhXyeuRQki3FE25M8pkU2bQYnR/dmiOKFwV25idsBY6NUNJ6RC8KbK3Z8EZEtW7WowvbgFLDhsftCgkx+BgUJWgzxzAYiie8FyUixmWJOiY28lFRks9aFvPGTPCW2xo7VUkm5ydieM3K/IRroKdhw9h0hYQYfJgOm0RpFgd2QklRFc9YIcNmpU7d4MGlYjl8bTqYn5/oY+J7dAY3fLFUI5MSCtmKdXc4AAAgAElEQVRw23FQjJmRXJqUmA5/Ri+YaP28CUQJr3mymJ8fV0iIwceY4YMxkStN/FCSyaRVMXKMVFebmN1kYXMdYH6TRY8b2/HowjSLa8pSsJEDCp98FhHj3ssgzU/Yh/S0ynJLVyWJyJtF5BkROSci727z/k+LyFdF5HER+bSI3N7yXl1EHvN/Htr7WVcYM76woyfPWj9vRpYaKYHhkA/cGbcUIoi6yexUakXz6EwmwqMUKXiyZClWzLWiqNYblKqNUP2JWmUxa1hFDK9ZyHdEXbM2QrGuEJsYRCQN/A7wFuBO4J0icueeyx4FziilXgN8BPj3Le9tK6Xu8n/eGleeqDBtkTZd4ZCbLJ0ShnNp4wt7NB/eIrUT1orqlput8IhTIgqmFU/48A14+RGlvJYapuSA8GOiPzMIHp0d4g7XmbhVllu5XPVe4JxS6rxSqgJ8ELiv9QKl1GeUUvrZmV8Ajhn4XqPQ9c+mytwGbZNFUcam2y7o/vRRx8TkAar4oYr+E4Pp8x1RQ476MzaqkkZCNKCEHUPMfI4hwpj4UQhTHt2NjRL/fG6RTQdkY4IYjgKXWv6e91/rhB8FPt7yd0FEzorIF0TkbZ0+JCIP+NedXVhYiCdxG4wWMtQNtnaO6grrz2wYVIJRShC1HGDuRKtWHIOQY1iPOD8Zv7mgSZKKUqkF5nsURR0T0KE+c2NS9HNRYVqig62qvujzY7K54OefX+IH/vCLXF0rGblfN5gghnYz15YiReQHgTPAb7S8fJtS6gzw/cBvicgd7T6rlHpQKXVGKXVmbm4ursz7YDpsErUOG/xNZviAW9jqKLAxJuEb6GmYbi4YpYuohh2PLlo5JJgj7qglvGAn1BdrTAzJUqs32KzUI5ergjnvJU4UIixMEMM8cLzl72PAlb0XicgbgZ8H3qqUKuvXlVJX/H/PA58F7jYgU2iYDhFE6cioMSiKx3S8Nl54zWxzwY2IFqkni+lEeDziNjc/cbxc0+XE0das6VYUUctmvc+Y1SlRHhUcFSaI4WHgtIicFJEc8A5gV3WRiNwN/D4eKdxoeX1KRPL+77PANwNfNSBTaJiuwInSvrhVFuMlopFipGZPYcfZZM18hzHrOJpF6slizjpWKlpLdDDfiiIOcY+aNmYitNxuymKwFUWsMTFubFZJp4ShkCewoyA2MSilasC7gE8ATwMfVko9JSK/LCK6yug3gFHgL/eUpb4COCsiXwE+A7xXKdUXYhiz4PaN5NKRLFIvlNT/EtFCNkUmJQNCluaVYFRiMNnqYLNSpxGhJTqYb58StVILPGOmUm9QMpSj0zmGKDA5Pzvhz+g5OlOy6DEJW10YBUZ8EqXUx4CP7XntF1p+f2OHz30OeLUJGeLCdAVO1IQVmA0l6fbFUWTxWlGYk8WE9WVMCcaySDNcM5QAjBu+ab1HfFnih02K5VronkKdZDk+PRzpsya9F1NjYkqWqGQZFsnJZx+m2T2ORTpayLBlqBVFudag1lCRZTF5SCeORWq6AifO/AwKWQ77rSiMKZ5yjXwmRS4TXi3YKN6IY1gNQsLX9JiEfTxvHCTE4MNG8jmOMgYz1nGcEkQw26hNJwSjViWB2QqcODFsc4on+vzoVhRmvdyIxkzefI4uSvjGk8VgKGk7eiWdyX0M8eYnLBJi8GEjURTH4gEzidamxRPRBR0rZIxWJWXTQj6GRWqyAieOx1As14x4dHHKmr3PmWtFEddKBzNKsFJrUK41IodNxgrZgci7DGfTiBjMMZSTUJJzpFPCSC5tNh4YdbMbJKm4tc8m69Ojtubw5DCbY4h66K9VFhMHl+KcHdCymDz5HHdMTBB33LJMs15U9BxDKiVGe3zF8XLDIiGGFpiMp6+XorWX1nKAGWIoxljY3ufMueXFWOEbc2RZrtWp1MK3RNcwGTuOE0rSsgxCKGnMYCgp7piMFzLGmgvGybsARp+yV4xhbIZFQgwtMJu0ilH1YvBEaxxXWH/OpBcVVQ6TrShMeFFgxnuJK4vRHFDEZnFgtgIn9pj4zQXNeHTR9zGYLt5Iks99gakyNx0jjRPXh0EJJXljYqIVRdyFbco6NqGMvfuYIe44h5YGRfGYLPeOWzBh8mE9GzE8fzD3FLdyzX8UbZJjcA8vkdf/GKlJYmhussiWYNZYc8H1UjVSr3+NUUOtKDZij4lZ4o5zaMl8KCnamGR9j84EScXp8gpmw45xcoVgbn7ihoTDIiGGFozlMxQHIEa6E681Zx1HXdwmD5ZFfaynhqlEuAkvCjBEUvG8KFNkWW8oNiO2RG+VZRBCSSaJO2q7Eg1TT3GL0/QxChJiaMGghCpMtqIolqO35gCz1SaxQ0mG6tPjJ3zNJlrjWIHjhawfuozn0cW10vVnzayT+An51vvElSWOMh6UKERYJMTQAlOJvLgxUhExaH3FVzz6PnGw0ywuTiLPjPW1vm0mVGHmAGJ8i9SELHFaomuY9ujinGMAMzmGOJV0nixmHjCl5yepSuoDxgpZtip1avV4rZ1N9E036b3ElUPfJw50i484C9tUfXpcJTjse2Cm5ideeM3M/MQNOYLBUGzMElHjOYY4HkM+Q6naoBpTp2jij0PcYZAQQwuaB5fK8dzyOM9iaMqSN1NtYiKGDfGtLxOusH78alzoe0RVgroVhRlZ4pdDeveJSwzxypr1Z80ZM/GsdIjvRTUaimLFDHHHlSXJMfQRo814ejyFbGKTmSqd3ShVGTWgeOKSVNy4MewkN+O2oojTEr0pSz5jpBVFbOLOm1mzcZ6VoWGOGKL3SQLvOdEmWlFsVmooFc+LGi2YKSRJcgx9hKme/ybc8nFTmyxmVYWpUEXcnkCwMz9xDy6ZOEFqIt+x0xK9//NjIvxpqrlg3DFJpYTRXPxqLVNkCfH7nsV5lkkUGCEGEXmziDwjIudE5N1t3s+LyIf8978oIida3nuP//ozIvJdJuSJiuaDTwwsqEI2RTYdfXhNhSrixrBHcmaqkuI28wNzseM4z2LQMGEdl6oN6g01EGETU6EkMx5d/Pkx4XGbCN+Y6nu2Ua6Ry6TIZ+w/vQ0MEIOIpIHfAd4C3Am8U0Tu3HPZjwIrSqmXAL8J/Lr/2TvxHgX6SuDNwH/w79cXmCpzM7GwTZ1ojStL2m8EFlfxmDigY6oVhYkHnowVsuaswJjlkK33iixLjJboO7KYyUeZmZ/4a9ZMEYmZUJLXrsSNtwBmPIZ7gXNKqfNKqQrwQeC+PdfcB7zf//0jwBvEO+p5H/BBpVRZKfUCcM6/X19gKtFqYhJNtKKo1huUqvGP0ZtopGcq79J6r+iyxO85Y4IsNwzEjXeeUxF/zWZS0Vqia5gqbY57qAz8CjZDxG2meKP/azYMTBDDUeBSy9/z/mttr/GfEb0GzAT8rDOYOsxlRPEUMtQainItepmbiVwHmAmbmCrhhfgnjk0onkEZEy+8kDIgi9dZNc7zhI0aVgY87therskcQ2yPu+osvwBmiKHdStpr5na6JshnvRuIPCAiZ0Xk7MLCQkgRg2HcWKjCTCgJ4lWbmOqvYsr6EtnJWUSBsXhtqRq5T5KGiVYUJiq19OfjGjNxD3J5csSfn3pDmfEYBiTHYCwvFqPzbRSYIIZ54HjL38eAK52uEZEMMAEsB/wsAEqpB5VSZ5RSZ+bm5gyIvR/5jJlWFCaetKQTxnFIat2AK+x9Pr71tVGuMZrLkIpRImoqxxDngTQaJlpRFA0oHjDTGt1U3sW7VwxjxlBZ5rgB4jbRJqSQTZNLx/foTJBlGJgghoeB0yJyUkRyeMnkh/Zc8xBwv//724G/V17w/CHgHX7V0kngNPAlAzJFgogYCxGYiJHqe8WRA0wQw4CMiYEcg24WZyK8BvFIyuz89D+GbcJjMEUMJp77bMLLBX3+Jv78uAwlxf4mpVRNRN4FfAJIA+9TSj0lIr8MnFVKPQT8EfBnInIOz1N4h//Zp0Tkw8BXgRrwk0qp+P2dY8BEjyIzdfLxqxn0xoh7jN5EqMJEeG0k5z1DN878mAyvgTc/M6P5SPfYMBDD9j5vJhF+dLIQWw6Il6MzGV7TrSiilo2b8HI9WUwYVlVn7TDAADEAKKU+Bnxsz2u/0PJ7Cfi+Dp/9NeDXTMhhAnFbUehj9OYUT3y33ESoYhAsUt2KIhZZlvWzGAaHuOPOz2g+w+LGVqx7eMQ9FuseJh7vacqLam0uODWSiyyLifBNXOLWDShdtcOA5OTzPozGbB2sj9HHOVQGZipwjIUq8hnKtQaVGBVSxbIZV3g8ZgdP04onTlLeRGsOMHPmxUQMW7eLjxdeM+UxxK+QMvWM5bjGzGalTkO5a4cBCTHsw3hMdjfV7GrcpEVqKp4ek6RMPH0qbuw4bgM9DSPxdEOKJ26oQrfmiLtmTeToTO2fnbBWDOI2cELekyU7ENWFYZAQwx7EPdG6YSqGbSDRauoYvYlqE10nHxdxK3BMWaRGiNug4ilWajQitqLYrtZjt+ZolcVEKCm+xx2/gq1ogCwB/wFT/TfwwiAhhj2Ia/Ho6oO4k5hOCSO5dGzrK+4GAzPWcdwH0mjErU83bZHGVYKmFI9SUIzYXNBEWaZG7ByQheKAOLKYMmbiGnj6Pq6QEMMexG1FYaKL6I4s8a0vI+GbmG55uVanUovfmgPiP5PBlEU6aoAsTSqeOLKYXbPxreN0Sihk46kmE+FPU8aMXrNRdYqJBpRhkRDDHowVstQbiu1qtKrZoiHF48kS03uJ+bxajbgnwk3GSAclx5BNpxjKpmPLYqIEMW6oz2Svfy8UG+8cQ9zWHGAmFBv3IUoaY4UMDeUlkSPJkeQY+o+41tdOqMKMpf4vwSI1VQkE8Z9ToS3SoWz8Jr4mrGMjoaTY82Mm7wJ6fuJ6ueaMmagkpRtQmpmfeMRtoplfWCTEsAeDNIlmQkmDsLDNegxxSmd1O4y4FimYIAbTxB1vfkyRVHyyjL9Odtrb9D/vYszYTIihf4h7erNYrpES74HxJmSJu8lMKWPvfhEXtk7IG7WOoytBUweF4pQh1huKLQOtObQcEH1+jCaf/aqxOPF0E3Lo0tmo4U/TZOndM+Ka1eHPmK05wiAhhj2I+3hPrXhMWKRxG4FtGDot2WztHFEWk6GkuM/QXTeUkId4xG0ybjwe05hpNls0YKnrHN1WxHi6qUo60KHYqMrYXHhtp1Ny9FDfaD5+a44wSIhhD0yETcwpnuihpIbfvtjUJhuLceJ4pxLInBKMrJDLVWPVHXHKEE215vDkiFccYDJUEbcayNTZDtDtbfpvzMRes4a8qDBIiGEPTCTyTE3iWD7TbAQWFro1h6lNFiepaDrv0nrP8LKY22SDonh0K4o4VUkmWnOAmfkxF+oz4dENxpp12ScJEmLYh4FSPDFIylRZpkacCqmiBYs0Tg7I1JjEUjyGOqtC/FYUGwafDqY9oCjzo5Qyah2PD0XPAdnJi0VfK4nH0GeM5NKkJN4kmkxuQjSSMmmR6vvESZ4VsqnI7Y9bEfe5wmaJO8t2tR7JozPd5iBO2NFTPObyLhBt/5SqDWqGWnNoWQYhBzTse2NxPG6XZxggIYZ9iNva2eQkxtlkJmvTIW7YZDDGBHT/G7OyRInt2yHu6GEtk2Tp3TOKMWOWLMdjVI2ZPA0eW6cY9HKDIiGGNohThmg2VBH9uc8my+08WaI3rzPVWgDitaIoVetU6o2BCPWZbnMQhxjWDcf1IRpZamVsrmDCW7NRmgsWyzWyaSGfMaMi4xK3qTEJilj/axGZFpFPishz/r9Tba65S0Q+LyJPicjjIvJvWt77ExF5QUQe83/uiiOPKcTdZIOkeMyW/sWoqjCkeLLpFMO5dCSyNNn6wbtPfOI2WcEW2Zgx+HQwE3kxk/tHKa8QIyx0iaiJsnNPluihPlMn5MMgLh2+G/i0Uuo08Gn/773YAn5YKfVK4M3Ab4nIZMv7P6uUusv/eSymPEYQNZ6um8UNwiazoXiK5Rr1CNaX6Rhp1PkxHb6JU4ZYLJtpFqcRp7WzyVDSSC6DSMxQkqFQ33gh+vmBosGyc/DWXBQ5dGuOmy3HcB/wfv/39wNv23uBUupZpdRz/u9XgBvAXMzvtYqoNftF4+Gb+PFa00owSjjJpOKB+PNj4iCXlgOiFweYas3hyRK9OMBkwUQqJYzmoilB83mXePNj0kqP2uPLZNlsGMQlhoNKqasA/r8Hul0sIvcCOeD5lpd/zQ8x/aaIRHuqumFEDSXZSCi23jcMTLbm2C1L/zdZ9PkxXQkUz6MzTZZRWlHU6g22KnULHl3/lWCs+TFcIho1lGTque1h0ZMYRORTIvJkm5/7wnyRiBwG/gz4H5RSur7vPcDLga8DpoGf6/L5B0TkrIicXVhYCPPVoRHV+jI9idl0ikI2FdlKNx0jhWgeg8lySC1LJIKyEMOGOGRpVhlHae1s+ryLJ0u2+cCqMFg3XUkXc34GgSxNj0lQ9FwNSqk3dnpPRK6LyGGl1FVf8d/ocN048LfAv1VKfaHl3lf9X8si8sfAz3SR40HgQYAzZ85E69AVEDpUoZQKpVhtTGJUJbhuOK4ftZFe3W/NYdL6Gi9kmF/eCv25DWuhpGjei2mLVN83jGFi2svV94rjcZsPxUbLAY0VxozI4cmy01wwjE65WUNJDwH3+7/fD3x07wUikgP+GvhTpdRf7nnvsP+v4OUnnowpjxGMFTLUGopSNdzBJVubLHryzKwcEN76Ml1p4t0rWgWO6bxLnOaCxbK5Si2IHjYxXb0G0SvYNkrmWnMAjA9FP4VtI9QXpbmgDZ0SBHGJ4b3Am0TkOeBN/t+IyBkR+UP/mn8NfBvwI23KUv9CRJ4AngBmgV+NKY8RRE1a2WD3qIlWGwtb3zecHOYfMjIegyzBfNgkTvLZnBzxiNtsWCtqPN2sl9usStoOJ4tuzWE6LwYR9o/B1hxhEOvblFJLwBvavH4W+DH/9z8H/rzD518f5/ttobWN8YHx4J8zfdpYyxItnl7lwFjBqBwQfmGb7AmkMVbIUKk1KNfq5DPBk+smW3NoRCYpw6dZo7Z2tkHccUJJJsckn0mRTYd/WI/p1hyw29g8NBF8X/bjsZ6QnHxui6itg21UEMSp8DApR9QTx3bCa1G9F7MJX0+W8POjlDJ+tiMqcQ9ajsGkHF5zwfDeS9NKt+DRhSVuk605wiAhhjaIGkraKNXIZ1LkDB2jB92jqP+hiqFstEZgNryo6PH0qvHWAlEUT7nWoFpX1pLPYbBhoyopn6FS9zy6ULJYaBYX5fyAjbzLeIxQn8nWHEGREEMbRI8Hmi1xg3iJPJOyRG3tbLrSBKJ3WDUdvoFo1rHpPklajtZ7B5fFG0NTp/U9WSJ6dIYT8lqWqLlCs55/9BzdWCFrrOw8KBJiaIM4HoNpl2+skGGrUqcWorWz6WZxrbJEGRMwa3013fLt8Ap5MMbEvBcVtbXzRsm8RRqnQsrG/ITPu5iP649HJAbTIeGgSIihDeJYX+YXdviDZTZKRMGrXAmbd7GxyeJUjZk6w9AqS9SEvMkNH7W1s1Y8Ji3S6IaVjf0TnriLFiqBolaN2SDLIEiIoQ1G/UZgYS0NG+wehaRs1T5Hs76qZAw2i9NyePcOL4uNUFJYj87W/EQhBhtx/SjzY6tZXBTitpHw3fHo+u9FBUFCDG2gG4ENArtHqTZphioMW8dRE3mjBpvFeXJEa3e9XqoZjaVDNI/OdM+mHVkieAwGG+hpRDklbyOuD95aCU0M/rmHiWGzOTqPuMMXB5iupAuChBg6IFpSsWqhHDK8W27jIJeWJUrC17hlHIEsa/UGxXKNiaH+W8c7D6QxTdzh58fk80Na5YBwa9amlxu2Xfz6dhURL3JgWpYoJGV6zQZBQgwdEEUJrpcGS/HYiddGCFUYJst0Knw8vZkEH7Lj0YXxXmxYpBC9QmoQQkm2msVFeaLcesmrjkoZas2xI0s2/DmG7arxNRsECTF0QNhNZs8i9a2vEN0q1y2UIILn5odt7WzDIoXwSUVbYxKlDNGqRRqyq6kNizSKR2fLY4gSdvSUsXkrPeyarTcUG2Xz4c8gSIihA8ISgy2LNJL1Zc0i9RqBbVeDH1xa3zaf3PRkyYTc7Hp+BsE6tmeRRglVmF6z2XSKoWw6EnEPhsdtJ3wTNke3YWlMgiAhhg4IG0qyZ5FGIwZbFmlYWTYshNc8WcIpQXuKJ3w8fW27apy0PVk8xRPUo9MWqY35GfVj+0Gxtm1nfrQhEHZ+bFjpY4VsOM/fkjETBAkxdEBYj8HWws5n0uTSqVDWsV7Y5i3S8LXYa5aSZ2GtLz0/g+LR2VI8YTw6G6eed2QJNz/ayx0Mw6pmJa4fekya85PkGAYGoS1Si+wehaTsWOnhNpnOu9jZZCE9OuuKJ5x3OQjzY8uY8WQJ98wM7eXaOiAaSpaSLeIO59HtGDOJxzAwGCt4jcBKAa0vu5ssrKVhTxnr+weBlnkwxsTOJstn0uQyqdDEbUvxQHCSsmnMRPHo7ORdInp0lsgyjEe3blGn9EJCDB0Q9mDZjuLpv3Vsy2PQ9wz64BPbFmkY62t9u0Y6JYzkgj+/ISjCPpPBVqhiPCRx2zdmwljpNWveNgQny1q9wWalPhAenS1jJghiEYOITIvIJ0XkOf/fqQ7X1Vue3vZQy+snReSL/uc/5D8GdCAQNqloc5OND2VCPYXKVgxb/9/WQhKDLevYa+0crBWFZ6WbPYG9I0u4HlK2Q31B14pNY2ZiKMtaiCaHtsYkn0l7j18NrIzNN33UiKpTbsYcw7uBTyulTgOf9v9uh22l1F3+z1tbXv914Df9z68APxpTHmMI+2CN9W2vJ9BQ1rxFOjmUC6yMwb7HEFSWZiWQhQqcsAfL1kt2wgPgl84GHJNKrcF2tW6VuMN6DHZkybG2XQnh0dkxZiDcwbJ1i3H98DqlRkrcP9YT4hPDfcD7/d/fD7wt6AfFM91eD3wkyudtY3I4vBIcH7LTN318KBuKGGwpwULWi6eH9RjseFHhDpbZVDwTQ8ETrRsWwwOagNe2KoGutxnDnhjKUq0Hj6fbMmYg3ONxbZWdt94zjEdnS6f0QlxiOKiUugrg/3ugw3UFETkrIl8QEa38Z4BVpZTe2fPA0U5fJCIP+Pc4u7CwEFPs3tCLdDXgJlvbtlMPDh5JrW1XA1lf5VqdUrVhT5ahLGtb/Q8lhd9k9uZnIsKYDIJHt7ZdJZ0Shi3kXaLIYmt+xobCeAz2EvJRxqQfp54BevooIvIp4FCbt34+xPfcppS6IiKngL8XkSeA9TbXddR8SqkHgQcBzpw5E7wnQ0RMDHnpjsCKZ9v8YyN3ZNmxvoZ7HFqzHZecCOG96E1mRQn61vFqiE12cDxvXA7wiDuoHM0YtoW4fj6TZiibZjUESU1Yski1x726VeXwxFDP6z3r2M6anRwKMz/2iDt0FKJPfZIgADEopd7Y6T0RuS4ih5VSV0XkMHCjwz2u+P+eF5HPAncD/wmYFJGM7zUcA65E+D9YwY7HEILdLVrpWpZexGD7tGQYYljbrpJLp4w+i0FDj0lQS912KEl7dL2UrE2PoVWWIPDakNszICCYErTu5Q5nubi0GehaWwchIbxOsenl9kLcHfsQcL//+/3AR/deICJTIpL3f58Fvhn4qvLiIp8B3t7t8/1CLpNiOJcOZWnYVMYQbJPZVjxhrOM13+KxY5F6Hl3QUJ/N+ZkcylFvqECVSbYO2jVlCTk/NglKf0cv2PQsIaTHYHF+sukUo/lMOGOzT6GkuMTwXuBNIvIc8Cb/b0TkjIj8oX/NK4CzIvIVPCJ4r1Lqq/57Pwf8tIicw8s5/FFMeYxiMlTYxOImGw5uadiufR4fyoZOnlmRw7d0g2x42xbpoM1PuFCFZWIIMCa2T/hODHtVfY0Az2RYL9nLu0DYUOwA5xi6QSm1BLyhzetngR/zf/8c8OoOnz8P3BtHBpuYGM4F2uxKKe/QksVQBQS1vgYoVGGRLDPpFGOFYNZXM7zmIGxyvMe11j26oSwXl7YCXbu+XeXoVO/4fxRMhIinWyeGoSxK+Q0de5ROe/vYjperZVnbDu7l2ij1DoLk5HMXTAQ8WFaqNqjUG9YSRTuKp/eCsh2qmBjyDnNVAzzj2LYrrKu1esG2lT4ZMmySS6fIZ+xsvXA5BnvzM5rLkJKAY2K5vXQzRxdk/1jqY9WUZTgbyJjRXm4/DrdBQgxdMTmUC7yYwGZc34unh7O+7FV4QLBqLduPJZwcygXKMdiP6+t8RzAlaCvv4smSDbRmlVJWcwyplDAxFEwWF3kXCDY/a5aeH9IqS5DwZz9bbkNCDF0RlN1tL+yRXJp0SgIv7EI2RT5jKUYaMkRgkxiCWsfWyXI4uEVqs3oNvDEpVXs3fyxVG1Trymo5ZNC2GLbDn5MhSptXt6rN621gYihYeNrmQbsgSIihC8IrHjuTKCKBE+HrFg/aQfB8h1LKWpfXpiwBrS+9EbVlb1yOkDkgm5t9YjjY+Rubhw+bsgwIceszSUG8y9WtirV14sniFW/0Oqyq1+zNWq76LxoTw1nKtd7WV1PxWFbIQZSg7bh+UCW4WalTbyjLoaRgJ45XfIUwZWnDF7Jeo7agskxZtUiDzY/t8Cf41UABlLFtLzfMwbLV7arV+ZkczlKpN3q2CtEkNjXSn76iCTF0QdADKVrxTFucxInhYGWiq9sV664w9N5kLixSHa/tZX2tbHkPgRmEpOLKZtUaQUFrorXHmt20S5YQ3GNw5eX2mp96w8u72DTwJgPrFO99myTVDQkxdMFkQCWoicGuQg6meFa37CqeoBbpqoMxCXqwbHWrwnghS9rwQ2BaETTRajtU0bSOAxozducnGDF4XpS9MQl6sMwL8dgLOULwRPjO/kk8hoHDziR23/ArW17LbZvtcYNusuVNu5ss6MGl1abFY9eLav2uTlOEcKcAABrzSURBVFjZshsegGCt0Ss17yEwLkJJPT0Gf8ysern+mu11sGxly66Xq2XpRdzNkOOIPVnGm/PTW5Z0SpJy1UFE0E2mrUCb7XE9j6H7YlJKeVUVFhe2bhXSSwkub9oPrwU9P2DbSgc/ER7UCrQ6JsG83GVHoaSGgmKlu0e3slW1uk7AP/MSMHxj1WMI2JxzZcsLafWj5TYkxNAVQcMmXtzYssUznGOjXOtqfW1W6lTqDaZtK8EAifAVB8mzoOcHbJcgQrB4uou48Vghg0jvZzKsblXIZ1IMWWr9AK3Ph+id77BJUBDs/IA+QGo1xxAilGR7zXZDQgxdECZea3thtx7r7yiHAytQy9JrYWuL1MkmC+CWW1c8gYjB/vykUsJ4IRhJ2bbSgxhWjYZyND+9D0OubNoPfwY9U2G7SKEXEmLogtF8xjtYFkTxWAzfwI6CXemyuF1Y6QAzozmWN8tdr1nZrDBeyJBJ21tiQatNXHgMk8NZtip1yrXOZYiauN3E03tb6bbDa5p4tJHQDhulGg1lf81OBGif4oK4h7JpculUoErHfiWeISGGrhDxjvWvBEpuWt5ko979l7psMlclblPDua6bXcsyCBZppdagWK7Z96J025Iua8VFwhc8Jdt7fipMWzZmghDDcrPU275htbrVvbR5datKSnaezWwDIuJ3wO0V6rMfnu6GhBh6YGYkx3Kx8yR6CV/77D4TYJO5OhQzE1Dx2JajkE1TyKa6hgi0t2d7k+m8znIQj87yWpkdybHUZc2C9qLcrNluxkwz5Oggx1BrKDYrnT261e0KE0NZUhbLmsFbi4Owf7ohIYYemO6hBDcrdap1ZV/xNImhcwjHRaWJJ0ue9VL3Dqu2y2Y1ZkbyLG92ttJtt8NoyuF7dN2MiNWtCoVsikLWXsIXeq9Z8AjM9prVZ0e6rVlN6rYLJvT8d5sfF54/6FBsZzm2K3XKtUaSfB5kzI7mWeqysF0lfGdGvOcVL/ZY2LZP+MJOWGulW1jLFTGM5gLNj+1Npq3jxR6hPhdjMu0rnk5hE33C17YyTqWkZ9jRRVkzwGwzFNudpFwo45nRfI+QsBud0g2xiEFEpkXkkyLynP/vVJtrvkNEHmv5KYnI2/z3/kREXmh576448tjA9Egu0CTaXlBDuTTDuXTXTbay6bnCNk/4QrAQgZdjcLDJeoRNdvIutj0Gj7iXir0UjwsvKkel3uh4ItzFCd9WWbrNz45HZ5u49fx0lmWpWLFOUNA71LdDDDevx/Bu4NNKqdPAp/2/d0Ep9Rml1F1KqbuA1wNbwN+1XPKz+n2l1GMx5TGOmVGvTW6tQ9ikqXgcLKheIQIXZX9aDuic79iu1Nmu1p2MycxovqsyXvTfmxvLW5VjcihLSnokWjftJ3yhtxLcSfj2f80ub1XIpu12DYCdUF83j2GxWGF21O46AS8Uu7ZdpVJrr1P0vE2P2JelE+ISw33A+/3f3w+8rcf1bwc+rpQK9uzBAUAz6dshqbi44S00FwtqJoD34sQV7uExuHSFZ0ZzLHYJm2hisK0EUylheiTXNdTnTPH0qGBrhj9dEEOPeLoOOdo+4avHvdP81BuK5c2ym32sQ7GddIojY6Yb4hLDQaXUVQD/3wM9rn8H8IE9r/2aiDwuIr8pIh1HQkQeEJGzInJ2YWEhntQhoEMEnRa3y0n0rK8uFs9GhTkHC1srlE45BldJcIDZkXyzJLUdFotlpoazZC2ep9DwEuGd52dhw5Hi6eHRLfjGjIu10suYWdp0E74pZNOM5jMdvaiVrQoNtZOLsIlmvqODLFqnuJClE3ruFhH5lIg82ebnvjBfJCKHgVcDn2h5+T3Ay4GvA6aBn+v0eaXUg0qpM0qpM3Nzc2G+Ohb0ou00iQsbZQrZFCMWWwvsyJLvWlWxUCw7ISjPwutskS44JMuZXptsw42VDn4+qoMcm+Ua29W6ozHpnu9wOT/TI15zwU4VbAsbbtYsdC9U0PM263J+OsiysFEmn0lZD691Q89vVkq9sdN7InJdRA4rpa76iv9Gl1v9a+CvlVLN2kLtbQBlEflj4GcCyu0Msz3c8kVfGbtoduUtbC9ssvf7qvUGy5sVJ5ssnfKeKNfJOl5Y914/4HiTnZgd2ff+YtGNle7JkuOpK+tt31twHHKELsS9USYlbnIMWpaVrQoHxgptZTk1t3/ebMmy2IEsd6x0h/PT0WPwjJl+NdCD+KGkh4D7/d/vBz7a5dp3sieM5JMJ4o3A24AnY8pjHNMjva0vFy45eAuq7Ldu3gu9yFxZX92Sik49Bl0m2sUtd2EFgl/a3EPxuBiTQrZ7Bdtiscz0SN569Rrs7J92siilHHsM+QDhGxdrVuc7Oq+VfuYXID4xvBd4k4g8B7zJ/xsROSMif6gvEpETwHHgH/Z8/i9E5AngCWAW+NWY8hhHr2oT16EKaH9Ix2XcGLwNpL+znSxjhYz1g1xaDuhlfbmJ1U6P5Fgv1dpWm7iOG3c7ROU6fAPt52d9u0al3mjrSdjA7Gjn4gCX+2d8KEMmJV09Olc6pRNiBbGUUkvAG9q8fhb4sZa/LwBH21z3+jjf7wJetUm+I7svFMucObHv+IYV6M18Y6PEbTPDe+Qo7brGNg5NFHj0xdW2793YKDkJI0FrDmj//JSqdYrlmtNQEnhGxKGJ3cquqXgcei/diNuVHAfHvXG4tlbaL4fjNauLAxoNta/txWLRK5sdH7If1xcRj7i7eC933zZpXY5uSE4+B8CBsTzX1/dvMh3Xd6V4Dk8MAXBtvc0mc6x4Do0XuLZealsm6lLx5DIpxguZtsTt2ovS39NOIS8UK4jYb/2goeenHRY23IU/D/nEcH1jvyw31l17uTkaqn3L68VimZkRd3H9mZF8M+TaCq9s1k11YTckxBAAhycKbS0e7aq7VMbQwfpymNwEODBeoFJrtO1semOj7Cw8oGVpR9zN8M2YG2WsifvK2va+9xY2ykwP56y2IW/FwfH2a1Yp5ax6DbwT++OFDNfbegx+kcK4I+L21+SNNiS1VCw7Wyfg6ZSrHXRKQ7mpjuqGhBgC4NBEe+vLtZU+PpShkE11JIZxR3F9aCGpDuPiMnnmbbL9yljHk12RpQ4ftZsf1wnFwxMFiuUaG6XdxO2Vjiqnshzs4L0493L9+bm62p6kXMb1D08WuNZ2zbo18DohIYYAODxRYHmzQqm6uxpowfEkikjHEIFLKxDgoG/l7VWCxXKNrUrdWY4B4MjEEFfaKGM9TjrObRszIzly6VRbS9B1QlErwet71oprZaxludbGo7vh1+uPOarXPzLpjUk7j+7aWonDE+683MMTQ6xsVdneU2F4ow/z0w4JMQTAIT9EsHeTaaV40JEr7H1XYZ8c4MVrXVuB+ntb0Q/Fc3iywGKxvK8a6OrqNpmUOFPIqZRwcCLf1hJ0rXh2kr675+eG47yLlqVtKGmjzIFxd3H9A2MF0inZ5zGUa3UWixUOjQ85kQM6e9xXV72143KttENCDAGgJ3GvJXhldZuU7LzvRJYOYa0rq9scmXC3sHVceK8smrRc5hgOTxRQaj9xX10rcXC84KRevynL+NC+dVKtN7i+UeLwpLv50Yplb4jtsq94jjqU5dB4gYVimXpjd6HC1bVtDjpcJ+mUcHAsv89juO6T5+FJh2t2sv38aJ3iysvthIQYAqBT7Pjy6jaHxgvOEopalutr5V3VQNV6g2vrJY5Oudvs+Uya6ZHcPmV8ecVb6MccyqKTvu2I+4jDzQ7tifvaWgml4KhDWbRiaTc/Iuwrp7Uqy0SBekPtKym+vLrtdM0CHJ4c2ucxaKJwaVg11+w+WUocGCs46e3VDQkxBECTGNbbKR63C/vQeIGKXyarcW2tREO5tQKhfVhr3icGl9bXkQ7W19W1UnMDuoKuNmkl7iu+le5yrRSyaaaGs23X7IGxPLmMu61/cGy/d1lvKK6ulpyv2XaFCtrgc0mWnUJJV1a3ne6dTkiIIQBG8xnGCpl9HsOV1ZJzYtgJEezIcrkPigc8C1gTgcb8yhYHx/PkM26qo2AnB9Q6Jo2G8uL6ffAYKrVG8zkdsGORuiapg+OFfRbp5T4YM/r7LreslevrJWoNxbGp4U4fsybLPuJecx/XH8p5xK2NBo2ra+51SjskxBAQRyeHuLS88xiJRkNxdc09ux+f9jbSiy2y6MXl2i2/fWaEi0tbuzbZ/Mq2882uiftqyya7vlGiUm84l0Ur/1YlqH93HdY6Pj3MpZXdjz65vLrt3ErXp/QvtqzZy31as4cnCpRruz3ui4tbzI7mGHHczfTQxNA+Y6Yf89MOCTEExImZEV5Y2mz+fXl1m2pdcWLGTWdIjdv973thcUeWS8vuE4qeLMNsV+u7Tvr2a2EfnxreRZYvLHjjc6pNx1WbODHrKcHzi8UdWRY9L2o451bxnJwd4cLSFg0/6VurN7iy6p64xwtZZkZyXGzZP/M+YbleK/r/fqmFuF9Y2nS+jwFunx7mwuJunVKpNTjpeM22Q0IMAXFqboQXl7aaj/h8fsHb+HfMjTqVYzSfYW4sv2tBnVsocmxqyNnhNg1NUheWvE1eqtaZX9nixIxbxQNwx4FRnl/YGRNN4q432YmZEUR2E/f5xSKnZt2uEy1LpdZohkpeXN6iWle85IB7WW6fGd41Js/f2CSdEm6bdrtWdIvv52/sEPeFxc22Ldtt444DI/6ceDpFj49rY6YdEmIIiJOzI9QaqhlT10roDke95HfJMjPChRbr69yNYl82uyYATVLnFzZpKDh9cMy5LHfMjXBpZat5CPGFhU0K2ZTTUmLwkr5HJ4c4768PpRTnFzadPXOgFZoULyx6xH3uhjZm3Mtywg87ajx3Y4MTM8NOk+DgWenZtHDON+yK5Ro3Nsp9sdJPzY5Sa6imp3vel+lkH+ZnLxJiCAi9sXWI4PxCkYmhrJOHnbST5dyNIkop6g3F+YUiL3HsuYAXBihkUzxzfQPwNjvAS/tCDKMotWN1Pb9Q5MTMyL4umi5wcnakuU6WNyusbVc51Yf5aVrHvsJpGjN9MCLuODDK1bVSs7fWczeKnD7gfp1k0ilOzIw0SbKpjPtBDHu8lxcWNxnLZ/reQA8SYgiMlx4cQwSemPee0PX01XVOHxjty1OWXnl0gpWtKvMr21xa3qJca/TFY8ikU7zi8DhPXF4D4NnrG2RS0pdNpv//z1zbQCnFE5fXuPPIuHM5AF5+aIxnrxep1Bp89aq3Xl7WB7I8MJZndjTH4/M783NwPM94IetcllcfnQDgqctrlGt1Li5tcfqg+zULcPrgKM/6xoweGy2fS7zkwCgi8PRVT5YnLq/xskNjfX1ym0YsYhCR7xORp0SkISJnulz3ZhF5RkTOici7W14/KSJfFJHnRORDItK/p1/3wFghy+kDozx6aYVStc6Tl9d53e1unsOwF6895i3ix+fXOHtxBYC7+tS//dVHJ/jqlXUaDcWjL65y+uCY8/AAeMQ9ms/w8IVl5le2WSxWuPt4f8bkntumqNQaPHVljUcurpASeO1x94pHRLjr+CSPXfLWyMMXlrnntv6sWa14H7+8xuPza9Qbilf1QRkD3HV8kotLW9zYKPGVS6tMj+ScHsjUGCtkeemBMR55cYVyrc6TV9a5p086ZS/i7uAnge8B/rHTBSKSBn4HeAtwJ/BOEbnTf/vXgd9USp0GVoAfjSmPVdx9fIpHX1zl8fk1KvVG34jh5YfGyaVTPPriCl88v8TksLfA+oHXHJukWK7x6KVVzl5c4ZvumOmLHOmUcM/tUzx8YZnHLnkPELrreH/mR2/uRy6u8MjFFV56cIyxPljpAK89NsnzC5t87do68yvbfN2J6b7IMTWS4/j0EGcvLPO5c0uIwDec7M9audf/3odfWOGxS6u85thE36z0152Y4tGLKzx5eZ1KrcE9fX5Aj0YsYlBKPa2UeqbHZfcC55RS55VSFeCDwH3+c55fD3zEv+79eM99Hlh800tmWNuu8huf+BrplHCmT5ssl0nx9aem+S9PXeMfn1vg609O9yWWDvD6lx8gkxJ+/q+foFJr9I0YAL7h1DTPXi/yJ5+7wORwlpcf7g9ZHhwvcGp2hL957DJfemGZbzjVvzH5ppd43/3v/sZ7nPrXn+rPmgX4zjsP8Q/PLvDRr1zmzsPjTAz3hyxfeWSckVyaP/7nF3juRpFveclsX+QA+IZTM2yUa7z340+TEvqmU/bChc9/FLjU8ve8/9oMsKqUqu15fWDxxlccZHI4y8MXVnjjKw70JfGs8X1njjO/ss319TLfe8+xvskxPZLj2182x9eubTA7mueb+7jJ7rvrKOmU8MjFFb777qN97TfzfWeO8+Tldcq1Bm9/Xf/m5+7jU7zkwCgPX1jhFYfHufNwf/IuAN9zz1Gqda9K67vv7t9Wz6ZTvO3uo5y9uEIuneJfvfZI32R54ysOMOXrlNe//GDfn8Og0fPEjYh8CjjU5q2fV0p9NMB3tDNlVZfXO8nxAPAAwG233Rbga81jJJ/h93/wdXzyq9d54NtO9UUGjf/+1Ye5sV5CRHjTnQf7KsuvvO1VHJs6z313HXF+lqIVRyeH+N0fuIdHL63yE99+R9/kAPgfv+UE5Vqd26aH+xZLB68V+O98/z188OEX+cFvuL2vic1XHpngN97+Gi6vbnP/N53omxwA/9ebX85IPsM33jHT106mw7kMD/7wGT7x5DX+pz7rlFZIu2f2hr6JyGeBn1FKnW3z3jcCv6SU+i7/7/f4b70XWAAOKaVqe6/rhjNnzqizZ/d9VYIECRIk6AIReUQp1bFQSMOFr/0wcNqvQMoB7wAeUh4jfQZ4u3/d/UAQDyRBggQJElhE3HLV7xaReeAbgb8VkU/4rx8RkY8B+DmEdwGfAJ4GPqyUesq/xc8BPy0i5/ByDn8UR54ECRIkSBAfRkJJrpGEkhIkSJAgPAYplJQgQYIECW4iJMSQIEGCBAl2ISGGBAkSJEiwCwkxJEiQIEGCXUiIIUGCBAkS7MJNWZUkIgvAxYgfnwUWDYpjColc4ZDIFQ6DKhcMrmz/EuW6XSk11+uim5IY4kBEzgYp13KNRK5wSOQKh0GVCwZXtltZriSUlCBBggQJdiEhhgQJEiRIsAu3IjE82G8BOiCRKxwSucJhUOWCwZXtlpXrlssxJEiQIEGC7rgVPYYECRIkSNANSqmb6gd4M/AMcA54d5v388CH/Pe/CJxoee89/uvPAN/V657ASf8ez/n3zLmSCziO15b8aeAp4Kdarv8l4DLwmP/z3/ZhzC4AT/jff7bl9Wngk/6YfRKYcjhmL2sZk8eAdeB/DztmUeXC6xD8GaAI/H97PvM6f7zOAb/Njrdufbw6yQUMA38LfM1fY+9tee9H8J6XosfrxxyP12f9e+rvP9BrTTgYr7E962sR+C2H4/Um4BF/HT0CvN7k+tolQ5CLBuUHSAPPA6eAHPAV4M491/wE8Hv+7+8APuT/fqd/fR5P4T/v36/jPYEPA+/wf/894H9xKNdh4J6WBflsi1y/hPdgpL6Mmf/eBWC2zff9e73YgXcDv+5Srj33v4ZXtx14zGLKNQJ8C/Dj7Fd0X8JrTy/Ax4G3OByvtnLhEcN3+L/ngH9qketH9v4fHI/XZ4Ezbb6v7b1cybXn848A3+ZwvO4Gjvi/vwq4bGp97f252UJJ9wLnlFLnlVIV4IPAfXuuuQ94v//7R4A3iPc8w/uADyqlykqpF/CY9d5O9/Q/83r/Hvj3fJsruZRSV5VSXwZQSm3geQ5RHpRrY8y6ofVeTsdsz2ffADyvlAp7EDKyXEqpTaXUfwVKrReLyGFgXCn1eeXt0D9lZ1ysj1cnuZRSW0qpz/i/V4AvA2EfUG1crh7otCacyiUip4EDeGQaBnHkelQpdcV//SmgICJ5Q+trF242YjgKXGr5e579yrJ5jfIeErSG5xp2+myn12eAVf8enb7LplxNiMgJPGvhiy0vv0tEHheR94nIVAe5bMqmgL8TkUf853FrHFRKXfXvdRVv87iUS+MdwAf2vBZkzOLI1QlH/fu0u6eL8eoJEZkE/hXw6ZaXv9cfr4+IyPE+yPXHIvKYiPy7FuUf9F5Wxwt4J54l31q943K8vhd4VClVxsz62oWbjRjaWQZ7y6o6XWPqdVdyeR8SGQX+E16sfN1/+XeBO4C7gKvA/9NBLpuyfbNS6h7gLcBPisi3dZHBpVz4j5B9K/CXLe8HHbM4cnVC2Ouj3iPS94hIBo9Ef1spdd5/+T/jxbZfA3yKHavTlVw/oJR6NfCt/s8PhbyXtfHysdfwcDZeIvJK4NeB/znEPUPhZiOGebykrMYx4Eqna/wFPwEsd/lsp9cXgUn/Hp2+y6ZciEgWjxT+Qin1V/oCpdR1pVRdKdUA/oDu4R0rsmmXVil1A/jrFhmu+66tDqHccCmXj7cAX1ZKXdcvhBizOHJ1wjy7QzSt93QxXr3wIPCcUuq39AtKqSXfGgVvvF7nUi6l1GX/3w3gP7IzX0HvZW28ROS1QEYp9UiLvE7GS0SO4e23H1ZKPd9yfdz1tQs3GzH8/+2bu0pDQRCGvwFFMYigNpYGAj6AhU9gISIIFrHRV7AVrUSwEwRbK5/AgK1WVlYJiogXfAErtT4WOyFnjzleiFkU/g+WLCdnJ8O/y85eJpdAzcymfVVYBxqFdxrAutdXgDPf7jWAup/JTQM1woVNV5ve5txt4DZPUvnlW+cj4CbLsv28oXZHO8vAVYlf/fKtYmaj7ksFmM/5kLeVVLNcu1UKx0g/0KwXv7riW/gXM5vzfl2jo0sKvUoxs13CxLNReJ7Xa4lwx5XELzMbMLNJrw8Ci3QfX5/Z6otezlfjqy96+XHfKbCZZdlF++VfGl8xZbfSf7UAC4QMnQdgy5/tAEteHyYcIdwTJotqru2Wt7vFb+3LbPrzqtu4d5tDqfwiZEVkQItCiiVwTEhNa3nHT6XUzHVperkuaDZBOKe+88/xxH05AjwDY4Xf+rZmPfr1RFjdvRJWcu1MslnC5PYAHNJJJ0yl1we/CCvLjDCJRWmWwJ73bZOwQJpJ6FeFkPHTch8O6GTDldpK0Y/+3WNRjxR6AdvAG3HKbDuNt+fxlS/657MQQoiI/3aUJIQQos8oMAghhIhQYBBCCBGhwCCEECJCgUEIIUSEAoMQQogIBQYhhBARCgxCCCEi3gF8B0BeJTUv2AAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# 另一种生成正弦信号的方法，生成时长为t的序列\n",
    "def generate_sinusoid_2(t, A, f0, fs, phi):\n",
    "    '''\n",
    "    t  (float) : 生成序列的时长\n",
    "    A  (float) : amplitude\n",
    "    f0 (float) : frequency\n",
    "    fs (float) : sample rate\n",
    "    phi(float) : initial phase\n",
    "    \n",
    "    returns\n",
    "    x (numpy array): sinusoid signal sequence\n",
    "    '''\n",
    "    \n",
    "    T = 1.0/fs\n",
    "    N = t / T\n",
    "    \n",
    "    return generate_sinusoid(N, A, f0, fs, phi)\n",
    "\n",
    "A = 1.0\n",
    "f0 = 440\n",
    "fs = 44100\n",
    "phi = 0\n",
    "t = 0.02\n",
    "\n",
    "x = generate_sinusoid_2(t, A, f0, fs, phi)\n",
    "\n",
    "n = np.arange(0, 0.02, 1/fs)\n",
    "plt.plot(n, x)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Scipy FFT\n",
    "\n",
    "介绍如何Scipy的FFT模块计算DFT\n",
    "\n",
    "注意，理论上输入信号的长度必须是$2^n$才能做FFT，而scipy中FFT却没有这样的限制\n",
    "\n",
    "这是因为当长度不等于$2^n$时，scipy fft默认做DFT"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD8CAYAAAB5Pm/hAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3Xl0XOd53/HvM1gJLgCBAVcQBLGQImXJFAWJCyiJWi0rjhXHUrzVVlIlTGylsVrHrp30ZGlPW/vkHG9tjms29rHdOpbt2q4UxamjyLYWUqJIURS1UCQArhApYgdIYp2Zt3/cO+AQBEAsg1kufp9z5szMncu57zsaPXjned77XnPOISIiwRVKdwNERGR2KdCLiAScAr2ISMAp0IuIBJwCvYhIwCnQi4gEnAK9iEjAKdCLiAScAr2ISMDlprsBAOFw2FVVVaW7GSIiWeXll19ud86VX22/qwZ6M1sFfA9YBsSAXc65r5lZKfBDoAo4AfyOc67LzAz4GnAf0Af8rnPuwETHqKqqYv/+/VdrioiIJDCzk5PZbzKpmwjwGefcemAL8IiZbQA+DzztnKsDnvafA7wXqPNvO4FvTLHtIiKSRFcN9M65s/ERuXPuPHAYWAncD3zX3+27wG/5j+8Hvuc8LwIlZrY86S0XEZFJmVIx1syqgBuAvcBS59xZ8P4YAEv83VYCpxP+WYu/bVb902tnuevLzxCNaTVOEclc/3joLHd/+RliKYxVky7GmtkC4CfAo865Xi8VP/auY2y7okdmthMvtUNlZeVkmzGuP/3xq1wcitI/HGVBQUbUmEVErvDvfnSQwUiMoWiMwlBOSo45qRG9meXhBfnvO+d+6m8+F0/J+Pet/vYWYFXCP68Azox+T+fcLudcvXOuvrz8qkXjSdP6+iKSydIRoa4a6P1ZNN8CDjvnvpzw0hPAQ/7jh4DHE7Z/wjxbgJ54imc2TfALQ0RkTptMjqMB+Djwmpkd9Lf9GfBF4Edm9jBwCnjQf+3neFMrm/CmV/5eUls8jvhIXuN5EckGsRRmH64a6J1zzzN23h3gzjH2d8AjM2zXtClzIyLZIJXzRoK3BIICvYhksPioOZUj+sAF+lR+eCIi0+ViqTuWAr2ISBpoRD8DOl9KRDJZPEQp0E9DfHql5tGLSDZIZaQKTKCP04heRLKBRvQzoBy9iGSDVIYqBXoRkTTQiH4GFOdFJBvohKkZ0IheRLJBKpcpDmCgT3cLREQm4Mco5einIR2nFYuITJdL4QTLwAT6+EemefQikg2Uo58BpW5EJBto1s0MaEAvIhnNzzOnMvsQuECvHL2IZAOlbmZAgV5EMpofopS6mQHFeRHJBjGtRz91ml4pItlE0ytnQLNuRCQb6ISpGdCIXkSygXL0M6ATpkQkG2jWzQwodSMi2UAj+hlI5YpwIiJTFS/C6oSpGVCcF5FsoNTNdKThtGIRkenSrJvpGDnbLL3NEBGZDOXoZ0DTK0UkGyjQz4DCvIhkMvPzzErdzIBG9CKSDTSinwEVY0Ukk8WnV2rWzXT4s25SuSKciMh0aR79DCh1IyLZQDn6GdD0ShHJBsrRz4By9CKSDZSjnwGN6EUkG2TUiN7Mvm1mrWb2esK2UjN7yswa/fvF/nYzs6+bWZOZHTKzTbPZ+LEoRy8i2SDTirHfAe4dte3zwNPOuTrgaf85wHuBOv+2E/hGcpo5eQr0IpLJXBqWa7lqoHfOPQt0jtp8P/Bd//F3gd9K2P4953kRKDGz5clq7ETi14xVnBeRbJANs26WOufOAvj3S/ztK4HTCfu1+NtSRiN6EckGGZWjnyIbY9uYvTGznWa238z2t7W1zfjA8YOoGCsi2SAbAv25eErGv2/1t7cAqxL2qwDOjPUGzrldzrl651x9eXn5NJtxJY3oRSQbZEPq5gngIf/xQ8DjCds/4c++2QL0xFM8KaM4LyIZzOLLtaQw0udebQcz+wGwAwibWQvwl8AXgR+Z2cPAKeBBf/efA/cBTUAf8Huz0OYJaUQvIpksHbNurhronXMfGeelO8fY1wGPzLRRM6EcvYhkg2zI0WeceBVYI3oRyQpZkKPPWFrrRkSygUb0M6DUjYhkg4w6MzbbKHUjItlAI/oZ0IheRLJBpi1qllWUoxeRTJaOs/gDF+iVuhGRbKDUzTSYf7qZUjcikg2yYQmEjBNP2WhELyLZQCP6GVCcF5FsoBH9DMSUuxGRLKAR/QwozItIJru0XEvqjhm4QK8cvYhkskvTKzWinzJdYUpEsolOmJoO/zPTCVMikg1UjJ2GmKZXikgWUY5+GuIfmlI3IpLJ0nHOT4ACvUb0IpL54hFKOfppcO7yexGRTOOcS8s1YwMT6EdG9MrdiEiGShyIKnUzDZpeKSKZLjG4pzJUBSbQK0cvIpkuphH99CXmvTSPXkQy1WUjeuXop+byvFf62iEiMpHLYlUKg1UgAn3iX0mlbkQkU10eq1J33IAE+kuPFeZFJFMlxifl6Kfo8ryXQr2IZKZ0xapABPpEsVi6WyAiMjaXEJ80vXKKlKMXkWyQrlgVkEB/6XHfUDR9DRERmUDf8KX4pGLsFCX+ZTzb05/GloiIjO9s96X4pBz9FCXmvc72DKSvISIiEziTEJ9SWU8MRKAfjHg/h0qK8jjXO0AkqoqsiGSe+Ih+cVEeA5HUpZkDEehPdPQBsHlNKTEHrecH09wiEZErne0ZYGFBLtcsW8RJP26lQiAC/fH2CwDctnYJAEfeOZ/O5oiIjOnIO+epKC2iunw+x9oupCxPH4hAf6ztIvm5IT5ww0oWFeby+MG3090kEZHLnO3p58XjHdy9YSnV5QvoHYjQeXEoJceelUBvZvea2REzazKzz8/GMRLtO9FJTfkC5uXn8GD9Kh5/9Qz/8OqZ2T6siMiknB8Y5rM/PkReTogHNlWwbulCAPad6ErJ8XOT/YZmlgP8LXA30ALsM7MnnHNvJvtYAHua2jlwqpv/8BvrAfjMPWs5cKqLf/ODV/jOnhPcs2Ep9VWlrAnPZ3FRHmY2G80QERkxHI3xdlc/b71znhea2/m/B89wfmCYL37weirLilhRUsjy4kK++Wwzd65fQl7O7CZXkh7ogZuBJufcMQAzewy4H0h6oP+7547xX35+mOrwfD500yoAivJzeWznFh576TT/68WT/Nd/emtk/4UFuZQuyKd4Xt7IbX5+LgV5IQpyQxTk5nj3ed7jvJwQOSEImZETSriZEfLvc3L8+5CN7Bcy8P6eGGZg3ufg34P52xn1/LLHXPkeIbOx//2o9x7LWH/fxv2TN84LY733eH83x9o83h/Z8doxdpuncMBx32O8fSceBIzOp7rLXhu98+ino/7tBAvxTeU4o9+XCd93/H975Wvjt+lqaeXL+zZ+v696nAneN75HzHnn0URj3jUpYu7SNhd/HPPuXcJrI7fYpX8z+nXnIBJzDEViDEaiDEVi3i0aY9B/POi/dn4gQk//MD19w/T0D9N2YZCof0bUvLwcbltbziO313JdRTEAuTkhPnfvOv7tD1/l7547zid31Ez8oc7QbAT6lcDphOctwOZZOA71VaV8ckcND22rYmFh3sj2gtwcHtpWxUPbqmjtHeBQSw8nO/s43dlHV98Q3f5/jJaufvqHogxEogwOxxiIRHVxcRGZkBkU5IbIzwmR7w8OFxbmUjwvj9VlRRTPy2NZcSGVpUVUly/gupXF5OdeOWL/wA0VhMy4Z8OyWW/zbAT6sYZFV4RPM9sJ7ASorKyc1oE2riph46qSCfdZsqiQuzYUTur9nHNEYo7BSIyB4SiRqCPqHLGYN2KI+iOH+C026nk0YYTg/PdzAM4b2TjnjUoSX3P+Dpe2J+zLpRFOfHssxvjvPW6/xtg2zt7j/aEbc/M4O4+1ddz3TcZ7jL15SjMaxvuMRv+CmGjQP/oXwehdR/9bu+y1iY9jE7w40XGu1n6b8LXxd77ymFNo01WOM5XPOMfiv6Av/ZqO//INmfk37/WQQSiU+PzS/iEzQiFG7e/9Ws/P9X7x5/u33JAlLQV8/8aVSXmfq5mNQN8CrEp4XgFcURl1zu0CdgHU19dnxDjazMjLMfJyQiwomI2PRkQk9WajArAPqDOzNWaWD3wYeGIWjiMiIpOQ9GGrcy5iZn8M/ALIAb7tnHsj2ccREZHJsUy4IpOZtQEnp/nPw0B7EpuT6dTf4JpLfQX1NxlWO+fKr7ZTRgT6mTCz/c65+nS3I1XU3+CaS30F9TeVArEEgoiIjE+BXkQk4IIQ6HeluwEppv4G11zqK6i/KZP1OXoREZlYEEb0IiIyAQV6EZGAy+pAn+p171PBzL5tZq1m9nrCtlIze8rMGv37xf52M7Ov+/0/ZGab0tfyqTOzVWb2KzM7bGZvmNmn/e1B7W+hmb1kZq/6/f1rf/saM9vr9/eH/hnlmFmB/7zJf70qne2fDjPLMbNXzOxJ/3mQ+3rCzF4zs4Nmtt/flhHf5awN9Anr3r8X2AB8xMw2pLdVSfEd4N5R2z4PPO2cqwOe9p+D1/c6/7YT+EaK2pgsEeAzzrn1wBbgEf+/YVD7Owjc4Zx7N7ARuNfMtgBfAr7i97cLeNjf/2GgyzlXC3zF3y/bfBo4nPA8yH0FuN05tzFhvnxmfJedv25ztt2ArcAvEp5/AfhCutuVpL5VAa8nPD8CLPcfLweO+I+/CXxkrP2y8QY8jnfBmsD3FygCDuAt4d0O5PrbR77XeMuIbPUf5/r7WbrbPoU+VuAFtzuAJ/EWtQxkX/12nwDCo7ZlxHc5a0f0jL3ufWrW/Ey9pc65swD+/RJ/e2A+A/+n+g3AXgLcXz+VcRBoBZ4CmoFu51zE3yWxTyP99V/vAcpS2+IZ+SrwOSDmPy8juH0Fb8Xsfzazl/1l2CFDvsvZvBbvpNa9D7hAfAZmtgD4CfCoc653grW+s76/zrkosNHMSoCfAevH2s2/z9r+mtn7gFbn3MtmtiO+eYxds76vCRqcc2fMbAnwlJm9NcG+Ke1vNo/oJ7XufUCcM7PlAP59q7896z8DM8vDC/Lfd8791N8c2P7GOee6gV/j1SZKzCw+6Ers00h//deLgc7UtnTaGoD3m9kJ4DG89M1XCWZfAXDOnfHvW/H+iN9MhnyXsznQz6V1758AHvIfP4SXy45v/4Rfwd8C9MR/JmYD84bu3wIOO+e+nPBSUPtb7o/kMbN5wF14hcpfAQ/4u43ub/xzeAD4pfMTupnOOfcF51yFc64K7//NXzrnPkYA+wpgZvPNbGH8MXAP8DqZ8l1OdwFjhsWP+4CjeHnOP093e5LUpx8AZ4FhvL/6D+PlKp8GGv37Un9fw5t51Ay8BtSnu/1T7Ot2vJ+rh4CD/u2+APf3euAVv7+vA3/hb68GXgKagB8DBf72Qv95k/96dbr7MM1+7wCeDHJf/X696t/eiMejTPkuawkEEZGAy+bUjYiITIICvYhIwCnQi4gEXEbMow+Hw66qqirdzRARySovv/xyu5vENWMzItBXVVWxf//+dDdDRCSrmNnJyeyn1I2ISMAp0IuIpFBT63kGI9GUHlOBXkQkRZ5rbOOuLz/L/37xVEqPq0AvIpICvQPDPPrYQQAGhjWiFxEJlIHhKH/yg1fouDgEQGFeTkqPr0AvIjKLznT38/Fv7eWZo2382X3XAJDqpWcU6EVEZkEkGuP7e0/ynq88yxtnevn6h2/go5tXA5DqJcYyYh69iEhQRKIxnjx0lq893cjx9otsXlPK3zzwbirLiugb8i6uFUtxpFegFxFJgtbeAR7bd5q/33uKd3oHuGbZQnZ9/Ebu3rCU+FXTQv59TCN6EZHscHEwwlNvnuPxg2/zXGM7kZjjlrow//H+a7lr/VJCocuvGBi/SqZG9CIiGax/KMpzjW38w6GzPPXmOwwMx1hRXMjD29fwoZtWUV2+YNx/Gx/Rp7oYq0AvInIV7RcG+eXhVp46fI7nGtsYGI6xuCiPD26q4P6NK6lfvfiK0ftYlLoREckgTa0X+JfD53jqzXMcONWFc7CiuJAP1a/izvVL2VpTRl7O1CYuhpS6mZlYzHF+MELxvLx0N0VEstDAcJQXj3Xw7NF2fn2klWPtFwF418pFPHrnWu7asIQNyxeNFFanwzSin5nP/p9D/ORAC03/+b3kTvGvrIjMPc45mtsu8szRNp452sbeYx0MRmIU5IbYXF3G7zVUcef6pawomZfU44ZMOfpp+8mBFgAiMUduas8uFpEscX5gmN1NHTzb2MYzR9p4u7sfgJry+Xxs82puW1fO5jWls7pEQchMqZvpaO0dGHmc6jPORCRzxWKON8/2jozaD5zsIhJzLCjIpaG2jEdur+XWtWEqFhelrE1eoE/Z4YCABPqDp7tHHqf6L6WIZJZzvQPsbmrn+cZ2nm1sp/3CIADXrljEzluruW1tOZtWL55yITVZzFSMnZZDLT0jjxXoReaWnv5hXjzWwZ6mdnY3d9DUegGA0vn53FIX5ra15dxSV075woI0t9QTMgvGWjdm9m3gfUCrc+5ds3GMRIfP9o48TvVPIhFJrYHhKAdOdrG7uZ3nmzp4raWbmIN5eTncvKaUD9WvYlttGeuXLZrU3PZUC5mXUkql2RrRfwf478D3Zun9L9PS1T/yONXVbBGZXdGY440zPTzf1M6epg72nehkMBIjJ2RsXFXCH99RR0NNGTdULiY/N/Nn3AUmR++ce9bMqmbjvcc4Fqe7+sgJGdGY04heJMs55zjWftFLxTR18MKxDnr6hwG4ZtlCPrZ5Ndvryrh5TRkLCrIv+zyncvRmthPYCVBZWTnt9+nqG6ZvKMqa8HyOt19Ujl4kC7X2DnipmMYO9jS3c7bHm0m3smQe9167jG21ZWyrCWdMnn0mQiGbO/PonXO7gF0A9fX10+716c4+ACpLixToRbJE78AwLzZ3sKe5g+eb2kcKqIuL8thWE2ZbbRnba8NUlhbN6EzUTBSY1E0qtZ33pk4tW1QIaB69SCYaGI5y4FQXe5q8wH4ooYB605pSfqe+gm01YTYsz8wCajKF5lLqJlm6/dxd6YJ8QNMrRTJBNOZ480yvV0Btbuel46MKqLfX0lAbZmNlCQVz7FR2C8qI3sx+AOwAwmbWAvylc+5bs3GseJFmcZG3mJmKsSKp55zjePtFdjd3sLux/bIC6rqlC/no5kq214a5eU0pCwvn9sKDgVnrxjn3kdl437H09A1hxsiqlamenyoyV7X2Dozk2Pc0tXMmoYD6nmuX0lAbZmtNGUsWFqa5pZlFa91MQ0//MAsLcskNefNnlbkRmR29A8PsPdbJ7qZ2dje10+gXUEuK8thWU8anasJsrw2zuix4BdRkUjF2Grr7hykpyseP88rRiyTJYCTKgZPd7Glu9wuoPURjjsK8EDdVlfLAjRU01M6NAmoyzal59MnS3TdMSVFewiW6FOhFpiO+0uPz/oh934lOBoa9Auq7K4r51I4attWE2bR67hVQkykwa92kUk//MMXz8tJ25RaRbOWc40RH30gq5oVjHXT3eQXUtUsX8OGbvALq5moVUJNJ0yunobd/mJWL541ci1Fr3YiMr/X8AHuaOtjd1M6e5o6RC2+sKC7k7vVeAXVbTRlLFqmAOluUo5+Gi0MRFuTnpu3q6iKZ7Hy8gNrsjdqPnvMKqMXzvALqH+2oYXttmCoVUFNGOfpp6BuKMi8/J21XVxfJJIORKK+c6h5Jx7w6qoD625sqaKgJs2HFInJUQE0LL0evQD8lA8NeoDcVY2UOihdQd/sX3XjpeAcDwzFCBu9eVcInb6uhoVYF1EwSMiMWS+0xszrQD0djDEcd8/JyRlI3ivMSZM45Tnb0jaRiXmjuoMsvoNYt8QqoDX4BdZEKqBlJqZsp6h+OAlCk1I0EWNv5Qfb4gX1306UC6vLiQu5cv5QGfwnfpSqgZoXArHWTKv1DXqAvTBjRqxgr2e7CYIS9xzrY7c+OOXLuPOAVULdWewXUhpoy1oTnq4CahbxBqUb0kxYP9EX5OZhG9JKlhiIxXjnVNZJnf/V0N5GYoyDXK6D+1g0raagt49oVxSqgBoCmV05RPHUzL+9SMVbz6CXTxWKOw+/0jqRiXjreSf9wlJDB9RUl/OFt1V4BtXIxhXkqoAaNTpiaoj5/RH/59Mo0NkhkDM45TnX2jaRiXjjWQefFIQBqlyzgQzetYltNGZury0ZWYZXgUo5+igYSRvQR/5PTgF4yQbyAGr+iUryAumxRIbevW0JDbRkNtSqgzkWBWY8+VfpGcvS5nB/0ppgpRy/pcGEwwkvHLxVQ33rHK6AuKsxla00Zf3RbNdtqw1SrgDrnBWY9ejO7F/gakAP8nXPui7NxnJEcfX6Ii0M6YUpSZygS4+Dp7pGLbhz0C6j5uSFuqlrM5+5dR0NNmHetVAFVLheIE6bMLAf4W+BuoAXYZ2ZPOOfeTPax+ociwOXTKxXnZTbEC6jxVMy+E530DXkF1OsqSth5azXba8NsWq0CqkwsKCdM3Qw0OeeOAZjZY8D9wCwE+kupG50wJcl2yj8D9Xn/DNR4AbWmfD4P3ljBttowW1RAlSkKmRFNcTV2NgL9SuB0wvMWYPPoncxsJ7AToLKycloH6h/2fv8kTq/UrBuZrvYLg+xp7mBPkxfcW7ouFVB3rCunoSZMQ22YZcUqoMr0hUIwHM3+QD9WQvKKXjnndgG7AOrr66fV63+1pZL7rltGYV5II3qZsouDEV467l0D9fmEAurCwly2Vpex89ZqttWEqSlXAVWSJyjF2BZgVcLzCuDMLByHhYV5I1e+CemEKbmK4ahfQG1sZ09zO6+curyA+tn3rKOhNsx1KqDKLArKPPp9QJ2ZrQHeBj4MfHQWjnOZkbVuUlzNlswVizneeuf8yIJge48nFFBXFvMHfgH1RhVQJYUCMY/eORcxsz8GfoE3vfLbzrk3kn2c0bTWjQCc7uwbubj1C80ddPgF1Ory+TxwYwXbasJsrS6juEgFVEmPwKx145z7OfDz2Xjv8Wj1yrmpI15A9WfHnO70CqhLFxVw29pyttWGaagtY3nxvDS3VMSjtW5mIBTy7pWjD7aLgxFeOtHJ7kZvpcfDZ3sBr4C6pbqM399eTUNtGTXlC1RAlYwUlBx9WmhEH0zxAuruJm/dmAOnurwCak6Ier+Auq2mjOtWFpObE0p3c0WuKhA5+nTR9MpgiMUcR86dH7m49UvHO7k4FMX8Aurv3+IVUOurVECV7BSU6ZVpoYuDZ6/TnX0jF914obmd9gt+ATU8n9/eVEFDbRlbqssoKcpPc0tFZi4wxdh00Fo32aPz4pA/5dFb6fFUZx8A5QsLuKWunG013hK+K0pUQJXgCcpaN2mh1E3m6hu6dAbq7qYO3owXUAty2Vxdxr9uqKKhNkztEhVQJfhCZikfkAYo0KsYmymGozEOtXTzfGMHu5vbeeVUF8NRr4B64+rF/Ok9a0fOQFUBVeYaTa+cAZ0wlT7OxQuoXipm77GOkQLqu1YU87A/5bF+dSnz8lVAlblNxdgZ0Fo3qfV2d78/l91Lx7RfGARgTXg+H9i0koaaMFtrVEAVGc38C49cGIzwqe8f4He3reaOa5bO6jEDF+iVupkd3X1DvNDsXXRjT3MHx9svAhBeUMD22jL/DNQwK1VAFZlQfB5918Uhnj3axm9ev3zWjxmgQO/dK3WTHAPDUfaf6BpZN+b1Mz04B/Pzc9hSXcbHt6ymoTbM2qUqoIpMRXx6ZU+/d53rVFy4JjCBXhcemZlozPHa2z0jJyrtP9nFUCRGXo5xw6rFPHrnWrbXlXF9RQl5KqCKTFso5A1Iu/u8QJ+K9GZgAn18RK8c/eQ45zjWftG76EZjOy8e66B3wLsG7/rli/jEltU01IW5uaqU+QWB+ZqIpF18rZvufu/EQI3op+DSevQK9ONp7R3wroHa6M2Oead3AICVJfO477rlbKsNs62mjPCCgjS3VCS44jn6eOqmJAVLZgcv0CvOj+gdGGbvsc6RdExj6wUAFhflsc2//mlDbRmVpUXKs4ukSHx6ZTx1oxH9FJifNp7LxdjBSJRXTnWPBPZXW3qIxhyFeSFuqirlgRsraKgNs2H5IkK6VJ5IWhiMFGMLckMpWZwvqYHezB4E/gpYD9zsnNufzPefyFxc6yYWcxx+p9e/uHUH+4530j/sXSrv3atK+ORtNTTUhtm0uoSCXJ2oJJIJzB/R9/QNp2Q0D8kf0b8O/DbwzSS/71XNlemVpzr6vDy7f6m8Tv9SebVLFvChm1bRUBtmc3Upiwp1qTyRTBQyAwddfUMpyc9DkgO9c+4wkJZ8b1Bz9F0Xh/yzT6+8VN6OdeVsrw2zrSbMsuLCNLdURCYjvtbN2Z4BlqXoEpeBydHHZfuIfjAS5eWTXTzf6AX21972TlRaWJjL1pFL5YWpKZ+vAqpIFgqFvOmVLV19XFdRnJJjTjnQm9m/AMvGeOnPnXOPT+F9dgI7ASorK6fajCtk61o3zjmOnrvAc41tPN/Uzt5jXp49N2TcUFnCo3eu5Za1Ya7XSo8igWAG/cNR+oejrFpclJJjTjnQO+fuSsaBnXO7gF0A9fX1M47Ol06Ymuk7zb7W8wPsbmrnuUbvZKXW896CYNXl8/md+gq215WzpbqUhcqziwROKOGXeMVipW6mJJNz9APDUV463slzjW0819jOW++cB7z57A21YW6pC7O9rlwLgonMAYkzmytLM3REPxEz+wDw34By4B/N7KBz7j3JPMb4x/buMyFHH4s53jzby/NN7TzX2Ma+E966Mfk5IeqrFvO5e9dxS205167QfHaRuSY+KM0NGeuWLUzJMZM96+ZnwM+S+Z6TZWaYpS9H33lxiOca23jmSBvPNraNXOB63dKFfGLLarbXhbl5TSlF+YH5ESUi0xCfRHHN8oUpOVkKApS6gdReXT0ac7za0s0zR9r49dE2DrV045yXjrl1bTm31pVzS12YJYs07VFELhmKxAC4YdXilB0zYIF+dlM3becHefZoG88cbeO5xja6+oYxg42rvNkxt60r57qVxeQoHSMi43jhWAcAd6xfkrJjBirQW5JH9M55ufan3jzH04eixQlUAAAF7ElEQVRbee3tHgDCC/K5/Zol7Fi3hFtqwyyer8vlicjkfPrOWv7Tk4dpqAmn7JiBCvShJOTohyIx9h7v4Kk3z/Evb57jTM8AZrCpcjGffc86bltbrkXBRGTa7rhm6axfI3a0gAX66V1dfWA4yq/eauUfXzvLM0faOD8YoTAvxK115Tx691ruuGaJ1mgXkawVwEA/uX2HIjGeb2rjH149yz+/8Q4Xh6KEF+TzG9cv5+4NS2moDaesIi4iMpsCFehtEsXYptYL/P3eU/z0lRa6/WVC379xBb95/Qo2V5epkCoigROoQB8yG3MJBOcczxxt438808yLxzrJyzHuuXYZH9y0ku215eTnag0ZEQmugAX6K0f0u5va+dL/e4tDLT2sKC7k3997DQ/WVyjnLiJzRsAC/aVibHffEH/5xBs8fvAMq0rn8aUPXscHbqjQ6F1E5pxABfr4PPqj587z8Hf38U7PAH9yZx2f2lGjwqqIzFmBCvQhg+bWC3z0f75IyIwf/eFWbqhM3WnGIiKZKFB5jEjMsfd4J87B3//BFgV5ERECFujjF8r+mwevp3bJgjS3RkQkMwQqdXPj6sUU5IZSfnqxiEgmC1Sg//EfbkXXyxYRuVygAr0WGhMRuVKgcvQiInIlBXoRkYCzdF1j9bJGmLUBJ6f5z8NAexKbk+nU3+CaS30F9TcZVjvnyq+2U0YE+pkws/3Oufp0tyNV1N/gmkt9BfU3lZS6EREJOAV6EZGAC0Kg35XuBqSY+htcc6mvoP6mTNbn6EVEZGJBGNGLiMgEsjrQm9m9ZnbEzJrM7PPpbk8ymNm3zazVzF5P2FZqZk+ZWaN/v9jfbmb2db//h8xsU/paPnVmtsrMfmVmh83sDTP7tL89qP0tNLOXzOxVv79/7W9fY2Z7/f7+0Mzy/e0F/vMm//WqdLZ/Oswsx8xeMbMn/edB7usJM3vNzA6a2X5/W0Z8l7M20JtZDvC3wHuBDcBHzGxDeluVFN8B7h217fPA0865OuBp/zl4fa/zbzuBb6SojckSAT7jnFsPbAEe8f8bBrW/g8Adzrl3AxuBe81sC/Al4Ct+f7uAh/39Hwa6nHO1wFf8/bLNp4HDCc+D3FeA251zGxOmUWbGd9k5l5U3YCvwi4TnXwC+kO52JalvVcDrCc+PAMv9x8uBI/7jbwIfGWu/bLwBjwN3z4X+AkXAAWAz3kk0uf72ke818Atgq/8419/P0t32KfSxAi+43QE8CVhQ++q3+wQQHrUtI77LWTuiB1YCpxOet/jbgmipc+4sgH+/xN8emM/A/6l+A7CXAPfXT2UcBFqBp4BmoNs5F/F3SezTSH/913uAstS2eEa+CnwOiPnPywhuXwEc8M9m9rKZ7fS3ZcR3OZtXrxxrqcq5NoUoEJ+BmS0AfgI86pzrtfHXms76/jrnosBGMysBfgasH2s3/z5r+2tm7wNanXMvm9mO+OYxds36viZocM6dMbMlwFNm9tYE+6a0v9k8om8BViU8rwDOpKkts+2cmS0H8O9b/e1Z/xmYWR5ekP++c+6n/ubA9jfOOdcN/BqvNlFiZvFBV2KfRvrrv14MdKa2pdPWALzfzE4Aj+Glb75KMPsKgHPujH/fivdH/GYy5LuczYF+H1DnV/HzgQ8DT6S5TbPlCeAh//FDeLns+PZP+BX8LUBP/GdiNjBv6P4t4LBz7ssJLwW1v+X+SB4zmwfchVeo/BXwgL/b6P7GP4cHgF86P6Gb6ZxzX3DOVTjnqvD+3/ylc+5jBLCvAGY238wWxh8D9wCvkynf5XQXMGZY/LgPOIqX5/zzdLcnSX36AXAWGMb7q/8wXq7yaaDRvy/19zW8mUfNwGtAfbrbP8W+bsf7uXoIOOjf7gtwf68HXvH7+zrwF/72auAloAn4MVDgby/0nzf5r1enuw/T7PcO4Mkg99Xv16v+7Y14PMqU77LOjBURCbhsTt2IiMgkKNCLiAScAr2ISMAp0IuIBJwCvYhIwCnQi4gEnAK9iEjAKdCLiATc/wcaHs17oJKj3gAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "from scipy.fftpack import fft\n",
    "\n",
    "# generate sinusoid\n",
    "N = 511\n",
    "A = 0.8\n",
    "f0 = 440\n",
    "fs = 44100\n",
    "phi = 1.0\n",
    "x = generate_sinusoid(N, A, f0, fs, phi)\n",
    "\n",
    "# fft is\n",
    "X = fft(x)\n",
    "mX = np.abs(X)   # magnitude\n",
    "pX = np.angle(X) # phase\n",
    "\n",
    "# plot the magnitude and phase\n",
    "plt.subplot(2,1,1)\n",
    "plt.plot(mX)\n",
    "\n",
    "plt.subplot(2,1,2)\n",
    "plt.plot(pX)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 自己实现DFT\n",
    "\n",
    "自己实现DFT的关键就是构造出$S$，有两种方式：\n",
    "\n",
    "+ 我们可以单独构建一个$S_k$，然后计算 $sum(S_k * x[n])$，得到 $X[k]$\n",
    "+ 或者，构建整个矩阵$S$，一次性计算 $S*x[n]$，得到$X$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def generate_complex_sinusoid(k, N):\n",
    "    '''\n",
    "    k (int): frequency index\n",
    "    N (int): length of complex sinusoid in samples\n",
    "    \n",
    "    returns\n",
    "    c_sin (numpy array): the generated complex sinusoid (length N)\n",
    "    '''\n",
    "    \n",
    "    n = np.arange(N)\n",
    "    \n",
    "    c_sin = np.exp(1j * 2 * np.pi * k * n / N)\n",
    "    \n",
    "    return np.conjugate(c_sin)\n",
    "\n",
    "def generate_complex_sinusoid_matrix(N):\n",
    "    '''\n",
    "    N (int): length of complex sinusoid in samples\n",
    "    \n",
    "    returns\n",
    "    c_sin_matrix (numpy array): the generated complex sinusoid (length N)\n",
    "    '''\n",
    "    \n",
    "    n = np.arange(N)\n",
    "    n = np.expand_dims(n, axis=1)      # 扩充维度，将1D向量，转为2D矩阵，方便后面的矩阵相乘\n",
    "    \n",
    "    k = n\n",
    "    \n",
    "    m = n.T * k / N                    # [N,1] * [1, N] = [N,N]\n",
    "    \n",
    "    S = np.exp(1j * 2 * np.pi * m)     # 计算矩阵 S\n",
    "    \n",
    "    return np.conjugate(S)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD8CAYAAAB5Pm/hAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3Xl0XOd53/HvM1gJLgCBAVcQBLGQImXJFAWJCyiJWi0rjhXHUrzVVlIlTGylsVrHrp30ZGlPW/vkHG9tjms29rHdOpbt2q4UxamjyLYWUqJIURS1UCQArhApYgdIYp2Zt3/cO+AQBEAsg1kufp9z5szMncu57zsaPXjned77XnPOISIiwRVKdwNERGR2KdCLiAScAr2ISMAp0IuIBJwCvYhIwCnQi4gEnAK9iEjAKdCLiAScAr2ISMDlprsBAOFw2FVVVaW7GSIiWeXll19ud86VX22/qwZ6M1sFfA9YBsSAXc65r5lZKfBDoAo4AfyOc67LzAz4GnAf0Af8rnPuwETHqKqqYv/+/VdrioiIJDCzk5PZbzKpmwjwGefcemAL8IiZbQA+DzztnKsDnvafA7wXqPNvO4FvTLHtIiKSRFcN9M65s/ERuXPuPHAYWAncD3zX3+27wG/5j+8Hvuc8LwIlZrY86S0XEZFJmVIx1syqgBuAvcBS59xZ8P4YAEv83VYCpxP+WYu/bVb902tnuevLzxCNaTVOEclc/3joLHd/+RliKYxVky7GmtkC4CfAo865Xi8VP/auY2y7okdmthMvtUNlZeVkmzGuP/3xq1wcitI/HGVBQUbUmEVErvDvfnSQwUiMoWiMwlBOSo45qRG9meXhBfnvO+d+6m8+F0/J+Pet/vYWYFXCP68Azox+T+fcLudcvXOuvrz8qkXjSdP6+iKSydIRoa4a6P1ZNN8CDjvnvpzw0hPAQ/7jh4DHE7Z/wjxbgJ54imc2TfALQ0RkTptMjqMB+Djwmpkd9Lf9GfBF4Edm9jBwCnjQf+3neFMrm/CmV/5eUls8jvhIXuN5EckGsRRmH64a6J1zzzN23h3gzjH2d8AjM2zXtClzIyLZIJXzRoK3BIICvYhksPioOZUj+sAF+lR+eCIi0+ViqTuWAr2ISBpoRD8DOl9KRDJZPEQp0E9DfHql5tGLSDZIZaQKTKCP04heRLKBRvQzoBy9iGSDVIYqBXoRkTTQiH4GFOdFJBvohKkZ0IheRLJBKpcpDmCgT3cLREQm4Mco5einIR2nFYuITJdL4QTLwAT6+EemefQikg2Uo58BpW5EJBto1s0MaEAvIhnNzzOnMvsQuECvHL2IZAOlbmZAgV5EMpofopS6mQHFeRHJBjGtRz91ml4pItlE0ytnQLNuRCQb6ISpGdCIXkSygXL0M6ATpkQkG2jWzQwodSMi2UAj+hlI5YpwIiJTFS/C6oSpGVCcF5FsoNTNdKThtGIRkenSrJvpGDnbLL3NEBGZDOXoZ0DTK0UkGyjQz4DCvIhkMvPzzErdzIBG9CKSDTSinwEVY0Ukk8WnV2rWzXT4s25SuSKciMh0aR79DCh1IyLZQDn6GdD0ShHJBsrRz4By9CKSDZSjnwGN6EUkG2TUiN7Mvm1mrWb2esK2UjN7yswa/fvF/nYzs6+bWZOZHTKzTbPZ+LEoRy8i2SDTirHfAe4dte3zwNPOuTrgaf85wHuBOv+2E/hGcpo5eQr0IpLJXBqWa7lqoHfOPQt0jtp8P/Bd//F3gd9K2P4953kRKDGz5clq7ETi14xVnBeRbJANs26WOufOAvj3S/ztK4HTCfu1+NtSRiN6EckGGZWjnyIbY9uYvTGznWa238z2t7W1zfjA8YOoGCsi2SAbAv25eErGv2/1t7cAqxL2qwDOjPUGzrldzrl651x9eXn5NJtxJY3oRSQbZEPq5gngIf/xQ8DjCds/4c++2QL0xFM8KaM4LyIZzOLLtaQw0udebQcz+wGwAwibWQvwl8AXgR+Z2cPAKeBBf/efA/cBTUAf8Huz0OYJaUQvIpksHbNurhronXMfGeelO8fY1wGPzLRRM6EcvYhkg2zI0WeceBVYI3oRyQpZkKPPWFrrRkSygUb0M6DUjYhkg4w6MzbbKHUjItlAI/oZ0IheRLJBpi1qllWUoxeRTJaOs/gDF+iVuhGRbKDUzTSYf7qZUjcikg2yYQmEjBNP2WhELyLZQCP6GVCcF5FsoBH9DMSUuxGRLKAR/QwozItIJru0XEvqjhm4QK8cvYhkskvTKzWinzJdYUpEsolOmJoO/zPTCVMikg1UjJ2GmKZXikgWUY5+GuIfmlI3IpLJ0nHOT4ACvUb0IpL54hFKOfppcO7yexGRTOOcS8s1YwMT6EdG9MrdiEiGShyIKnUzDZpeKSKZLjG4pzJUBSbQK0cvIpkuphH99CXmvTSPXkQy1WUjeuXop+byvFf62iEiMpHLYlUKg1UgAn3iX0mlbkQkU10eq1J33IAE+kuPFeZFJFMlxifl6Kfo8ryXQr2IZKZ0xapABPpEsVi6WyAiMjaXEJ80vXKKlKMXkWyQrlgVkEB/6XHfUDR9DRERmUDf8KX4pGLsFCX+ZTzb05/GloiIjO9s96X4pBz9FCXmvc72DKSvISIiEziTEJ9SWU8MRKAfjHg/h0qK8jjXO0AkqoqsiGSe+Ih+cVEeA5HUpZkDEehPdPQBsHlNKTEHrecH09wiEZErne0ZYGFBLtcsW8RJP26lQiAC/fH2CwDctnYJAEfeOZ/O5oiIjOnIO+epKC2iunw+x9oupCxPH4hAf6ztIvm5IT5ww0oWFeby+MG3090kEZHLnO3p58XjHdy9YSnV5QvoHYjQeXEoJceelUBvZvea2REzazKzz8/GMRLtO9FJTfkC5uXn8GD9Kh5/9Qz/8OqZ2T6siMiknB8Y5rM/PkReTogHNlWwbulCAPad6ErJ8XOT/YZmlgP8LXA30ALsM7MnnHNvJvtYAHua2jlwqpv/8BvrAfjMPWs5cKqLf/ODV/jOnhPcs2Ep9VWlrAnPZ3FRHmY2G80QERkxHI3xdlc/b71znhea2/m/B89wfmCYL37weirLilhRUsjy4kK++Wwzd65fQl7O7CZXkh7ogZuBJufcMQAzewy4H0h6oP+7547xX35+mOrwfD500yoAivJzeWznFh576TT/68WT/Nd/emtk/4UFuZQuyKd4Xt7IbX5+LgV5IQpyQxTk5nj3ed7jvJwQOSEImZETSriZEfLvc3L8+5CN7Bcy8P6eGGZg3ufg34P52xn1/LLHXPkeIbOx//2o9x7LWH/fxv2TN84LY733eH83x9o83h/Z8doxdpuncMBx32O8fSceBIzOp7rLXhu98+ino/7tBAvxTeU4o9+XCd93/H975Wvjt+lqaeXL+zZ+v696nAneN75HzHnn0URj3jUpYu7SNhd/HPPuXcJrI7fYpX8z+nXnIBJzDEViDEaiDEVi3i0aY9B/POi/dn4gQk//MD19w/T0D9N2YZCof0bUvLwcbltbziO313JdRTEAuTkhPnfvOv7tD1/l7547zid31Ez8oc7QbAT6lcDphOctwOZZOA71VaV8ckcND22rYmFh3sj2gtwcHtpWxUPbqmjtHeBQSw8nO/s43dlHV98Q3f5/jJaufvqHogxEogwOxxiIRHVxcRGZkBkU5IbIzwmR7w8OFxbmUjwvj9VlRRTPy2NZcSGVpUVUly/gupXF5OdeOWL/wA0VhMy4Z8OyWW/zbAT6sYZFV4RPM9sJ7ASorKyc1oE2riph46qSCfdZsqiQuzYUTur9nHNEYo7BSIyB4SiRqCPqHLGYN2KI+iOH+C026nk0YYTg/PdzAM4b2TjnjUoSX3P+Dpe2J+zLpRFOfHssxvjvPW6/xtg2zt7j/aEbc/M4O4+1ddz3TcZ7jL15SjMaxvuMRv+CmGjQP/oXwehdR/9bu+y1iY9jE7w40XGu1n6b8LXxd77ymFNo01WOM5XPOMfiv6Av/ZqO//INmfk37/WQQSiU+PzS/iEzQiFG7e/9Ws/P9X7x5/u33JAlLQV8/8aVSXmfq5mNQN8CrEp4XgFcURl1zu0CdgHU19dnxDjazMjLMfJyQiwomI2PRkQk9WajArAPqDOzNWaWD3wYeGIWjiMiIpOQ9GGrcy5iZn8M/ALIAb7tnHsj2ccREZHJsUy4IpOZtQEnp/nPw0B7EpuT6dTf4JpLfQX1NxlWO+fKr7ZTRgT6mTCz/c65+nS3I1XU3+CaS30F9TeVArEEgoiIjE+BXkQk4IIQ6HeluwEppv4G11zqK6i/KZP1OXoREZlYEEb0IiIyAQV6EZGAy+pAn+p171PBzL5tZq1m9nrCtlIze8rMGv37xf52M7Ov+/0/ZGab0tfyqTOzVWb2KzM7bGZvmNmn/e1B7W+hmb1kZq/6/f1rf/saM9vr9/eH/hnlmFmB/7zJf70qne2fDjPLMbNXzOxJ/3mQ+3rCzF4zs4Nmtt/flhHf5awN9Anr3r8X2AB8xMw2pLdVSfEd4N5R2z4PPO2cqwOe9p+D1/c6/7YT+EaK2pgsEeAzzrn1wBbgEf+/YVD7Owjc4Zx7N7ARuNfMtgBfAr7i97cLeNjf/2GgyzlXC3zF3y/bfBo4nPA8yH0FuN05tzFhvnxmfJedv25ztt2ArcAvEp5/AfhCutuVpL5VAa8nPD8CLPcfLweO+I+/CXxkrP2y8QY8jnfBmsD3FygCDuAt4d0O5PrbR77XeMuIbPUf5/r7WbrbPoU+VuAFtzuAJ/EWtQxkX/12nwDCo7ZlxHc5a0f0jL3ufWrW/Ey9pc65swD+/RJ/e2A+A/+n+g3AXgLcXz+VcRBoBZ4CmoFu51zE3yWxTyP99V/vAcpS2+IZ+SrwOSDmPy8juH0Fb8Xsfzazl/1l2CFDvsvZvBbvpNa9D7hAfAZmtgD4CfCoc653grW+s76/zrkosNHMSoCfAevH2s2/z9r+mtn7gFbn3MtmtiO+eYxds76vCRqcc2fMbAnwlJm9NcG+Ke1vNo/oJ7XufUCcM7PlAP59q7896z8DM8vDC/Lfd8791N8c2P7GOee6gV/j1SZKzCw+6Ers00h//deLgc7UtnTaGoD3m9kJ4DG89M1XCWZfAXDOnfHvW/H+iN9MhnyXsznQz6V1758AHvIfP4SXy45v/4Rfwd8C9MR/JmYD84bu3wIOO+e+nPBSUPtb7o/kMbN5wF14hcpfAQ/4u43ub/xzeAD4pfMTupnOOfcF51yFc64K7//NXzrnPkYA+wpgZvPNbGH8MXAP8DqZ8l1OdwFjhsWP+4CjeHnOP093e5LUpx8AZ4FhvL/6D+PlKp8GGv37Un9fw5t51Ay8BtSnu/1T7Ot2vJ+rh4CD/u2+APf3euAVv7+vA3/hb68GXgKagB8DBf72Qv95k/96dbr7MM1+7wCeDHJf/X696t/eiMejTPkuawkEEZGAy+bUjYiITIICvYhIwCnQi4gEXEbMow+Hw66qqirdzRARySovv/xyu5vENWMzItBXVVWxf//+dDdDRCSrmNnJyeyn1I2ISMAp0IuIpFBT63kGI9GUHlOBXkQkRZ5rbOOuLz/L/37xVEqPq0AvIpICvQPDPPrYQQAGhjWiFxEJlIHhKH/yg1fouDgEQGFeTkqPr0AvIjKLznT38/Fv7eWZo2382X3XAJDqpWcU6EVEZkEkGuP7e0/ynq88yxtnevn6h2/go5tXA5DqJcYyYh69iEhQRKIxnjx0lq893cjx9otsXlPK3zzwbirLiugb8i6uFUtxpFegFxFJgtbeAR7bd5q/33uKd3oHuGbZQnZ9/Ebu3rCU+FXTQv59TCN6EZHscHEwwlNvnuPxg2/zXGM7kZjjlrow//H+a7lr/VJCocuvGBi/SqZG9CIiGax/KMpzjW38w6GzPPXmOwwMx1hRXMjD29fwoZtWUV2+YNx/Gx/Rp7oYq0AvInIV7RcG+eXhVp46fI7nGtsYGI6xuCiPD26q4P6NK6lfvfiK0ftYlLoREckgTa0X+JfD53jqzXMcONWFc7CiuJAP1a/izvVL2VpTRl7O1CYuhpS6mZlYzHF+MELxvLx0N0VEstDAcJQXj3Xw7NF2fn2klWPtFwF418pFPHrnWu7asIQNyxeNFFanwzSin5nP/p9D/ORAC03/+b3kTvGvrIjMPc45mtsu8szRNp452sbeYx0MRmIU5IbYXF3G7zVUcef6pawomZfU44ZMOfpp+8mBFgAiMUduas8uFpEscX5gmN1NHTzb2MYzR9p4u7sfgJry+Xxs82puW1fO5jWls7pEQchMqZvpaO0dGHmc6jPORCRzxWKON8/2jozaD5zsIhJzLCjIpaG2jEdur+XWtWEqFhelrE1eoE/Z4YCABPqDp7tHHqf6L6WIZJZzvQPsbmrn+cZ2nm1sp/3CIADXrljEzluruW1tOZtWL55yITVZzFSMnZZDLT0jjxXoReaWnv5hXjzWwZ6mdnY3d9DUegGA0vn53FIX5ra15dxSV075woI0t9QTMgvGWjdm9m3gfUCrc+5ds3GMRIfP9o48TvVPIhFJrYHhKAdOdrG7uZ3nmzp4raWbmIN5eTncvKaUD9WvYlttGeuXLZrU3PZUC5mXUkql2RrRfwf478D3Zun9L9PS1T/yONXVbBGZXdGY440zPTzf1M6epg72nehkMBIjJ2RsXFXCH99RR0NNGTdULiY/N/Nn3AUmR++ce9bMqmbjvcc4Fqe7+sgJGdGY04heJMs55zjWftFLxTR18MKxDnr6hwG4ZtlCPrZ5Ndvryrh5TRkLCrIv+zyncvRmthPYCVBZWTnt9+nqG6ZvKMqa8HyOt19Ujl4kC7X2DnipmMYO9jS3c7bHm0m3smQe9167jG21ZWyrCWdMnn0mQiGbO/PonXO7gF0A9fX10+716c4+ACpLixToRbJE78AwLzZ3sKe5g+eb2kcKqIuL8thWE2ZbbRnba8NUlhbN6EzUTBSY1E0qtZ33pk4tW1QIaB69SCYaGI5y4FQXe5q8wH4ooYB605pSfqe+gm01YTYsz8wCajKF5lLqJlm6/dxd6YJ8QNMrRTJBNOZ480yvV0Btbuel46MKqLfX0lAbZmNlCQVz7FR2C8qI3sx+AOwAwmbWAvylc+5bs3GseJFmcZG3mJmKsSKp55zjePtFdjd3sLux/bIC6rqlC/no5kq214a5eU0pCwvn9sKDgVnrxjn3kdl437H09A1hxsiqlamenyoyV7X2Dozk2Pc0tXMmoYD6nmuX0lAbZmtNGUsWFqa5pZlFa91MQ0//MAsLcskNefNnlbkRmR29A8PsPdbJ7qZ2dje10+gXUEuK8thWU8anasJsrw2zuix4BdRkUjF2Grr7hykpyseP88rRiyTJYCTKgZPd7Glu9wuoPURjjsK8EDdVlfLAjRU01M6NAmoyzal59MnS3TdMSVFewiW6FOhFpiO+0uPz/oh934lOBoa9Auq7K4r51I4attWE2bR67hVQkykwa92kUk//MMXz8tJ25RaRbOWc40RH30gq5oVjHXT3eQXUtUsX8OGbvALq5moVUJNJ0yunobd/mJWL541ci1Fr3YiMr/X8AHuaOtjd1M6e5o6RC2+sKC7k7vVeAXVbTRlLFqmAOluUo5+Gi0MRFuTnpu3q6iKZ7Hy8gNrsjdqPnvMKqMXzvALqH+2oYXttmCoVUFNGOfpp6BuKMi8/J21XVxfJJIORKK+c6h5Jx7w6qoD625sqaKgJs2HFInJUQE0LL0evQD8lA8NeoDcVY2UOihdQd/sX3XjpeAcDwzFCBu9eVcInb6uhoVYF1EwSMiMWS+0xszrQD0djDEcd8/JyRlI3ivMSZM45Tnb0jaRiXmjuoMsvoNYt8QqoDX4BdZEKqBlJqZsp6h+OAlCk1I0EWNv5Qfb4gX1306UC6vLiQu5cv5QGfwnfpSqgZoXArHWTKv1DXqAvTBjRqxgr2e7CYIS9xzrY7c+OOXLuPOAVULdWewXUhpoy1oTnq4CahbxBqUb0kxYP9EX5OZhG9JKlhiIxXjnVNZJnf/V0N5GYoyDXK6D+1g0raagt49oVxSqgBoCmV05RPHUzL+9SMVbz6CXTxWKOw+/0jqRiXjreSf9wlJDB9RUl/OFt1V4BtXIxhXkqoAaNTpiaoj5/RH/59Mo0NkhkDM45TnX2jaRiXjjWQefFIQBqlyzgQzetYltNGZury0ZWYZXgUo5+igYSRvQR/5PTgF4yQbyAGr+iUryAumxRIbevW0JDbRkNtSqgzkWBWY8+VfpGcvS5nB/0ppgpRy/pcGEwwkvHLxVQ33rHK6AuKsxla00Zf3RbNdtqw1SrgDrnBWY9ejO7F/gakAP8nXPui7NxnJEcfX6Ii0M6YUpSZygS4+Dp7pGLbhz0C6j5uSFuqlrM5+5dR0NNmHetVAFVLheIE6bMLAf4W+BuoAXYZ2ZPOOfeTPax+ociwOXTKxXnZTbEC6jxVMy+E530DXkF1OsqSth5azXba8NsWq0CqkwsKCdM3Qw0OeeOAZjZY8D9wCwE+kupG50wJcl2yj8D9Xn/DNR4AbWmfD4P3ljBttowW1RAlSkKmRFNcTV2NgL9SuB0wvMWYPPoncxsJ7AToLKycloH6h/2fv8kTq/UrBuZrvYLg+xp7mBPkxfcW7ouFVB3rCunoSZMQ22YZcUqoMr0hUIwHM3+QD9WQvKKXjnndgG7AOrr66fV63+1pZL7rltGYV5II3qZsouDEV467l0D9fmEAurCwly2Vpex89ZqttWEqSlXAVWSJyjF2BZgVcLzCuDMLByHhYV5I1e+CemEKbmK4ahfQG1sZ09zO6+curyA+tn3rKOhNsx1KqDKLArKPPp9QJ2ZrQHeBj4MfHQWjnOZkbVuUlzNlswVizneeuf8yIJge48nFFBXFvMHfgH1RhVQJYUCMY/eORcxsz8GfoE3vfLbzrk3kn2c0bTWjQCc7uwbubj1C80ddPgF1Ory+TxwYwXbasJsrS6juEgFVEmPwKx145z7OfDz2Xjv8Wj1yrmpI15A9WfHnO70CqhLFxVw29pyttWGaagtY3nxvDS3VMSjtW5mIBTy7pWjD7aLgxFeOtHJ7kZvpcfDZ3sBr4C6pbqM399eTUNtGTXlC1RAlYwUlBx9WmhEH0zxAuruJm/dmAOnurwCak6Ier+Auq2mjOtWFpObE0p3c0WuKhA5+nTR9MpgiMUcR86dH7m49UvHO7k4FMX8Aurv3+IVUOurVECV7BSU6ZVpoYuDZ6/TnX0jF914obmd9gt+ATU8n9/eVEFDbRlbqssoKcpPc0tFZi4wxdh00Fo32aPz4pA/5dFb6fFUZx8A5QsLuKWunG013hK+K0pUQJXgCcpaN2mh1E3m6hu6dAbq7qYO3owXUAty2Vxdxr9uqKKhNkztEhVQJfhCZikfkAYo0KsYmymGozEOtXTzfGMHu5vbeeVUF8NRr4B64+rF/Ok9a0fOQFUBVeYaTa+cAZ0wlT7OxQuoXipm77GOkQLqu1YU87A/5bF+dSnz8lVAlblNxdgZ0Fo3qfV2d78/l91Lx7RfGARgTXg+H9i0koaaMFtrVEAVGc38C49cGIzwqe8f4He3reaOa5bO6jEDF+iVupkd3X1DvNDsXXRjT3MHx9svAhBeUMD22jL/DNQwK1VAFZlQfB5918Uhnj3axm9ev3zWjxmgQO/dK3WTHAPDUfaf6BpZN+b1Mz04B/Pzc9hSXcbHt6ymoTbM2qUqoIpMRXx6ZU+/d53rVFy4JjCBXhcemZlozPHa2z0jJyrtP9nFUCRGXo5xw6rFPHrnWrbXlXF9RQl5KqCKTFso5A1Iu/u8QJ+K9GZgAn18RK8c/eQ45zjWftG76EZjOy8e66B3wLsG7/rli/jEltU01IW5uaqU+QWB+ZqIpF18rZvufu/EQI3op+DSevQK9ONp7R3wroHa6M2Oead3AICVJfO477rlbKsNs62mjPCCgjS3VCS44jn6eOqmJAVLZgcv0CvOj+gdGGbvsc6RdExj6wUAFhflsc2//mlDbRmVpUXKs4ukSHx6ZTx1oxH9FJifNp7LxdjBSJRXTnWPBPZXW3qIxhyFeSFuqirlgRsraKgNs2H5IkK6VJ5IWhiMFGMLckMpWZwvqYHezB4E/gpYD9zsnNufzPefyFxc6yYWcxx+p9e/uHUH+4530j/sXSrv3atK+ORtNTTUhtm0uoSCXJ2oJJIJzB/R9/QNp2Q0D8kf0b8O/DbwzSS/71XNlemVpzr6vDy7f6m8Tv9SebVLFvChm1bRUBtmc3Upiwp1qTyRTBQyAwddfUMpyc9DkgO9c+4wkJZ8b1Bz9F0Xh/yzT6+8VN6OdeVsrw2zrSbMsuLCNLdURCYjvtbN2Z4BlqXoEpeBydHHZfuIfjAS5eWTXTzf6AX21972TlRaWJjL1pFL5YWpKZ+vAqpIFgqFvOmVLV19XFdRnJJjTjnQm9m/AMvGeOnPnXOPT+F9dgI7ASorK6fajCtk61o3zjmOnrvAc41tPN/Uzt5jXp49N2TcUFnCo3eu5Za1Ya7XSo8igWAG/cNR+oejrFpclJJjTjnQO+fuSsaBnXO7gF0A9fX1M47Ol06Ymuk7zb7W8wPsbmrnuUbvZKXW896CYNXl8/md+gq215WzpbqUhcqziwROKOGXeMVipW6mJJNz9APDUV463slzjW0819jOW++cB7z57A21YW6pC7O9rlwLgonMAYkzmytLM3REPxEz+wDw34By4B/N7KBz7j3JPMb4x/buMyFHH4s53jzby/NN7TzX2Ma+E966Mfk5IeqrFvO5e9dxS205167QfHaRuSY+KM0NGeuWLUzJMZM96+ZnwM+S+Z6TZWaYpS9H33lxiOca23jmSBvPNraNXOB63dKFfGLLarbXhbl5TSlF+YH5ESUi0xCfRHHN8oUpOVkKApS6gdReXT0ac7za0s0zR9r49dE2DrV045yXjrl1bTm31pVzS12YJYs07VFELhmKxAC4YdXilB0zYIF+dlM3becHefZoG88cbeO5xja6+oYxg42rvNkxt60r57qVxeQoHSMi43jhWAcAd6xfkrJjBirQW5JH9M55ufan3jzH04eixQlUAAAF7ElEQVRbee3tHgDCC/K5/Zol7Fi3hFtqwyyer8vlicjkfPrOWv7Tk4dpqAmn7JiBCvShJOTohyIx9h7v4Kk3z/Evb57jTM8AZrCpcjGffc86bltbrkXBRGTa7rhm6axfI3a0gAX66V1dfWA4yq/eauUfXzvLM0faOD8YoTAvxK115Tx691ruuGaJ1mgXkawVwEA/uX2HIjGeb2rjH149yz+/8Q4Xh6KEF+TzG9cv5+4NS2moDaesIi4iMpsCFehtEsXYptYL/P3eU/z0lRa6/WVC379xBb95/Qo2V5epkCoigROoQB8yG3MJBOcczxxt438808yLxzrJyzHuuXYZH9y0ku215eTnag0ZEQmugAX6K0f0u5va+dL/e4tDLT2sKC7k3997DQ/WVyjnLiJzRsAC/aVibHffEH/5xBs8fvAMq0rn8aUPXscHbqjQ6F1E5pxABfr4PPqj587z8Hf38U7PAH9yZx2f2lGjwqqIzFmBCvQhg+bWC3z0f75IyIwf/eFWbqhM3WnGIiKZKFB5jEjMsfd4J87B3//BFgV5ERECFujjF8r+mwevp3bJgjS3RkQkMwQqdXPj6sUU5IZSfnqxiEgmC1Sg//EfbkXXyxYRuVygAr0WGhMRuVKgcvQiInIlBXoRkYCzdF1j9bJGmLUBJ6f5z8NAexKbk+nU3+CaS30F9TcZVjvnyq+2U0YE+pkws/3Oufp0tyNV1N/gmkt9BfU3lZS6EREJOAV6EZGAC0Kg35XuBqSY+htcc6mvoP6mTNbn6EVEZGJBGNGLiMgEsjrQm9m9ZnbEzJrM7PPpbk8ymNm3zazVzF5P2FZqZk+ZWaN/v9jfbmb2db//h8xsU/paPnVmtsrMfmVmh83sDTP7tL89qP0tNLOXzOxVv79/7W9fY2Z7/f7+0Mzy/e0F/vMm//WqdLZ/Oswsx8xeMbMn/edB7usJM3vNzA6a2X5/W0Z8l7M20JtZDvC3wHuBDcBHzGxDeluVFN8B7h217fPA0865OuBp/zl4fa/zbzuBb6SojckSAT7jnFsPbAEe8f8bBrW/g8Adzrl3AxuBe81sC/Al4Ct+f7uAh/39Hwa6nHO1wFf8/bLNp4HDCc+D3FeA251zGxOmUWbGd9k5l5U3YCvwi4TnXwC+kO52JalvVcDrCc+PAMv9x8uBI/7jbwIfGWu/bLwBjwN3z4X+AkXAAWAz3kk0uf72ke818Atgq/8419/P0t32KfSxAi+43QE8CVhQ++q3+wQQHrUtI77LWTuiB1YCpxOet/jbgmipc+4sgH+/xN8emM/A/6l+A7CXAPfXT2UcBFqBp4BmoNs5F/F3SezTSH/913uAstS2eEa+CnwOiPnPywhuXwEc8M9m9rKZ7fS3ZcR3OZtXrxxrqcq5NoUoEJ+BmS0AfgI86pzrtfHXms76/jrnosBGMysBfgasH2s3/z5r+2tm7wNanXMvm9mO+OYxds36viZocM6dMbMlwFNm9tYE+6a0v9k8om8BViU8rwDOpKkts+2cmS0H8O9b/e1Z/xmYWR5ekP++c+6n/ubA9jfOOdcN/BqvNlFiZvFBV2KfRvrrv14MdKa2pdPWALzfzE4Aj+Glb75KMPsKgHPujH/fivdH/GYy5LuczYF+H1DnV/HzgQ8DT6S5TbPlCeAh//FDeLns+PZP+BX8LUBP/GdiNjBv6P4t4LBz7ssJLwW1v+X+SB4zmwfchVeo/BXwgL/b6P7GP4cHgF86P6Gb6ZxzX3DOVTjnqvD+3/ylc+5jBLCvAGY238wWxh8D9wCvkynf5XQXMGZY/LgPOIqX5/zzdLcnSX36AXAWGMb7q/8wXq7yaaDRvy/19zW8mUfNwGtAfbrbP8W+bsf7uXoIOOjf7gtwf68HXvH7+zrwF/72auAloAn4MVDgby/0nzf5r1enuw/T7PcO4Mkg99Xv16v+7Y14PMqU77LOjBURCbhsTt2IiMgkKNCLiAScAr2ISMAp0IuIBJwCvYhIwCnQi4gEnAK9iEjAKdCLiATc/wcaHs17oJKj3gAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# 生成信号，用于测试\n",
    "N = 511\n",
    "A = 0.8\n",
    "f0 = 440\n",
    "fs = 44100\n",
    "phi = 1.0\n",
    "x = generate_sinusoid(N, A, f0, fs, phi)\n",
    "\n",
    "# 第一种方式计算DFT\n",
    "X_1 = np.array([])\n",
    "for k in range(N):\n",
    "    s = generate_complex_sinusoid(k, N)\n",
    "    X_1 = np.append(X_1, np.sum(x * s))\n",
    "    \n",
    "mX = np.abs(X_1)\n",
    "pX = np.angle(X_1)\n",
    "\n",
    "# plot the magnitude and phase\n",
    "plt.subplot(2,1,1)\n",
    "plt.plot(mX)\n",
    "\n",
    "plt.subplot(2,1,2)\n",
    "plt.plot(pX)\n",
    "plt.show()\n",
    "\n",
    "# 结果和scipy的结果基本相同"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD8CAYAAAB5Pm/hAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3Xl0XOd53/HvM1gJLgCBAVcQBLGQImXJFAWJCyiJWi0rjhXHUrzVVlIlTGylsVrHrp30ZGlPW/vkHG9tjms29rHdOpbt2q4UxamjyLYWUqJIURS1UCQArhApYgdIYp2Zt3/cO+AQBEAsg1kufp9z5szMncu57zsaPXjned77XnPOISIiwRVKdwNERGR2KdCLiAScAr2ISMAp0IuIBJwCvYhIwCnQi4gEnAK9iEjAKdCLiAScAr2ISMDlprsBAOFw2FVVVaW7GSIiWeXll19ud86VX22/qwZ6M1sFfA9YBsSAXc65r5lZKfBDoAo4AfyOc67LzAz4GnAf0Af8rnPuwETHqKqqYv/+/VdrioiIJDCzk5PZbzKpmwjwGefcemAL8IiZbQA+DzztnKsDnvafA7wXqPNvO4FvTLHtIiKSRFcN9M65s/ERuXPuPHAYWAncD3zX3+27wG/5j+8Hvuc8LwIlZrY86S0XEZFJmVIx1syqgBuAvcBS59xZ8P4YAEv83VYCpxP+WYu/bVb902tnuevLzxCNaTVOEclc/3joLHd/+RliKYxVky7GmtkC4CfAo865Xi8VP/auY2y7okdmthMvtUNlZeVkmzGuP/3xq1wcitI/HGVBQUbUmEVErvDvfnSQwUiMoWiMwlBOSo45qRG9meXhBfnvO+d+6m8+F0/J+Pet/vYWYFXCP68Azox+T+fcLudcvXOuvrz8qkXjSdP6+iKSydIRoa4a6P1ZNN8CDjvnvpzw0hPAQ/7jh4DHE7Z/wjxbgJ54imc2TfALQ0RkTptMjqMB+Djwmpkd9Lf9GfBF4Edm9jBwCnjQf+3neFMrm/CmV/5eUls8jvhIXuN5EckGsRRmH64a6J1zzzN23h3gzjH2d8AjM2zXtClzIyLZIJXzRoK3BIICvYhksPioOZUj+sAF+lR+eCIi0+ViqTuWAr2ISBpoRD8DOl9KRDJZPEQp0E9DfHql5tGLSDZIZaQKTKCP04heRLKBRvQzoBy9iGSDVIYqBXoRkTTQiH4GFOdFJBvohKkZ0IheRLJBKpcpDmCgT3cLREQm4Mco5einIR2nFYuITJdL4QTLwAT6+EemefQikg2Uo58BpW5EJBto1s0MaEAvIhnNzzOnMvsQuECvHL2IZAOlbmZAgV5EMpofopS6mQHFeRHJBjGtRz91ml4pItlE0ytnQLNuRCQb6ISpGdCIXkSygXL0M6ATpkQkG2jWzQwodSMi2UAj+hlI5YpwIiJTFS/C6oSpGVCcF5FsoNTNdKThtGIRkenSrJvpGDnbLL3NEBGZDOXoZ0DTK0UkGyjQz4DCvIhkMvPzzErdzIBG9CKSDTSinwEVY0Ukk8WnV2rWzXT4s25SuSKciMh0aR79DCh1IyLZQDn6GdD0ShHJBsrRz4By9CKSDZSjnwGN6EUkG2TUiN7Mvm1mrWb2esK2UjN7yswa/fvF/nYzs6+bWZOZHTKzTbPZ+LEoRy8i2SDTirHfAe4dte3zwNPOuTrgaf85wHuBOv+2E/hGcpo5eQr0IpLJXBqWa7lqoHfOPQt0jtp8P/Bd//F3gd9K2P4953kRKDGz5clq7ETi14xVnBeRbJANs26WOufOAvj3S/ztK4HTCfu1+NtSRiN6EckGGZWjnyIbY9uYvTGznWa238z2t7W1zfjA8YOoGCsi2SAbAv25eErGv2/1t7cAqxL2qwDOjPUGzrldzrl651x9eXn5NJtxJY3oRSQbZEPq5gngIf/xQ8DjCds/4c++2QL0xFM8KaM4LyIZzOLLtaQw0udebQcz+wGwAwibWQvwl8AXgR+Z2cPAKeBBf/efA/cBTUAf8Huz0OYJaUQvIpksHbNurhronXMfGeelO8fY1wGPzLRRM6EcvYhkg2zI0WeceBVYI3oRyQpZkKPPWFrrRkSygUb0M6DUjYhkg4w6MzbbKHUjItlAI/oZ0IheRLJBpi1qllWUoxeRTJaOs/gDF+iVuhGRbKDUzTSYf7qZUjcikg2yYQmEjBNP2WhELyLZQCP6GVCcF5FsoBH9DMSUuxGRLKAR/QwozItIJru0XEvqjhm4QK8cvYhkskvTKzWinzJdYUpEsolOmJoO/zPTCVMikg1UjJ2GmKZXikgWUY5+GuIfmlI3IpLJ0nHOT4ACvUb0IpL54hFKOfppcO7yexGRTOOcS8s1YwMT6EdG9MrdiEiGShyIKnUzDZpeKSKZLjG4pzJUBSbQK0cvIpkuphH99CXmvTSPXkQy1WUjeuXop+byvFf62iEiMpHLYlUKg1UgAn3iX0mlbkQkU10eq1J33IAE+kuPFeZFJFMlxifl6Kfo8ryXQr2IZKZ0xapABPpEsVi6WyAiMjaXEJ80vXKKlKMXkWyQrlgVkEB/6XHfUDR9DRERmUDf8KX4pGLsFCX+ZTzb05/GloiIjO9s96X4pBz9FCXmvc72DKSvISIiEziTEJ9SWU8MRKAfjHg/h0qK8jjXO0AkqoqsiGSe+Ih+cVEeA5HUpZkDEehPdPQBsHlNKTEHrecH09wiEZErne0ZYGFBLtcsW8RJP26lQiAC/fH2CwDctnYJAEfeOZ/O5oiIjOnIO+epKC2iunw+x9oupCxPH4hAf6ztIvm5IT5ww0oWFeby+MG3090kEZHLnO3p58XjHdy9YSnV5QvoHYjQeXEoJceelUBvZvea2REzazKzz8/GMRLtO9FJTfkC5uXn8GD9Kh5/9Qz/8OqZ2T6siMiknB8Y5rM/PkReTogHNlWwbulCAPad6ErJ8XOT/YZmlgP8LXA30ALsM7MnnHNvJvtYAHua2jlwqpv/8BvrAfjMPWs5cKqLf/ODV/jOnhPcs2Ep9VWlrAnPZ3FRHmY2G80QERkxHI3xdlc/b71znhea2/m/B89wfmCYL37weirLilhRUsjy4kK++Wwzd65fQl7O7CZXkh7ogZuBJufcMQAzewy4H0h6oP+7547xX35+mOrwfD500yoAivJzeWznFh576TT/68WT/Nd/emtk/4UFuZQuyKd4Xt7IbX5+LgV5IQpyQxTk5nj3ed7jvJwQOSEImZETSriZEfLvc3L8+5CN7Bcy8P6eGGZg3ufg34P52xn1/LLHXPkeIbOx//2o9x7LWH/fxv2TN84LY733eH83x9o83h/Z8doxdpuncMBx32O8fSceBIzOp7rLXhu98+ino/7tBAvxTeU4o9+XCd93/H975Wvjt+lqaeXL+zZ+v696nAneN75HzHnn0URj3jUpYu7SNhd/HPPuXcJrI7fYpX8z+nXnIBJzDEViDEaiDEVi3i0aY9B/POi/dn4gQk//MD19w/T0D9N2YZCof0bUvLwcbltbziO313JdRTEAuTkhPnfvOv7tD1/l7547zid31Ez8oc7QbAT6lcDphOctwOZZOA71VaV8ckcND22rYmFh3sj2gtwcHtpWxUPbqmjtHeBQSw8nO/s43dlHV98Q3f5/jJaufvqHogxEogwOxxiIRHVxcRGZkBkU5IbIzwmR7w8OFxbmUjwvj9VlRRTPy2NZcSGVpUVUly/gupXF5OdeOWL/wA0VhMy4Z8OyWW/zbAT6sYZFV4RPM9sJ7ASorKyc1oE2riph46qSCfdZsqiQuzYUTur9nHNEYo7BSIyB4SiRqCPqHLGYN2KI+iOH+C026nk0YYTg/PdzAM4b2TjnjUoSX3P+Dpe2J+zLpRFOfHssxvjvPW6/xtg2zt7j/aEbc/M4O4+1ddz3TcZ7jL15SjMaxvuMRv+CmGjQP/oXwehdR/9bu+y1iY9jE7w40XGu1n6b8LXxd77ymFNo01WOM5XPOMfiv6Av/ZqO//INmfk37/WQQSiU+PzS/iEzQiFG7e/9Ws/P9X7x5/u33JAlLQV8/8aVSXmfq5mNQN8CrEp4XgFcURl1zu0CdgHU19dnxDjazMjLMfJyQiwomI2PRkQk9WajArAPqDOzNWaWD3wYeGIWjiMiIpOQ9GGrcy5iZn8M/ALIAb7tnHsj2ccREZHJsUy4IpOZtQEnp/nPw0B7EpuT6dTf4JpLfQX1NxlWO+fKr7ZTRgT6mTCz/c65+nS3I1XU3+CaS30F9TeVArEEgoiIjE+BXkQk4IIQ6HeluwEppv4G11zqK6i/KZP1OXoREZlYEEb0IiIyAQV6EZGAy+pAn+p171PBzL5tZq1m9nrCtlIze8rMGv37xf52M7Ov+/0/ZGab0tfyqTOzVWb2KzM7bGZvmNmn/e1B7W+hmb1kZq/6/f1rf/saM9vr9/eH/hnlmFmB/7zJf70qne2fDjPLMbNXzOxJ/3mQ+3rCzF4zs4Nmtt/flhHf5awN9Anr3r8X2AB8xMw2pLdVSfEd4N5R2z4PPO2cqwOe9p+D1/c6/7YT+EaK2pgsEeAzzrn1wBbgEf+/YVD7Owjc4Zx7N7ARuNfMtgBfAr7i97cLeNjf/2GgyzlXC3zF3y/bfBo4nPA8yH0FuN05tzFhvnxmfJedv25ztt2ArcAvEp5/AfhCutuVpL5VAa8nPD8CLPcfLweO+I+/CXxkrP2y8QY8jnfBmsD3FygCDuAt4d0O5PrbR77XeMuIbPUf5/r7WbrbPoU+VuAFtzuAJ/EWtQxkX/12nwDCo7ZlxHc5a0f0jL3ufWrW/Ey9pc65swD+/RJ/e2A+A/+n+g3AXgLcXz+VcRBoBZ4CmoFu51zE3yWxTyP99V/vAcpS2+IZ+SrwOSDmPy8juH0Fb8Xsfzazl/1l2CFDvsvZvBbvpNa9D7hAfAZmtgD4CfCoc653grW+s76/zrkosNHMSoCfAevH2s2/z9r+mtn7gFbn3MtmtiO+eYxds76vCRqcc2fMbAnwlJm9NcG+Ke1vNo/oJ7XufUCcM7PlAP59q7896z8DM8vDC/Lfd8791N8c2P7GOee6gV/j1SZKzCw+6Ers00h//deLgc7UtnTaGoD3m9kJ4DG89M1XCWZfAXDOnfHvW/H+iN9MhnyXsznQz6V1758AHvIfP4SXy45v/4Rfwd8C9MR/JmYD84bu3wIOO+e+nPBSUPtb7o/kMbN5wF14hcpfAQ/4u43ub/xzeAD4pfMTupnOOfcF51yFc64K7//NXzrnPkYA+wpgZvPNbGH8MXAP8DqZ8l1OdwFjhsWP+4CjeHnOP093e5LUpx8AZ4FhvL/6D+PlKp8GGv37Un9fw5t51Ay8BtSnu/1T7Ot2vJ+rh4CD/u2+APf3euAVv7+vA3/hb68GXgKagB8DBf72Qv95k/96dbr7MM1+7wCeDHJf/X696t/eiMejTPkuawkEEZGAy+bUjYiITIICvYhIwCnQi4gEXEbMow+Hw66qqirdzRARySovv/xyu5vENWMzItBXVVWxf//+dDdDRCSrmNnJyeyn1I2ISMAp0IuIpFBT63kGI9GUHlOBXkQkRZ5rbOOuLz/L/37xVEqPq0AvIpICvQPDPPrYQQAGhjWiFxEJlIHhKH/yg1fouDgEQGFeTkqPr0AvIjKLznT38/Fv7eWZo2382X3XAJDqpWcU6EVEZkEkGuP7e0/ynq88yxtnevn6h2/go5tXA5DqJcYyYh69iEhQRKIxnjx0lq893cjx9otsXlPK3zzwbirLiugb8i6uFUtxpFegFxFJgtbeAR7bd5q/33uKd3oHuGbZQnZ9/Ebu3rCU+FXTQv59TCN6EZHscHEwwlNvnuPxg2/zXGM7kZjjlrow//H+a7lr/VJCocuvGBi/SqZG9CIiGax/KMpzjW38w6GzPPXmOwwMx1hRXMjD29fwoZtWUV2+YNx/Gx/Rp7oYq0AvInIV7RcG+eXhVp46fI7nGtsYGI6xuCiPD26q4P6NK6lfvfiK0ftYlLoREckgTa0X+JfD53jqzXMcONWFc7CiuJAP1a/izvVL2VpTRl7O1CYuhpS6mZlYzHF+MELxvLx0N0VEstDAcJQXj3Xw7NF2fn2klWPtFwF418pFPHrnWu7asIQNyxeNFFanwzSin5nP/p9D/ORAC03/+b3kTvGvrIjMPc45mtsu8szRNp452sbeYx0MRmIU5IbYXF3G7zVUcef6pawomZfU44ZMOfpp+8mBFgAiMUduas8uFpEscX5gmN1NHTzb2MYzR9p4u7sfgJry+Xxs82puW1fO5jWls7pEQchMqZvpaO0dGHmc6jPORCRzxWKON8/2jozaD5zsIhJzLCjIpaG2jEdur+XWtWEqFhelrE1eoE/Z4YCABPqDp7tHHqf6L6WIZJZzvQPsbmrn+cZ2nm1sp/3CIADXrljEzluruW1tOZtWL55yITVZzFSMnZZDLT0jjxXoReaWnv5hXjzWwZ6mdnY3d9DUegGA0vn53FIX5ra15dxSV075woI0t9QTMgvGWjdm9m3gfUCrc+5ds3GMRIfP9o48TvVPIhFJrYHhKAdOdrG7uZ3nmzp4raWbmIN5eTncvKaUD9WvYlttGeuXLZrU3PZUC5mXUkql2RrRfwf478D3Zun9L9PS1T/yONXVbBGZXdGY440zPTzf1M6epg72nehkMBIjJ2RsXFXCH99RR0NNGTdULiY/N/Nn3AUmR++ce9bMqmbjvcc4Fqe7+sgJGdGY04heJMs55zjWftFLxTR18MKxDnr6hwG4ZtlCPrZ5Ndvryrh5TRkLCrIv+zyncvRmthPYCVBZWTnt9+nqG6ZvKMqa8HyOt19Ujl4kC7X2DnipmMYO9jS3c7bHm0m3smQe9167jG21ZWyrCWdMnn0mQiGbO/PonXO7gF0A9fX10+716c4+ACpLixToRbJE78AwLzZ3sKe5g+eb2kcKqIuL8thWE2ZbbRnba8NUlhbN6EzUTBSY1E0qtZ33pk4tW1QIaB69SCYaGI5y4FQXe5q8wH4ooYB605pSfqe+gm01YTYsz8wCajKF5lLqJlm6/dxd6YJ8QNMrRTJBNOZ480yvV0Btbuel46MKqLfX0lAbZmNlCQVz7FR2C8qI3sx+AOwAwmbWAvylc+5bs3GseJFmcZG3mJmKsSKp55zjePtFdjd3sLux/bIC6rqlC/no5kq214a5eU0pCwvn9sKDgVnrxjn3kdl437H09A1hxsiqlamenyoyV7X2Dozk2Pc0tXMmoYD6nmuX0lAbZmtNGUsWFqa5pZlFa91MQ0//MAsLcskNefNnlbkRmR29A8PsPdbJ7qZ2dje10+gXUEuK8thWU8anasJsrw2zuix4BdRkUjF2Grr7hykpyseP88rRiyTJYCTKgZPd7Glu9wuoPURjjsK8EDdVlfLAjRU01M6NAmoyzal59MnS3TdMSVFewiW6FOhFpiO+0uPz/oh934lOBoa9Auq7K4r51I4attWE2bR67hVQkykwa92kUk//MMXz8tJ25RaRbOWc40RH30gq5oVjHXT3eQXUtUsX8OGbvALq5moVUJNJ0yunobd/mJWL541ci1Fr3YiMr/X8AHuaOtjd1M6e5o6RC2+sKC7k7vVeAXVbTRlLFqmAOluUo5+Gi0MRFuTnpu3q6iKZ7Hy8gNrsjdqPnvMKqMXzvALqH+2oYXttmCoVUFNGOfpp6BuKMi8/J21XVxfJJIORKK+c6h5Jx7w6qoD625sqaKgJs2HFInJUQE0LL0evQD8lA8NeoDcVY2UOihdQd/sX3XjpeAcDwzFCBu9eVcInb6uhoVYF1EwSMiMWS+0xszrQD0djDEcd8/JyRlI3ivMSZM45Tnb0jaRiXmjuoMsvoNYt8QqoDX4BdZEKqBlJqZsp6h+OAlCk1I0EWNv5Qfb4gX1306UC6vLiQu5cv5QGfwnfpSqgZoXArHWTKv1DXqAvTBjRqxgr2e7CYIS9xzrY7c+OOXLuPOAVULdWewXUhpoy1oTnq4CahbxBqUb0kxYP9EX5OZhG9JKlhiIxXjnVNZJnf/V0N5GYoyDXK6D+1g0raagt49oVxSqgBoCmV05RPHUzL+9SMVbz6CXTxWKOw+/0jqRiXjreSf9wlJDB9RUl/OFt1V4BtXIxhXkqoAaNTpiaoj5/RH/59Mo0NkhkDM45TnX2jaRiXjjWQefFIQBqlyzgQzetYltNGZury0ZWYZXgUo5+igYSRvQR/5PTgF4yQbyAGr+iUryAumxRIbevW0JDbRkNtSqgzkWBWY8+VfpGcvS5nB/0ppgpRy/pcGEwwkvHLxVQ33rHK6AuKsxla00Zf3RbNdtqw1SrgDrnBWY9ejO7F/gakAP8nXPui7NxnJEcfX6Ii0M6YUpSZygS4+Dp7pGLbhz0C6j5uSFuqlrM5+5dR0NNmHetVAFVLheIE6bMLAf4W+BuoAXYZ2ZPOOfeTPax+ociwOXTKxXnZTbEC6jxVMy+E530DXkF1OsqSth5azXba8NsWq0CqkwsKCdM3Qw0OeeOAZjZY8D9wCwE+kupG50wJcl2yj8D9Xn/DNR4AbWmfD4P3ljBttowW1RAlSkKmRFNcTV2NgL9SuB0wvMWYPPoncxsJ7AToLKycloH6h/2fv8kTq/UrBuZrvYLg+xp7mBPkxfcW7ouFVB3rCunoSZMQ22YZcUqoMr0hUIwHM3+QD9WQvKKXjnndgG7AOrr66fV63+1pZL7rltGYV5II3qZsouDEV467l0D9fmEAurCwly2Vpex89ZqttWEqSlXAVWSJyjF2BZgVcLzCuDMLByHhYV5I1e+CemEKbmK4ahfQG1sZ09zO6+curyA+tn3rKOhNsx1KqDKLArKPPp9QJ2ZrQHeBj4MfHQWjnOZkbVuUlzNlswVizneeuf8yIJge48nFFBXFvMHfgH1RhVQJYUCMY/eORcxsz8GfoE3vfLbzrk3kn2c0bTWjQCc7uwbubj1C80ddPgF1Ory+TxwYwXbasJsrS6juEgFVEmPwKx145z7OfDz2Xjv8Wj1yrmpI15A9WfHnO70CqhLFxVw29pyttWGaagtY3nxvDS3VMSjtW5mIBTy7pWjD7aLgxFeOtHJ7kZvpcfDZ3sBr4C6pbqM399eTUNtGTXlC1RAlYwUlBx9WmhEH0zxAuruJm/dmAOnurwCak6Ier+Auq2mjOtWFpObE0p3c0WuKhA5+nTR9MpgiMUcR86dH7m49UvHO7k4FMX8Aurv3+IVUOurVECV7BSU6ZVpoYuDZ6/TnX0jF914obmd9gt+ATU8n9/eVEFDbRlbqssoKcpPc0tFZi4wxdh00Fo32aPz4pA/5dFb6fFUZx8A5QsLuKWunG013hK+K0pUQJXgCcpaN2mh1E3m6hu6dAbq7qYO3owXUAty2Vxdxr9uqKKhNkztEhVQJfhCZikfkAYo0KsYmymGozEOtXTzfGMHu5vbeeVUF8NRr4B64+rF/Ok9a0fOQFUBVeYaTa+cAZ0wlT7OxQuoXipm77GOkQLqu1YU87A/5bF+dSnz8lVAlblNxdgZ0Fo3qfV2d78/l91Lx7RfGARgTXg+H9i0koaaMFtrVEAVGc38C49cGIzwqe8f4He3reaOa5bO6jEDF+iVupkd3X1DvNDsXXRjT3MHx9svAhBeUMD22jL/DNQwK1VAFZlQfB5918Uhnj3axm9ev3zWjxmgQO/dK3WTHAPDUfaf6BpZN+b1Mz04B/Pzc9hSXcbHt6ymoTbM2qUqoIpMRXx6ZU+/d53rVFy4JjCBXhcemZlozPHa2z0jJyrtP9nFUCRGXo5xw6rFPHrnWrbXlXF9RQl5KqCKTFso5A1Iu/u8QJ+K9GZgAn18RK8c/eQ45zjWftG76EZjOy8e66B3wLsG7/rli/jEltU01IW5uaqU+QWB+ZqIpF18rZvufu/EQI3op+DSevQK9ONp7R3wroHa6M2Oead3AICVJfO477rlbKsNs62mjPCCgjS3VCS44jn6eOqmJAVLZgcv0CvOj+gdGGbvsc6RdExj6wUAFhflsc2//mlDbRmVpUXKs4ukSHx6ZTx1oxH9FJifNp7LxdjBSJRXTnWPBPZXW3qIxhyFeSFuqirlgRsraKgNs2H5IkK6VJ5IWhiMFGMLckMpWZwvqYHezB4E/gpYD9zsnNufzPefyFxc6yYWcxx+p9e/uHUH+4530j/sXSrv3atK+ORtNTTUhtm0uoSCXJ2oJJIJzB/R9/QNp2Q0D8kf0b8O/DbwzSS/71XNlemVpzr6vDy7f6m8Tv9SebVLFvChm1bRUBtmc3Upiwp1qTyRTBQyAwddfUMpyc9DkgO9c+4wkJZ8b1Bz9F0Xh/yzT6+8VN6OdeVsrw2zrSbMsuLCNLdURCYjvtbN2Z4BlqXoEpeBydHHZfuIfjAS5eWTXTzf6AX21972TlRaWJjL1pFL5YWpKZ+vAqpIFgqFvOmVLV19XFdRnJJjTjnQm9m/AMvGeOnPnXOPT+F9dgI7ASorK6fajCtk61o3zjmOnrvAc41tPN/Uzt5jXp49N2TcUFnCo3eu5Za1Ya7XSo8igWAG/cNR+oejrFpclJJjTjnQO+fuSsaBnXO7gF0A9fX1M47Ol06Ymuk7zb7W8wPsbmrnuUbvZKXW896CYNXl8/md+gq215WzpbqUhcqziwROKOGXeMVipW6mJJNz9APDUV463slzjW0819jOW++cB7z57A21YW6pC7O9rlwLgonMAYkzmytLM3REPxEz+wDw34By4B/N7KBz7j3JPMb4x/buMyFHH4s53jzby/NN7TzX2Ma+E966Mfk5IeqrFvO5e9dxS205167QfHaRuSY+KM0NGeuWLUzJMZM96+ZnwM+S+Z6TZWaYpS9H33lxiOca23jmSBvPNraNXOB63dKFfGLLarbXhbl5TSlF+YH5ESUi0xCfRHHN8oUpOVkKApS6gdReXT0ac7za0s0zR9r49dE2DrV045yXjrl1bTm31pVzS12YJYs07VFELhmKxAC4YdXilB0zYIF+dlM3becHefZoG88cbeO5xja6+oYxg42rvNkxt60r57qVxeQoHSMi43jhWAcAd6xfkrJjBirQW5JH9M55ufan3jzH04eixQlUAAAF7ElEQVRbee3tHgDCC/K5/Zol7Fi3hFtqwyyer8vlicjkfPrOWv7Tk4dpqAmn7JiBCvShJOTohyIx9h7v4Kk3z/Evb57jTM8AZrCpcjGffc86bltbrkXBRGTa7rhm6axfI3a0gAX66V1dfWA4yq/eauUfXzvLM0faOD8YoTAvxK115Tx691ruuGaJ1mgXkawVwEA/uX2HIjGeb2rjH149yz+/8Q4Xh6KEF+TzG9cv5+4NS2moDaesIi4iMpsCFehtEsXYptYL/P3eU/z0lRa6/WVC379xBb95/Qo2V5epkCoigROoQB8yG3MJBOcczxxt438808yLxzrJyzHuuXYZH9y0ku215eTnag0ZEQmugAX6K0f0u5va+dL/e4tDLT2sKC7k3997DQ/WVyjnLiJzRsAC/aVibHffEH/5xBs8fvAMq0rn8aUPXscHbqjQ6F1E5pxABfr4PPqj587z8Hf38U7PAH9yZx2f2lGjwqqIzFmBCvQhg+bWC3z0f75IyIwf/eFWbqhM3WnGIiKZKFB5jEjMsfd4J87B3//BFgV5ERECFujjF8r+mwevp3bJgjS3RkQkMwQqdXPj6sUU5IZSfnqxiEgmC1Sg//EfbkXXyxYRuVygAr0WGhMRuVKgcvQiInIlBXoRkYCzdF1j9bJGmLUBJ6f5z8NAexKbk+nU3+CaS30F9TcZVjvnyq+2U0YE+pkws/3Oufp0tyNV1N/gmkt9BfU3lZS6EREJOAV6EZGAC0Kg35XuBqSY+htcc6mvoP6mTNbn6EVEZGJBGNGLiMgEsjrQm9m9ZnbEzJrM7PPpbk8ymNm3zazVzF5P2FZqZk+ZWaN/v9jfbmb2db//h8xsU/paPnVmtsrMfmVmh83sDTP7tL89qP0tNLOXzOxVv79/7W9fY2Z7/f7+0Mzy/e0F/vMm//WqdLZ/Oswsx8xeMbMn/edB7usJM3vNzA6a2X5/W0Z8l7M20JtZDvC3wHuBDcBHzGxDeluVFN8B7h217fPA0865OuBp/zl4fa/zbzuBb6SojckSAT7jnFsPbAEe8f8bBrW/g8Adzrl3AxuBe81sC/Al4Ct+f7uAh/39Hwa6nHO1wFf8/bLNp4HDCc+D3FeA251zGxOmUWbGd9k5l5U3YCvwi4TnXwC+kO52JalvVcDrCc+PAMv9x8uBI/7jbwIfGWu/bLwBjwN3z4X+AkXAAWAz3kk0uf72ke818Atgq/8419/P0t32KfSxAi+43QE8CVhQ++q3+wQQHrUtI77LWTuiB1YCpxOet/jbgmipc+4sgH+/xN8emM/A/6l+A7CXAPfXT2UcBFqBp4BmoNs5F/F3SezTSH/913uAstS2eEa+CnwOiPnPywhuXwEc8M9m9rKZ7fS3ZcR3OZtXrxxrqcq5NoUoEJ+BmS0AfgI86pzrtfHXms76/jrnosBGMysBfgasH2s3/z5r+2tm7wNanXMvm9mO+OYxds36viZocM6dMbMlwFNm9tYE+6a0v9k8om8BViU8rwDOpKkts+2cmS0H8O9b/e1Z/xmYWR5ekP++c+6n/ubA9jfOOdcN/BqvNlFiZvFBV2KfRvrrv14MdKa2pdPWALzfzE4Aj+Glb75KMPsKgHPujH/fivdH/GYy5LuczYF+H1DnV/HzgQ8DT6S5TbPlCeAh//FDeLns+PZP+BX8LUBP/GdiNjBv6P4t4LBz7ssJLwW1v+X+SB4zmwfchVeo/BXwgL/b6P7GP4cHgF86P6Gb6ZxzX3DOVTjnqvD+3/ylc+5jBLCvAGY238wWxh8D9wCvkynf5XQXMGZY/LgPOIqX5/zzdLcnSX36AXAWGMb7q/8wXq7yaaDRvy/19zW8mUfNwGtAfbrbP8W+bsf7uXoIOOjf7gtwf68HXvH7+zrwF/72auAloAn4MVDgby/0nzf5r1enuw/T7PcO4Mkg99Xv16v+7Y14PMqU77LOjBURCbhsTt2IiMgkKNCLiAScAr2ISMAp0IuIBJwCvYhIwCnQi4gEnAK9iEjAKdCLiATc/wcaHs17oJKj3gAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# 第二种方法计算DFT\n",
    "S = generate_complex_sinusoid_matrix(N)\n",
    "X_2 = np.dot(S, x)\n",
    "\n",
    "mX = np.abs(X_2)\n",
    "pX = np.angle(X_2)\n",
    "\n",
    "# plot the magnitude and phase\n",
    "plt.subplot(2,1,1)\n",
    "plt.plot(mX)\n",
    "\n",
    "plt.subplot(2,1,2)\n",
    "plt.plot(pX)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 总结\n",
    "\n",
    "+ 回顾了DFT的计算公式，并尝试用矩阵相乘的角度来理解DFT\n",
    "+ 介绍了两种生成正弦信号的方法\n",
    "+ 实现了两种DFT的计算方法\n",
    "+ 完整代码在[这里](https://github.com/jiemojiemo/audio_signal_processing_demo)"
   ]
  }
 ],
 "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.6.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}