ericteunis/m-thesis-introduction
1
0
Fork 0
mirror of https://gitlab.science.ru.nl/mthesis-edeboone/m-thesis-introduction.git synced 2024-11-13 10:03:32 +01:00
Code Issues Projects Releases Packages Wiki Activity
m-thesis-introduction/fourier/03_phase_by_correlation.ipynb

441 lines
91 KiB
Text
Raw Permalink Normal View History

Some work on (auto) correlations
2022-03-07 18:31:31 +01:00
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Phase by Correlation"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"##### Continuous Cross Correlation\n",
"For signals $u(t)$ and $v(t)$ \n",
"and correlation parameter $\\tau$:\n",
"\n",
"$$\n",
" Corr(\\tau; u, v)\n",
" =\n",
" \\int_{-\\infty}^\\infty \\mathrm{d}t \\; u(t) \\, v^*(t - \\tau)\n",
"$$\n",
"\n",
"\n",
"##### Discrete Cross Correlation\n",
"For signals $u[n]$ and $v[n]$ \n",
"and correlation parameter $k$:\n",
"\n",
"$$\n",
" Corr(k; u, v)\n",
" = \\sum_n u[n] \\, v^*[n-k]\n",
"$$\n",
"\n",
"---\n",
"\n",
"Writing\n",
"\n",
"$$\n",
" u(t, \\phi_t) = A_1 \\exp(2i \\pi f t + \\phi_t)\n",
"$$\n",
"\n",
"\n",
"Correlation $u \\triangleq u(t_1, \\phi_{t1})$ and $ v \\triangleq u(t_2, \\phi_{t2}) $:\n",
"\n",
"$$\n",
" Corr(\\vec{\\tau}; u_1, u_2)\n",
" =\n",
" \\int_{-\\infty}^\\infty \\mathrm{d}\\vec{\\zeta} \\; u_1(\\vec{\\zeta}) \\, u_2^*(\\vec{\\zeta} - \\vec{\\tau})\n",
"$$\n"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"from scipy import signal\n",
"rng = np.random.default_rng()"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"class SineWave:\n",
" \"\"\"\n",
" A periodic sine wave, totally described by its frequencies, phases and amplitudes.\n",
" \"\"\"\n",
" def __init__(self, frequencies, phase = 0, amplitude = 1):\n",
" self.freqs = frequencies\n",
" self.phase = phase\n",
" self.amplitude = amplitude\n",
" \n",
" def __call__(self, time, complx=True):\n",
" if complx:\n",
" return self.amplitude * np.exp( 2j * self.freqs*time + 2j * self.phase)\n",
" else:\n",
" return self.amplitude * np.cos( self.freqs*time + self.phase)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Initial signal"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"mysignal = SineWave( 5e1, np.pi/2, 1)"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Frequencies: 50.0\n",
"Phase: 0.5π\n",
"Amplitude: 1\n"
]
}
],
"source": [
"# show the signals\n",
"print(\"Frequencies:\", signal.freqs)\n",
"print(\"Phase:\", \"{}π\".format(signal.phase/np.pi))\n",
"print(\"Amplitude:\", signal.amplitude)"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"[<matplotlib.lines.Line2D at 0x14591292bb10>]"
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAD4CAYAAADhNOGaAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOy9d3hUx9m/f8+qFyQQaiABoghUEFVUYxtTTTHgXmM7ceIkb5z+S+LkTU+cOMk3id/0OHYc3AvGBoMwohiMbZooQh0JUSTUJYSEunbn98esHBlLqO3u2ZHOfV177e7ZUz6I3fPMPPMUIaXExMTExGToYjFagImJiYmJsZiGwMTExGSIYxoCExMTkyGOaQhMTExMhjimITAxMTEZ4ngaLaA/hIaGypiYGKNlmJiYmGjFsWPHqqSUYVdv19IQxMTEkJaWZrQMExMTE60QQpzvarvpGjIxMTEZ4piGwMTExGSIYxoCExMTkyGOaQhMTExMhjimITAxMTEZ4jjEEAgh/i2EqBBCZHbzuRBC/EkIUSCEOCWEmNXps4eEEPn2x0OO0GNiYmJi0nscNSP4D3DzNT5fBcTaH48CfwcQQoQAPwHmAXOBnwghRjhIk4mJiYlJL3BIHoGU8n0hRMw1dlkPPC9VzetDQojhQohRwGJgl5SyBkAIsQtlUF5xhK6r+c2R35Bbk+uMU5t0g01KrjS309hmxWqTCMDXy4Nhvp54eWjgmbS1Q0s9tDWCzQYWD/DyA58g9drEpBsk0NDSTmOrlTarDQAfTwsBPp74efXvuxMXEsf35n7PgSoVrkooiwKKOr0vtm/rbvunEEI8ippNMHbsWOeoNHEYTa1WSmqbqGlsxWr7dM8LISDQx4tRwb6EBHgboLAHmmvhcgk0X4KuenYIAX4hEBQFvkGu12fitrRZbZRebqbqSgut7bYu9/Hz9iBimC/hQT5YhHCxwk/jKkPQ1b9UXmP7pzdK+TTwNEBycnK/uuk4w5KafJLLTW38OiWHV48WEeDtwboZUaxIiCBhdBAj/L1pbrdypuIK75+uYvOJYtJPN5IUFcyTtyeRODrYaPlw6Txs+wac2QuBkTD9XohdDuEJ4DMMmmqhIhvyU+Hky1CYA/G3wM1PQnC00epNDMRqkzz7QSF/3JVPS7uVZfERrJk2ihljhhMR5IuUcKGmkSNnq9lysoS0E5eoHubDz9cncvPUUYZqF47qUGZ3DW2TUk7t4rN/AvuklK/Y3+eh3EKLgcVSyi92tV93JCcnS7PEhPuRXlTLl148RnldM48sGs//LJ7EiGuM9tutNt45VcIT23OpbWzlB6vj+ex1MQijRkiZb8LWr6vXN30fkh8BL9/u929thMN/h/2/A08fuO1pmLzSNVpN3IqK+mb+58XjpJ2/xLL4CH6wOo4JYYHXPObI2Rp+ujWL7NI67pgdzS83TMW3ny6j3iKEOCalTP7UdhcZgjXAY8Bq1MLwn6SUc+2LxceAjiii48DsjjWD7jANgfux7VQJ33o9nbBAH/52/yymjxne62NrG1v57qZTpGaXc1dyNL+6NQlPV64fSAn7fwv7fgVj5qsb+ohxvT+++gy88RCUZaiZwfwvO0+riduRW1bHZ587Sm1jG7+6bSobZkT1ejDTZrXxpz35/HlvAdOjg3nus3Od6irtzhA4Knz0FeAgMEUIUSyEeEQI8SUhxJfsu6QAhUAB8C/gfwDsN/xfAEftj5/3ZARM3I9Nx4r52isnmB4dzDtfXdQnIwAw3N+bfzwwm68umcTracV8/dWTHy+uOR0pIfWHyghMvw8e2to3IwAwciI8sku5iN59HPY96RytJm5H5sXL3PP0IaSETV9ewK0zo/s0o/XysPDtFVP4xwOzySmr556nD1JR3+xExV3jsBmBKzFnBO5DSkYpX3n5ONdNDOVfDybj5z2wqe2/3i/kiZQcbpsZxe/vmu58N9GeX8CB/wdzvgCrf6cWgfuLzQpbvwonX4IVT8DCxxyn08TtKKio5/a/HyTQx5NXvjCfsSP9B3S+j85U8fmNaYwbGcBrX5xPkK+Xg5T+F6fOCEyGJkfO1vCN104ya+wInnlo4EYA4As3TODbyyez+cRFfrszzwEqr8Hx55URmPXgwI0AqHDSdX+GhA2Q+r9qzcFkUFJR38xD/z6Kl4dwiBEAWDgxlH9+Zjb55fU8+nxatxFHzsA0BCb9ovRyE19+8RjRw/145sFkhy5yPbZkEvfOHcvf951ha3qJw877Cc4fhG3fhIlLYc0fB24EOrB4qDWGsQvg7a9AyUnHnNfEbWiz2vifF49T09DKcw/PdYgR6OD62DB+d+c0DhXW8MT2bIedtydMQ2DSZ1rbbXz5xeM0t1l5+sHZ14wM6g9CCH62LpHkcSP43qZT5JXVO/T8XKmETZ+F4ePgzufAw8FR1J4+cNfz4B8Cr38Gmi879vwmhvKrlBzSzl/iN3dMIyna8SHPt86M5vOLxrPx4HneOlHs8PN3hWkITPrMU7tPc7Kolt/eMZ1J4cOccg1vTwt/u38WAT4efP3VEzS3WR1zYinhrS9C0yV1s/Z1Uu5CYDjcuREuX4R3vtF1UpqJduzNLee5D8/x8MIY1k0f7bTrPL4qjrnjQ/jhW5mcr25w2nU6MA2BSZ84craGv+8/w93JY1gzzblJMOFBvvzujunkltXzO0etF6Q9C2f2wIpfQuSnIp0dy5g5cNMPIGsznHrdudcycTrVV1r47qYM4iKH8f3VcU69lqeHhT/ePQOLRfCN107S7uQoOtMQmPSa5jYr39mUztgQf358S4JLrnlTXDgPzB/Lvz88y7HzlwZ2spqzkPojtS4w5/OOEdgTi74JY+bBju9CfblrrmniFH6yNYu6pjaeumcGPp7OrzMVNdyPX26YyokLtTz34TmnXss0BCa95qnd+ZyvbuTXtyUR4OOq6iTw+Kp4RgX58vibp/ofSSElbP8WCHtkj6uyly0esO4v0NYEO77jmmuaOJw9OeVsO1XKV5dMIi7SdbWl1k0fzbL4cH6/K48L1Y1Ou45pCEx6RV5ZPf86UMhdydEsnBjq0msH+njyy1unkl9xhX8dKOzfSTLfVPWDlv4Igrusa+g8wibDjd+B7C1QuM+11zYZME2tVn70diaTIwL54o0TXXptIQS/2DAVT4uFH23JxFl5X6YhMOkVT6TkEOjjyfdXxRty/SVxEaxIiOBv7xX0PfOytVFlD4+e6TqX0NUs+KqKUnr3+2BtN0aDSb/414FCSi4388sNSXh7uv6WOSrYj28un8z+05XszqlwyjVMQ2DSI+/lVfD+6Uq+tjTW4aGifeH7q+Nptdr4Q+rpvh148K9QX6rqABnVQ8DLVy1QV2TD8f8Yo8Gkz5TXNfP3fWdYNTWSueNDDNPx4IJxxIYH8ott2Y6LoOuEaQhMrkm71cYT23OIGenPZ+b3sQaPgxkfGsCDC2J4La2I7JK63h1UXw4f/BHi1sLY+c4V2BPxt0DM9bD3CWg0S2rpwO9T82i32Xh8lXOjhHrCy8PCT9cl0thq5UzlFYef3zQEJtfklaNFFFRc4fur4w2ZFl/N15bEEuznxS+3Z/fOX7rv12BtgWU/c764nhACbv61anqz/zdGqzHpgeySOt44VszDC2MYNzLAaDlcNymUA9+9ySl9O4z/ZZu4LXXNbfxx12nmjQ9hRUKE0XIACPb34htLY/noTDV7c3vwl1bkqnpCyZ+D0EmuEdgTkUmqttHRZ6D2gtFqTLpBSskTKdkE+3nx2E2xRsv5GEfU8+oK0xCYdMt/PjxHTUMrP1yTYFyzmC64f/44Ykb684ddp689K9jzM/AOgBvdrDPdDd8BYYEDvzdaiUk3HCys5sOCajUD9Xd8FVB3wzQEJl1S39zGsx+cZVl8uFPqqQwELw8LX7lpElkldd3PCkrTIS8FFn4VAlwb7tojwdFqVnDiRXNW4Kb8eU8B4cN8uG/e0OiPbhoCky554dB5Lje18dUl7jMt7syGmVGMCfHjT3vyu54VHPg9+ATB3EddL643LPqmOStwU46eq+F
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"timing = np.linspace( 0, 5/signal.freqs, int(signal.freqs*4))\n",
"raw = mysignal(timing)\n",
"\n",
"fig, ax = plt.subplots()\n",
"ax.plot(timing, raw.real)\n",
"ax.plot(timing, raw.imag)\n",
"ax.plot(timing, np.abs(raw))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Correlate signal with a sine wave\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$ u[n] = \\exp(2 \\pi i f \\frac{n}{f_s} + i\\phi_t)$$\n",
"sampled with $f_s$, such that $N = f_s T$\n",
"\n",
"$$ v(t; f, \\chi) = \\exp(2 \\pi i f t + i\\chi ) $$\n",
"\n",
"---\n",
"\n",
"$$\\begin{align}\n",
"Corr(\\eta; u[n], v(t)) \n",
" &= \\sum_{n=0}^{N-1} u[n] \\; v^*(t - \\eta; f, \\chi) \\quad \\quad \\quad \\quad \\quad \\quad \\quad \\quad \\quad \\\\\n",
" &= \\sum_{n=0}^{N-1} u[n] \\; \\exp(-2 \\pi i f ( t - \\eta ) -i\\chi ) \\\\\n",
" &= \\sum_{n=0}^{N-1} u[n] \\; \\exp(-2 \\pi i f ( \\frac{n}{f_s} - \\eta ) -i\\chi ) \\\\\n",
" &= \\sum_{n=0}^{N-1} \\exp(2 \\pi i f \\frac{n}{f_s} + i\\phi_t) \\; \\exp(-2 \\pi i f ( \\frac{n}{f_s} - \\eta ) -i\\chi ) \\\\\n",
" &= \\sum_{n=0}^{N-1} \\exp(2 \\pi i f \\frac{n}{f_s} + i\\phi_t -2 \\pi i f ( \\frac{n}{f_s} - \\eta ) -i\\chi ) \\\\\n",
" &= \\sum_{n=0}^{N-1} \\exp(2 \\pi i f ( \\frac{n}{f_s} - \\frac{n}{f_s} + \\eta ) + i\\phi_t -i\\chi ) \\\\\n",
" &= \\sum_{n=0}^{N-1} \\exp(2 \\pi i f \\eta + i\\phi_t -i\\chi ) \\\\\n",
" &= \\exp(2 \\pi i f \\eta + i\\phi_t -i\\chi ) \\sum_{n=0}^{N-1} 1 \\\\\n",
" &= (N-1) \\exp(2 \\pi i f \\eta + i\\phi_t -i\\chi ) \\\\\n",
"\\end{align}$$\n",
"\n",
"---\n",
"Max correlation if $ 2 \\pi i f \\eta + i\\phi_t -i\\chi = 0 + 2 \\pi i k$ :\n",
"\n",
"$$ \\eta = \\frac{ \\chi - \\phi_t }{ 2 \\pi f } $$\n"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"def correlate_function( eta, timing, raw, func, func_params=[] ):\n",
" N = len(raw)\n",
" \n",
" \n",
" np.sum( raw * func(timing - eta) )\n",
" \n",
" \n",
" "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"---\n",
"\n",
"$$\\begin{align}\n",
"Corr(\\eta; u[n], v(t)) \n",
" &= \\sum_{n=0}^{N-1} u[n] \\; v^*(t - \\eta; f', \\chi) \\quad \\quad \\quad \\quad \\quad \\quad \\quad \\quad \\quad \\\\\n",
" &= \\sum_{n=0}^{N-1} u[n] \\; \\exp(-2 \\pi i f' ( t - \\eta ) -i\\chi ) \\\\\n",
" &= \\sum_{n=0}^{N-1} u[n] \\; \\exp(-2 \\pi i f' ( \\frac{n}{f_s} - \\eta ) -i\\chi ) \\\\\n",
" &= \\sum_{n=0}^{N-1} \\exp(2 \\pi i f \\frac{n}{f_s} + i\\phi_t) \\; \\exp(-2 \\pi i f' ( \\frac{n}{f_s} - \\eta ) -i\\chi ) \\\\\n",
" &= \\sum_{n=0}^{N-1} \\exp(2 \\pi i f \\frac{n}{f_s} + i\\phi_t -2 \\pi i f' ( \\frac{n}{f_s} - \\eta ) -i\\chi ) \\\\\n",
" &= \\sum_{n=0}^{N-1} \\exp(2 \\pi i [\\frac{n}{f_s}(f-f') + f' \\eta ] + i\\phi_t -i\\chi ) \\\\\n",
" &= \\exp{(i (\\phi_t - \\chi))} \\; \\exp{(2 \\pi i f' \\eta)} \\; \\sum_{n=0}^{N-1} \\exp(2 \\pi i \\frac{(f-f')}{f_s} n) \\\\\n",
"\\end{align}$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"---\n",
"Write $\\alpha = \\frac{(f-f')}{f_s}$,\n",
"$$\\begin{align}\n",
"0\n",
" &= \\frac{\\partial}{\\partial \\alpha} \\sum_{n=0}^{N-1} \\exp(2 \\pi i \\alpha n)\n",
" \\\\\n",
" &= \\sum_{n=0}^{N-1} \\frac{\\partial}{\\partial \\alpha} \\exp(2 \\pi i \\alpha n)\n",
" \\\\\n",
" &= \\sum_{n=1}^{N-1} \\frac{1}{2 \\pi i n} \\exp(2 \\pi i \\alpha n)\n",
" \\\\\n",
" &= \\frac{1}{2 \\pi } \\sum_{n=1}^{N-1} \\frac{ \\cos(2 \\pi \\alpha n) - i \\sin(2 \\pi \\alpha n)}{ni}\n",
" \\\\\n",
" &= \\frac{-1}{2 \\pi } \\sum_{n=1}^{N-1} \\frac{ i \\cos(2 \\pi \\alpha n) + \\sin(2 \\pi \\alpha n)}{n}\n",
"\\end{align}$$\n",
"\n",
"Actually, we are looking at the amplitude (or power) so we take the absolute value:\n",
"$$\\begin{align}\n",
" \\\\\n",
" &= \\frac{1}{2 \\pi } \\sum_{n=1}^{N-1} {\\frac{ \\cos(2 \\pi \\alpha n) + \\sin(2 \\pi \\alpha n)}{n}}\n",
" \\\\\n",
" &= \\frac{1}{\\sqrt{2} \\pi } \\sum_{n=1}^{N-1} {\\frac{ \\sin( 2 \\pi \\alpha n + \\frac{\\pi}{4}) }{n}}\n",
" \\\\\n",
" &= \\sum_{n=1}^{N-1} \\frac{1}{\\sqrt{2} \\pi } {\\frac{ \\sin( 2 \\pi \\alpha n + \\frac{\\pi}{4}) }{n}}\n",
" \\\\\n",
" &= \\sum_{n=1}^{N-1} \\sqrt{2} \\alpha {\\frac{ \\sin( 2 \\pi \\alpha n + \\frac{\\pi}{4}) }{2 \\pi \\alpha n}}\n",
"\\end{align}$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"---\n",
"\n",
"We get a $\\sin(2\\pi\\alpha n + \\frac{\\pi}{4}) $ dependency, \\\n",
"the amplitude is zero when\n",
"$$\n",
" 2\\pi\\alpha n + \\frac{\\pi}{4} = \\pi k \\; \\text{with k integer}\n",
"$$\n",
"or\n",
"$$\n",
" \\alpha n = \\frac{4k-1}{8} = \\frac{k}{2} - \\frac{1}{8} \\;\\text{with k,n integer}\n",
"$$\n",
"\n",
"Likewise, the amplitude is maximal when\n",
"$$\n",
" 2\\pi\\alpha n + \\frac{\\pi}{4} = \\pi (k+\\frac{1}{2}) \\; \\text{with k integer}\n",
"$$\n",
"or\n",
"$$\n",
" \\alpha n = \\frac{4k+1}{8} = \\frac{k}{2} + \\frac{1}{8} \\; \\text{with k,n integer}\n",
"$$\n",
"\n",
"---"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\\begin{align}\n",
"Corr(\\eta; u[n], v(t)) \n",
"&= \\exp{(i (\\phi_t - \\chi))} \\; \\exp{(2 \\pi i f' \\eta)} \\; \\sum_{n=0}^{N-1} \\exp(2 \\pi i \\frac{n}{f_s}(f-f'))\n",
"\\\\\n",
"&= \\exp{(\\pi i ( 2f' \\eta + \\phi_t - \\chi))} \\; \\sum_{n=0}^{N-1} \\exp(2 \\pi i \\frac{n}{f_s}(f-f'))\n",
"\\end{align}$$\n",
"maximal with\n",
"$$\n",
"n = \\frac{f_s}{2(f - f')} {(k + \\frac{1}{4})} \\text{, and }\n",
"\\eta = \\frac{ \\chi - \\phi_t }{ 2 \\pi f' }\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": []
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"ename": "TypeError",
"evalue": "almost_sinc() missing 1 required positional argument: 'N'",
"output_type": "error",
"traceback": [
"\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
"\u001b[0;31mTypeError\u001b[0m Traceback (most recent call last)",
"\u001b[0;32m<ipython-input-4-a147ffca3697>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m 33\u001b[0m \u001b[0;31m#ax.plot(alphas, np.abs(misfit)**2, label=\"Power\")\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 34\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0;32mTrue\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 35\u001b[0;31m \u001b[0max\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mplot\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0malphas\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0mN\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0malmost_sinc\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msinc\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mpi\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0;36m1.5\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0malphas\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0mN\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mabs\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msinc\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mpi\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0;36m1.5\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0malphas\u001b[0m \u001b[0;34m-\u001b[0m\u001b[0;36m1\u001b[0m \u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 36\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0;32mTrue\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 37\u001b[0m \u001b[0mn\u001b[0m\u001b[0;34m=\u001b[0m \u001b[0;36m10\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;31mTypeError\u001b[0m: almost_sinc() missing 1 required positional argument: 'N'"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAs8AAAEGCAYAAACafXhWAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOzdeXycVb348c+ZfTLZl6Zt0jbdtyzdoJRFNgFZvQiugFdRREQBxQWvV7j4cwEuqChwBUWKgqCyyKIsZau00JY2TZOm+5Ku2ZNmm33m/P6YTNo0yzzJTJvt+3698qrO82Seb0r6PN8553u+R2mtEUIIIYQQQsRmGuoAhBBCCCGEGCkkeRZCCCGEEMIgSZ6FEEIIIYQwSJJnIYQQQgghDJLkWQghhBBCCIMsQx3AQGRnZ+uCgoKhDkMIIYQQQoxyGzZsaNBa5xz/+ohKngsKCli/fv1QhyGEEEIIIUY5pdS+3l6Xsg0hhBBCCCEMkuRZCCGEEEIIgyR5FkIIIYQQwiBJnoUQQgghhDBIkmchhBBCCCEMkuRZCCGEEEIIgyR5FkIIIYQQwqAR1ed5T30Hn330w5N6zQ59gBZTKS7ySdUlqJH1VyaEGKM0mna1lXb2YCGJdL0QKxlDHZYQQhjio5YWVU5Ya8Zz4VCH041kgjG41T6O2P7FEaWpD08kL/BlnHrKUIclhBB9CtDEQetjeMx7ul5r0M+TG7yKjNC5KNQQRieEEH3ThKmz/ING85ugwphD2eT6L0Cp4XPfUlrroY7BsCVLluih2GGwxeNhedlrvLDvIfxhD49d8BhFOUUnPQ4hhIjlQOsBvvzGl2n1t3N5/g18ueRqKmr38cTWB9naso6vl3ydmxfcPNRhCiFED6FwiB+t/hH/3PNPlmZ/gi/Pv5Flk6dhMg1NlbFSaoPWesnxr0vNswFpTiefnHkR35z9IKm2dL71zrdo9DQOdVhCCNGNN+jltvduwxv08q05v+TqWZ8mPz2Nc6bO52uzfsb5+Zfzu02/4+XdLw91qEII0cNjFY/xzz3/5NPTbuC6md9n6aShS5z7M/wiGqbyMpxkO3O5veQXtPnbuOuDuxhJo/ZCiNHvVxt+xc7mnXxnwd1McE5nUmYSAE6bmdxUJ1fk38LicYv5xdpfUNNRM8TRCiHEUZvqN/Hopke5uOBSzsj+LHnpDizm4ZmmDs+ohiGr2URWsg2HzuPWRbey8uBK3jvw3lCHJYQQAGxv2s6z25/ls7M/yyTnQlIcFpLtR5e1jE9zEA4pbl/4Y0I6xC/W/mIIoxVCiKPCOszP1vyMbGc2X5nzHbSG8WnOoQ6rT5I8D8C4FAf+YJhLp1zN1LSp/HLDLwmEA0MdlhBCcN9H95FmS+PrxTfT6gmQnWLvdjzTZcNkArsax1eLvso7B96hrK5siKIVQoijXt3zKlubtvLtxd+mw2clyW7u9uF/uJHkeQAyXTYAWr2a7yz+DlWtVby6+9UhjkoIMdZtqN3Aupp13FB8A6GgA60hq/N+FWU2KdKcNhrb/Vw791qyHFk8WPrgEEUshBARwXCQR8oeYX7WfC6a8gla3AGyXPbY3ziEJHkeAJvFhMtuodnt5+z8s5mdMZsnKp8grMNDHZoQYgz7ffnvyXRkcvWsq2lo92M2K1Id1h7nZblsdPiCmLDzlaKvsL52PeX15UMQsRBCRLxZ9SaH2g/xteKv0e4LEQprMpJ63r+GE0meByjTZaPFHUBr+HLhl9nbspeVB1YOdVhCiDFq95HdrD68mmvnXovT4qTZ7SczyYbJ1LMnalZyZDS6qcPPp2Z+imRrMn/e8ueTHbIQQgCgtWZ55XKmpk3lnEnn0OyOlMKmJ9lifOfQkuR5gDKSrITCmlZvgIsKLiI3KZe/7vjrUIclhBijntvxHBaThU/N/BTeQAiPP0RGHw+eZLsFi1nR4gngsrq4etbVrNi3QjpvCCGGRHlDOVubtnLt3GsxKRPNbj/JDgs2y/BOT4c0OqXUt5VSlUqpzUqpZ5RSjqGMx4jop6EWT6DrgfXBoQ841H5oiCMTQow1vpCPl3e/zPmTzyfLmUWrNzJqk+rsfaGNUooUh5VWT+S8z8z+DCEdkr7PQogh8fftfyfJksSl0y5Fa02LJ0D6MC/ZgCFMnpVSecAtwBKtdSFgBj43VPEYZbOYcFjNtHmDAFw540qUUjy/4/khjkwIMdas2LeCVn8rV8+6GoA2bxClIKWXeueoNKeFDn+QcFgzKWUSp44/lRd3vihrN4QQJ1Wbv403qt7gkmmX4LK6cPtDhEK61/Uaw81Qj4tbAKdSygIkAYeHOB5DUhyWrhGeCckTWDZhGf/a+y/ZNEUIcVI9v+P5rgQYoNUTwGW3YO6l3jkq1WElHIY2X+cAwMwrOdh+kA21G05KzEIIAfD2/rfxhrxcOeNKgK5ByRTH8G1RFzVkybPW+hBwP7AfqAZatNZvHn+eUuprSqn1Sqn19fX1JzvMXqU4LLh9IYKhyEjNxVMv5lD7ISoaKoY4MiHEWFHnrmND7QYun345JhW5lbd6gzFHbVKdkePR0o3zJ5+P0+Lktb2vndiAhRDiGK/vfZ285DyKsosAaPMGMJnAZZPkuU9KqQzgk8BUYCLgUkpde/x5WuvHtNZLtNZLcnJyTnaYvYpOibZ3jtycN/k8rCarPHyEECfNin0r0GgumnIRAB5/iEAwHHPUxmE1Y7WYumbPnBYnH8v/GG/vf5tgOHjC4xZCiCZvE2uq1/CJgk+gVGSmrNUbJNlu7bVT0HAzlGUbHwf2aq3rtdYB4AXg9CGMx7Dow6lrisGWwpl5Z/Jm1ZtSNyiEOCnerHqTGekzmJY+DYA2X+diQQP1gikOC+3eo4nyhVMupMnbJKUbQoiT4q19bxHSIS6eenHXa23ewIgo2YChTZ73A6cppZJU5GPH+cDWIYzHMIfVjO2YkRuIlG7UeeoorS0dwsiEEGNBnbuOjXUbubDgwq7XOnwhAJLs5pjfn2y34PaHutZpnJV/Fk6LkzerelTOCSFEwr1e9TpT06YyK2MWEJk5C4a0JM+xaK3XAs8BpUBFZyyPDVU8A+WyW7oeVgBn55+Nw+xgxb4VQxiVEGIsOL5kA6DDF4yUZJhj39ZddguhsMYTiNzDoqUbb+1/S0o3hBAnVJO3ifU167mo4KKuko0Of+S+k2yX5DkmrfVdWus5WutCrfV1WmvfUMYzEMl2S9d/bIAkaxKnTjiVlQdXStcNIcQJ9c7+d7qVbEBkDYbLwKgzQHLngpzoug2AC6ZcQJO3ibK6ssQGK4QQx1h1aBUazTmTzul6raPzXuSS5Hl0S7KZCYU03kD30edD7YfY27J3CCMTQoxmbf42SmtL+Vj+x7pe01rj9gcNj9pEk+xjZ89On3g6FmXh/UPvJzZgIYQ4xsoDK8lx5jA3c27Xax2+EDaLydDM2XAwMqIchqIPqY5jRm6iD7OVB1cOSUxCiNFvTfUagjrYLXn2BEKEw8ZHbSzmyGZPx96/UmwpLMpdJMmzEOKECYQCfHD4Az6W/7GuFpsQKdswOnM2HEjyPEhJvYzcjHeNZ3bGbEmehRAnzL8P/psUWwolOSVdr0U7ZwxkytNlN3cr2wA4K+8sdjbvpLq9OjHBCiHEMUrrSmkPtHf78A/RsrORUbIBkjwPmt0S6ZV6bN0zREafy+rKaPG1DFFkQojRKqzDrDq0ijMmnoHFdPRB0+GPfIh32YyP3LjsFtz+YLc1GtEHmow+CyFOhJUHV2Iz2Thtwmldr3kDkW25R8LmKFGSPMfBZes+7QmRh09Ih/jg8AdDFJUQYrTa1rSNBk8DZ+Wf1e11tz+IzWLCMoB6wSSbmXAYfMGjvemnpk0lLzlPkmchxAnx/sH3OWX8KSRZk7peG2mLBUGS57gk2SK9Uo9VlF1EijWFtdVrhygqIcRoterQKhSKMyae0e11jz9E0gBGnQGc1sj5x97DlFKclXcWa6vX4g/54w9YCCE61XTUUNVaxbKJy7q9Hr0HDfQeNpQkeY5Dks2MPxgmGDo6cmM2mTll/Cl8ePh
"text/plain": [
"<Figure size 864x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"# Frequency misfit term alpha = (f - f')/fs\n",
"\n",
"N = 10\n",
"\n",
"def misfit_term(alpha, N):\n",
" n = np.arange(0, N-1)\n",
" \n",
" alpha_n = np.outer(n, alpha)\n",
" return np.sum(np.exp( 2j*np.pi* alpha_n), axis=0)\n",
"\n",
"def almost_sinc(\n",
" x, phase=np.pi/4):\n",
" return np.sin(2*np.pi*alpha*n + phase) / (np.sqrt(2)*np.pi*n)\n",
"\n",
"def almost_sinc(alpha, N, phase=np.pi/4):\n",
" n = np.arange(0, N-1)\n",
" \n",
" alpha_n = np.outer(n, alpha)\n",
" return np.sum(almost_sinc_per_n, axis=0)\n",
"\n",
"####\n",
"alphas = np.linspace(-1, 1, 2*500)\n",
"misfit = misfit_term(alphas, N)\n",
"\n",
"\n",
"fig, ax = plt.subplots(figsize=(12,4))\n",
"ax.set_xlabel(\"alpha\")\n",
"ax.set_ylabel(\"misfit sum term\")\n",
"ax.axhline(N-1)\n",
"\n",
"ax.plot(alphas, np.real(misfit), alpha=0.3, label=\"Real\")\n",
"ax.plot(alphas, np.imag(misfit), alpha=0.3, label=\"Imag\")\n",
"ax.plot(alphas, np.abs(misfit), label=\"Abs\")\n",
"#ax.plot(alphas, np.abs(misfit)**2, label=\"Power\")\n",
"if True:\n",
" ax.plot(alphas, (N-1)*almost_sinc(np.sinc(2*np.pi*1.5*alphas)) + (N-1)*np.abs(np.sinc(2*np.pi*1.5*(alphas -1 ))))\n",
"if True:\n",
" n= 10\n",
" ax.plot(alphas, -1*(N-1)*np.sin(2*np.pi*alphas*n - np.pi/4))\n",
"\n",
"fig.legend(loc='center right')\n"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [],
"source": [
"# Header"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"s\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"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.7.6"
}
},
"nbformat": 4,
"nbformat_minor": 4
}
Reference in a new issue Copy permalink