283 lines
153 KiB
Plaintext
283 lines
153 KiB
Plaintext
|
{
|
||
|
|
"cells": [
|
||
|
|
{
|
||
|
|
"cell_type": "markdown",
|
||
|
|
"metadata": {},
|
||
|
|
"source": [
|
||
|
|
"\n",
|
||
|
|
"# Plotting and manipulating FFTs for filtering\n",
|
||
|
|
"\n",
|
||
|
|
"Plot the power of the FFT of a signal and inverse FFT back to reconstruct\n",
|
||
|
|
"a signal.\n",
|
||
|
|
"\n",
|
||
|
|
"This example demonstrate :func:`scipy.fft.fft`,\n",
|
||
|
|
":func:`scipy.fft.fftfreq` and :func:`scipy.fft.ifft`. It\n",
|
||
|
|
"implements a basic filter that is very suboptimal, and should not be\n",
|
||
|
|
"used.\n"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"execution_count": 1,
|
||
|
|
"metadata": {
|
||
|
|
"collapsed": false,
|
||
|
|
"jupyter": {
|
||
|
|
"outputs_hidden": false
|
||
|
|
}
|
||
|
|
},
|
||
|
|
"outputs": [],
|
||
|
|
"source": [
|
||
|
|
"import numpy as np\n",
|
||
|
|
"import scipy as sp\n",
|
||
|
|
"import matplotlib.pyplot as plt"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "markdown",
|
||
|
|
"metadata": {},
|
||
|
|
"source": [
|
||
|
|
"## Generate the signal\n",
|
||
|
|
"\n"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"execution_count": 2,
|
||
|
|
"metadata": {
|
||
|
|
"collapsed": false,
|
||
|
|
"jupyter": {
|
||
|
|
"outputs_hidden": false
|
||
|
|
}
|
||
|
|
},
|
||
|
|
"outputs": [
|
||
|
|
{
|
||
|
|
"data": {
|
||
|
|
"text/plain": [
|
||
|
|
"[<matplotlib.lines.Line2D at 0x10bbc9dc0>]"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
"execution_count": 2,
|
||
|
|
"metadata": {},
|
||
|
|
"output_type": "execute_result"
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"data": {
|
||
|
|
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAgMAAAGsCAYAAABAeaTxAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjYsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvq6yFwwAAAAlwSFlzAAAPYQAAD2EBqD+naQAAhHhJREFUeJztnXecFdX5/z9z73bK0svK0hEsiAgiRRQsKMYWjTUhajSJRmIs0WiaJn7zw6jRJNaYKJqYRGNEY2JF6aIgXZReF+ltd1nYes/vj+XePTNzZuZMn3vneb9e6N0pZ87Mac95nuc8R2GMMRAEQRAEEVsSYWeAIAiCIIhwIWGAIAiCIGIOCQMEQRAEEXNIGCAIgiCImEPCAEEQBEHEHBIGCIIgCCLmkDBAEARBEDEnL+wMmJFKpbB9+3a0adMGiqKEnR2CIAiCyBoYY6iurkZZWRkSCfO5f6SFge3bt6O8vDzsbBAEQRBE1lJRUYEePXqYXhNpYaBNmzYAml+kbdu2IeeGIAiCILKHqqoqlJeXZ8ZSMyItDKRNA23btiVhgCAIgiAcIGNmJwdCgiAIgog5JAwQBEEQRMwhYYAgCIIgYg4JAwRBEAQRc0gYIAiCIIiYQ8IAQRAEQcQcEgYIgiAIIuaQMEAQBEEQMYeEAYIgCIKIOSQMEARBEETMIWGAIAiCIGIOCQMEQRAEEXNIGCAIgiCImEPCAEEQWQtjDOt2VSOVYmFnhSCyGhIGCILIWn73wVqc+/gc/OadVWFnhSCyGhIGCILIWp6cuR4A8Py8TSHnhCCyGxIGCIIgiECpqWvE7urasLNBcJAwQBAEQfjG6p1VuPNfy1Cx/3Dm2NBfT8eI33yEPdV1IeaM4CFhgCAIgvCNrz81H9OWfIUbX/osc6y+KQUAWFZxMKRcEVpIGCAIgiB840hDEwBg7a5DunOM0SqQqEDCAEEQBBEKJApEBxIGCIIgCCLmkDBAEARBhAJZCaIDCQMEQRBEKJDPQHQgYYAgCIIgYg4JAwRBEEQokF4gOpAwQBAEQYQCWQmig6/CwJQpU3DqqaeiTZs26NKlCy699FKsWbPGz0cSBEEQWQIj3UBk8FUYmD17Nm699VZ8+umnmD59OhobGzFhwgTU1NT4+ViCIAgiCyDNQHTI8zPx9957T/X31KlT0aVLFyxevBhnnHGG7vq6ujrU1bXEqq6qqvIzewRBEARBIGCfgcrKSgBAhw4dhOenTJmC0tLSzL/y8vIgs0cQBEEECCkGokNgwgBjDHfeeSdOP/10nHjiicJr7rvvPlRWVmb+VVRUBJU9giAIImAozkB08NVMwDN58mSsWLEC8+bNM7ymsLAQhYWFQWWJIAiCIAgEJAz88Ic/xFtvvYU5c+agR48eQTySIIQwxqAoStjZIAgC5EAYJXw1EzDGMHnyZEybNg0zZsxAnz59/HwcQZiy7cBhnPqbj/D7D9eGnRXCIYwxrNh2EIfqGgEAJNdFn4RJGdHSwujgqzBw66234uWXX8Y//vEPtGnTBjt37sTOnTtx5MgRPx9LEEIem74Wew/V4fcfrgs7K4RD3v9iJy5+8mNc/ESzuZFkgeiTNJEGSDMQHXwVBp555hlUVlZi3Lhx6N69e+bfq6++6udjCUIMdTxZz5tLtwMANu6lWCXZQsJEfUPCQHTw1WeAPEWJSEHTyKwnpelTFEWhESXimAkDRHSgvQkCYvbaPbjy2U+wcc+hsLNCEFmLdtinYSb6mPsMEFGBhIGAuO6FhVi4eT9ue2Vp2FmJLQoNHVmPVttIk87okzCRBrSaHiI8SBgImL3V9WFnIbbQwJH9pGjsyDrMHAiJ6EDCQMDQUprwoC4p+9H5DFCpRp6kRgpXaXeoO4wMJAwEDGnFwoPvk+obU/jqIC1xzTZ0mgGSBSKPmZmAJkfRgYQBIjbws8grnp2PMQ/NwGeb94eYI8IuOp+BkPJByKPXDIh/E+FCwkDAUN2PBsu3Ne+g+doi2gwrm9AvLQwpI4Q0Wp8BZvCbCBcSBojYQANH9pNKhZ0Dwi4JzSjDa3dIMxAdSBgIGKr84UHCQPbzycZ9qr/JgTD6mEYgJN1AZCBhgIgNh+qaws4C4YJFAv8OEvCyD5WZgGSByEDCQOBQ7Q+DZ2ZtwH+Xb9cdp84oe6g4cFh3TEYWqNh/mEKjRwgqimhCwkDAUEMIh9++tzrsLBAucWISeGn+Zox9eCYe/N8qH3JEOIE3DVB3GB1IGAgYqvwE4QyRSUCxsBP8v3eahYAXPt7kR5YIt9DsKDKQMEAQRFYgGvitdAW0Y170oACE0YSEAYIgsgLhsG4x1lNc/PAxm/yTYiA6kDDgIyKnJXJkIgg5pi3Zhtv+uRR1jc2rQESzfGvNgA8ZI1yhjkBI/WFUIGHAJxZvOYCTfz0d/6IId5GGuqLocue/luOt5dvxysLmNuRkYDeLi08Eg5mlhtpfdCBhwCdufnkxKo804J5/r1AdZwAqDzfgnMdm47EP1oSTOYLIIg4cbt7224kDoTYuPhE+qtUEJA1EBhIGHMAYw5fbq3Ck3jiITYrbXo3XDjAG/H3hFqzffQh/nLHe13wSRG4hMBNYjPWkGYge5EAYTUgYcMCHq3bjgj/OxSVPzTO8hq/kvHaAMUYezgRhg/TgIRrXyWcg+yABIJqQMOCAN5ZuAwCs3XXI8Bozx5ji/KTneSKIXCXdkoRLC8lMkNWQA2F0IGHAATKR0IyqOANQlE+fPSpQX5Q9OBnWyUwQPUgAiCY0KvmEWX0v4jQDTSlqGAQhg3YrXICCDmUjtFFRNCFhwAkS/Yuh9MvUwsCRBtpJzy/mrduLtwSbExG5g9VYT0GHoofagZCkgaiQF3YGchWzKl6Q1yKDHalvQutCKgY/+NbzCwAAJ/doF25GCE8Qy9fmgz3JAuFjVgSkGYgOpBnwCRPFgAqz5YmEN+yurg07C4QHaC1qdY1N2HuozvQeMhNEEFpaGElIGHCATPdiZCZgjKlaAJkJ/IfcMrKco20ppWlTT83cYHkrbyZIUUXwncojDZi5Zjcam1LC84wxCjoUUUgYcIDVcibAeABiUNvJDtc3epQrwgjtIMJDNsvsQStgL9i4z/IeXjNQbzBAEd4x6fkFuGHqZ3h2tlhQY4wEgKhCwoADpDQDJoMM3xhqG6iD8hszYSDN4i37ccmT87B4y4EAckTYIV16Tib2vGaAhAH/WbGtEgDw2uJtwvPaIiRhPDqQMOAThj4DGsmY1tz6j8wnvvyZT7B8WyW+8ex8/zNEOEJGqNPCOxDWN5IwEBSNTeKy2ryvRlWO1P1FB3Jj9wmzOs4MfhP+YGcQoc4pujjRDDRyN9WRMBAYjamWb80X29m/m43zTugafIYIS0gzIElTimFH5REA1mubARiO8gxMpQ2gwcd/vPQbq21owt8XbEHF/sPeJUpI4USLxgf1Is2AN1QebsBjH6zBhj3G4diNNAMA8P4XuzK/STMaHUgYkOT7f1uEUVNmYOaa3VLXG9nCtHXfieqTsIbvZLz8xk/MWIefvbES5/1+jmdpEnJoy1GmVBtJGPCcX761En+csR7nm7SBBs4/w2zuRAs8ogMJA5J8uKpZCHhh3iYpB0Lz1QTqvwnvkfbLsFkAs9fuAQAcpvgQgZEuvpSDsZxf4kbCgDcs2tzsZNtgMvuvqm3EL95cia37zDVoNBeKDuQzYBOZZYWA+QBEDoT+8uX2KuQn+fXl3qVtpv4k/OXVRRXqAxJFodIMNJEA5wWymra/fboFc9ft8SQtwn9IGHCAjEBgXsWZ4BfhBQcP1+OCP85VHXPb4ew9VIc7Xl2Ga0f0VA0uRHDUNjRh4ab9tu9rIgdCz7HTnjbvO4xeHUsMz9NkKDqQMOAThnVcG3SD2oKn7KzShx52O37/v3dWYe66vZi7bi96m3RshDWMMaz8qgoDurZWbdhleg+YygZtB16VTWYCb7DbnraYmAqo+4sO5DNgEzOdwNx1e/Do+2sstyVWywItf/1jwVZ87Y9zKZa+CxRBCZmabCTS3FPdEv/ezE5KWPOPhVtx0ZPzcNNLiwJ5Hj+LJWHAOWt
|
||
|
|
"text/plain": [
|
||
|
|
"<Figure size 600x500 with 1 Axes>"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
"metadata": {},
|
||
|
|
"output_type": "display_data"
|
||
|
|
}
|
||
|
|
],
|
||
|
|
"source": [
|
||
|
|
"# Seed the random number generator\n",
|
||
|
|
"rng = np.random.default_rng(27446968)\n",
|
||
|
|
"\n",
|
||
|
|
"time_step = 0.02\n",
|
||
|
|
"period = 5.0\n",
|
||
|
|
"\n",
|
||
|
|
"time_vec = np.arange(0, 20, time_step)\n",
|
||
|
|
"sig = np.sin(2 * np.pi / period * time_vec) + 0.5 * rng.normal(size=time_vec.size)\n",
|
||
|
|
"\n",
|
||
|
|
"plt.figure(figsize=(6, 5))\n",
|
||
|
|
"plt.plot(time_vec, sig, label=\"Original signal\")"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "markdown",
|
||
|
|
"metadata": {},
|
||
|
|
"source": [
|
||
|
|
"## Compute and plot the power\n",
|
||
|
|
"\n"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"execution_count": 3,
|
||
|
|
"metadata": {
|
||
|
|
"collapsed": false,
|
||
|
|
"jupyter": {
|
||
|
|
"outputs_hidden": false
|
||
|
|
}
|
||
|
|
},
|
||
|
|
"outputs": [
|
||
|
|
{
|
||
|
|
"data": {
|
||
|
|
"text/plain": [
|
||
|
|
"[]"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
"execution_count": 3,
|
||
|
|
"metadata": {},
|
||
|
|
"output_type": "execute_result"
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"data": {
|
||
|
|
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAjcAAAHACAYAAABeV0mSAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjYsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvq6yFwwAAAAlwSFlzAAAPYQAAD2EBqD+naQAAYi5JREFUeJzt3Xl8k1W+BvAne9M0DS3doSyCsljUsSiLS0GggGzqjCDVSgeGEWUdYBwZxxG9Co6D4L0guFwHXBCcGcQFsdOCrMpa6ZUKIiLQQjeWNk23rOf+UfK2oQttSZM0eb6fTz62b07ynr4geXp+57xHJoQQICIiIvITcm93gIiIiMidGG6IiIjIrzDcEBERkV9huCEiIiK/wnBDREREfoXhhoiIiPwKww0RERH5FYYbIiIi8itKb3cg0DgcDuTn50Ov10Mmk3m7O0RERO2GEAImkwlxcXGQyxsfn2G48bD8/HzEx8d7uxtERETtVl5eHjp37tzo8ww3HqbX6wHU/MGEhoZ6uTdERETtR1lZGeLj46XP0sYw3HiYsxQVGhrKcENERNQK15rWwQnFRERE5FcYboiIiMiveDXcLF26FHfccQf0ej2ioqLwwAMP4MSJEy5t0tLSIJPJXB4DBw50aWM2mzF79mxERERAp9Nh/PjxOHfunEubkpISpKamwmAwwGAwIDU1FaWlpS5tcnNzMW7cOOh0OkRERGDOnDmwWCwubY4ePYqkpCRotVp06tQJL774IoQQ7rsoREREdF28Gm527dqFmTNnYv/+/cjMzITNZkNycjIqKipc2o0aNQoFBQXSY+vWrS7Pz5s3D5s3b8bGjRuxd+9elJeXY+zYsbDb7VKblJQUZGdnIz09Henp6cjOzkZqaqr0vN1ux5gxY1BRUYG9e/di48aN2LRpExYsWCC1KSsrw4gRIxAXF4dDhw5h5cqVWLZsGZYvX95GV4iIiIhaTPiQ4uJiAUDs2rVLOjZlyhQxYcKERl9TWloqVCqV2Lhxo3Ts/PnzQi6Xi/T0dCGEEMeOHRMAxP79+6U2+/btEwDEjz/+KIQQYuvWrUIul4vz589LbTZs2CA0Go0wGo1CCCFWr14tDAaDqK6ultosXbpUxMXFCYfD0ayf0Wg0CgDSexIREVHzNPcz1Kfm3BiNRgBAeHi4y/GdO3ciKioKN910E6ZPn47i4mLpuaysLFitViQnJ0vH4uLikJCQgG+//RYAsG/fPhgMBgwYMEBqM3DgQBgMBpc2CQkJiIuLk9qMHDkSZrMZWVlZUpukpCRoNBqXNvn5+Thz5kyDP5PZbEZZWZnLg4iIiNqOz4QbIQTmz5+Pu+++GwkJCdLx0aNHY/369fj666/x2muv4dChQ7jvvvtgNpsBAIWFhVCr1QgLC3N5v+joaBQWFkptoqKi6p0zKirKpU10dLTL82FhYVCr1U22cX7vbHO1pUuXSvN8DAYDb+BHRETUxnzmPjezZs3C999/j71797ocnzRpkvR1QkIC+vfvj65du+LLL7/EQw891Oj7CSFc1sE3tCbeHW3ElcnEja25X7RoEebPny9977wBEREREbUNnxi5mT17Nj7//HPs2LGjydspA0BsbCy6du2KkydPAgBiYmJgsVhQUlLi0q64uFgaVYmJiUFRUVG997pw4YJLm6tHX0pKSmC1Wpts4yyRXT2i46TRaKQb9vHGfURERG3Pq+FGCIFZs2bhk08+wddff43u3btf8zWXLl1CXl4eYmNjAQCJiYlQqVTIzMyU2hQUFCAnJweDBw8GAAwaNAhGoxEHDx6U2hw4cABGo9GlTU5ODgoKCqQ2GRkZ0Gg0SExMlNrs3r3bZXl4RkYG4uLi0K1bt9ZfCCIiInKfNp/a3IQnn3xSGAwGsXPnTlFQUCA9KisrhRBCmEwmsWDBAvHtt9+K06dPix07dohBgwaJTp06ibKyMul9ZsyYITp37iy2bdsmvvvuO3HfffeJW2+9VdhsNqnNqFGjxC233CL27dsn9u3bJ/r16yfGjh0rPW+z2URCQoIYNmyY+O6778S2bdtE586dxaxZs6Q2paWlIjo6WkyePFkcPXpUfPLJJyI0NFQsW7as2T8zV0sRERG1TnM/Q70abgA0+Fi7dq0QQojKykqRnJwsIiMjhUqlEl26dBFTpkwRubm5Lu9TVVUlZs2aJcLDw4VWqxVjx46t1+bSpUvi0UcfFXq9Xuj1evHoo4+KkpISlzZnz54VY8aMEVqtVoSHh4tZs2a5LPsWQojvv/9e3HPPPUKj0YiYmBixePHiZi8DF4LhhoiIqLWa+xkqE4K31/WksrIyGAwGGI1Gzr8hIiJqgeZ+hvrEhGIiIm86WWTC4/84iO9ySyCEwDObvsd/bzvp7W4RUSv5zFJwIiJvSVt7COdLq7D7pwv4cs7d2HgoDwAwd/iNXu4ZEbUGR26IKOCdL62Svq62OrzYEyJyB4YbIiIi8isMN0REfmbdunWQyWTSQ6lUonPnzvjtb3+L8+fPt9l5Fy9eDJlMhosXL7bq9R9//DFuvvlmaLVayGQyZGdnu7eDFDA454aIyE+tXbsWvXv3RlVVFXbv3o2lS5di165dOHr0KHQ6nbe75+LChQtITU3FqFGjsHr1amg0Gtx0003e7ha1Uww3RER+yrkfHwAMHToUdrsd//Vf/4VPP/0Ujz76qJd75+qnn36C1WrFY489hqSkpCbbVlZWIjg42EM9o/aIZSkiogAxcOBAAMDZs2cB1GyBs3r1atx2223QarUICwvDb37zG/zyyy8ur8vMzMSECRPQuXNnBAUFoWfPnnjiiSeaVX768ccfccMNN2DAgAHSXnxXS0tLw9133w2gZrNkmUyGIUOGSM+FhITg6NGjSE5Ohl6vx7BhwwAAFosFL730Enr37g2NRoPIyEj89re/xYULF1ze32q14umnn0ZMTAyCg4Nx99134+DBg+jWrRvS0tKkds6y2tWcZb4zZ864HP/4448xaNAg6HQ6hISEYOTIkThy5Ei9ny0kJAQ///wz7r//foSEhCA+Ph4LFiyA2Wx2aWs2m/Hiiy+iT58+CAoKQseOHTF06FB8++23AIBhw4ahd+/euPr2dEII9OzZE2PGjGnw+gYihhsiogDx888/AwAiIyMBAE888QTmzZuH4cOH49NPP8Xq1avxww8/YPDgwS6bDZ86dQqDBg3CmjVrkJGRgb/+9a84cOAA7r77blit1kbPt2vXLgwePBi33HILduzYgaioqAbbPffcc3jjjTcAAEuWLMG+ffuwevVq6XmLxYLx48fjvvvuw2effYYXXngBDocDEyZMwCuvvIKUlBR8+eWXeOWVV5CZmYkhQ4agqqp2Bdz06dOxbNkyPP744/jss8/w61//Gg899FC9DZdbYsmSJZg8eTL69u2Lf/7zn/jggw9gMplwzz334NixYy5trVYrxo8fj2HDhuGzzz7D1KlTsWLFCvztb3+T2thsNowePRr/9V//hbFjx2Lz5s1Yt24dBg8ejNzcXADA3LlzceLECWzfvt3l/b/66iucOnUKM2fObPXP43fa/F7J5ILbLxD5nq5/2iI9Dp+5LH3dXq1du1YAEPv37xdWq1WYTCaxZcsWERkZKfR6vSgsLBT79u0TAMRrr73m8tq8vDyh1WrF008/3eB7OxwOYbVaxdmzZwUA8dlnn0nPPf/88wKAuHDhgvjggw+EWq0Wc+bMEXa7/Zp93rFjhwAg/vWvf7kcnzJligAg/vGPf7gc37BhgwAgNm3a5HL80KFDAoBYvXq1EEKI48ePCwDiD3/4g0u79evXCwBiypQp9fp/Nef1PH36tBBCiNzcXKFUKsXs2bNd2plMJhETEyMmTpxYr////Oc/Xdref//9olevXtL377//vgAg3nnnnYYujxBCCLvdLm644QYxYcIEl+OjR48WPXr0aNFWQO1Vcz9DOXJDROSnBg4cCJVKBb1ej7FjxyImJgZfffUVoqOjsWXLFshkMjz22GOw2WzSIyYmBrfeeit27twpvU9xcTFmzJiB+Ph4KJVKqFQqdO3aFQBw/Pjxeud9+eWXkZaWhld
|
||
|
|
"text/plain": [
|
||
|
|
"<Figure size 600x500 with 2 Axes>"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
"metadata": {},
|
||
|
|
"output_type": "display_data"
|
||
|
|
}
|
||
|
|
],
|
||
|
|
"source": [
|
||
|
|
"# The FFT of the signal\n",
|
||
|
|
"sig_fft = sp.fft.fft(sig)\n",
|
||
|
|
"\n",
|
||
|
|
"# And the power (sig_fft is of complex dtype)\n",
|
||
|
|
"power = np.abs(sig_fft) ** 2\n",
|
||
|
|
"\n",
|
||
|
|
"# The corresponding frequencies\n",
|
||
|
|
"sample_freq = sp.fft.fftfreq(sig.size, d=time_step)\n",
|
||
|
|
"\n",
|
||
|
|
"# Plot the FFT power\n",
|
||
|
|
"plt.figure(figsize=(6, 5))\n",
|
||
|
|
"plt.plot(sample_freq, power)\n",
|
||
|
|
"plt.xlabel(\"Frequency [Hz]\")\n",
|
||
|
|
"plt.ylabel(\"plower\")\n",
|
||
|
|
"\n",
|
||
|
|
"# Find the peak frequency: we can focus on only the positive frequencies\n",
|
||
|
|
"pos_mask = np.where(sample_freq > 0)\n",
|
||
|
|
"freqs = sample_freq[pos_mask]\n",
|
||
|
|
"peak_freq = freqs[power[pos_mask].argmax()]\n",
|
||
|
|
"\n",
|
||
|
|
"# Check that it does indeed correspond to the frequency that we generate\n",
|
||
|
|
"# the signal with\n",
|
||
|
|
"np.allclose(peak_freq, 1.0 / period)\n",
|
||
|
|
"\n",
|
||
|
|
"# An inner plot to show the peak frequency\n",
|
||
|
|
"axes = plt.axes((0.55, 0.3, 0.3, 0.5))\n",
|
||
|
|
"plt.title(\"Peak frequency\")\n",
|
||
|
|
"plt.plot(freqs[:8], power[pos_mask][:8])\n",
|
||
|
|
"plt.setp(axes, yticks=[])\n",
|
||
|
|
"\n",
|
||
|
|
"# scipy.signal.find_peaks_cwt can also be used for more advanced\n",
|
||
|
|
"# peak detection"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "markdown",
|
||
|
|
"metadata": {},
|
||
|
|
"source": [
|
||
|
|
"## Remove all the high frequencies\n",
|
||
|
|
"\n",
|
||
|
|
" We now remove all the high frequencies and transform back from\n",
|
||
|
|
" frequencies to signal.\n",
|
||
|
|
"\n"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"execution_count": 4,
|
||
|
|
"metadata": {
|
||
|
|
"collapsed": false,
|
||
|
|
"jupyter": {
|
||
|
|
"outputs_hidden": false
|
||
|
|
}
|
||
|
|
},
|
||
|
|
"outputs": [
|
||
|
|
{
|
||
|
|
"name": "stderr",
|
||
|
|
"output_type": "stream",
|
||
|
|
"text": [
|
||
|
|
"/opt/miniconda3/envs/pyscixx/lib/python3.12/site-packages/matplotlib/cbook.py:1719: ComplexWarning: Casting complex values to real discards the imaginary part\n",
|
||
|
|
" return math.isfinite(val)\n",
|
||
|
|
"/opt/miniconda3/envs/pyscixx/lib/python3.12/site-packages/matplotlib/cbook.py:1355: ComplexWarning: Casting complex values to real discards the imaginary part\n",
|
||
|
|
" return np.asarray(x, float)\n"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"data": {
|
||
|
|
"text/plain": [
|
||
|
|
"<matplotlib.legend.Legend at 0x11cf1ae70>"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
"execution_count": 4,
|
||
|
|
"metadata": {},
|
||
|
|
"output_type": "execute_result"
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"data": {
|
||
|
|
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAhYAAAHACAYAAAD+yCF8AAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjYsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvq6yFwwAAAAlwSFlzAAAPYQAAD2EBqD+naQAAzjtJREFUeJzsnXV8G0fax3+zkjGxww46ccBhcJjaNNg0pZSZ+dper/BemfnapnxlvDJzmyYNtGkYHEaHmRM7dmxL2nn/UCTtzM6SwALP9z65WquF0e7OzDMPEkophUQikUgkEkkUUOLdAIlEIpFIJKmDFCwkEolEIpFEDSlYSCQSiUQiiRpSsJBIJBKJRBI1pGAhkUgkEokkakjBQiKRSCQSSdSQgoVEIpFIJJKoIQULiUQikUgkUcMd7wbUJKqqYseOHcjJyQEhJN7NkUgkEokkaaCUoqysDC1atICiGOslapVgsWPHDuTn58e7GRKJRCKRJC1bt25Fq1atDL+vVYJFTk4OAP9Nyc3NjXNrJBKJRCJJHkpLS5Gfnx+cS42oVYJFwPyRm5srBQuJRCKRSMLAypVAOm9KJBKJRCKJGlKwkEgkEolEEjWkYCGRSCQSiSRq1CofC4lEIkk1KKXwer3w+XzxbookyXG5XHC73RGnY5CChUQikSQp1dXV2LlzJyoqKuLdFEmKkJ2djebNmyM9PT3sc0jBQiKRSJIQVVWxceNGuFwutGjRAunp6TLxnyRsKKWorq7G3r17sXHjRhQWFpomwTJDChYSiUSShFRXV0NVVeTn5yM7OzvezZGkAFlZWUhLS8PmzZtRXV2NzMzMsM4jnTclEokkiQl3VSmRiIjG+yTfSIlEIpFIJFFDChYSiUQiSRo2bdoEQggWL15s+5gPPvgA9evXj3s7tAwfPhy33nprVNtkhyuuuAJnnHFGTK8hBQuJRCKR1Chbt27F1VdfHXQ6bdOmDf71r39h//79lsfm5+dj586d6N69u+3rnX/++Vi7dm0kTY463377LR577LF4NyMmSMFCIpFIJDXGhg0b0K9fP6xduxafffYZSkpK8MYbb2DKlCkYPHgwDhw4YHhsdXU1XC4XmjVrBrfbfuxBVlYW8vLyotH8qNGwYUPLYl7JihQsJBKJRFJj3HTTTUhPT8ekSZNwwgknoHXr1hg3bhz++OMPbN++Hffdd19w34KCAjz++OO44oorUK9ePVx77bVCE8SPP/6IwsJCZGVlYcSIEfjwww9BCMGhQ4cA6E0hDz/8MIqKivDRRx+hoKAA9erVwwUXXICysrLgPhMnTsRxxx2H+vXro1GjRjj11FOxfv16R7/1tddeQ2FhITIzM9G0aVOcc845we94U8jOnTtxyimnICsrC23btsWnn36KgoICvPjii8F9CCF45513cOaZZyI7OxuFhYX48ccfg9/7fD5cffXVaNu2LbKystCpUye89NJLjtocDaRgIZFIJMd46Y91uPubpaCUxrspjqGUoqLaG5d/du/XgQMH8Pvvv+PGG29EVlYW812zZs1w8cUX44svvmDO9+yzz6J79+5YuHAhHnjgAd05N23ahHPOOQdnnHEGFi9ejOuvv54RToxYv349vv/+e/z888/4+eef8eeff+Lpp58Ofl9eXo7bb78d8+fPx5QpU6AoCs4880yoqmrrty5YsAC33HILHn30UaxZswYTJ07EsGHDDPe/7LLLsGPHDkyfPh3ffPMN3nrrLezZs0e33yOPPILzzjsPS5cuxcknn4yLL744qOVRVRWtWrXCl19+iZUrV+LBBx/Evffeiy+//NJWm6OFzGMhkUgkx3jhD78d/pJBbdC9Zb04t8YZRz0+dH3w97hce+WjY5Gdbj2drFu3DpRSdOnSRfh9ly5dcPDgQezduzdouhg5ciT+7//+L7jPpk2bmGPeeOMNdOrUCc8++ywAoFOnTli+fDmeeOIJ07aoqooPPvggaI649NJLMWXKlOBxZ599NrP/u+++i7y8PKxcudKWf8eWLVtQp04dnHrqqcjJyUGbNm3Qu3dv4b6rV6/GH3/8gfnz56Nfv34AgHfeeQeFhYW6fa+44gpceOGFAIAnn3wSr7zyCubNm4eTTjoJaWlpeOSRR4L7tm3bFrNmzcKXX36J8847z7LN0UJqLCQSiYSj0iPrbsSDgKZCm0E0MNEasWbNGvTv35/ZNmDAAMtrFRQUMD4OzZs3ZzQE69evx0UXXYR27dohNzcXbdu2BeAXGOwwZswYtGnTBu3atcOll16KTz75xDD1+po1a+B2u9GnT5/gtg4dOqBBgwa6fXv27Bn8u06dOsjJyWHa/cYbb6Bfv35o0qQJ6tati7ffftt2m6OF1FhIJBJJCpCV5sLKR8fG7dp26NChAwghWLlypTDkcfXq1WjQoAEaN24c3FanTh1mH4+PNUX4VBXl1T6UVnqQm5kGALZMM2lpacxnQghj5jjttNOQn5+Pt99+Gy1atICqqujevTuqq6stzw0AOTk5WLRoEaZPn45JkybhwQcfxMMPP4z58+frQl+N2ivabtbuL7/8ErfddhsmTJiAwYMHIycnB88++yzmzp1rq83RQmosJBKJhCP5PCz8E0x2ujsu/+zWKGnUqBHGjBmD1157DUePHmW+27VrFz755BOcf/75huejlKJkzxEAgKr6n1Krgg4oXrgAm/aVB/dbsGBBOLcwyP79+7Fq1Srcf//9GDVqVNBE4xS3243Ro0fjmWeewdKlS7Fp0yZMnTpVt1/nzp3h9XpRXFwc3FZSUhJ0PrXLjBkzMGTIENx4443o3bs3OnTo4NjhNBpIwUIikUg4ktB3M2l49dVXUVVVhbFjx+Kvv/7C1q1bMXHiRIwZMwYtW7Y09Y1QNQ/Gd0ywOO/SK7Fx/Tq88ORDWLt2Lb788kt88MEHABB2UbYGDRqgUaNGeOutt1BSUoKpU6fi9ttvd3SOn3/+GS+//DIWL16MzZs343//+x9UVUWnTp10+3bu3BmjR4/Gddddh3nz5qG4uBjXXXcdsrKyHP2GDh06YMGCBfj999+xdu1aPPDAA5g/f76jdkcDKVhIJBIJRzJGhSQLhYWFWLBgAdq3b4/zzz8f7du3x3XXXYcRI0Zg9uzZaNiwocnR2knW/4zy2xRgwhsfYOpvP6Nnz554/fXXg1EhGRkZYbVRURR8/vnnWLhwIbp3747bbrst6Bxql/r16+Pbb7/FyJEj0aVLF7zxxhv47LPP0K1bN+H+//vf/9C0aVMMGzYMZ555Jq699lrk5OQ4KgR2ww034KyzzsL555+PgQMHYv/+/bjxxhsdtTsaEFqLelBpaSnq1auHw4cPIzc3N97NkUgkCUbB3b8AAD6/bhAGtWsU59aYU1lZiY0bN6Jt27ZhV6FMNnwqxYodhwEAnZrlIMPtwoa9R3CkygsA6NmqPgDgiSeewBtvvIGtW7fGq6kRs23bNuTn5+OPP/7AqFGjauy6Zu+V3TlUOm9KJBKJJPnQLIm/+PAddOvVB3WrCzBz5kw8++yzuPnmm+PXtjCYOnUqjhw5gh49emDnzp248847UVBQYJr7IlGRgoVEIpGANX/UHj1u8qJ9RFs2bsDbL09A6eGDaN26Ne644w7cc889cWtbOHg8Htx7773YsGEDcnJyMGTIEHzyySe6KJBkQAoWEolEAlaYoEkZF1IbED+Xfz/8JP798JNBU0gyMnbsWIwdG59w4WgjnTclEokE3JQl5YqER2qVEhcpWEgkEgk4U0gc2yExhpp8kiQOSSNYPPXUU+jfvz9ycnKQl5eHM844A2vWrIl3syQSSYogp6kkgAr/lCQYSSNY/Pnnn7jpppswZ84cTJ48GV6vFyeeeCLKy8utD5ZIJBILGB8LOWslPPIZJS5J47w5ceJE5vP777+PvLw8LFy4MCnDcSQSSWKhddiUzpvJQ7jZNSWxI2kEC57Dh/1JUsyytFVVVaGqqir4ubS0NObtkkgkyYnUWEgk0SFpTCFaKKW4/fbbcdxxx6F79+6G+z311FOoV69e8F9+fn4NtlIikSQTUpiIP8OHD8ett94a/FxQUIA
|
||
|
|
"text/plain": [
|
||
|
|
"<Figure size 600x500 with 1 Axes>"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
"metadata": {},
|
||
|
|
"output_type": "display_data"
|
||
|
|
}
|
||
|
|
],
|
||
|
|
"source": [
|
||
|
|
"high_freq_fft = sig_fft.copy()\n",
|
||
|
|
"high_freq_fft[np.abs(sample_freq) > peak_freq] = 0\n",
|
||
|
|
"filtered_sig = sp.fft.ifft(high_freq_fft)\n",
|
||
|
|
"\n",
|
||
|
|
"plt.figure(figsize=(6, 5))\n",
|
||
|
|
"plt.plot(time_vec, sig, label=\"Original signal\")\n",
|
||
|
|
"plt.plot(time_vec, filtered_sig, linewidth=3, label=\"Filtered signal\")\n",
|
||
|
|
"plt.xlabel(\"Time [s]\")\n",
|
||
|
|
"plt.ylabel(\"Amplitude\")\n",
|
||
|
|
"\n",
|
||
|
|
"plt.legend(loc=\"best\")"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "markdown",
|
||
|
|
"metadata": {},
|
||
|
|
"source": [
|
||
|
|
"**Note** This is actually a bad way of creating a filter: such brutal\n",
|
||
|
|
"cut-off in frequency space does not control distortion on the signal.\n",
|
||
|
|
"\n",
|
||
|
|
"Filters should be created using the SciPy filter design code\n",
|
||
|
|
"\n"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"execution_count": 6,
|
||
|
|
"metadata": {
|
||
|
|
"collapsed": false,
|
||
|
|
"jupyter": {
|
||
|
|
"outputs_hidden": false
|
||
|
|
}
|
||
|
|
},
|
||
|
|
"outputs": [],
|
||
|
|
"source": [
|
||
|
|
"plt.show()"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"execution_count": null,
|
||
|
|
"metadata": {},
|
||
|
|
"outputs": [],
|
||
|
|
"source": []
|
||
|
|
}
|
||
|
|
],
|
||
|
|
"metadata": {
|
||
|
|
"kernelspec": {
|
||
|
|
"display_name": "Python 3 (ipykernel)",
|
||
|
|
"language": "python",
|
||
|
|
"name": "python3"
|
||
|
|
},
|
||
|
|
"language_info": {
|
||
|
|
"codemirror_mode": {
|
||
|
|
"name": "ipython",
|
||
|
|
"version": 3
|
||
|
|
},
|
||
|
|
"file_extension": ".py",
|
||
|
|
"mimetype": "text/x-python",
|
||
|
|
"name": "python",
|
||
|
|
"nbconvert_exporter": "python",
|
||
|
|
"pygments_lexer": "ipython3",
|
||
|
|
"version": "3.12.11"
|
||
|
|
}
|
||
|
|
},
|
||
|
|
"nbformat": 4,
|
||
|
|
"nbformat_minor": 4
|
||
|
|
}
|