{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# MATH7502 Group project Topic3: Signal processing¶"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## VMLS 7.4 Convolution"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Vectors Convolution¶"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"vecConv (generic function with 1 method)"
]
},
"execution_count": 2,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# convolution of an n-vector a and a m-vector b\n",
"function vecConv(a,b)\n",
" n = length(a)\n",
" m = length(b)\n",
" c = zeros(n+m-1)\n",
" for k in 1:n+m-1\n",
" for i in 1:n\n",
" for j in 1:m\n",
" if i+j == k+1\n",
" c[k] = c[k]+a[i]*b[j]\n",
" end\n",
" end\n",
" end\n",
" end\n",
" c\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"5-element Array{Float64,1}:\n",
" 2.0\n",
" 1.0\n",
" -3.0\n",
" -1.0\n",
" 1.0"
]
},
"execution_count": 3,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"veca = (1,0,-1)\n",
"vecb = (2,1,-1)\n",
"vecConv(veca,vecb)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Properties of convolution"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"1.6653345369377348e-15"
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Symmetric\n",
"using Random\n",
"Random.seed!(1)\n",
"coma = rand(20)\n",
"comb = rand(15)\n",
"convab = vecConv(coma,comb)\n",
"convba = vecConv(comb,coma)\n",
"sum(convab - convba)"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"-2.8033131371785203e-14"
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Associative\n",
"comc = rand(18)\n",
"convbc = vecConv(comb,comc)\n",
"convabnc = vecConv(convab,comc)\n",
"convanbc = vecConv(coma,convbc)\n",
"sum(convabnc - convanbc)"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"-7.494005416219807e-16"
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# For fixed a, the convolution a ∗ b is a linear function of b\n",
"n = length(coma)\n",
"m = length(comb)\n",
"matrixb = zeros(n+m-1,n)\n",
"for p in 1:n\n",
" matrixb[p:p+m-1,p] = comb\n",
"end\n",
"linearab = matrixb*coma\n",
"sum(convab - linearab)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Application: timeseries smoothing"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"using Plots\n",
"using DelimitedFiles\n",
"ECGsample = readdlm(\"ECG.txt\")\n",
"ECGsample = ECGsample'\n",
"plot(ECGsample,title=\"Original ECG signal and smoothing signal\",legend = false)\n",
"kernal = (1/8,1/8,1/8,1/8,1/8,1/8,1/8,1/8)\n",
"smooth = vecConv(ECGsample,kernal)\n",
"plot!(smooth)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"However, this is just an example. The average filter is not suitable to most of biosignals. Bandpass filter is commonly used on biosignals smoothing."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### 2-D convolution"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"matCov (generic function with 1 method)"
]
},
"execution_count": 9,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# convolution of a m*n matrix A and a p*q matrix B\n",
"function matCov(A,B)\n",
" mn = size(A)\n",
" pq = size(B)\n",
" m = mn[1]\n",
" n = mn[2]\n",
" p = pq[1]\n",
" q = pq[2]\n",
" C = zeros(m+p-1,n+q-1)\n",
" for r in 1:m+p-1\n",
" for s in 1:n+q-1\n",
" for i in 1:m\n",
" for j in 1:n\n",
" for k in 1:p\n",
" for l in 1:q\n",
" if i+k == r+1 && j+l == s+1\n",
" C[r,s] = C[r,s]+A[i,j]*B[k,l]\n",
" end\n",
" end\n",
" end\n",
" end\n",
" end\n",
" end\n",
" end\n",
" C\n",
"end "
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"5×6 Array{Float64,2}:\n",
" 14.0 -9.0 43.0 16.0 11.0 21.0\n",
" 60.0 -15.0 179.0 25.0 166.0 204.0\n",
" 98.0 -16.0 414.0 186.0 214.0 394.0\n",
" 68.0 -23.0 287.0 309.0 190.0 420.0\n",
" 60.0 13.0 -91.0 104.0 -245.0 49.0"
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"matA = [1 -1 3;4 -2 6;5 -1 7]\n",
"matB = [14 5 6 7;4 -3 66 54;12 5 -34 7]\n",
"matCov(matA,matB)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Application: image blurring"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAQAAAAEACAAAAAB5Gfe6AAAABGdBTUEAALGPC/xhBQAAACBjSFJNAAB6JgAAgIQAAPoAAACA6AAAdTAAAOpgAAA6mAAAF3CculE8AAAAAmJLR0QA/4ePzL8AAHV3SURBVHjarP1bk6xJlh2GrbW3+/dFRGaeW926q3umB5eZaQxAEBjACBMMBpGE+CBKZqTM9KifwD+mN8lESXyATAbSDAbqhQRxGRCDy6Dn1piq6jqXzIz4Pve9Fx/cI8+pRrcazamwsq7qk3ky4/Nw37732mutzf9qvdyen91nSZkp+iKBGXGXPUQnCUASaGYGxmJdXrgD7iZ1goJkkCBgjUhBMCgBwCNZmKHKSGPKCFr1pFldj0vP1oK1yJgZwCMJRU9IoLn1JmPSIfmyFOuHgy+nY+lYn536l19vH706HW5vCiAAADPpQBKEIJkkACiCFC0AXyqg0P2Pv+qHtXgs53rxDkACCIA0oQsGBRzjK6ZMQGAyVJoIZFISCBIm0ITErgRBJCQAbAICUqcLJJsRYFgHRa+lKBNe64FeSbNPaYpt73tEpKi7G39339HbspZ8CJXtQTDEsXcWU+rNnxxuP3p5fE4zSFnMlNnjAJLgeCIAF9DdF7xTNlnGEbWQvhTrx4fD43Esk8Dxtg1hTmWEQxj/CJLkyjChrwKNKQdAMyYNFJUwjO+GAWAzIwSSkCSZGRKRQWTKSBL04q+9urnbsdRaWGKBSDc7flr/6A/Osd8cltisrj+JSHP114fToq1FPH7J5eam/Mp6e3dzdDt7Ja0AgjKhMnfAYawD7ua/RYPgS4GnNQQEQqJJonH8M7eyQCDHM2BsqYyQaMhcU4KbAXMrkIhMEOP0wAsjaSxUBKKYAQG5USlwN0J9t7O5iTSUZS0GFLNaS/GXz9bMB9izX/kzt2/+4N8+fNIfHsNCL9AvgbqGGeLt6/yjcjgeFuer5bDU9bAsgI3tO7b2O3MjyZRAAw3r6RQ3dyWX7Xg+NQeAJCwkiAQRgOiEABA295JJmh8nEpmpSARZlElzk5siAIOR0DxYtCABigkgUqBICnAzZWauRmWKuT2QJCpobs5X5VDO5/X5J3fVXj7/DwA8vHncHh/P7XLZt60pzSCx9e2tkYji7stxfWl1ORxWvyNpBhx9vn9AANGqpIxWgnl4KLuLlGCmxHirGnsBGofaIaRSNEsIToDE0zOSmQlL24ojZL7RSEnsCdDUYAmDySgaCEkhUQJp4vxpXDMjaXwASQJfX/wWl+Xlv/n8z3/v2VLs7XK4AbAvAPav/ujLP9m2vWE5pDJTgvcmeCn/jF6WtforK3Vdq39UF5PMwogUrYJW4Czqaxx6lThiQMmEEUgAHCd07maKsLFIshErAQhGCArQoYggMlAsmIQy1xDcxU5PkWkEDWadFAnrBAByN1MKFOikMUkjiWIZmdvX7d/+j1ZPN4eXt8/u7p4fDArWz55v/+r8cNn9+ctPom/nyxZ/Ej0i9vOzjP5wT/yI7rUU03o6HI53hxdlMZqj4+FhI62gPXt7ejxuACSCRoAGyEBAKEmCtAYzAtiNCpr18b6p8d3YzU0RsRgl5biTlGkZMFCEMoV0QmDh+LFuIAVBByNCI9YoQ4RCkqrzkqb7VVC8e5usDDudCm8++uRFadu77s/r7Xc/DSOgHofYLpeHt+evYt9aT1RldAFuiCzH43k9rrXUcsi3X9wfXj3jf5XwaoqSKQPE7VaJy91mCjnjcGnl6NvZXR0FSWfIkTGiHgokZarOeze8IMQZ/TNNCUOO0CqhziuBKSuO9J6OsGLIVCZWIKUEiyNEf+vFxv38/sRB0DqjeSmHyuN3v/sX3SlxRCoRmX2/3N9fvtzOe2QP9QZ36uJOJABjhoqXkiAyQ0qKyqLMpDKEJGnpxoh0N7kMWiBzjyjj+bGTZAHKvOZoFBKwsbXN0kBYjv8LdGicrRQVkLqETFBjCzYaKUihkYDNDEcAOP+fIMRcgDi/Y+hHt//1evv81Scf3X5MQATk5fDso4gj9vPj+bx/0bZt3y5thVJWylsohbBCGwfaZDQkzZBZiBwpEVuxbEHnyCmYSjhiu5lJxvXDSQBKiSMBGrfGNXoANI27t2PcAEqNKBEyJGfIAXdz58jojFS4GwViZGHizPeQcwHWnuaxXW7O938oVj/4sh5Op/WHy2Hx4jiXZXkOyXK/PN6/uzz27fH+3LZDNJDI0scR14hmQAkDvdAMIoDMsT4pQaK6mRKh5IhYCyQFGLzmkgSc4MhlZ4KKmT+QNo9ABxQEogMIIGesVShIwJRwKGaKwRxvB9RTcgcAOCdrKT3uS6lQ6jVJK0v5b9e1sKzrd0/P7o6HQ+12OLwA0BCXt6/vt/vLw6X11ssIZOPgEZmuDCC7EZCUVKA4mSBoVJQakX68Pphr5MjjxoBCAj94h3MX6/pHIz0lbRQY8xOFZr4Fz54wI1OZzHBpBEYQwtMGeFqEcOZG6OCOiMwCRcL98VIUqPWf17WY1VLr6fb2dKjrcjgcP019fX73cL5cWnGjMsGSmYB0abjQtlalDiFKpJkhx8dHbAvUYa75LFKKJPt8vpohSXJwXqZPL2ouBwCRUMKYhF2/LYFRc6QxAQaUYZRE85nbc+6CmeRhqexhznc0sqzlnVmFu5szeyK0bz0SDC+l1qWUsqzr8bgco2mpN1mup9YyQxR6x2ZlU3UFCqNzprSJNDOpK9PwtBEliTQLSJB0gERqlmOSX2MHACqfdoJgAmAJ8oMPFp10SJlGhMgQJFnqw4AyFwIA9m4dlfkM2fu225GIdgkcOqWUFsmcEql27mQBQV+KvLi7e2kFAOmt7V0gFJnIyDAGxezFFGnmhWm0XEyoBREQoRnfpSyCUqnoQTMbZwjKYglD2tjukmM8Sjc9xYeZhhMENy/G6KI5AmYCIYKd78u7DxdgcdFc7VyXZVHi3t2sFBhgRESmzBkR7k6nhTsyz+ilIEQrDhFqeSEw7vOH1E52kOiSBczcqQZoA2C1lE4cyYgQLnVh9qYadAZXp7oW5NpZtat0wQ3RC0dQ6e4K0U0dMKgiJRLlKar3BSp1U1hhpgF0AE5ERyl9YQsS+8JgUaMbeyOrRUpEVYQp02ik3JCBgLhKoGXzzFHYRmSyWjlkbz2k5LjGQqO8NY1w3YUIAGXkH3DsZm58vh4r+h7q+5a2LucQfRV3XkNcF4J0CUmrPk5LFiRA5IIc0VLQhztBYo5MAsA3Iz5JMtNobmqax0+JFDgyVrFIpMigRo4+f0bO90+IIlCJTBgLYr+0BPfreYZAOqlMmNko/VJ9piCWKdLsXV2Px9PtskNtu7TLqiazyChlxj4laCYiE245YAaOzFApKSHC9Y3nv6IPSbyPC3Nh6HQaQbMBQyE14zDYZuUSAohMt/GZjguIGjH4WtMANi5slft22XrSHudnkQBgo/YHSbPi4/fRrqloAlz2xy+5ntZXh9Mpt0u7mEUz0Iupx7jJQLPdzEJpHSNjjEyCUrsu6Deff4R3cj76h18USWNxpBJBCfNuGnmTBKNgEqnMg3IE4blGRgHvr11qZIJZvo597ym0p5Rz1P9pZlRKVouTqjQzIzthyMzdSunnB31d19Pt3avDv65t77JZPgIxPqbei5HK8eNJWEIiafMT69/4nAHjAGuuO+B6kZow604lILlEJJgwwGAlAiTSJPJaYWVoEUzMdEIB87z+tgHsoLyL1nqX6sDCni635gXMyLBmRmEh3dxY6UakHt2Ni9njGeVw9/zmE/aHh9Y7EkmnLpUJyw6lkZmLpEzCTADoSdKuHzTHNTmeFAE+fRj6MOsDpJBkDngHkbQGUKDWTplpLIBIaFyDo2aIDEMGfD4l4EoSyPIuew9hbkVDzntf2UWlqEgQCJJm5MHMi5m9yrZ3mVdKj49f2G8tpTyLy6MQYWS6IUWgIGWGHEsuGEsq3RrHR2nfyOvGuZAhbT76UwapgcMaATPBxmdlbgZByAjSxZynRoSJoyyXMgzMhKHPS1QpM2aWc/aAu6ckUNxoNnA99qdcbkA4TAIXmpel+uOyHrO1aJJU3P9/dy9evrrbl8dISOqLRxdNSw/RyJAEI2VQgs6xneXiN3JnDXg9jd8MgqQEKQ+IJDI5YhmrUUxw790sCR95GK51R5sH34xKmG9zATLlzgz/KEEq+ibEtgcMGT2UxwHh506IxdSUe1CRymj73rZ9bwHzgkhj3w7b/bvX7/zlaW91zU5JUHLpMjdlHLMuCZwLALhb2KjCjMi9i31Ae95Bg2wZBx7UwFx5wQAGCDeF3LOnFe2ekYK0eHEIVMgQKJ7y3OYutlLvV2RlC+vuG0r2gfyZv4xRpiVtJK6plDIUojKCxhQz3Ai6QRzPtkXvrbW2eq2UkD0yLpcHOzxfLo92TM7PIHuCmdpSYZYHQoreqnIijwP7WKqb0d1roTKzJWgAFhIZEcWg1LiXIKiCXtx95sjGCeUOIBsWASOtalxqyiU7yGKolvT0UcqjAKMKU/KKQFAJ0ZyZiXFUETL10f+ghBkd3elLtQwyYNou5zfnT1/5loJlglSvKeDaYWLJEMwhjY2VvbTwansokK0HrBRnAXeboWgExvE33W2fJ6IJStHbKKMJUfOKn0dp/LlN1LZ6D6PVTqMp46CRORaTlCJSQirkAAmxOxVJS44UKpQdcI1cII0KmLm1Za0kLZzI3Pt2eXdzvH14wwGcpLEAMisxljcEFIea4NX2iMdHHkvb7ols29ZjPd2eDtX6qHzGcc4ZQyRN9Nw6iYDZhOuvwf1pBdyIhHoZfzJaMAb0SJOUidGYKH2kf0wzCsqcJUdXZKogQY09ZgajdjOKo3CgueVjWY+LhagIAbn95O7zTw5bjEyaF3dF1AIhg55LKtsuIaEe948Pjxdb7PK4eC2I1s9Lhxl8gQAad6O6NG6FCByUIXohkTkyOkB6Ko8/SJskIHldQ8GL9ciBw3szNzCzJAQz5+P45JUjQWEVaHTvpqRB3SVERvdRns0MKxU9kFQllEB63J+3x1e3D2QA7gmlSM+ZszUviEjs0p6XNw/vUub9/vHoy/G4HHTwpSA21pGO8bGuXmB+BqGBWUAkkSPMjhplINrQ+x3QnUl7X0coQaO6w8xYsEukAeW6bQYaniKUtBFUnAgICTAkCqk+ajrmACiTAvPSyMPiMKY82H+ybZ+AVNDd2pZlzaYO7NF3nar1fdeuzP38cLkAiN7R27a124OD2S5JOInsPS7Hm9Xp/kB3M2OQRioVQbMcTTiAI6P4AIAhlBi13QQ0UgRcUBpUoCDoJTITNPpsgtRMCcRZ8IIeKyBBsbolmAaOsBjz2LFU7DLfT4tgiuwocc/+qTIaFuvv3ubptrNclGfwcrcizl/dt2fKvl+amPvWuJ5s27bWT/WBfjgdHcUs9/Pjtmx7yZCl1WVdq3d3V6pXyapnD4pGMOYOuD5/KYpE7scZHQzKJEgqIhlrRsDc+Z1aXZlagQxxhvlUgREDmBy5Fc0m4pwBs3mW2JdCY7ZlLfCVvYlmZV0Px9UMl4fLw7ut3D07Zdv2HhmvAv3ta944Te1y2R8xu64SWJZyUTneVjUsi7ZdNMW+hfmh3Dx/cXKgVGTClHTtYVaYdMRsHAykOnM0+CQUjRwoTR2FyUg6oyfohoyC7CmBwXENDjjcFGOPj/Y3iJgVxMDNwdkkpZRIT0ULuTx9/LXctbNwO+8BZ3ts50F60FdQP3fbKEXb99jcM+mlI8EMrOZ4G6GlLswAOpGg4lJsOa5L6Sb1QXsYWcGgIvSDkBPaJmHWR58LIaRBqVFjZDGo0WofS6ViVBeNzQ2AcqLB8BETMXulA7RVDgSUOS8dCAlRUAhyVVZkBjt2gNX3S3RR+44tRSNyM+1h2DKj9Z5qGLGiKWXBRCmx7yzn/eJU6lwKQbBie5vbTb0TQLpvAysawNs3i4mcGZFIYwgC8oochJv1QK2joSsVG/09xnhKVUmpVJnbKXglm0hUpqBxv0xKjGHcnIhUGl1Ezwxj76we+/grEV2gIVMlLjygNWWEBkBCdyATUmSOXollZiB7NK/VzXjT2znavnpd3MZ+Bw3eAUNhJkd9RIUShCsFiH6tDs0Iy9TIMEeRSEMJcyIjNVv5HYMg0GeKmBxNHUjJkS1JAwWGIK2zzMsMshF7ydYdWVuzwkxESsrcRgqVh+iNrsvOcYA8wpaifYJystbVQjxKoejdAqC7SZkh+sPBndlnlgeP0RHk++A/dqwLqYRsoK0TUiWRsoq4VMlIY+mOkRENzgP3gQPxiqbNh5dDCQ3GiMSB4gJIQkgmQo7ozBW99dzXnt5Sht5SlHaNUmJDdPXsFy/OjCiZkHrzCDMo2QFQ8UggM0XBMjMfs2UH42aBE5FlpLt0SglmmTHAB3aI0W7QE6QgRV7hFFePksjRfh99ETohifMMiQ69b+5hgpkA0zi2Kg0UGByPX2DVU9qdGWjnnqbYWdj3LkBMUU5ckoZNBuXYeW6O6KEIdym3mXWe3dSTDCCDPanouJwPH3eMhx/dXMPE7q47IGfswuzbT9SMkMYpJ9GFQkpJqQxY2gfMA7BmRghYR+eKRqOudSmJ7k4oRq0hQUGJikVl8Q7b0iVES9boG2W99ZS0RKCZgVGKOhaRFI10p0IWPZmJTJMiB3iZ3TxIqUdPDlj0xzeXZ4dS3IRkWkg0VlObMaDPa3DBtUKCGa/oAU3O7PI64AWJn4+WF/vIMJlGZGgURTRebLSBRvA1s8vup7Kf37cm6MWUe6nr6RBbyRpfb9oPUEbvmREwRBwSZWFvXpwJM/VAdKu+7aqVCIlOqdGKI2Jte4DSoVjS1JcJ1LTT3bGa2e3NkVYtrXerrhifFngxN8Ds0RFwhM/sPma/0YjoSatmyBQKB2iOffz8cRU8wctza2mixAntadybrEpJkpapDmk1ZiSdih4y71JG9LTx160lGSYEzdzcLhnKkGUPZlJt5Nfa3RGSeo8goQaTTU6RCFZdwt3t8XC6uTtZLyPrHTl9AisHiFcchIOm+dmP9GTQL92sjyqYZRZ2o/RSZrm2uOfRiplfpsZ1m6WqyUcT1sz2FIsbqFRZajWNS3eDMiJUlYKBmYxm2ZUqlmkZXRli9qSCgzBFCVaYom3RRQINg6jWSCs0fzx3mRmf+eF5ltJ9pGdPQcCJxAA8A0D42L3iIA7KIM06RgJYrnDIgYQy0mZV6DN4BGbjeu550pKi9flDGujuhvMIIsZROrQmZWSKEQmaOaWOiKKMBGEZkjJ3DVyITziGQenWJMnIlMREBgbVdN0jU6E38sfQs2Uf/VYGRyuiG0e2eg2KkgYcOboxJmUqWAiBbqXboEVVKVNmCdIU80xNRiAxUkwCuzroiEElUy6gmVrKirtatp3RO/YuZeZotGi0YyQM9rEjB6uAUIJXNGowY6KrZw6alciRyCuAVGZ3D9Ziyj2RsHoaSOt7+DznnT8aile8a7Bm8HTCr+zXZNH8GwO1NbYcPeLBCr0mxkAffQwhRaO6zIiMAMxoJl+rWY/2sDvPyRYaDyE4RM6eoaNxpoQ9zebNoiQzR25q47bKnGRMAkKoL8K4HVhKdeCm9/4WjhfmhsxccaX22TfaCU957ABYkEqaQX0g41E4v7tNQkNpPdLEMlormIuFhIGSfDQ0rSmhTIWQWBaezRi5b20jLvCGTElEH9cykxRpJmSmeqDH4EhASjNKSShVEAQVI99gIpGQso3LimtZitOYpXF7bc3rodrExISnjGBU7EAYQRMY1ywxxxkfJz9YbBSig7kakc96ZJax9VPAMn/uAkLKpBtCPqpkOo6tdwDenMhsfWnZG9QtR4uimylgaRgciwNYPBWKHLkFZr9k5G/qEjl4qrN/KoHS2Wrx4hTUN9GaSFd7fH0QKjh6jwY22bjA5y2wD3Bn3OQcNC9kZrHJS/uIRmTmgSQyVQzR5P5Y3AGyGzPljkyCmSOr0CiSSG5mNDPekLHvLXPfRMTgPxAcRIlRT0hSTa4rLpenGusavcdJY8zQswmEkb2ulX3vMDMrtRyuuEXA10p8cnN3coBwh+Aes71Ze6A65Nf9nzBnD4bolvmkIyiD0CjOaB8xiuo7IlsMUvBoV8821qj3BwOD0JpJN8W9W2x70zKw6fdncOw66bpDqxeF/BvsoQ9eA9rVExXmTnkBbB0c1uwbTZqXWPfqyL5ZAuluoLsNXIQypxGK+Xt3jC65ACgy6jwyZRD8J60NYJdoxnw0Jxf3h3RaKsqAhTSRMw2QMeKQSSuRSy1BOFpKcX1eTU64Jr4PIN3Res4s9X33+wPyFJDXzJXIjIS7DdQ6sbvbCNoZKXMfCR12MwKDcm9Gurv18UjjbJnjSlKnhsJFgwo8FntcxqSFaJZ5I/W+Dy4C37O7aJOg04ZeIPaEgiyZGb0He0QCmdcbiLMWu+IWQTDjCb79YAXmmZk4y0zf+OBldUAti2E0+gkANricg6gs0olMNcZcAE4OUJ1NXborkiPcCDax5AGIGMHBnKEGjXlQ5DkY4YSl3gPPs+liMoAW5th7rYaMzBycPtGuFdn1+Z6IED3cPK8KjQ93vwaa/02Ih1YWH1REjb0xxDjmRtZ1TQ2Qo0zO/kSyGJ0KeGmjnmNEQYrWjQQLtslZKkkqr41ICsbxu85m5uZ+UQ4wab6fnLtqaInco3iPKFgRKUKkIYfcYjzU3Mw+8+0p2biezW8cganY0futQVW0ZiTrzO4SSc2eRb1eahBiQpQziZcNNmtu8wh0L24ku5uN1TPCpblwIwxQUtKVCTsQ0ffEIZU2++yYhQdIyhgdpdLgZktlu2yttV40gMOxU8TrMfQZ8dOIDqOobzQyrkv8TS7wyBpobmVCnomkYLToDDJcxY1AKG2GpgGDK+lQm5koguaFGGIIq6Wbm0kodYaJNnl8vRgk2GsvRi+UgvqAyFMGMZKXwmhiQcLWtfrN2bz31u6RkaPp8vRcBOmZpKAwiwYuP+MC0FP8n+INUeyTmrPRQKeHKND8brukMvCTG2AxMwsQyrBJum4Jc/S+zuaoZ0joUQcoVh/N3QD58+wJxU4aJdpe410/4F0ddUAu7gRoBTM3iUxJQoZ8cfZL88XNWgsp9h0Qvbhdj8BMQhWRE8dLGZmO0ZG5rmwMitNMiEZT7z1XCsqZeneIy8Fa7kGqXdal0pYyMvaMcFyr1sUVNiOnBHf0bstGInsGMtreei8onoKVpkSk4KlitPXKX9+uFYmut8HIKq5Eo6WSl0dfP2GkecV5nvyfdcM/MZ/exzl98OVJZhuXUT4lBLhWs/OzReypVgcpgG3faNVoACizOgqtay45K1ghu0RTKwqYK2wS8srupi6n22QrNLmlagDQlbykmS+Pq4CCVDDTkbJw1+Fu5JvmnFyWn5HgjA7+0wp8s58/v2lSpoxX9eHkio1agaQ7ewtmTKx4M5Otq1s6CKUpdSVF4omuTdQeViyj9qBb9PFEYslRkMTA2g22aeWuumSmkqiQFAEXJ9Vw7oMBFzNkxuV097C3Fr31nyZ/fnDdP8kdftYOuDKmk8ar6OP65VGMATDQfDx8zKZI9p0lZBx0cvR5Xz/SYBosMQHweZeJNhJfaVz7S6YVZu7myBSdmkTLUYwus1Ficz8OclnK5uZc8gwZ1XqLiLZTP/3p80pb+JkH4xufvoh0cw7FKqGnNdLglAFDYefoZWSXhtjKZUUdJd5ghJBYQTBTZeAgCnjsSduNaml2/RhKRg6R8NAx0EqSTuzzduiYUjIhjWDMlYwZBTwC0DlK9Os7/enHf7+Zv8H7/GY//7oIQ+Z4bci+/4pN5KvBkobMTCbAokDr3Q8cHCLDk1jaqMgyPj/J2BNWLqQiNTJIAGUH1SEVMySETubQi46cs9FGgjFwpal8McTgieAiL5RZRkjysuuaMv70HfdTG+Cn+vnXI7CMVMlmWvxElvN5CjIEozpGl85C8Hq4PR6ZfeuzvTVkMAN5pzi3axnE58HVvP7sglIQstqVipCcSnMLf8/CAAT54G3SBJHGPimOW1kWNSWzJ8EnKubPWoAPU9yfsQMoEDUHNhQU32cEQ32UY4ubeeRgsxg2mR3uXtw6YpMCHaSZWZNMyJlxEuhBt4ybHnBExEjEUE57KwxYL5U0Wo4ukiqRisyCIIzWxs/Fo5kTqWVu0erb/XrHs8/6qnaJvFLhJR/eA9cMdZSeADuSk+Eze7lN5sUvIA1EbYHi2dNckWYSfFHfPC1bBxcO7lJd+fhYXj67e3ybya0VSZ2GbsbMnmXWLo0LIViyACgOgNn3VnKA5ameNliG45MAQBhYhkoaBkUHRk9pVo6ctIllNYiZEZEOJJm6aubtgxz3KZrPX/F+iwOTJ5HX6D8ARQGzozUpFHYlwOVkVcQG9vriO5fiikjFlc8uxSBTz7J8no34ILaYF5aYuNMSgrO3mMqcfEIIRjrvGIpYU0IgO4xGYi+SkiZk9BzRJT94pikdym/eAxOd+nABaOCsOaYGcWwPJWxUK8pRUiDdGDSK0sF6bPcP53fnx8uljdb4QDAT+qCmuHo96AqYCKLDizi4lYMfaCUmqhjzJMfckjlLwLhaW8wPNQS0rWT23mUO2iRvTn5/pz0pvcbij/vTvsmRxjgO74GK2dIzxZRfzeZkSjLAml1bt27tq39z2aLtPYSYamRTQqTFtbc1q5IrTT9zSDgKRsBXM8Mu2uy8XrsgzImG95EnTdxEKjPfWShYBiJSoDPN1YfO7/oyQIqZ9dmA3n5WRjDw/dEeGCgpjJnXXICDjzEPBCdU8Qgvcf/6mGUp+6XlYMxetfC0uBYT82adKcnEqqWSPsiENE/N61HASHUnKDLEVOPnuAKK0T6XoAIpd7cBLYKCuQYGTgJZnzChkfraZJjMvPh9NTyT4Gt67MAQbQCpIStF6iouG6xhCbCe3M5b3h4OeY+eA78czNcJKfNpU80qZihTJUpZchTOWRB6Ithj9kM4BHdD7DQLoUmbmOecISG6+zLazUi7VnFOZuY6tfSc2MD1hvx3t8BsT86UwKSE4Ab2kBnoA5XnyPA5KQILZDo/PD6r61C+jR8zeRZiWo5uyJU2NqNxSiTTBkFDUEaHW49y1eCkcrxtNxq4v89lSNrgFF2zFwm9DzpHDMxtovRKjMOB5YPnnG/iG7XRAFCNMbcAZ8g1Z3aNqD8aM+MImMZN0cKra9/uYegJXk96jCVHwnIygIeGaL6J0dJiloUsgJx1lVTsIpF0nhfPRHZDgnRsNBIa1ZjqEqOm4Johq4tfoqczw9TBCqBiS3c0sJLETvNJJyYBDLx1yFSIDKvosNzKwtZp6kAZ8ru0BQAjwALuXo0itxDdDTL2zsv9wW3/8suNp1qy7TR3BkzhuSdLWTcbbNnd3IgRw0IDEhtR/2ogM2hlsw0E89FOepJDqQFuzkgqRW4zMuXoyc3bA1CMNqRPNvsq0WA+mpzSLPkIKQ2jAKbDzRykqYECAja2yuzoKX3Av1wiQIRcFDP6fbw7v36wdfTCzMOGy4uPhmiEhgbZrnY+c+MVPVFhrozhkfJKGXLDEDym6hULyGRx9LHRYQEws0t0RGqcG3DS62A+g86aAdr45rF3MFjoQ1dAc3MQNZ6ugxHJnCSVI6iOlhYFwgI0yxiFc9+82/kh6oKZomUODHW0fCJniNKQL0E241nBEHHPjuHTKdXwObInvvVcFxjsA+LRLKpCPWGR8UEzQCMpv8Z3v2aFpOXo7pAkSqbMqLRCiNaLz7J9tOoEM2hyszSgRiPZe/go+UEq9gUup6FMEte18ZJzOxrJHD4EMyEajYgyr6K5pafzxftscfgGTHIcBKaX6HLfBz40RE9DzzQZuiPRMNEIK4NRTMDMqJ6SMMgSU5gwIH6kuSvMuHCzGPQbTf8mKqftEucbIYfhSsZIDKJHYt9hzRWDJfrNmDv5fdDoEg8aoAgUXcnm5NX6YDSKDJYcgeP97kjuy8CV2+BVjxa/uW/RWoh6D+mO42sYqC5Gw2Bk+ZJiAo0p0JwqLCVpFUV9bC6CoI1sY3504yfl4AAZIzJdKSk6+r5FWbibjUL2WlLMHCty6CzzynMY+xgl56W/jILDJsGeskEIuOYxM5+QKWBAN9oUpKZAB2LvU74xGIbXzjdBc+MigfJKZURkwtwHnclLccDhJVkOzQc5ZRJE0J+wkbEV6iyT0xwRMCS7K/qWYbZUywSc0jUG2Lh23ldDQ24wW/NkufrOBZmZnGdzGBnh2qwcqSmk1KKkKZosIYhrZAaYyN7tahkyaOpQUmmDqXWIEOGq2VsbTG5nJpxeq0OEO1FXOGvA+E64qstHPgK6pFxIFNAeq4lmHhq9yLdel2LoBUAi05I5VWl4YvuA1GiyqlyD08txP459EflUP3sPKxY9q2dPLzmj4rgVh6YgIRTR3KStpSNguw0iPWjmZhxODcVryUsrS4m+7y0ShUT2iLCyVFNux5OHvHRl9n1rOzCYfPMBwgDBvI1aRacUhn7q0JajfV6a3xzIMoo1ms2UsE+84alRMz/wdFemWZlAP5JmdFPOrHnYHyV8gtRGWg5WuSDCwGQCPoTwMEVySq5GpTJ5h+bFze18WA6uaIxIeGUgM1M80GuheFprAc0fRu+TdzloqKbMSMFnCTEJvBP/IzGgv34aeUIfn/OsBb5JVbgWxrPLP7ZqmffjE9FuQiGEYzgIIIcD20yRx64ECaeYSUrpwBpBuzafwSvlBealurlpWdIzeqRgNO6jjYAD6M4kFSzFbAnIrMDRA0r55OePiJLXdzffM6cYBp2Fkw/43mBjgOBXotT7W2GGBAjKorQJeGBYYo01ptwiB6uApBT+5HTzBFqMdtr4fiMAu7bPzTsNMGOpS3XSDiV7p9vFaRFdzayac3yXWWbuciCXPhpZrbcEBv5g9Em3kyB7gheGlef4jwAHr3ACIrOsfOJdfNCSGv2J+cmX8SXOkzWsJiYkgAz4OIRKzE7TjIwcrDJyaEgMewu4wKRArzXpXtxsqUsxEIXRuhloRb3tWcpSnZO8ShpTpVSOPM4RHEJvPPFWv7mTZ3DMLDO4p9vA9ofu8+mbJsT3YfGpD/x8VCZGaTEzs3YFH5iZyikkAcsowu2KO2PuwQQnya0AgnXRfDlUmwtw8GLjLlwtkM3MWUvpVmpl9LCh5UJNuBvV5oJUEwOpOn0CecXPZv05N5ueGBvuGRi3Nwa8f82kdGXFPZ0fSDYYLrpWv3yyt3mPB8BoTDOFzMsMK1dgr2Rmjn6UURFWXT1ooFlZj+vt4JbzaBwiq8Iqxr7ToLIk3Q1uYTE4lwU0hwA6LC0Dw2fRr3Do+3D2QVmOqzxo2BzSJt5xVVRN2ssHgdBAIaWS4x4u1ySjamRo12tQM53dzTNFt1QG33czazQIORyeIrqB2Ucl46Wuh8NoKPFIZB94TlgppUgjtBgNZmIfLhjLLMcSMGVEh5XB9MVIJCbMfD2sY99azgwne3VYcbXxPXIOF7L+zdT+CSAzwNwyC6IJ5kOtrVBN+WLRPdMKejNGoHLcy4QuNDPtMi/eI9IEs2EmqerRiwK1HG5QrMNdUaqsiYNKQ9KZvQ+OgGjmpqRRPosNSxj7dlHNC2/rHm5967AWXTQiMFaKAqneF11M/fDwHWZxtNxt8b7bgZ47HKGD7VrZRqhPcbGuZbHwWrwutfi1fhxuOTl2SPQy0a+nG2i+FnKIZCDRaN0NTRP/kuhgXdd1sUJPpxYki0UsZSBEHMY0aKDB3IJOM2iZ10sIBErJkbuZ3KnhwRmTOzgwn8yAWdmtmvxw92iL0UdZIPVzSeTgcslJlt2LsaxVWA/U8ZkvFenroVzvzSs9DTBKPXxCmPGUGly3ETkOkgQza87sWd1kffgoWlnWpaKQ6dTaA16irxTBaLPwMpoTdNucAyfANVtz0iOFrPTCUgrWntR+2VODUTqBlME9NkTI19Vce0auirKs4G4gq2VPO1WzelxW53p78ra+4La+ssWj8XgsMUnSPqyCJupxrWo5NLkf5FMxfFFnLiY5h/lsDuNZdJKlGJs5vbjMB9xjEGma7Wvz0yBD06+U0ys3gnLzSGQmanXWYtkT2IwtJRtmJ4iRzLFmpHw5jJLLGDXZBd4UQ1mLvh/ry5PVwwsvZmWtvmO1PetnUOustVx7UVcge/yvX6Hdoau097fwQIigmikRsWTAjdlFUwJh8FJ8YCFumWaZUeplOMDIJmBtozoecc3sSYKUlLuFImJlLTCvhki5Z6Knnqw0h50PlfJyXIukSBbX4bgsh5sljrfPXn38rLxo9fnxfQ7QCy6b+cPrH23nXcgo08tn8grIGMWOixzRaaB3ulLSD8qrnAGzlwenMqbAyOhlWaqXxTMNkhsEK24apg2AA2bNhwB0cql9Sr9VLd1dPKO6udHMnSWzqHXN9gsIyMwGUmBGK+vz5bC41/Lq2Yub2xcvjp+amhyxBhzouWSH9cfzW92/8eO7f/v9h/sduW1l+gc9Eft8YsijLekTW8ss7zlsQ7GcMCLVOHzkDcoeNNK8FKIMUS6NMKOpFPaA1QbzhHmam0SUQVsro3cjLD1YAMuoA0YjzZAKrzVSiokkmSElM1KB21//4bPnH7+6Pa712aECDgRYSBQ4Qnl+fN33iP2ylTifWfb7L7Zzp/a9aAJhPt9AlTJghgEoMnOwFa97aMNYrybQ2bUvRS2sFPVhHDTSjoygaIJHWrWMFGIXcnzkbnBjhnCEQrbU4UaccIFOYE34NOswg/RIcx+VlgFAGZ+BSRmH7/+t/5jPXt0BIRWmAduq9rj13t5632Dnt4x2OW89+tYVXb5HyqAsQ4UrPdZqSloApRD01lUMeRqSo2WnmRE8n3DhiWe61CAsTBb2S1HaMqqFpZC1KEtcuNKlHUCV/ABBxuyZLKUgnOYdlSn1iVRiYaZKLQKoHgU0U0Qx7TTrh4FysZtXKPc1qkNsH92sgEA/3z/urTUosmfo3C6Pl5Z67CFmb3sxwbUpUub0YtctpY55QAEA6xWOeQJVBqCoxuTOqJM7N9MFC1uYbY8y2doJA918usUyJju+GzNEFAznh+vN6hxdQCctxXEfycOnCfsQtUWfWN8gNGfqQjfL+y9aXdCaYL1H7Jf9K0XvPcSMHpHygS8Miy8pKpAIDZU8CK6T4j3r5VFQTsoZNM1GJHnQ2Ae7/JopAYZ4aga7uQ0Eh8XKwBjtiec1daxmSPeBJzzdubC0MiyHp9DIA8PicMr8RgwAUpUGjLvV1oNt94rL4y7eR0/0vVXFGNEwgY599AqsMpVQ7jGlMxMMYc306siH2SSJkX3jSjX1+e/SUdBR+qzF5r1OaM9M+RRmku6E+ei8KIYkZgKRnKinGWQGp2zkI4SclOnaGtGow3NEf5pZTtf+abWMNRV+c8A59vO5gdofd0HZ8EH5nsphnA+7WOxpyG34DkbBrLOGnFGyKymFT4XvKL451cgmk4HXbnFOBFC61vbmZoCXAhEzx59EEwmBoUdYcl6bDqfcRTPARgOI82qBCzE6Uz46gQPjMmRM9DkD9NJ+0hS9J3j00WFuV1KWIjOB5swW5MWyyRmhgSkWv5IAnMw9dUDmsM8evbSp35vWGWAH+qTaXLvuBJTdCrPtfR3okLx4YrKlBt8QcAlLzqkVQ+8IEjSZTYYB+jRu8hyTKMbn4BT9SlMaqaQEGC3o5bSe/9hnOv/Q05itn4ZwU+fhqIIzoQghB24SFpI0yuHZGfHr5IcBYAy+g+kKAk5QsC1oXNCeCMQDYdI1AkwERzRLml0tyicjCPBMGCbamHgKpdPiYYIsRg7rONISoKEr5pSPAXpN6w3rNF8Odrm7kkjmVnw7QFfdX7GNIbW6MnGTw1XHyuxmoZUKJfj49FYGK2uwDUGfyFQaBGfzGQNmYwVVESmrV9xw3A02iA7XUD8CPYBIw7CWVsIkKqZoGROHuJ5+zSN0UUtazcexAcY0m5JuuUe3EpYhQdrP95cAZHp6D4aI3FN0ZkgMMXuNgLmh+Kz1KhpYaMNnJ8NoBcpYew+asRuhiDhGSA+0R5q5G/vesrihs57yvJM8rL1XLoWAkeeJjU/tnnzMJjDnOAa7uTmZbj7MZgcbMqN4sLrHBqoFqyIfHtIp0CFVY+9Zypax3HIvOxT7vvcCL2btfOTUNrj36JnHkUNbUYJWDvuqxqq9cBCkM7w6M1rVbO0rBIH7TAIv5m5mtcMcIJ4hB6tkWQD1VO977kEavC4HdLMYIzUmr3A42Q6YEpglMK2MNjCKIolrWSRcFRA0Axwo6XXZW5bIhJF96CuMa7k74oKY3dmu7RwF1gjQTRfTFQUbhC2aAdRqaSipMozaphUrh7FiiiyZgJvt02TrTpk9U8sAybWImZHahu+gblAXPTw0wbwUA42WiSzITJjiSZ12xWOf6FpDcI9hiTIP0PS9He7yA7VNlNJ6QWL6oY9DLS/cH2XZIwHeaztHXX2mF+I1eg0KI4boO9WpZGSWTkcmzZE76aVPKslFAHMMUJrsZXiZDUNJZ5BeCJTqXtdlrXen7Q//9RdGM2TPVWPQEaEImAeGRduTmmAmPG1CXMO9JmYnEyo2pB/qAWVSfW/JARmSGrYK2WGx9+hNJfZ976kLYle6X/OWyShBRSqHZzulSMscLLG+uHrWcsVZc8h4WIbAUosPNwEbqRh2kKUWdyu11mK3y2HxUutxuTu8/WeX/d68OJXbiLm+IyPhOShbnOzgD4UTQ8Afw8Y4bazTSDafuvgc1BQrEXBSEbPrm0vAq857V2+t9VhoC0tdrk0xt2vLHnPVqew9biMREVk4hQKbl8LsOe1IYDQoMlspiqRVEebOI305HI/1zrzW6laWldv9u4eveDic//CLnyzmpRbLHkJm4SaFrHidmMrMf55eJTNHRgkoI53mIrEHBdd0giS9Lj2j5wYzKhAgncaOsujd/jBtRhioixPDyva6sxyzRZaMTEX0fL217JemUtGsPEF/hAPmBpxJM3NkXQysdSnr6eZ0qAeWuqzVB0+a3Nyj76/f3PTtknb30aO7IdWst2jN1wsk+VLKJEpwDlGxeQSYGZkYiVfmADgIZHWZkMMoSrQlkW0UQRwlG2lEL677y77fgLTi1kI09esUGo0bH0T2zIxI9R49tW0dsXWUuvdyKL0tiKBZUcKrI29pdVmW8nw5VKuH9XbWCD4foGCkvk746SVO54y9b/vuZoiu/bA/bo9nrgUQSqugubuZD0UIpnTHYjw1lJkpoI3+Mg5TyBWRYgiyJYp6M2kg44PpqlqrzmZLGQpLHnd6jdH/HxzNjBZi73vb295LH5ziZc/CkJfNVl1Ylk3rXen1Jsu6rMtaPrKyLMtiizRIdzYobjYzt+5IOKKCWl4+zy/Kl3+4ubVYnQD668vDOaQ28J4L7VBZV+zmxUhyoxuJ1N6JEFsPOnprSpgXb1luDgF6RXaoewB+iLpxSBBRSt9UV+R+KFVpKI7oYyKK0/KSi23n3rPvbe/5SERP4MGRtrjekbJClmOttSxrfaXjs0XrzQsvbqX4cGSYZKIcDLs5fo4EsRckpr8WnPZrebv8OE7HR6Naaw/7fr70HlFSLMVYL8XKwcNqreZuebUVaGGxN2VvolqXBPPi91oillKGkzhF0uvShjcPyUrBWGy/fXm0WNY3vlJhtsLZWtrrhwu07QmNnnznwPrWYcwUy2SNlP+yruuyLuVZlAO7r9s8mzejltBwoMhMGiAqZoOSV3KhRMF52F7h3buXjz9eivbHy7voe++9RWayOeHFra4FNyCLaViDZCh7sj+ew9uewIju7sWIBTqukFcbs3m8ZgR7RNANqUi6gfHwGOleF6Vy4S5if9iFh8fIKY/LlIZc2higw0LLrFDK/9G9GIm1J2NLv4OyhwZTVELn5P1icp4m89xs9N9dY2zifrk9IQ93b1bP7Xy5RETIFg1PFCm7WevVDjQqSySkiMjeQu3xEhZNHNaokGQ1+8V6O6iqzD6El1oDCo0LM71Uw2l78Jvl/O4kq856QPrBLg/tNtVzcAlgENyGLiAAWHFfrzvgVzWuHxrgi7hP88iKAX17RsTMvzUEWKOH5zL5aKHCyKr1HC25FrXL1tUiuuCegx+fnmGZvbwtay2ljEKj92yXPfrWkBlJjuZ2htNSjPN6WBbIJ4bibhx8jYxCr0vJ9MPNi+8+u//yi/t7VPf7LQ8fPT+mVA5o1/wbg8orkC8LZNXySLqbsRwSoHI8Zk2tgy3Ge13BrtEsykmVjIGhI005UtlZAC/93Z7pxZk9YNX7mOKSKSjTuhTRS9RYF8ExF6C3rfeeY9Dj9PNIdkKWu63bXRLVbdjGuZHukKbZETLbx9/5+OVx+/jmyy83NwVnr3OvBw/aAGjcxmRL0I6mzmJ5GOQlloxBiD/23oXsV3HP3TCCwgWspGF/Eh1ceUIkp9QGEvbl4csLgsuylmXPKFTvaWbDUjPZ0tT19NqBHOn7tVIq0MhDBAIBpgg+tIy1Vs8QvJRSkerg8H/kyueffXbUFv4p87731EFrFRbHsbQ++o9e3Ew23OHpDDqiTuiuHHrIbbh5uskMyJ5izCGbJ2J6jI2kdLkCGA6HY1j/pXI53T907b2ux3JoLcJ9iaQNmnLEuqflpmIGxTCYGEtMM4Nk5k7a1b8TAS9I6rFLmSpMGUvxig7Bi0tWj6fl8+eHEGm/cnN683DpfFPqstR1iRX3DwmtNDMzHK0YzX2phebqRprRWMYIogh4gaDSAFsgNGP0gFFK+YI2rCdhitTQfM7xP2Nhsr1+e2k8nQ6Fx11etj2Jx8sycFBsgOyY2Xq6d1sbTYjombTiKasWKpWX5VjVW8Zw4Ow1zn0PYF2UKEvj4yNr53IpbsePXx2P1fT4Nbw/+81ti8uDn8td/+rBT9DtaVhsjTkxHDn6cXCN1zQDfKmFzqlwGlvZB+imhYyesAuGDdtxqPEQsvdcGejKWIOdt0uHs2UiMvOxtVS21qfGpoI2mIhUWuwD2izcB4SYXo2qzuJK+KAZG6Qlpb5ZFmYaS10KIUVrvh5ubypd2+Obt5f8QTm+sJp7XLhsd7p5N2XsR7rXWnyjuw9PHDc3vgEl5F7cBYATU3wvcs4R+VhpRES4Mn3ifRpFvZ5wC4C8f9xUaOlMQbH18fDKACF1mntxQ0L0YsUcoZiVAEZeqBQC3aw4Z01gSnld1wGmjMzDXMH1+bPbU1GP7eH+3PTF+uxYK8qayv3Wlu/O/vZPOILAcXoS1Mnd/n721rZLK5XDpe1DHpWmVblROJBjOFDEpCRxGJvgyUNAU7MQhZ3rWps5kx4KeVUfStI+ZcUppB0OZTFqz13TIgcZamlDckQvXkblYEQZHnxGKiOSpgLhePv8+almO2+P516O5b6hFEkPjr0bY5Kc8SszyBbSvRTTiH1wnpbFei8rkJZAn1hFuWp8Z7S2q76kd9GgjYDl4Odd6TMicGmdBa3X6g02EW9f187R6IIZ1aEkQ2VdDmZZsg+mqaie2ZLZzAygWSm1lLIUN1Nc4vJscctwW45KQ2RdajFknNs+rhhbbo8k3i247MZ8PltCl5nKHyErGK30ZSm2sKyLRS+rUj4/xOlRMByjc7hGF02nvzFfISqEJH0MW4Um7RwbYNgv7U6ZMkfPUKnHwyUUdKEUR2aApTBR6sFdRcm32XtAMuScQTjQdi9L9Qr0uhjgpXlhsUQeth7RuldGi76/LbW2vV3Kdt6K01ZXTdH7yBrw0dzyj33ftwv4WTb15v5WEbHvvSwRqopwSMNgeh7tzCmyoF+10s5AVYpy+Si1rwpj3BzXusfWTkn5Gi26GZbDij4OU6lVvSMPx4pSDzdH92TbYlbBDInugY4BDDeMKa9gcffD2ijIFjNk29re7ebose8X+YGW8tzfvekGO4OUrYePh9h9BBwjPo+27XvPP8oM0GwzM7NipSjMSRUoYwhVr6I+GKGYw8QPNHM22FV8+bRLxva5PR2WPhGNUrKXtjXVaodWg5WiFynFsq5aDjd3i1mEERVJkXgrDkNdKUX3MdXU7WDq6WZrGSbbS7W27b37adXeL48vqrfWUbwa3clDltNhffHiSgef3X+hHp5Tmb+RPWDFVErx5bgULlXSUje6l4xEDsbEFaDHGObiMGbQy1a8i4peCpDR3NBl6Ofv/uCPHgsjZB442H2pXTSLG3942M1vchcWs9NSEimiOOpN4A1ttTiffSnlZnn7dsBjEcXdqctmkaUU9BVWIJUaPVhoNR4BsezLofSmGu30Kx8fl3CWUqqbJgCzTC8xbTisgNVIVsd72dxVrb3yisqMLTk1O09GPhrzzlRNGBpbEbQ67GpEsK4MofhwkF378PJz0lb5+D10w7Kg3t6c0qaYIBiQ1SSUqDdfFYfg5SbaJjPr6Cb6lbZwNZwbPaSKejgsrEtAPH76eUU70L0Up5NuTh5mb9JH3DYUDdbqdQFmQ0rzGk5NM1RkCMgxA1c5lAiSzKU0DfoK6dsTb+jm2VfCoFUVAI2aRifF6BZE9IFilXWtLA7SSn1xedgzBVUnsdylad+xHNY9YU400ZxWrp2AIVI39Uy4r6iF5aAtjeurzzzjdhY/PvwNnhDYaS40eg54Lwccky4nckxLszFqSlSkZMP+IoeO7WmmOIA1ZUyzZkOSVvzuox/TffTZ6SujWQZiEowEtc51OZBWKnMxDTVJWZhprsXhZsh27lpPN3XpLGptcMhKubK4xh/Q3vSWVpey5eIH5HZYj7d3hXZFevBE/pwLMHgFaVCmvwemy9xYMw0afaSh47cYA4WN4CzDlZLJaNLSZBKs+Ci63Z99fFNrGAE6k6WUBg1PqxzTE1gPhxv35XRcplDV7BJWh722u5H1xX7k4fZUL/ueSOOxLMt6WJeOeQ6phMwU+5b1sDB9qdo23n788Qszzyu34fr98/kMyL7v+ycczkL1ugDfkLEPiwn4GKHNGOYRo4M/EmWly0TYgiyZ8MXQZYLp+PKm1r04nQUR9FIiB4Q5FrceTqdaBKv1SsYRnme07bFvyKgFhpsVh9OibkTaemIptZRafTBSiMrskbQl0fY9+uKyUovvt9/9eMV7cRR/qgHTiExt9/vh5jh4+R/ugA9eIgjPIbK1MR5RBpXZ0wsDZKLRa6YqfTLyez0dquNJSAOvkep79DSSq/xwc3QgR3CZo7rzVrF7dqi3jDA22eUerR1S9MNpkRcHFNMGYRDnaH7ibuq9bFy725rts+89Q1h+c+LU1UOEqABwyMsfr8fT4Xh4eu7yXlyiq5oIo/8Gc5nntdwd2jCkzakVoFlXgSlzqJtrGQTKZGYfJzdUBvfWbuCnm9r3OsYZZCYFKd/SVG9rfzxfAmnmyA6zcvFaqynW4lTE03luRtCWpehSdxnPWFbHYodf+Q6Ro7SbO+DJC+KJhbu9+RLney3PP3q2/vQOuPbP/P0OYsG0556GR0qpgIsA9KOkorBpjhhZbl+9TL+xvudhPT2QYNXIHM3py3FlT7WN1fZ6gkK9Y91phuK9727WN7oPhUoezerNamuoAl62UYQnIhPV7O5u7ZcGxuph1OHmP/r+Za3x4cyxJ5b7fKWWX/3Vf/j6Ifcfvfx8+/Q7rs3rv3MEfvr1ZCoCQJaYxsVXk7I5zxdisfXj529Xyya3DNJoZSSUXosVLy6ywLCOMa45jZzNQNnL5TECxjbsb1iW0+FwqnDQ4LKhcgcJI8zQudzcntGONXYalx9+DAPkP+85HlUXXN48LtuPN7fHL5//Xtx977vP9AsXgNfdNKulJx35KBSv9vXkmsfv/eEXUXZYcQQJKyhVpJda3Efnz/pw4I86NlQWOqVit8tyvgQsSJiX4jd3h3VhTuH2sMEYrq0jpBvL0lPQVsve7/7GZ8ME4Oc9xwn9YT2sX9+/uD2+vInLY55f/9P/cT38wgX4YDeBnNIWCcEpcZmOEKhb/fSzf7WnoC6FmWCOdJCD8D20ojsGjBpTRHnlmB+8VGNXgWi1lno6Vkd2H6s8kxMgJ8nJ1LvckN1u1of1N354eLoAfubrJ69K28uf/JOLtfP6ne8+/zWIePvHX/7CBdA3F5VXU5qZmU3jACnd7ObjF48PqdzT3E2CWa7jLw1r1Aa4L4s71bLHkM6bNExyl+W8d6Xgtdiy+uh/DvecDKUoKWRkBkLyg0e4HU7JX/tbLyaS83OPsvj4xU9+9KPL558c+o//RDy9+s5Hz579+r/HAgyqwQeDCkgIPuGzmTsrrMo//ZW3aQoESnPKzKZt/NTzQ6rLWpi91SkHNLeUmTe3erNcLsyg+QAKRbNuSGC25JCDJWeWF/eT75fLAbXdv/jLf3EofO3nLsGL3A/2xeVXfqPsr+9btsKvfrydvv/df78jwCfl5rwmr/UIQMmUQpE1fvpnfyxs5gH1UYTBpht3KHvvPZd1dTfImVMKOTptDVSpTo8YelRpRHwh0zKGTkCB4arFgFgp2XIr/PCvmrKgm8XPi4IR6/pi3S9/7/Ts7tmx2vndV73if/6nv3gB+OF1Mjxf584AgAGZUUAW5fPvPP86ty2GWoBmk58EoUYqBdb1UJajR4kgh+9QdDCribVUsoeoDgEKZSxUT8t+tfR6asrkvsMq8vaj7fiX/jxAIIn8eQuQy+O/+uc/+sM/2bNv7fSdT//MR792G/eP+7/PNYj3C/Ce3/Kk0J5/klu6591ffveu4u0FcWO+FCNy2j3n5dKjw2m1emQSXgYveUzEGaQ4+WlriFZPKIOf1Rqgniw9nZ2lWGwquizwpbWul7EdPv+baNWBBag/7znq63/5b76Wl/b2/ny239vt9uPPv/+9z+7+PW+BX/wyS8EOzw43j2EHi2shGrOIVGt7aD3eHFebBFm6RzKbuV0d7qa5PVkbrlNhoExJMUhk3SXzYorEcjRu3//68f9wfO/g9/Ne97d/7a/94//Ovvdlefbw1aVvp3evfxellm9tAehk2PHVs1PJNBt8BEwzoEy17dJVD3eHpUgaXy2lZWTS3bhysu+NJoLGjHQFpNRw6RKEnCRbuJJ1cbN8/O2/grRf9P5uc//yX/zO77aWazneHfIPgFDkLz4Cv/DBn/6jpFBvPvvDBt8ih/SI8/1H6NK4nG5Oyq40w/BRoWDsDV6KOOa/FTEzO6FIIeNKM3UzUnTCEKlMlOro67/89P/U7edv/afX9q9+5w9e/uqfvDvfx8WW/mlG7JfLn34B3r/cE/DPai+HLUpKGUZmhJS9c9fp9uboMeZ5q+TIBK1yv9+9Wh1AQwzsIWfiGe1q5jykJmOid8Iuvh4Wtv4u/5MfvjvpF+6Avn4vX20Pb99+8fWbnzyeH95ZOZw+MnxrCzBBQz67PWjLyVo3khlCRqSW47EiQPUm5QIqtj06F7Nse6wLjaYoQ2sXNpLvGGwUBpW0DA+jwcu6rCX3bW9/+7+43GH7hc9R8vnztz9+3d9+8tjO2/2XXzzeP3zde3x7MWCQPOurP/v1v3nnFmXITKg+QORcT0dc+ip52zNzd8u275Gk3/Bhi2HqLog5B+MMp6LZe+yjsxY+7FPqshT19MNf+i+ff334+SXA+9cWp+e3rx9/H7dx6ftHf/PhzdvHx8fLt3cExgIspx8+vvviyIfrhLPH1mHZ9uqH0vbql16iSdGK57Z3LWZHP5gex7wILQWRibw6IPrT0K1BVvJhFs5zLMhlXf/OX3h8/ui/+DHaEWdbT+8+vzmf34nPzzc332N7OFPDLjxEo/49FvIX/qL2P/3ff+/8J8Vz9CdT0dve4vT8uV92X3goewf6KZXR9va8FDdsjxu8LMVHMx8YUpdMnH1Y+y9Xz4hY+Vhfrq/39bTxP/jf/Ke/5PuLaOfH+/N+7pdt70LZzQ3DYBBQ6BeH0///r4rjsXitbZpwdUl0+KXY0reQHWvCTVtG7y10WefbQo8ozqMiOW7+J6cIABgWHtJxy3Kwh9ff/erM0+e//Td/2ffnvtwAmV/1/Xw+770sUIBWesrqtxET/XQyFt8AZEqbMjNSSzxcsrfI80rWwlRcLh12Hyc6rZgQwJjrNIR2o9Kcdcdg/ovI3Hg4HOzNx+3Ns7/1Nw+/7NsbTpbGz6Bo297LmdWhzis+/qc/BesRyUwie08NKpT0Iva9J5hI417c2c8PGxw3+/FQes9BILY2uUJ9QjEzoWQDSDf6Mc8JnBbHZ3/9bz7/+uUv+e4WAKlUMZZyBMoRQO62nIsr5f5L/ryf8TrcONi7D46AUYOxb4xzgxeki1675fndOYnz5bhW9SRAz6T55B1NG71p2DOnzONx6XEq7fLqx8/+k//spX7pHSACBr+agsjznzzc3hZHdY65K3/qV/7+j/a27aQkL5cegpeFXgqzNfXIvqf2vm/b3juN0aO3ORoBCwnkiAMxRcvDjXEcjqhvfhL7Q77+zt/5zz6O/fjLvrt9tD4H0YfGsv7xj/7u4fn3vvfXe9pi2Zdf9if+O0tcD0u1RihamEcOV+AqHMxx7shU9LRsDcUAZmct5m7ZBXXzYYal6ZUwzUuufha/+ecuv/vl13r5+d/5j489yi8uAn7qdYXDcwh+UX6sg12+/tHp/827H/zwB3d/6ufH0P23QvWtsSwxKkPfd7DUDamMFixta6Jj6yZnLXCLyMCD1eIGMGKoPsZPZXAwG3/t//z55Z/+jv7CXzXsWK8DrX6Jz+dqPDW801Cs+KXL2p4PX/8jO9394Dd+zUVclrG02zo8jNIBtFFzTP3Vz3616h8dHrCeIRyWFntvcqPpsOw7lhtlBCv2LdONGbekr6e7w8N0UDgMQakFkYEySgarqI8nwvTsL75A/yu/vS8Y8eyXjlkfCHZIOMqWooGqPS6Kt1//7j9Y7dVv/aUX0UL0um4YEN3F3CuURhi0x887fAaV42EfoiY4Tr0FmH0FzQHsmIaBHEaFxZbDaUEUjkGXmqJ8zJZrB01DviQry4vvLPjFlc8v8SoHM0p9XyKM6nHsrf34d/6fh48++c3fOKKPM5MRB2QY6ReA5j//pJi03h4fdPVG2QSnGbsoumR9WHMM1bTR6/HmsLgWRSZonVe9MWC0SdSA0pjA6XvfL/1byFc/WIA6sKgcwTrzq6ynNSL+4H/+e6fT8x/84C/QSjWrD+6uDB3et9x/9ittWcyWbVpaXoYWlC3ZBeDQh0NhDN2ol+V0OhTCIjITkQNE5fQAHzPpkfJa6H76/gto+cDY6k+/AJE5PH+Q2dHbJ4/JftF+qEs9f/UH//j/dvroe7/y+atyA+SG1e5LdWCMrPxZL4GLyw77tDm4JbJdZU4EF5RRz5AwK74sZbREkBk9sjLf+xgjlEOlYayGaofvUOZPDmnfxgKYjXRjJ5Vm9jZrZW670ntani8PP/m9tRJ/+bu/9vERqdsh9vr5xzB9Xd0WH1wzxqDNFwyjRKO7+zBaoA06z9CQASSi2/B/tq6hJRs+1q5kptvh02Ho8+1FgfLlsjpgti0GFlvY+lsrz5YHErbkrsxopvz7ZuvHf+63vje8IOLnhl/rWG4OroJh93GBeudikUardKImmJK70YczTnYR9AKmco6U6JK5jTEKMlpxqdS771oA6f5L3/8/dwH+6e2r27VWjm5EEUqJzPP9Ea3ZumQd3pp5LvX85l/9f/zwnd/49c9c+8+9f4R6Olr4MDzASS0C5jkSOholWFd6sUHmJJXSQqvGvBopqw/f6tGdBq2WRF1vn4PpcsS3tgD/UFyff/zx7UdZzHJn6Yo0Xx+8HJF7hJXqShyVQVq/f/3Pg8++9+nfWg4FyF6HNftIFEJ0w93na3z0Y8CM0RH7Hmg4IpFR/MHKzam1/nzc78W09+Oal4OUPVD3HY5ddVFEJ2kjDdAet0z+pV6GUdSftmp/vwCe8bi9/v1S7z766Nmh8FyXNff9/EzZEvRnfb8kyZie8k7Efn74w/9vvfv0V//crx7g6HuWYi3ghUDIji/4BWbOEfseRG8XL2Nqt5XiyzLMvVl8I9CZfbtyZgeliQEbpTGN7kavOHB99S2CuHMBmpSx38P/SH7z8uO7Tw83ByuMM0Rztw1WNQOaSGxZrPS3757H11/9s1LLd37tt35wAgAzDSzf33316Md3Yy6cAX3PuFwWK7VehRu1XME/M1ffWvaMHvIxBlMgO51EqtMJRRr2Q54+/bafH6XT6hjXgnz79vfNjs9fvfroxSnMqNhGh2fOFDESF6CsS8SPy1JxuY+3v/t37eZ7f/57n3ldBwHx3e/+4z/c69O0KhLZkzX3zdy5cnBYR6NvEGL3VLLvHWUZoDcJIWJMWxr6aVl3vfrk218Av/pl0jwTOsQXP+ZyWH9Yj6fjuhiy78FaJzSDY0YaI58rA6xrydgev/4X9fLy8+9+9snHR3zxO//wd798+4g5Q/ccLFxL8XbZWzvHYbC4h52EuS+rtQZqa3tDzcqhSR7KbJAVNPNiMGP9/reZBM8FwIQtPC8AluUNSWvo/1093L189eJgh+MzRaJMqczaxgD2LsGsWLSWhgch3vxOlNub//D17/3u77+T9fH9SJixHNaoS8bl4XIZPR6EGWHub9ciuePcWyCiHQYxr4c0hmoOhRpb4PLye8M9/ltdgOFQarYvz0uct92M6u0M3++//NGylI+fffR8deMyXYJawNyNF7ohL3lgNUdE3/bIkvt/nfc//uPz7YsxgSKxBM1KQW1GL+v+NjMzyT6cyK1dfDT6lbIM9DJUczH4hplmoCKJuh+/n7+wCfpLL0CRlBlgawkW38ZkXCaN3Dd95XW5+eizFycrxc3Qh3RVa0aTeW2jBVaSLLQ4v3v91f2eD+3aA3ZhPa7M5eHNRcb10/182RKMIT4DxKCh3wyP8TnpSWN0caRWK7UQOB7Xy7Pv2p8er/ipl0cpQvEsY85wv/pbGUc3wtS3hy9/9M//2Z9s0TucRqXGJG1SWTJCinZY0fwUP3n95uHy9vVmGRH03OjLWg16SC+EcDgd3RfXCtZCX21ZCuTHKi4uVgeLQ1jUuryU1hPRUXxf+l/9G9/24wNlKcpIoQFXUdigxU9q4HUy0/ZHP4bfvrr75Hh7qsp+sVpckftw87av1kqcH97tl8eH89bLnaPviP04pEWD8iJoK2KtWjvlMPpK45Hk0lnQgK1aoUHnTCJDNwGl0WG9ftL9Wy2FAaDcMj2p2OcffDCrGFeCJABYNuDxzfpPj89evHh2u34cve0iDsMsxo9ePc5Y317uH867P26LQyxe6MUtmZlmltiTcEONkCNYUizVkFXm2htqcapnisUdIc8ecEt5PP9z8S2A1j+9ABVYaNCDMiPn3TO+xg8okugsztw3PHzpy/G0fuf0/O5gUDOI5bCc9rY/vHnMtl8uHbnf3x5d866rVPqYfiMFpU5UwugoIt1JmZWSZgEbY1iG1wiZOUyV4PbsB/atbwCUyOEndqvo0WNo4jWVUSByUoE8e6dRN8bs717jX/rx7tWnnz6LZa3mpWS/XN69PefjpfWE+uUAoqfc4SKw7hFEdgWhhEcASGbxokBZrBYmsV8SMKM35nDIZzGzWgA/nAq+jfblNxfgdtuzS1xFy5rjTgglrlMZxsRplOFrxIfx19yg88OG5bQcS27bYyrS6uP564etJ6W3XmpK4aCZD2vezOgQHRg0GIK21cyg8aZadmRWIhKQUUpz3zTmEKMt/avPvvUNgPLdN29aih5TroxhYDlnc72fzLEPBJYlU3Y4HT558dmrUzEriN638/5YtO2Xd1++aT2FxL5tcCPc3YdiOkUzRhJqc1IZgEiar8e1FCbLMHOQhNVGw/0ikxu9bOubf/HZtweFPS3Ai/7YJC/tOh6ENIfUkTmHpwGAKgQ/3Bxw+/Kj57eHujK3bSv18XLurO5L2fc94nyeUfTuWChFK8yIzKggrURGArEjGegUYwn4ze1peTSEnHrQsIrY3ZQplZii6b2//X18+1Gw3IGvH3tebOFll2ExRENdH2hJkZspy2m1H56e3R7XQ21FDSW3++w9EijvzvXQL2mXvV++/vLSbMyi8o8Pb/ZSnrkvSzWiH2pru62Hvm3d3PcOMswf1uNhtYyTULLtrbaeli0Z5mTm6pmktu3zt5e3P3m1xQlQ/xbxAJaTjvu2Z+stvK6vl3XxnpuVMZz4cPP8xasXB+xmAPamyD3Y96bhfvga62Kt6cjob79+t03fex4Oh0siIwoJJ1Is9DX3XmvrmYcGp1lZ62Gty1K9t55yXmBmQ3kDCLqHFS9reczbr/+nv/aqrDjbWi/2bW2Fklxc7fH8tvFGl/N+g7h0lnJ3evHi2elQbks11/nB2t5zDAcImeJMEplJL9lws75Tf/j66/MUU1s5HI67do1ZX8NcX1WJ6C2ibbv3JGTlzksptTqnxJUmTOOgkYkRY9zNuzv8zn//55/vQEX+8h2hn7sAEmvh8ZKxU9vF3Je7m+fPTy+Ww1rdzGO79HZ+uNv3tEI9+jD/DDOo94BayMVzv399v9OHe5HXpdY1siWn77ANojG8d1hu+zBOMDtxJAIN7rHvuwOZmuMhABrdieT6ENj+h1//21xir/zWTgAKlOR6uFvrVygvl9vTzatXd6tzRez3LRWIboXr1pssMguh6CNDidbi4Ovi2t9uj1+/bUZGpgFWT9WXvm/THk52ter2VlgsuisTtDnImbRERuY0LgCvKq7mZnWpjuNjnL78+8e/AX/k+u2VxSWBDC3l7vbzsi6H411dCvre7tVbyG0vQNv6bplz1JdTokOZPWEn5LbUmnG5f+zIVMoo86PD1w29zsl+17mLMpFGXxzBwmwZIpRlO29dtlxAEm5Th4HqtRgS/ZPLTxb+3t/94rd/cIefTwr/5RfguN7I6lqXF74uzIh4EAn0bHtL+nbwtkeGEYQVf0R1yEqmxLow47yXU4ltb60NMRwtYZ4N6+myD1e7RAx/dZnGnIkwh1X0Q4aYHUD2fnW5lhwa43rX4U9I3Pd8eNl/9/d/76/9xqeIb+8W+MhMYimwbLvFbmYW+2VvGlaYCl0u6XYhMm2pVYa+MT0i6aVoPTV4tK1lRg6VmInsO33NU7k3QUqrUmaGDl7N6LaB5o6cuL96siAibdAHfGRO4LkcilW3/kd22i8vz8t//w//9n/+nfzWGiP+v63uy2Fx0GtxqxmtJcB07VELcvHQUgDGtmM5ojgEpGX0QK1NXYbL2/3rr9pScw9DQsfPDodVjdi17eXAtquLZsNFIRVdOaapwozZY2QVitint1IXrdR1qdURkZn7u+20hgq//J0/4EdH9DS0cBG9i9PL4pevFPxve12WWuzEdj7viZ6C+rYJGVaX0q2fL314GIpjMB5oxTPS3DTkfJfznzw29V01qntROR2rJ5idsmqCbTO93gVFH2OtMft/vfc471trbe8to/ce8TiMH1HdvXix2Lsfj4Vfp/d//Y/arwct6a6dZm7Ryhij+EunikWCMaWGcjO8vSBF2x7Zd62r49LPG4tfp7oRIg143JqfDgMYUuyX89bHQA6aOYw+OA6le27NSulWihuxl6pAKWY1ww2u7PveIoasQPswuscBYKmL50wjtg5TtPjsJ1+hfvHf/N7/5c849orHmwGflS74/5rsqERvph6MshTF3iNCiB4GQM1tywxkmrJtvcILRpGw753miEh47ufz+nCpN7jfkDB6XerouNNL7MGlMrOTQF+6di2LW+/FjamM1vqwnwohicwc436zi8uQ36OhVEfq4ZX/uN28+fs//t/9p58U4iaChUjAh/XmL0sbK7kHYktfL70nbLC7ET2yZ1cgOoptnY2iUb1mZmTk0btb77vLGJfHS0vz4aynMFkpNoZvVUOE4smxvUXJHZFuvrsbu5SRoYui9cjhhwcg6e7T0InISPOlEvziWC9xe9P+5eu/d7r7s3/9L7gP5l9zeNgvnyGV/dLZL6qXMemL6iGi93eZGe6LgTUuTbYuDHlepid+uhcTfDdnnM/7T3C0y65jeg76VZLMHLNMs+1jYg65dc8YfCCbBrWSoEdkbx1o01cfZqUUssDMkH2zmt1kxzfvltv18fLm3x6/8/JH/+D2/Plf+e3P0BPAuTv73S+7AK83LWriQ1mqZU+ElO2yMVNKGTqFjLLWMX1zy6RR+bAcGGJ5qGb9con0otitgs6Sxk53krYTdOzb+GjIbMMgrE5x+9WXoI/EmmOUtFmL3spSrQpA9r1VKRqMfHHbvuo3z98pme/ePvsX/+L/cXz2vc//cwRv/lfcjeXr5rIuO/SHlixGCfH47vHUMxv2WHceIBbPx3334zEj4cPnWtvOEtWYrefd3lmW9mhmpWbxzuEAZsa0CF0zu2VOmh59I2GZsEuneoshCoCZ3WbA61pCssjWal2ql1Lu8PZcT1G/bHyT61JeZ8l3f/hvXv03681nv/Vbr/r6yy7AVvOLvtTWjjcLZN6jxOPxxbvt8nDhacW7x9rizfn48WMu/e3jqwXmOp97rtz3zIdS8fXXzdrrsjhL1zGBS677K1NLlr5eLp1+Wh5b6wJQAUU0CEYz8muvzuy57peOUuxxjLvD/VId6qoN5tvDw8GHJ+pPKn1336MlKCtwBg/2Vpd3X/7z/1e1z3/7rywZhQ9YKjoKUo4eVrOp+hx49gFbsLxeqh9vDnwwbNFTtztiM7WAl9zi3RanSji+xrosxe1+a4oWvGnn1kzp6PveE26KzGEnIQCX6mYIbKoFrfVj3VpEZJu/2JEj742ElLoI7sw42hjee3D3YWVtuDyc2zNnRhrXOpibi1pu5l58Js9A28z5x//o/1rvfvAf/uYNeoeIzRwoRV8fDoOjAgCKvLaYyg+ogB7a2nrrrfU/2Q1xKGeplswMlngdW1y+o/5lt3VHTyvFa8bDzhLtNtr5srVwQ/YcNuY0N9CQ2TaYqcPK4lZ6b/39JKOhQl8y5yRic2fmGJRptDGaJ+AWj/cblutYG7FAiRIRj006VBQgaVOniId89vYf/P3jwye/9Vf/7NLePgN6X4wvgexppoTVD5iuZemtZ+z7m95las1Uc1uW3XwQGhz7w46KtyU3u7n1t1yrgqXGFgYSfd/avjfQxuQDM9KLRQyqZMneWpYaVj17i46racCIgjEtNtYMpbm3YThgklhM2hfr5zPWOWBNe3fMoVUZbTqsDvCEBLHc/+TtspwO9ua//fvPXt4d/uLnn5SCfLc0HX9WY7H82+0S2Xv0xLJ6b9Xr+fJw6UFf1uNabVnLTal4c3eQ1Xx8JsXljOVV9CAKXKOcFQw2uEKZjtwV9LKWbcjnLWgG92g5hl0EMEZvsQwqTSZG+jNKwpIqpUAbcr+0uoTDa0FqVIoSHMq22QJlcZsz0/hwuEnh9R+VUvOrn0j/5KbYZ3/lt57jCKCh3hdHvrfmQvmjFoRoH8vWg0drXv1hM/TYz2U91OHkIL5auO2P2+Pzvl8e2uH2dYQCS6l7ikbvk9CWpkzP/swQEYgOW3xvrQweQjqRw6vyKkg3IDLDvBZkrAMxF0hkkhX7eYvFx92QeYhhelYkIpiPPfu6LMWe3D7bgz1/9eLR0cIKvz57+4N/UmW//h/9sKTqLYAP2yvlDdbTwdfDvZQ7My2isVrWiGixma/H0ndWZE8W5+8C2fy4pLkkWy0ul72LbsgQXe4CoTFiL3uJBL3lMPOCzSHCTy5XRlIGFA7eDGmYEx7Uw43t8thZa3VmCrAuo2AlFGNAVrT9EMtxiA3lwTUfHrCmAhbtee6hLL7/D79zKHd//td/s4dXvtdFlB+g1tz6fhs9JfHwsJ0b4kwYMkP++NZbO9xekqa95cvlUFGev3woHruOt1s7n7cWckNE1qLqYCl+X0xe0P2yd8LVh/XDHJphvM6pG90WoPbezI3d3EC3lYwQrG2PW5S6LKYxMquBTqOjNwiZXWN63RCV0Ha/sUsuLTNraW1vKp6bL4+PC//wd//eS97+4Dd/9eYpHPB/T6rvWz90Ves9he1deF6yxLnXgkczSOBRMDfpdKjym5vypizOlHR+d3/uMTRSIZ5gFnz2alksxdES662nZQ6rUfXLJcwLhu3aKTMSwBrp1aLXyZ5MmJsibLt/6Ifnt7541+G0RLcDz4/hRLu05Jo9D6eaH6/Hgp5aAMHMxiwjiNN7RGOyEfuyQOXm9Jt/8ddtZ0VBZvT2vzR1NkuO5Ep2PsfdgSCZWV19e+7odyszLeb9X0ALbbXRC4xkGo2uVXffrsokGQG4Hy0QbJtcVa2YTIIIwP3495Uyl+IiKdKt/RGaJXFbQGUa1mOqXXzP+ah7zB5mxM7w51jzfjQmVGU1PVFLEK6XbuwcUxKDWJwdkk+Q3d3up9kcIM2Idn5RPvc9vfdmxuab55e93K4d3wmPnJlBzkP4TF5jo+3mqJzcKudckNjTH7uQWM6qmrX/t/9xwdt/+af/HD0H4I5x8vZvwDNz7rlYtwyyZOATFk6Vd5/s9jQC3tyOTLPMTtWcRTul9MiTs6NTjbU2fQmSNa9aISmzKZVyAfAk8FRwcMXEKx/HXq0FCVhrGveh8koR8E2Zz0tjHbDPqexmqGbIkUqI68xNlCTDKrUrS4nx+Pyjuf3zf99i8QnBaR6gybGECg1lUp5bC/EmWDhRsIg3zIVkd+rYIfltiZFIsQTiFQFabeZzIgxSmcSo+cIwslVmSngvoaBaU2OLaVo5ZqboViMMyoFKi+7zOZbzAn4s3LEOcFiO/NqvzTeznctPP885bAoqGgZpQOUvZTa//8b4eLkgSSmLdxzJ5vzMVcxOUqJ55EptKClDDl9KVihHDXmnKSX3CRSaZpxwvpduTn8KPNf/bcEVkJAtm6kA0Qu2JB5GVNXMosdyQdVRzWRI5TxaCRYbo3Z0KL0mcDyO/xXb9fZ+7W/umKNkhGQRg5LEZVoCcJcF+9sWeba4IGAeZeCqeOQxzWGxnKFheyZLqlsN1eN4vC+TpFI1Dlihlp20uab17liyrD9jbWerr7A0ciegk+tcKGECBgt72hLtTAiVmbsswiPgkQmYaU9yzu5aceuaUs15qUKY9X/dB/q123+63C4tWpijZlVKa42aO0rg3oM5n9/jBtQpMs05xSZqVj7bBECPPJNCXfLWnH3ispldblkeqPyYmgqvmXOmgM3LLHyBqVF/HlyXQEtajRha2UIr2anZiXUG9Dj7J1lVNY45I3qYuQQ3w7gxJ821hzNL2qymqmaNElvTP/3+64ceP+a3y/vt+vZ+wXb1yqVEFYrlUJb71Y4nb9eKSyVSQJQKoj05Ega/ZNHxurMARQc1dSvb3rnBR1I5xo+9xBYEaS7jUTWsH3XRsoj4iezmaxvEcr8W1pDIaoy9JqOqtjU1KasclSC9Bdd+6GHjqeZiPj5wuRoJ692OPQuVGHlr9ks9ZlWmPfdv5ZfNvv7154sHy5cmY9CqRP+4xKW9jY+Io6xVQubkMWbTnHDVbyCL+aj3vlfjEVk1ZuvTbtf5kV13vPWpSwGJuFqMFah/k8Lc2gMWKnZXqm3zCSegzCPLNs29c92CPI+Et7C1KRbDRe17L9v2x/Rtq+PagUPV29ZC3VLmPgY8uo8SomiOkOEoxu0/3P74tO4CMWi6/2+096+3f983K5nfNFNGXTXruVvEbwNhgtWcSm0dNUKA+pkhBh53vG1dtAi29ga777tb7/K8PwGLzkbNmQLNDkhGezHJbVAqox8vlyQlnsGvxQixkwsDWritgeHSY0pjFC89vTeGGWsAeY3lX/d6ojUdmllJs63KPIwPtusYEwNmxDwGDcfv3+1/Xr+8v3/56ba7R9W+B9jCpKCvA1qbKfPWDyyQB6u03A9jLnO2udP9/vU9LYB7NrPgPHKU1Zg5SzBbpy87xQKg1XpPCwVa+ZK1LogYVHPNkylWkJLJ0gL4p8o8wuhhsq15TRpRY5cdxyNaC9uf20whzCQVp0G9vWvOKpBWGulgzqzrc//GuGz/brtdL1uLJ6gcxxHNaLVGloxzDs4pALZXyo02LW4eY9/sxOr2yx7OcUf3VhqPHLPhFAMDaGfdhacSLIJjlrJJqhzTuaSBdR4LXw+IUC0EakGwNVpY5k0lV9bcms2iWyoLqHx0lRq3GFkG8pAoL3WhvefImRIdBpbKmoXmmPjwf4nt+vbly/Vn72bm19jFkMzEBfJ/l+zVsoIlCxZuMjM3SPrLZf/183LlxkZBedkP9Sv883xDtYYdtGm5qdZTXTDl4hCa1ZLzrLqcgxINuU4Ga6RQqqpxpHF6gmsQZREUDguald88f+j61i0JqZiimRs4q29fy+xxTIPZeT7zB+lOt+cc32E9fr59uV1uX27RC8ipZXNXtF5RlZmNmKgESyJ7T3cTLOop9VYjXYlSUaI35ZkD4Qm7WNm+5Zgts1oSXHMuNjf+zKPGykNw96C9npcLRiytC160KVcqoaxJ86Apum8Dj/tPRdWBGjC4tN7LFd2sZq2hHGWOY8077Hozao79+Xe/Rly+3OKSlZA4igB8lemrNKvWuCegkR7FKsFi7vEW+HyMbmFUPp4HoTyyRJpZo7BYwCsTY1XlREWpJENyHUSWX5gvByTdABkxQaDmMDNj7yxaBBCsTKtZOIUcAy26ng+HKlNxwAsigtSM7rTaq2oEYSrgiTWq8t3CYcZb5eDnb4xHCWW+fbNwzZxtP0oS7xLdPHYL5E6n1gzwX/enSvvh5qY8HmNWmM2Tjxl2xotPwKmqrRPjccNJ7idl9WdhlABdVfV2wiqXHzmnuUB31gJJmRcA1rw5ZOY+P+56//nWvztzauVaVVRz58yIr/OwOY5as8thvUoE8KZKWUQewy6Rw//ri/REQinTmDJWpkdzVn45ngiT3ByI6+WPvz8d90cV3ZEJHvPyxmOWhBzyixnNIy4tjABTdBZieaPJUfOYMjuWsT0Xdcu8sUZiPmfOnPvQ0NjRuz/VulkoVEW//vREpWp//r5PAAzvQd9sPoRxf+TZvh/3ircL55gXgmaVI8uMJ5kLNWjNNYvxXc5KRMpCyly2eHDLORktjiJB5z+4q2JrPSpy3/W1b1FVKJkhx7Gp6DAbp1SyzHx1BghV6SDjrH2cmtMVTD0vCkyzACvPg7Nds/r1zWcTcrrTEA7TvGGU4lbzeIz549N+qZz0ULd1nDyUFi1ocLO2PQdKok4J/AvK8aLzAfHpm47DgCyjYbY10mV1kjOPMgoeP1qAzPnFxN22MDONfWrCw0tCpTmYbh5mZmYeazJ7NcXSHKiZXC5vNpx84PNPUGak5voAYGSVb30eZR4G1Qm5n7JMlvCXvd+PeZR6M8ecN9CDmXJUBoywuFwvj8dcHanXKMDrcnYi/xFXEFCOqCo38w6lQDvoTmRqFRVIM/rluv2tzPeHt1Nl7nU2AibMBHr3Fmbuvg5v6wqExbtYOnXQCmg6cWnnEnhtHidjFUPWGypnb5dFoC0YUR/dzTR2l29emffp9GCC0RDuzZVEViu4X+Tb8TxKwJGLOqk/NUJrm46f96GAvFAJGcx92Z9ey9NDRdVb4xQ0JkxibNVC5gJNmplQ65jVruEGmEc3X9NFKlot/DLc8AGeq0tamOTFFXCrCeO6RRIs366t0JaKIzwzEUGIBrNoH5lwd1ntT1wv/kc0Z2wbYMvIWpVB78oljvazOHXa9P5kcMc4dlBElgpVmFk5xsioWWwtnkaByM+G0bfiF5lyV6Gych4JA+YYc3osGeSJUXaPMKieEqqE47TCLlsBkMvisnqVp29d/ppRAEZEaCQvjlG2bcecpAnBnA7bHji/3BEXMfeSkHa5jK1F60ZAmsTySkmy03P7Z0NgvRLjxyyovA1FUwkrrbNG2nkOENCZNHh//9l/S+T9k/PqgzOPMat1l2tMUdQ0evMIM/M4n9GpXPKS9fJ/ij9M/+aTWEGJkx12BrFRSaMZ8mAVyKoUOcs4fdVWavhG5TMbBWLUo0f0WxNEFYSoXshUgb6+RqfY4LUHTHYcRab1XsllaTfnD7Nwq0NFuSO3mMqcuqzz3/VwhUXk7ta7YyxtYUBsW1voAwhV1zwIQZsyxxhr9VMvar1e9aKFWa46M8Xa4BECHtd2wfO5v4FGzQPRaTaPfc4USl/0eOJyuT6wBCfIOka1VRBNvVpta/sp/ts/gCAgLsnoY4w+92l2bR8I45Iw61AEwdCIi93r+n5847G1mr1vnVAq7JCbnn88bzXhHm1zJMKwtR6oktI2TPOaWWjdHyq3VUAmVHVdu0FZFSsLVmw6hm522SqTW89HuyCHc4y+Xd70PHgRnKhk+AG23IcTKkS3kHl9wp/b1nqbu8yOB5Dx4hK9atUx6qaHAkphJNK2rT5/5KXmnMc+Q8UwopIUk0/ouX3xvbvGIddhRuR+7E73drEUzQn5+YOzGjxJ0cRmc6Jqfc0pbDkLLaxOGVStHRB5fgE+mgWsRzWvAaMH3YK5m19ae56qO1Wuf9QSYxK3WSX07VFPu9xaf5iHl7jbopGcDiXBiIS3qEwRblV7olp3JEqAZllHZYRTdLcYGe/t+8Gak+SXqqzj437r18g5c9Yq5nm4R3jYWQ3ONYODtRqVqFoaYy6RwDi19/vZLFz+R1j0ruOhdnHLaRFpdI+ZOyXzHKtRpiLpqxJxVvAFrKjBUrEfLmMJw2MVxAlQomhZaO+hhdGJI+ckoo8z4A8DC5CZsUheyrcrRq2dVPoUvL3H23652DHGQDFlbiv4Ygxfjw/QZA013EDVnOdd4aA3krn4G3XOJglbgqLbZevz2FHR1oVwzTMptzyGkKcStlBGmsBPizBICW/KPLqzrO4HAlVT5hH80wuwRrWDimu8FHfFHtofjzZnZqnCWCWLNd4ntTtv8WOUU1I9yy0Cs6aWsJgDNcybvzSMiwUBwYoODW+cZXlaaFeilOfgRR5zthfnbDm5rFSJ5k10K9OSvsLsljSn+VwrQMMMEllcJoYFSQMvYSIGNFfTo5svkDmghFIycrIHsKSyngnJvf48lZ12pRNldxyXa39ONwGac1wbcnz8ePzjJbKwGOqnZPr8OU18tDLWXGmgmn6usKgqs/BhJhj0siOs+NgKC8iDzQ1kVtBAWBzeSVatnrCSEArgmy2d5zbqKLSunIqtXz4XxShW6xMy8PWwJWERM0FmVa+U5G3IvFCrekcUzVAg7f39ullIKnPQf7jVTP+peY2R82hueD1aJGjInLaoaeREFlL+KpiYfKm2uAFJb/a6plhB6d6vnlXIeaFAp569kSrdmwYKbV+VxWpEEeZjPVV0BWhYShODsY1Ronkt+6kVAXMRSBjocSyFOSqubXx8nz9lznHsxyXlLVATS/Fxub27OneWeUfEf3z+/rff1W/INpNtgsY6oaQkMFhu5yQyaTRj8MLItAXXaC2YM1slYF2DJI1IQUJsxZxwK9RkuNcnw5RzMueD3NqPRd2UpCwKjeaksCMaihWmtDom86gVpBcI94eZwZymhESPT/esEtg0Kv7ieBxq7j7CWcPDQYYZ55AfT14eyOfd3jf/8f9+CP2Ljam5H9kj9JyF0TdNGG+oAxRYpVHVm3SlkhGVCfNwKu3qnEfWLPUVm5CNYpjZ23MiMI+3Uh2HFSvDRZ/ynwxZlhm9CiRDqqqop1+8RtdHxTY1KjZFS2XrlTbWFoIvlVkyR73b43n8Pbo7VMJcCvZ8CBzHEFXEeo4TZmzb5twEwjr8sn3w/cj29sbj83EgwmpaeCORq/tD2WlJsajSLp3rX74Eji7JTMhYTTJaOJRRlTJzVsmt/31rFqBF3Y/eqAmUsSTVRM06zdciLKicW6OmE1uIpsn78N6nAF8mvH2Fr2s7ZkSrR6ziANfgQsKefcsUGjTNDTJbPbvrT+/I0LirSnXIvvh+9ItPtqw58zolF8tqwqCaZjqzcLQVFaGZlr/Qwp8nMXrJA4wGRjh0eOUT1jLhpLBtHRZUHElh9RqCqCnVUSuZt9qsNC6kf0nSb+5GMxsySu5pROVc6LysamPQWn7Ga8fmK0HhVgDZhsqkWr+n+Xbx+9iwZxWoMX5+Y9u+fvEyozBFB7PorIlYQWD7s0xeWSFa8Hwlc29mrLMrbWYWgDsLRhZo7XEoXNNuPapAVPQWi0LobtCs5e84s1aEihDskJVFMC4LazT7MSeieY6jxNDrsUaCyBl2QstzZVN80z6mqrqXTphflTmbzfshPv+aBTPJLXX7h69Hm/djz96nuXIuaFpVIc1WHmZ9QNUFM1Y1GkSzGyHNWY5C0V2gGyutZhnNNJ8etPYoL7Vmat1U5vupgixYgK/u+6pomVlfVZ2Gy7XlhGZslrNyNLlEt+eJKJrWmAVE1sLZnKPjdjlG0qsKJhQX69ra5esNF+tWv2SKJO/5HNH16L5ljbXEJNigoZJMFE2qBogmhEQD2VdY2WyRkhbc18JAC2ctNo3c+vd7Rrtt1mJkFYI1k1uvE/IalbBXlWOZPCjRF2QixcqUZ1U05bE/jqXtrL54Jsp1gmZM2mrFADXliHOA/jiVkLEsobcvN/oMbaqxT5l95K+ft6M+3j36BuEKpiF3b4SmWWIFGhepyuyFIlglBjvbhKZRIo1sS7qbxxhprbf73HPsWV96j6MqF8sixxTdTLrMo0Kll1aSgFOTCYNJ1SLoPR/TzNDseRoDFTRfA3wquCmcJtBsFjJLOiprztQRYYJ5oWYGaPBroXg/Hh/3Sfumz2oZ3ub9+yfeb7djlHHUm0xVZQkrqGqeKd7S0kAVSZB0CXbOHcCcHTRUzsfzEFvvzyHD/HHsX94tRR3BzTU/Pxl9s0C3LHFl01kGwpyjINDa5kkdh43j4wqVb4xtPPch40dEoGZ6VblpxuVEVoGO/QFjCKPgHcv6+vTbTzMvv+QfCdZx/Doej+eRy1c/2naHorO1MkpBs3yCqGrBHW6LjnQKpcC2qsGSQtQcCPc1XRxjbrf69du/9LcLxnP5PE35w3ovtcsvc4I0b2WRT277cckf2/XHu6qkia3XPa56RPn1avJ3plkwFS7CLQ+9afDC/dNvF9Zj1eci6uxp3oQyg4/KWVojsADQ87GnHeN5PPb9ed8fVQU4p0Qz53Wo0aJcsCbgiTX9n1x110k3IxBnPS7JV6WYpHRw2eCf12v++re//ci3TiB42lKF/VPdsH8PR8riOtGyxMs9vZksX8nbEmsHvhx5yCN2J+o4RjeauXtqLjZKszombavMylK8AK1fT6euzyyaEufpgJqUPf74Np77vh8vjaUIqOakH5PXQFUWQpVzhLEmQRXPq4CdagwRjDM0mAI9C7d1WYDrt1//z7eD71sHetPESTJ6iFubx/O2YU9rv2bcvLfbw3YLDcxzLNsKocF4DIHGOLZgzlFPpzw0e4xntUpFH8/d2+bGoVSsh3XpU0QK1QU2G6qTB5jyhnp+u1dm1hpik0oBaNzJ4zjMAjUyYajJQ9o0ZjunGLUk5qTVqyW04vFaqOk4aUWPb//31wcYP39pg60+7+OsD9ejXZu5e8gKuhwYx+F7D6nNipXCMFBsmXZ5oF0292W3GBVXVFqAsxnYtA+ah4x3dwGmiJOdNWVulumGdVU/JdeiO+c4lpGZNmnmWAosjN0PuBMW40XietzTm2RpIKTJ1R820ERynFrkFemV9mhWz8fxz398lMF44Tja/Pz1+VwbplD7p18uVgfpmBsXk8mvj2yc8QCUJiCXbvwL2uYiXGPu07sDBdIPk0jBNeWoqRItUvE0M3fjTBhTrFKpznIqiKjaue4RmALairQANAZFRqMKFBwG8fkcaNfw8/OTvVYSYK/hySVClmjS1Xx+fPv243fRE/1rVBHH5+e+IqYQOD4teh7phqoa1gx+/ex2l+aaFZE4pEGzDE26bdtz5DzKfGemGfPC3CbgmFlCzgvKjMg4uGzlnpmSWZYW3WJ910lTwgxwwxRjVVAgB81yaCOgFA+QYmvD8Pw+bu84C/yOdRfOc+m/lvyKoZja/sfvv/3+fZpUwO0fCWfe72cXFQJQx77HquoSVwPGve6Xt8vf6akVtiPKMaPZsbObCWDNBMgrMi2YmgrsafKzJ5czuteMEIpJYR8i3JfXCnZGflC2ylxbNzsGsEZgjBSpkW45DZU5zCkL733Pz8GrnWC6BMsWHHvJQs+HwOte8Psf//q3+7L6uflG0O348TGwklar9DuPuDlzpPkf6Rc2bNvt4zePYp3bbEVNNAxvtxvnUQEYyyyRQ64Re0U+sy1qeWv7qB5O/n/wsav4xlTw1AAAAABJRU5ErkJggg==",
"text/plain": [
"256×256 Array{Gray{N0f8},2} with eltype Gray{Normed{UInt8,8}}:\n",
" Gray{N0f8}(0.424) Gray{N0f8}(0.451) … Gray{N0f8}(0.808)\n",
" Gray{N0f8}(0.435) Gray{N0f8}(0.439) Gray{N0f8}(0.831)\n",
" Gray{N0f8}(0.443) Gray{N0f8}(0.427) Gray{N0f8}(0.855)\n",
" Gray{N0f8}(0.443) Gray{N0f8}(0.439) Gray{N0f8}(0.906)\n",
" Gray{N0f8}(0.435) Gray{N0f8}(0.459) Gray{N0f8}(0.886)\n",
" Gray{N0f8}(0.424) Gray{N0f8}(0.463) … Gray{N0f8}(0.89) \n",
" Gray{N0f8}(0.42) Gray{N0f8}(0.447) Gray{N0f8}(0.898)\n",
" Gray{N0f8}(0.42) Gray{N0f8}(0.427) Gray{N0f8}(0.757)\n",
" Gray{N0f8}(0.424) Gray{N0f8}(0.416) Gray{N0f8}(0.565)\n",
" Gray{N0f8}(0.439) Gray{N0f8}(0.431) Gray{N0f8}(0.533)\n",
" Gray{N0f8}(0.443) Gray{N0f8}(0.439) … Gray{N0f8}(0.525)\n",
" Gray{N0f8}(0.427) Gray{N0f8}(0.431) Gray{N0f8}(0.522)\n",
" Gray{N0f8}(0.427) Gray{N0f8}(0.435) Gray{N0f8}(0.643)\n",
" ⋮ ⋱ ⋮ \n",
" Gray{N0f8}(0.173) Gray{N0f8}(0.18) Gray{N0f8}(0.02) \n",
" Gray{N0f8}(0.161) Gray{N0f8}(0.165) … Gray{N0f8}(0.035)\n",
" Gray{N0f8}(0.145) Gray{N0f8}(0.149) Gray{N0f8}(0.043)\n",
" Gray{N0f8}(0.157) Gray{N0f8}(0.161) Gray{N0f8}(0.047)\n",
" Gray{N0f8}(0.188) Gray{N0f8}(0.184) Gray{N0f8}(0.063)\n",
" Gray{N0f8}(0.18) Gray{N0f8}(0.18) Gray{N0f8}(0.047)\n",
" Gray{N0f8}(0.165) Gray{N0f8}(0.165) … Gray{N0f8}(0.035)\n",
" Gray{N0f8}(0.173) Gray{N0f8}(0.173) Gray{N0f8}(0.043)\n",
" Gray{N0f8}(0.188) Gray{N0f8}(0.188) Gray{N0f8}(0.055)\n",
" Gray{N0f8}(0.157) Gray{N0f8}(0.165) Gray{N0f8}(0.059)\n",
" Gray{N0f8}(0.133) Gray{N0f8}(0.137) Gray{N0f8}(0.055)\n",
" Gray{N0f8}(0.149) Gray{N0f8}(0.145) … Gray{N0f8}(0.059)"
]
},
"execution_count": 11,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"using Images\n",
"image = load(\"man.jpg\")"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAQQAAAEECAAAAADrYxXWAAAABGdBTUEAALGPC/xhBQAAACBjSFJNAAB6JgAAgIQAAPoAAACA6AAAdTAAAOpgAAA6mAAAF3CculE8AAAAAmJLR0QA/4ePzL8AAFTzSURBVHjatb3pmiTHjSUKwMzcI7MWbi211DP3zp/7/k81X/dMLxLJqsqMCDfg3B8AzMwjixIlsVOkqpgZGeEOx3JwsBjXVltlIiIQE4HefDH5T4nz7/mH/+Urv/LV7z28JxMxsfg7gDg+ff1nvl1+Lvs/LOx/MmDEUkSY+Zc+DOMzQQQCx4eZmeJ6vV5re3p6ejrdC4jjPxgPQuD4z/XfX/hU5OXHe2C+2O+YiZiEhZkAYub1NgCDEeD3DcCvndkF57/GLASApbZWWxGhfJjzth+fCMabErQfXf/8408/1vr04cPHoQlY7yMVg+d/8/IU471CR9Z7PP97VgDkIyUmSiEQsbCMx0tEZgYDEbMUgmlXSjExx0uFhWHGUre97VurpcjX9WHVhSEE67fbvf/r/5VrbU8fv/uBCONZYb1mDE0cmhB3HJczXsmrWmDR56EJ848wBSISZqFVCEwsrhwGELFIqwK935Vrq0LWDUCYhTDDwFL3bdv21lotUlJRmE/PJT8YLgQmwvH6ej0u5fip1qcP3//erQ5nITAxThqxPHEG84MQVhMBxtuMDycmMGGqVDx21wSaQmARJgBGxMJ1a5X6cbd6uVwqo5uZkfnPCaZGpW572/fWtlKLsIhIkeKOQlJ/HoRATLi9fH694/7z/6nD5s96tMoDNAWUpg0hIlmsIATCaT1pSP7nEAnmn1OAxOkemJlFhCku1oS7MJmBuD69v1SomppZN5jBtBv04GupW9taK1X8K4QgUkIY05nOixGRUooI19Wax0Oj8dL0ijyU2V0djIlsSitUYdgMEPJASGh6iPlWq7DSpIbFExGLgJiNe1cyM9XCIkRggADA+nG73o8OBYHMRPzN/e5FWFIYIZkiEgpHTARmYWGqj0+fTvfyoCCLmTGFIzlriPs8/wZPVxBeF6f/Xt86VQNhWmFtIqWUImZgOu6tuDCZiViImUxv2+1+mJTW/KcgAqmm52YXQdpIWQRB99eXmx5qqAQA55vjvM5VleN+MqDz0ImUBVbxkAshjIn5LD7XFU5lGbeOiCjpisJrCINZ7v3+ctlqKSIsQizMQmZmZmBpT0+XJoC5mQAAQBbX5UYW1pFCQL/e7vbl9aY1o1+oK+PBOua3RvSFGwv4ZNfpP0JWYJwBxnCcCbuwCOSNO/J3cL33ZyrXL621WoqUUkspRUoV2NG7dkJRA4uwCBOZmampKdQFYhafJpy2QUz9dnR8ern3OsId5jUOdzgufw0Z+TLGQHkDSM6QGk8aiRfY4Q54vu30JOM5pA9ebDJcMqzfmFOpQwiF2bSbEko/qoBFCguxCBoTk5mqqiksVEYX4As9utKnl1uv6bEH6uFF+YdU+IxjH/T/DNBSW9IDMgZ04hE6w/Oyu0wXxxqKR2AKDxv6ZcYHpX6XUsThptTjuL1utdRSSxFJWwIJF5chUiBmMP+z964U5rA+SVpSCPBJS8d3zh7NdYKG8NICkDq/xAUaFsKp6A9IAxNHLW/n75RuKB/IeBmzSGlta62UUmsJ6y8l3tGjgbBUwGCAWdd+dCFiJZhVCqCKxZzz5nDKHiZk9uebeQgDUwCcqk9Eixh4jTJInQBogbl8+phpDsQRlXnRO/B8SGAjvjKL1JpCKKVUR0sOBfz+wibFiunW+3Z0vWxV6iJQl7Bx2AdCpznD4Qqhl5C/fHd1cjMa4NG1rsqFkV3hQQAjik4bZU5ldCkMDSUwSPUg8rAh4kIQFnEvmkmX/7LBDCAREmauS86EyNaMiBARGcyLpg9svHpEXuxmWn9eXBrH8hTTO55Eks/37HpO0IUHOg/vtcrbJc0E6LxYZhJxN5FQQYqIEJH7CFMDUClNAIAB/rep4eD0WcwISTLJOfBPKTA4NH2Ia4KNYeE8g9G0o6HlU1yrVwr8NHAvZw6bahrXgTBkjBQ33QynEFwmRFAYuhrVSHgQPgNwNDzNmJdnEP+ycCQ85Ph1aMaDrq8PleLqTrqfqr24nq/8sosn5Mocf+Nhp2AM1D1NKcOvP9lhbQ6VRZzP4de7WnXNMNgQAh4eHi+CSHjPDsvFEz4nQDKMrXc4FWLYw1cytek15jfHE8CUwlk8GJrgd87hMRjhalgGVOf89sqLMDHxl+tdqzMWDiVCE86PJVKvlPOQgohGasIEmDkq4yEBxno3i7WcHzUebj7RNHg8+0U5eKT74BDc+M6JQnArASeLcrorjmQHRHR0Q4Xpce+qOggtLASK6wB7tjLoNQ4Ds3Q2LMVMTcN4bT43XtzieLSLjmXWNC9wWiHxqglTJyPBwKJGKaMRaOGhJN98wBUMq0wqE6AK7ffXo5vZtJ1FZVwXJJicZPn8m56lFo3EtZiaYU0owY+WzYv/T5muIliN5CuOgR1pRsKXGWyoHzK6uBLQJFN4uiM6PZOw4Gr99toO1cUbjKtI+BK0ZjIT4aa1M3PxMKyem5mp4fHK1/9er+KUMCyYaVrOw/uMVGOElYEZlmxtJiuJAdmVk1aTGNgGAFU9bi+iajPVWVTTRZm3PtiqhcI7WKTUUmsthZmLWHoHGgnCSRQJIJY/TncwbWH6wTNMH0zdQM3+f5nYnkAbwn0sck/oHS8DUO241WI6X7RIysEHhRIMKQT6cv7HlJilttaqFBEuYjpCxZoQhW4+uMvBS0wHuRrIQ9gd0fCt4SQceIiznqGMp5uk2XwUAFE1Pe4FWMQ98SFx5JjDGkISEhymM9wA93utrdZaSmEh1ge+PS9glUAmp+stnG7rbW46nv/bF6dKZWo+JRk8OoYn5GH2SXdVaD/q+gRmVFgulhI1piSCtArTgGkXqa1utZUizKY2INCjafPA/o6AVnUHLdwkvnqnNPw+zezyJLbzB6Up0IIpRkbmQabCLBk5OlF/vLhgnpaX5N+gaFhYCGRK99LuW2utijCxLTLF2Ts4CI9EFLCJf2CA8818Yq6WO+MMEY/2sDiWQO80DfGNPFNZPETCTE+gLgoBkycdAJRmCjMUwsOCy1uP0tq27bVyYdbFCIISmIF7+GFo1yhTMAOmYCmtlpN1L88yI+KM4Et6cHr9AFHDME5ZF00IVQkwGxjP482aOqSDykgxvIVTvkEHR12VtN/aft+3JiREFveZH4j0S5GtEAG93zsIRiAWJhi4VCIWGmTmG4WYdzXMfv3pKfQQY3mWmRRzYg0mIlSkB6V4UCMNXgCu/zbGr62QKZjsgCbK2u/3y7ZVzmIDBkMW8Yop8hWQ9eN+6+YJLBcRIuLSQFyEJuYdDpiShXtEnm/thtKPPjz7lMYUI1cCoOf3Cd0NoSTeTEIjma2EEiwiToR7Uqn9fj8ufW8Z4kyNuBADagRIYSE3AqD32+12OCdMLC7M0gzWilckDSTiFReKSv6w9rO2PqhCGOD46SqL4bE9nlQCbNBoWOxueI/MyvPzcQ5RzCLmZIV7E+XeezerTMQw671DShNCVxhJEbHeu4Fg/X693Q718iuXUmthLgrr1UlUU5C4UgmXUovQyKtC0cfFxFOfFzciTEY8TojN8WD9JzWMc/pWLAg81YADfE58Hgbk8Qowdr6CESklYFslAqwf925SVQhdzUhE2PrRDWx63G+3wOwEllJbLQJA7y5TMzXP38Bc6wbm7EIYarn6gPW/z6SMfwuZfpxUwy/U01Y+4TwMH+auH+musxCQusIMgFhFSmEWBpH1G0GrEMGO+/0wrlXIVA1gEdLeOwjo/bgdqhYXo6rWq0jp96ijQNWCCWap227YqNDIjB7MPymNU7Y+stk1Pzy729SEMP/hyN/4mHS14wUjKYuYKixaRIoQMZTJeilM0H6/H8alsFdAiJmh2g1Epr0fPb7JBDOz0oOrYWFmM7VQYy5172a2tSIL6h4wYEkZT1nImpesvMyS4dbhWLDY0KnkkRaAuF0eBKrjck6yT0WLoAiB+DA9RIRg/X6oETMB6oECnmw619ldO4hAMAZg5kQoszADGv0ERYoFhF/LNG8A6QK6T9npWQppy36vdXBwJCPrfsNyToUbkD84rGE5xETGVkL1SD3RBFS7GcLRExGRM3nkxVM1gJiEsjLM5pwve0uSAeyx03ryR5NiZR4uCguH8kYhRmMEDSZvshrsFaj0chhqPiLiahFsPJ3D1BywBdfBCTTNhm2acy1mMHfLgPnNhSDgdggAJMZkwoktvGGFCyoJAOuhOFTe5lcP/UJ8UvjpIfIBrj+o43Z48SOu+18pmky3Ml870IsnAkoMWAcs/OrgLy1B9+T3h4eCR2oTAytnCZIM5ipKCgJxV4MZ9lom7xUBfHaVDO+4JgPp9h6kwsMnZLLPU789xUtsFvzV6Cyg6RGIEwTScJKGfpjBe++chAUARL8MEYFSBpgdPC6yRZMRGsGs8VlGnm60UjgK7DNcr7a/GPHXvpIgmJqQYlgVId8vgfl0kZxM7ZBYfGfpPFPTDoDNyxlm3i3BgfiG80I6F2aCgVhslmlADqG81CBEROwa1fsWbQqVhMm+0lTxYCu//EMQefEl/f0bCxol5xQNn6sLFF4w0xIWKd555jeOLGh4tiwB9jLRiA60+PzsmIlatdeDiFhgQRIygwC9X7dWa6lt27Za5XS3g5MfFMBDHBlas1DZlUZj5YTiC4v9Rj6PHEYAdQ6/K54JAqaw/H9YlvVATCRexgORRetC3jTZyAhAiHZOIg4bYWEo7LjVKrXWdnm6XCgZoahCrXo7bx7n2zrn4lTXlHs8/yTfUqqLT5nJ/OicGDjWI5grgnpnSPYPMYjYBOwxxs3PPIgQSISzDjqrZxbWAGPzYAhjMWPy4nu73Hvg6HMTwZI7TCWOJzXhQdzEqErPestE5OcwO3IUJD5/Q1A478pkdvQe7YaOELIAD4+uwgxvLIqfEsOYEkYwsyCE4A8DJu44XFfc7ljK0RVMm5zvnQPD8NDxIRGkMk99gVelgXNtPSDlOfPMxgCsUCMhhWRLMsMA1aPr0U0jQ6blzUZSx2RqZmAPCiYjToRwo2eV3UgsxAOIEDHEwKZe5cBWfsEDYjquAFVL5pByQeCEVUdoaVmZefcq50XfZhqV16ikqkfX7vXN8fO8Cud2CCAzTQMgIovcP8AfcdRF/SMsQrgEA8DCREamcCdbkhyYFk35UUMKy72ECWdWVZG5NKWqDNUZzgzj15cK4ypRBjw+mlnvvaupmhlsMUtOJBJAXXWULBPKBWj3QLxc1XiIACxaXokA6HGtpQihROhdu2aQij8yzEGOjFf5g147WgdGwEQcGGIZFU0etN2Il5HkQKG9d9UOdcps+cQFhgCWgHl85ICz0SCyPJjxkpTJHI3Q+621ItnrOBoppvNO66ME/POa8laRWeSCu5L7fWNfmVdxPtfFuCICaO9du6q5SQ8wteJ4/xgzGHE6gnwFVnEjoW/G4lH9iesEmR73W62oUwiDSg0S6GS+J41fccJ8FaYQTln446McLx9gJLs7dI0K9PhaSqwZYHAC08W2lgRl7eWbrceOqkAgCBNM+0EQZIAKZZjZZtCR68c8ZF11+fsQwykfYzpfZDiIB3Xxn5qZ64HNpOAhTIHICDwSa6yfg+XSou3EJyAGOgGMEX6SmIOyYMjogMwuEp4SWfscZo1rqnKdH3tGBitRk7/+KI/lZ+y9yIgOhfNbPSDy1IHTvNHyNa94ZFujaYuMCDnkk1YPYxsc80kIqxQo4VKGjYF1g1RZGzlOnMNUhpMhrFnpaN+JpDEwYrINy72l9hsHPDhLa6T0D2YBjtjowZQybWXP0EYdZ34Sz4GpVKa8UCxSSWnXtz7jTV76tWc11SkmdkqRyIJcHZagnyiCZoLOhMePWBDdqnyJ+h1dExHBmJkkL0K1F2E59eMQL+5k/FcGFSaXKhETqSH5hDdfwzueTHpi55RjAlYRkYxtiSim2+NI07CwYV/52K/LHAwTC/6OiZ19dL8AU+29FKfyloR/fPpy/Y9CYDBTV0Pl9AdLukUccXLJyc5oYd5lxiXvbVyHrvIasvo+SXF8/X7X702akxFWzJnqMkOoEEhA3KMozImKHlyNz0wkcspwkeNQzlXV9XLnlZyb02lK89SonRJjYSnCQA63Rf8TEl7MBCfk4Ep7gp9v7G0S3wQWv7t4KmJsQgwy5ej6Z6HoSjkRrok9Zj9jvn3giaO7ENJl5igqTfS2UI88pTDx1viO+0WRgWRm1jSpGE58P3po3hoBJ8pbh4ryr0I+hhdZHxNBmVlKLYV53iu9TZZO4W/55qGKmq8/Oaj54Eb0pGydfAQ1BJCxCslAwqfOKJooc4rjl77CiYcrnLMf8VvmfCrD3SMYMGU5jqNa9G6Nu1mQH2WP/hmLE4jILBu8Bz0yypKebixdh2nZWON4YhzM4goeHgCdXrlewFfS35QBpyBXhovWv3pK6fi7y3HUNvuVMO5k2iWWXO3RWDyLtMQS5+62mVHELyzJI5aXWI71qnob46qWQ6V4bZn6C2owutdB/AinHKGHUzYxEiJTkq4gLsLZqzYUIVMkJLZPH7lmPUmvpc5hhPJxWW+h0tQdf7UQAQaGLo2xDzeXOr1mXQ+qMLHNSJ35Qe/iESdxABhHrGdpTcjH3pKz8DtankoKYYlSRORE65TBCgxm5+EjTjija3ZQCx9SnJNeQ5TLr8/Roscem0nyDkg8QPtoughmkgc0Mg/PZiTbvgt6n+HQ7SpGN5BCWGxlaIhrQigYP1jMFP8CkmYb+wJxPTJi4Qge1WaY6/mnkx7Nj1h/j6YvTjmsA5xkgX8Vsl3e7aQcPeqnC+Ahdpyit9uHcI4JT5jPK5/Dyweu+fdkn+LxeDaHhVB7VPJhPA8/nAnIhF6PCdiinuNiM+SxByWu+9OOm/Iph1gLyDwS/6XHHTARpvrYCEvntziHgfn3WaJzKOMo0xBLMTDF9SCFr38xTTczrhLRRYjFIvwF0Tvjxu2dMbwAlNWFT+xOq9bFNJ5zlHUlEtbM6k2kA52cVIZrv9wMSIaApOBhlwvgenvrE4sMqLf6rLUTIYWPrAQmhWh2HP3oNct6I8I99Hqv9xHAFhBeepZopSGxPvy1CWQx7ox3+UN/T6+W5ajWCPrk42UPdO0Ah2+sA+O58dD7IfsHtSRov9+utyLgaJAdBrb0tZ6lOdrRKH1CjljNR/QW8SxavQD/iMqmROwdBxSp/JypPDUO5Jt91S6GSx8OgiZviuV9lpUETGDrt9eXS2nEIiOjWW8B498lAc4PIQ7YPLHAaQuOv1uq6QmVDxN0DEtEg0ZgZk99B+cZdrpQ4skq8FlZecLV+C6zwAgnGaxfICLocX152kWkyPiwR6GfbSIbRYnINSHnRlMPlwCPyPuyf3GFLWvPLYI7dZKNmOQ83xppw4lFGAh9quEAzBl/x9KL8y1l/3W+k/b7/d5B3vOX9OTygsd0LSfmYvLlrVzfwIQUxoxn42UYpZ8Jy0YbKefNLHe2/vPWHIbmjAslJiYxivI/olNpPC9/BKb9OI4OvR9HTxm8nbyLv73xEGkOEXDGew+1WK5/cbbTWYaWAZGHBxjh4TEHwTff0kH2m158fhOtMZsKySPmaG9eboOJrN+v11L0frurd4KcLjej4lDisxTqwNjjpjBfdnqf6dKx3NNIL4Y0FrSMeW+jFTQaEWgwPSMGMXPOJ3ufjgPesKfhUeJ3liovYNrvryTa70plHSd67D5ZnHr8LgXlHnowWLYEg/6TefVneJfLIQayn8zvcPI0h2Xy96IBxUL/h9/g6NkLBXdSeWrAsMyZ7Hhni6ugHrdXsHaFMM8kJW2PwUu4eCSGfEyYQgTTVod68FSkkVdhWQIT0C1xk9feshAeA8WzEBL1OstnwDGW5QPY3rjAcUUxGzTBYljnCr28DuWqcAOpQdhkCQQDzxFWe6BVRFRpeZyPXRczjRhJ0xlQrhlo9urR+hY0ZuhmnxJMzJRt+MxRY/SiirgkbYnM3qZy5u4TqIBAMLN+EBl4YP0zcD7bBRbfnj5hatrqbHj5nXlb5zJlSggWrcnRk5I9aDwEUXIeAgYzZWXL+/G+kzFtKVBjIh4k1QrWz+00cZeuCYfHZQJsFukn9BzfWGr0UwjnW8abb6wue3RbLIlkQuk1h4yHGA/Zd+KkEGLoyusfiMEvEYqZMmYSpUCfC8yZfmn0/6SBG5updgGJtzjxVOyRUDFONzQCAYjqYh/LhDUW8xixgJfJqPmO2emJYQwJdjhrErnjzFXVeKESsEy85P4jFwyip5EC2y0RFbMZMxM3MzWnWJx55/Pt0kLbnufIiRawFL3+0/Ro4UzysT5ohnN6EShGly6vcnNj93Vf8RA5uxLSW3nxhHwAnxgmJhJTGBl3VqtdC+JxR/EEJHbIvHl+tED2ky1jOsb5szcEYTzuRTve/sgR8RTAEBJnaaSUIiXCKIGZBZJKRFKKpMZIxrxoJ1swf9r06iOGzSaTIed7XDHNzBvWdwii9asINjHnanuruEbgnIhuBK0E/aHgpdZSfTYI00xk9Enm8kGOBAgkwizmadjSozrYvZGITDvPkW5aBRTAOtT7a8wORvdainkdfks1oxWBD+GOGH1+9MmtMOVoeREppVbfFenteAQWYtPYcMCUS0pJSpHKpGAWIWIzX3Y4wiFNYMJTwwM55bDgvP98XOdCDj14Bc8dhoEticvU/vmbD5ry9gd5Fez7c4U5ZJBrjkDGZiwGUzGinCmNNkgfOoagABlePAuxVdWWJzr/uvT7jVgwpnfiv09MzHrxdbkDPpkap8yme8RDpj7g3JgsDtURqbWKsA+O1jB5Bpij47uoVy6JooFZmLlIKQyCSCElwMQrj7MhcZo1rRTT5OtTOKuJJ7BbkpqpQ5KkysPkzHy8XysdOnfNY05mNFVPfpS4tK1VGVPEMVJIgES/t5bsBg+fFtugvOWjeAQvAfRgnPt71obNNW5iAOrJ25xMOKuIK+ElzLy1yg98wuRJvhYKziJacRGlywMRsUht+95KOHwRil5jGtucisFbJXwcgFhKgEpjAWoHRCxmA0zMeKRmj/GOaE5yR0a1pEk5HDV0VnxNTg79H69faqXhnjwi08I2TtCIh+A3DCShL2aazySlbvtTK8PQacxLI+aibfBvZGaAe5Fkk0SJyUQ4BuUcjNraFYaHa+QxbpA4ikamnnsexuqxMecugvuXrZ4QIpaR7vD9Tua8Cdmrnc3x9bTMUtvWthK5YXqTsDzjGALM5hMXggjHdAuLMTPMChDte2bmu9ImbDrpZCyz8HcItyDMLDEVUHIxaSlFirCUWtxe6fZpC5+wSGI8XFeJSdc98hEPjmfZT8BSaq21kMR6ovn8YCQmzpmkEMTANFb0gkHmysOuIcSI5qRDLTbgLKSOW2DJTb1wI2QpRWqtrZbSauwzlbGh1ZMZYvnp/aWOfYwz0NGiB+lS8m8ZOJKJm2XKdCVgFim11laFhbyNYIA9MoourJgKYx8pyRSSiHxgUISJxJfVMzMZrPc7H+rRJGFfPHmWWtu2eSz2W2+1lVprabXUVmrJddYJq+Nd0apwPRnYDI4LTF0c7KL0j5z2+l3xlUPie3bHC+CDG5yb6kfLWahzPFwhAlSYpcOF4EGkH248/JDcMBGX7fnjx3dba6XUVrfWWm21FBdC8T0vuR/Fh1IJpuhm93vXmHdYbu505+ekY7IJwx0sbF4+7IwJvhh6CTNIoktcg4MqDA3OCqbnVkWYpcKboUSIoF0wk/Vk6/zdy/7xh9/98HHft1pqba222motUkSqeGriI1Vqar2rqpr1jq7258+vPUcCIw86cSkjXzihw5OCLGE5QDODeO5SXx8XUUxBCaYmJGoJbs2FwCATkRqVzVieKFBVtZyx4WG5TGX/8MO//OHbp8tWSq2lliq1xPYfh5Oq2rv/39G7abd+HKpm//HnL/e63N4Zay0qMZY98vnVi1bA0bv4jHNogswlD0MTMuDH0FBKNmA1Mr1GkaIRVx1sCuk4uyBIudQklu3529/98YenvYX3j3B8mJqZqnXtvff4oytUYV1NDT+9xLbekzcYXne1ihM3OX9jWbiQk1eMGNPy1whNKQAQk1xhnKxjYAOZ2xOZCGQlJkoDNBtZtGzObG94orI9f/zuh+8vrTAR4fDNoxoLSIcQ1O0gJrdVDcDr9dAYBPO3ng2soyaYdzu+N7NSzpxxbfbiSI1lVBt4tYflbXNfV0J4mXmws9JGo82ICMhF9LNOEU6KmUvbL097sTugZl3jftW6qlkIxLGGjQ50h16xy/1BtmPXJyUhOW9wjQ8ZO3nRFGYOc5BZb+BhFFmncp/KsRrEB0YHITJCtNPIXoFA1i/GxVJmx8Gl6O0L2EzNtPfjCOfnk+ougWzBp8jozQiIjtbBBK5J01uEPBmtxS5O/RSOssIcYnmh0KjTE40GhnBqC153UPuYsICWQetM1dexNQplK4Lj5acX+H2HJpgvt49SsYVaGUWViFwqR1fMoY/ZET8e8ZJ74vz401OM/mWaZEuUUYJHJ5mOkSIfH2pDznhIoPu1dcjVQcz9wSkzXlvXPekshez+pUDVHFr2MZ5sWSLF6CJJIRgArA3eCwLIQLeI2jepjf06i6Awc4+oAKRys/NDMotoUaG3sQ19+IQc8BoxM4VgiagGOE9pZkXXOSkm61fvK4Ye/eipAmQYKpBVQhhR7LEgaLT1jtGQER1X+sGXgkzM7IgnE6pZ8SEiFhMDTX/Ag4dcLEmmU6Ch5+df4MxbyHOuRQtmvXfmuMxMdtzYx2pdE9IOhhCGLmTNEBbNW9nMOZTrIUNdQuVAaqnPq2d8cCBLEXbVsTAan5AeE3s5MHFyEiMRYTEZzjLbHwa5ErljLYx+j44hH800f6n3W9OYvlvrrUPlRs8S8gGc7je/M6q6a60xJ7OHInj/R4C8RApYG6+YYtA5VWw81IcUmdbrmGXR9FUr++tkNbT7sKX1fqgqZdaaW/CWzvPkRjjM6cwnLJ9/JgnGiRXn584ZNuIMDehx7xEcEh+lF6EEu0holSs4cg7ATv0co60olDi7ck/5fCBMJujhRziY6tG7IZMRO4s2ccipya/Omx73nPd6hghzeiQfY/wjwdjUWuz+8uUVklslM/MfWn92Pacrm/xC/oG8zaEkPLMnHprqqZj122GBCFUBWs8G89n8GeRmIuxSHd1rq+NbGKpp0bNXZhJ0nHspiYilFr3ScVevJPnCTSRpN7HTBOFDs5LVO02tMY1o5t4ho+4pXDihgv6qrLm6wogy3M4AhKnf01WaHUe34RM4MO4wlolOOectJVIjKUVq5khSRNjUILX2F7vfKPRgjoBE6ZkDKfMaVsd1DdsfhwYw8QjnNIYv33Q2gcj6Ve40Dv4hEvLZvMeKW1yRUwpmZvblergQeOhg+uFQ3VzNuxyOUIpkcXGcMSWk/VCIHHTdqhURjlNFfJ0WQMJ5zsQEv0tL5oiio5zGSRmORQpZ0aF0wZnhEOy4USFQLBDP1hhe7z6rtgYADiehZi+eQC2W7yzIODRpnDAmwk5Vxllj8iAE9Pvt3s0gqbIRqjrgM31xPpmM5p3V5wyXF8guG1ck1nMGJHfCSiSB0oRa1u+QeBkxc5TAiSJg5LkFY3H//O71dlgl5pKnynkXQRnnh5VxNIYLYUhB5olrgVSqyP0w611jFwz0OI57j/1jMs4tkuCVJ/INjQ9EZ3EwXrjaIuNOnb2tPbo3YvMSkwiTdQRI9TlVJqgyWSw7skgkNA9xQMoFx9Gt1v1yeZJSqpTqfMw8XW49T22csjY+bLhpIiXiut1vurVaqTAp9Dhux/3QaHd2EdQQZJxS5OEztNPz/FypIyKlttZa5ZGVlVJbO6plC3VGHbJ+R2huVFy6o8d48u4xkR+SWQQIMSa8ffj4zbelliqlxpMv8wSLJQykFj+y/uyRYev9/qKvrVQqjI7juN0XIUiRIjW2vtdaJAGStzB5cIuHREQiRUrdtt20lJQCl1pb6352AEdxgJlgna0Ezy5SfBePqVlXf+BqeTTYghu8eE6A1cs3//z7f3ZvXwr7+6QjzqMMJshYm0aIxhZfEgJwl/6ytcZC1vV2u98P3zyaTGEppZZWa2ta44wmwNRUe9eufZzpxiRSpdTj6H2rtTpDw64cvau7cHcAIsJSSmEcxH4eXiniCy/99A6k5mfGG1m7V4buDK2Xb//w//6vIqWwSCHPauOuJyKeE25v0vms/jNR6du+t8bM1o/b7X4czmB50PGqsz/L6lu+EWspez/8gnOZt2ti7UfX3jZzxRmecZUBs5S2X7aCfu8K921VmKwf2rvaCMCBc7wDhIm9U0xeBb1evv39//P/JSc4Uy0AHp1CfA9jdQsYm1yM+SUxMVnvakbEMoQJkHHp1WtpKMJxDNHR+3F01/KgoNjYWMzUtDfV1rKfp5RSypLchQyedrG7dVUTZq6tFEJXoyh8cIZ8V50QAkthkSZ6r+3pw3e/c0UPIiaYGErfkRNk8U8Gs0FuzZwMpsZSA5ByLDdfElgy6u7S4JU2Ne1HP+6Hqo5SarYuKxNMN1W1ViWdYxERiwEpKaVtl8u7501wL0RXdBCr7VJIHAfQIGyi7B/H1KTDa9SvlUvdNqIoTgw4uQwwzLsYA3WjNLtk+qzQfnQIM7FKMWK1IcU5J6U5+S5EaQ1qhgSJnGVBJphSRrNS3F/Oagaz1Nr2y/Pzu6eNcdTt9XpXUyrbFrSzmRFl44Mj1li6PoHI8VorM3Nx2DWJB++lHGciDUWYk4VZK1q4NtPjuHcqJAQqxQ9EtqXfzjMemIp61Z0SwhExl+T4sirjzw4qBxFhgzBxVFQD2ZW67U/Pz09PmxjkUp+OQ1UhtZRC5KvwIjIVyePSwvGnaVy/bGX0J+T05EK/JEc3ncHyXGkwSllw7tqPQ7lQIVDG7OFm0pER8wKFU51yfdtcSDG6w4VMj+TnuZRSYpOI1Hq5XJ6fnjbp1iGyP5OpryITPzQyQpNMsONCkAmqv+xNKomU4p5tZFg01SFv2KYMVnnQ8JtkbNo7hLiQIYZF1eaWGSbJqpH35wgx+0lu3pEQDyjyrGyBEz/xwEw4zyFjIiapbX9+vjztjbseneqlPm0C81XYQwhEM+kdSYCT4cRE9OWnrVRmKRXTvk/hwKIEMhRk6sRgGZNmULLejXzJizITYikzxnyYiYCd6JRSq7B1mBJ7ncrbeaKRze9jHO5npgXEDMkCLUvZLs/Pl62i9+Pe0eTSPuySOQHRHArLJDxPWI5MTJiJfv7TpVaWUutw+mnla9Y9HIJhqEFipkncpSNxjt3VoHefUTSLHgQiYiEQiZStFbLC0KKDMkVmkeKNFpZ1HInWLJFS26asIGnbvu+VzXo/7t20dSqXyuR5UuyLPalW+gcPtA42W2GqLLW21d5n2rakdzNkULLWiKtO3YFMBkKGb2MmI8uBFsqm3lK3fausB1tnnnqVlsMsUsyKO6ykqvyI0N3k6AbZLpe9CVR77wZ4Nw979iDuHbMWyPloyN0PYJHDiY8JS6kty1OekE+CORfT+JyCLJogmfkO6zFPixzClFJqL5UhzHY4Ek7WzzOh/Wmv3Av6UUQs69qxizJLKtUAX9ucLaulbpDWu4Lb5dIYpqphQ8mhMReW2gJuR7a3llkXjoTfP19qZZEQgt9b+IBgWBYhDJYxQykxTbMg+FFQ6fNKbQYqxMLqW0U9F4GnELVu+6VyF+v3IiwRhITY1KK7MzO+pdwkxMRS9967cWmVYaogjjggdbt4z6RIq7M/RXLONsltT5CEmfn50qSylFKHoWe1IsdqKLuCmSDx1MUmdp6+BNi2rRUykBCXugEkSixUGNyNY3qVS6tba1vb9kqF+r36pTKIpDA0im4MGPfwEZIVUyrMpXXtXYlrIZiBxQc/iKXuTzWTjJbHp2Y9cEWu+XbMrRapzCJlYkLYuMXpFH0ea/hCTh5rOAUCuLR936rzukyNmKWoEZMRuIAkjKXW1rZt37etgrXVVmv1wClSBdEL7/u5yYg7C0twZsQkbFK0VAWzMEAifrQGlbo9vfuwMXPxnp2MLYPTY5ATkcmdsvA37/ZayaesEwfaAg4GgsxQQcN75cqKAFjEQGn7ZW/eg8pSSVikqxFRqeSTOSCQJ5KttVYLUEptrVtnZQNLE1IGOoiLnybFBNNia2leiIgCCQHMpbDBjMr2/OH7bzfJxvoyzuoaha3sE5glA/7Xj0+1uuMY0Cffe5IvJ7Xw6IHougncE9WUbb88XdotvFT2qAAshWJ3GoRLa621rdVaBFJKre2IcyapNqEizAVcavXTDrw7cbZ9pHcYBwIzsxtau7z/9p++v4hwESmns/yynyOT1GzxIKJWhL0g60MBUWdchJApFNFYwh+psZcUiaKwSACVfd/3rSC8kdFAjEHGgYRLba22rS3ZcdsgrqtRyKlGUquwmYFYSmsteYRccsAkIHQDF3K7MGqXd998++0+uN5ZHqKsWkVGNToh4msum1qbExZawcIdjqNV3XLCUeYcIkCE/ely2aolB4lSteg6JcdSSov20iIu+1Kb8SHM3YgBFpZiLKVIrLiXUlrbtirCiBoUg4nFoN24NCk+VHN5fv/hw/stW0NoCSs0xdANPDvrViGsRBEtm2gMkDCDAOIeyz16gg3ZPAJC3fb9stce/EUQzIFiQexHHLUafGuWNqVUnx3sMUQgXJwz95JNKaV6uKNTlGMwVEmNmEUqlXcfPn5491SW+evHWwOhH4eS1Jo05yqE2a46C5Gh9FkWHUKAeM/PDOEBItt+ebpsmkPXzJ68ud90nrTWOtjbbFqR4j22pDBEwixxNli0JRYpywBgJNJsx2HGkFLBhS8fv/n48akWolUB3tQ9rd8VUltrW82C6Wlb7/n3OY5mPxFL6TUczBtIkjsybtt+ueyHpemKlFJBfiyBiPgD9Y02gR5BJNJCO1nNz3oh5wG4lLa1Fv3BpxIms7Cg9x4daST1/TfffvNuk9MtLFX4/H07Xm+dZNu2rdVSq+Er5nCSxdwOMUhGWxwGSWYTBDCwPz09P92PcfQhS6kUDZ6lSK1FGHHgEScYJhY2kVLKcXQvvXjOU1pp+7ZV8akwDFjsFS0h6kc383kJefr4/XffPJWvW0G6O7AI23H3tRGHiMjPp2bOXxTEJNRAft8jrR7sC0AM2p6e3z/fzFesMIilgGPGJceADIM9GumoD4yV434/bFCNUtu+7w55ArSnjrl8Gdu96uE5Q3337ffff7zIL9+HPyjTftxuKEZShPzQNvtLQojfx2gfOlnFhFJB4IAv795/+PB6QC2FX6iEr5MiwoDSXMYeq8yj3buU1uq9d82647Zt+1brqGNnIkHe48pehQiWff/4/Q/fvW/8i3cCAhTH65cvL9c7CazfW6tVQIS/IoRUp6kVAylM+tQxJhtvz+8/vP907erEAvJUvxj7ZEBhwlxL4UCmcSIQmFCt1+Y1K+ZS6rbtWysypJ8JSzxYGEltnUlEuL779vvvPl4K/6UbgfXr5y8vdxNi7aT3Wmq9Hka/QgjnJh1eHESqhMCDKT+9//Dx04tPZjARkSA7Gr1TncI0nOK3JJmDXC5S6tHNydHatlZTBsgz1SbdYwKSujGkSLl8/O6H795vf0EGDGbS2+unl7tJq7VUZrPOL6/3/teFgIf3Ik4/MQsy5tssBO8/fHj/8ytBNarEHNMmPHJZ5lJblaRShtMSplLq0T3fYKm1iC8bBOWCimWM3g88LZVQart8/O6H7795kr94G7Dj9ef//NPnu0mpdb9ctq0WqkL8KzThjXEskCLiowCAiW3PHz9//OnlSqYgT9mLDFLDUTSX1lp9VFwHeVKKdVXnn8fuWULAztlfBh+AkVKZpG3vvvvhh+/fVf6L123H9ec//dd/fbp7+G375fL0dHnp+BU+4XE8lJf/D7rBSSkx5su7jx8/fHq5RruSzI4MAKZdjaS0rdblQKhlXpOFUbvBKM/HC2LDHlY9JhVYqLDUp4/fff/9h13WC/zanRC0H9dr72qQ2rat7fv//rc/X3+NY/wLP8BMOBjMsr97/+H9U01gnlMrnniqdqUizbOnyOcTGHlCwEQmUY/JgOoGZ8DDfDZAXJhE6vvvvv/+m+fyV+6DYMfr5z/915fjMGaupbTa2r//x8+3v9kc3gg3WsOcc6mX53fPlypOvWRHg9uwqnYScc4ne8vYT/1hjsOhhNjJNhevYUaGyf8SZZVXhFkuH777/vsP25wb+AUZHNcvn376+edP1+sBsCf05U8/LuM/f78YxicLWdment89bdlEGEDCe6XUzHybgmdfU81AICtAFOCRS2WyLhyBaOkpBHmBjZmpPH/z7XeJk/6CWwCRtP3puUPtfnQQ1dral0+vvyI6/EpBEBGxoO6Xy2XzJw0LlxHTNiCW2loV9uWVGBVgA1irFET7ocPSOBgm+BxXHdDcV8MiLODLx+9/+PZdY+K/fJFle3r/8fVAe3q63u+3fvTjBry+vvyj5nD+GCYu29Pz89MW3H7U4+FpDon4GoFsRIDlqSCmxFZQsDZy0mB+MRq1KFoeo5GDier77//plwDzenEgSHv+2Hl/er4e9/txHL33jl7+9hD5lz6GibE9v//44enznclnXQMUmQJ+CDes+9ADMaz3rmAqZAorVit7Do3RBZqnbYbjnDVgb2tEff7+9//03btfcx9cnq0+f/Pz58+vt+u9H4f23hv131YTiFja5cPHj+8/3VRB5tOwPovlk1xQQoGTRkSmRweJkDBBQT4MEZTgqHfZWHiI7Mwe0bI9ff/HP/7+m7+qCEQkm2xPH15fvnz58vL6+nq73Xo/lO3+kgPjD+D475UCS7l8+Pjxw5+/HB7TvCfHVI3IDEpdSqE8G1O1K6T4UKSZAd4+IKPSlSR3IIrZDeBhpzx/+4f/8Yfv3zf+69fO0mrbn57fX6/3++12u15vt9vdbp/b2Ns8tOwflELZ3334+OH5p9yl4/1ppkYkehjBZzejP793UCEnmmBgWCUSDvL/lDQx8VijEbiKytM3f/if//JPHzb+NVfu/T11e7r3fvR+HMf9dq/68mNFVh9pbYb/uwTg71Avz8/Pe5Pg430KUxXMcijDuMSx0UTWewdXb0r3/Mofx9SEbG3ImwxK3EwAKvvHf/rjH3/4cHnkk75+fWAiqdI05oW79t6//Pif+9SEUW77+5Uh+zr3fSuc8S2m+sBQAAoqfhK5EOw4DnDd+l5IzWItAKKrJ5uo8t2XnW9uF+Xy4Yc//vGfPl7+KlZcLw+xAs9Utat+9+Gp1uPl86dPWy3/sEMYH1Vqa80BUabABhImg2pXGs0z6Mf9MG5H742hCkisoFo0Yd5B9hGbGRlY6uXbP/yPf/ndx/1XyoDG85Zs1lOzb95dau2vP/7Hv7V9a60WOc0W/d1SqLW1QtkUaFCD7xMCtHcbPcPox/1Qbl3VgQWy4Q48zQGnAgJLYTMoob77/p//5//4/a/Wg8erjBbL50sr9f75z//+v7d92+rWttZa/RXB5q+pQqm1MJaGau+ecrkovHOGGf04unLvqlqLMIiMlKO970ETsiO3Pl1awQ3t3bf//Ic//O7bp79RBstkt7c+CXO9f/nx3z+0ffNy+bZvW2tN/gFJgEurrcTWnCjlecsgrCgTlJWZmU27dot+lsicjSiEMEd/CIGUmEna+2/eP8vN6tO3v//dD988/z0P7VHb6/3zn95t275vrW2t7du27W3btlbrWqg6z8XPfu+va0Kttch07lEt8XqCBfpxBGE5JstShOIcSfU6N83McThF5v3j73/4tty0XL754dv3T0X+YXjDVMmOl5+4lK2lKuzb1ra2bVk1fdixBPprAUScP6PVpL0wWqyBSPPUNC/VluKTDUJmbJ52ZQtklj6ipsxc333/x3/5Xbt32d69v2yF/mYZfKU6FQff2HElltq2fWu7W8XWWtu2YDsfVwv9JTEw+0BDfiRyyo2IfE+KUuwDAIlwqdu2bVvlLP0SME6ajTHc/DzZ3n3zwx/+uKlJKW3O4/+dKhCCqEGjM4G0X0nikrZt29pW/Y9wE2NANDP7X7gAqbW2MnJiIiALiSxlbAAg8oJEbftlv9RCBsuSuOVknpMqziAwc9nevX//fGmYJfd/7IuZl/MiyS9M++2FS2lx79u2tXQSYR1//cv3LJVBmUWFOxvL/JN8ebXPcezb3oo4ks6t1jFlM/Il+DvvT09Pl98O1PhXpWw8Hi2tIBw3FilbS0GkcWzb1qbD/KUr4WybsuHXQHFUj5B4bzP5gEsptW6teTkmQAw5bk7SiYdbFCmX56fLRX4LOLMKwfr9fosOxLQ/IyLV40bkp9q3bW+ttbZve3jMOtpGviqEUltugRk9PdFaGpQkqfmQXK01VvnOZkQikrFyMWG0tyE8PV0uf5VE+puFcLx+atxaK1VilSkHp8FE1DXcRHNNCCFsW922VuvKMMYXwFxb27ZaRCxkYAaOLImZCoFi/2qUqr0HyrITLPVsbtDxmizLdrlc9sr/cLr7IIR+/cTa2hYtJCIDE6Q1qh6vzFJcE9w0tpZ/1hIVknFR4NLa1mJfTpYbGINBZK5UNHuVSxEy0/VsXQC5s3ldKsEs7XLZ81T1304Ktb+iX2vbamsxiFhzndGQtxvIK2foaG2ruweOtrVtaz6XNKc7Ss1uUpeCqsXW2uASKmq0abIzbhbnSC3dwlnRiGVzTMyyXy6XTZY62G8kBPTbl7j72vwvrXqbKQ144O5Z++2FpLaWuuBBZMusQ5xddjaNB2cKaNc5fhQNx0QxmcgU5LORzgbSoVfLQ+eyX/b91zEof5MQDP0m3gbrC/xCHt7blJEqW0QYZHojFonwGSHD/6lbq4Xt9vLly+uhFgxInrXNMMrZd8nTOEyCXPeZXu/PSKg4MsisUZTLZd/bb3n/LgQCw0wPFpFaUhW827ZFl5c/keXUGhjdibm0zVOv9BHb1raC64//8Z8/fjnMk6FIIikrUZ4nlJxf4dxbTBZTCjEfFR0rY2cDE3Hd9m37h9Pct0JA9mcZ8Z25lOKL/FptNa3Eh8NTRfNPI72/sJQ6LKO1vW0VL3/6t3/9z8/3ZTtSDpx6OcUHGHJ7qv/NGPBt7KWS5Cgh+5bZzDSo7dv+3yGEtesPoJ6T2bW1OjWi+ZJJOvk/Mgab3r8wV2cjtn3bRL/86f/8+4+fD4t+CpbCxMJQ6wrmUPoxlFI4i9sA+9FJMUmeWXaUm6RcLvu+/eYyiHOgxpfTwz2wXKttGkZrtZU61sRxem4CgM4vLKVsW2t8fP7zf/745W5JtRJXP9ZTezlEldnG3DYzxUCYt39zYRbJTuJs4IhEg2Tf90v7bYGSC2GGQ0p0RiDo4QR1rSmJjBuuEkj233IyhA96ZWa7ffn04+e7SS1MZAbylbkUk8rFVE0pp8PIMl4EPjQzCZg5Ni9HTl63bdvKbwuZQwjjicXcj++LJiLSxThqG3rRag0wEd3fo2ODYP368unnT6+dK8grr+KLVmEmfnxDURWNHcwZ7sYGiHEMfVbiKBuUiNt+2fdflcP9jULgXC6dHP+Kg0EgPdg34NYZNVqrreVEWeTIHvD1fv3y5cuXa5eNqEVLVilFCArk8AozW8ySpeuPla4OrEflKWE3mYBlv+yX/wa/SD754gWNMe7x8AUYHatKuChcLC4Jb0kx6P12fX15vV5V4DZAMeHnGNmgBOZCHGNGHBubHT2QnxU20ln/M9eW8rb/d6AE14RSQd63naM/S2ZNlHU/6zQma9xZbvt+uXBLhAyyftxer9fr7a6Fayl+kJdTzfDzWnN/faTvvt6aY3ENmXUlIsu69GmBurRa629ZQR5CkFK9y5nG+sKlcZJGJptZdr978+m2Xy4F4lkzZ1g5btfX1+vt6AqupUTy54ewx3703GnB0d3ie4yKU3LWWcwDw+yhN+boURD5b/AIRFVqbbVH3mKOYoMopxyLC1mMkUJHlrX6fvQEM6bH/fr65eXleu9qJF21RLyPcJ+JQgiYKc+P8CmAEpRh9uqdT1ak6Pj57xBCqVtr1AOj8NAJGZVxXy5disRC9Bi0raW1oMnZQW8/btfXl9fX6xHjMkmisQ+uiRSNDRAjLRImKiRSaquFvQBJmTOAmPNQJwB6v9/vx/YbZ09EVKXt20asMesmviFLvDTqPK+Ib8QpNZfil/yqTaBGIJiGMbzejh6z0cEG5wyiiFQQkelwv27yzjRWJpAUk3jas7GbRQCCvr68vFyfSirHbyaLWup2uXA5dCyEdWbJEY4vWCpBgZUxdV5GWgztnvOZ9fv15cvL6+0w87psLclQERMLCjGX4ltks/HAg862bVUYRmNQdE7VUE4F6Ounn3/+8p5/a2KJKpe2X6j04+g9uJ1x537rUqQULsJSuJTYZuQezncYgZgLwfr9+vp6vXcDEZe2bbvYsGpfR1uqqsbqH1+M4icAbK2VOJZxriSbjTnOMej1pz//+NOHUs6ngv4GQiCWukvr/XbzsScJEOCpQpyK4bspSkR8EDRW3pn5Xh8j67fXl5fX2xErs+q2b3J4lSWGvQEzqI1F2mFrJbgLZ+WDeY8CJEB5IAbs9ulP33z3sV2E8JvUz6cQiKVyNd2vt6ObSNm21nwrfM0NbMIg337tDwlzC7a3HnAowvV2qJ8+WWrbdhrbz4LApTy0IzYpE08xw1JcY3Qz/FLUtHF8+dO7D+8q9sq/qT1UYim1EKNfeldwKbVt/mhiL74QYN1UdZDQRDnnbzEyLOjH7fX1erikWGqtrRG66ShJJqJaO7LGMLgvsYSrkYOTCNKpfga8/Hh52un+4Xn7TfFCZanbhbiw9a4KX1URPDgAJSKYafdCe06cgrz31puzmWB63K7X273HiZE+lEOljkV4STWT5ING0pgOw3xwwoiIbbQxJiHHTKZA+fOl6esP375/+i1phcql7U8kpRDUQD5uLUJw04WzIKa5H2/kW0SEMchNZv3mMsiaYykiXH1pakwAStAlHP3cMh6xu+ToU0qmdWloji44u/68yfH58+fvv32/y+Sc8Y9ZR5XaLs+c59eNGQTz7ccgYoFCNbYNSXYOuAzydB6gH/fbvafbGKeAlqpGZspCAhk7WMU4DsvkWFqPFHvujQx3IDkuYKZajtefcP355x9//v0PH58aLYASw4v87eKoXLbLs7fPxGCBO/2eezWFYaYKIjPMleTRli0xx6PH7XY/uqUGB4copRQ1UyYUkZI7nmLxU4AlHsWJ2B+Za9iSg2QiIuvdUPrt8/3TT3/+86eX3333fpu1+cej0P5GIUhp2zBO5L5b9e5nYi6MbmZgn+nxYhrFTJKEk1D0437vgzRLVp1FipjBzCsv4k7BcyYxwM8O1K5H76px2JMGYFmWBbhnITM99Pr66efPX75cf/hwqVMXkM0Qf3tdpnrA83TVu838gfhKNBL2NSdGTAYjnn1S0XpgRAQ9jtvtrkbegx1LKKMZnckADYyJ3GEhpRat3qLlp5D0ww8lGKvIwsEyc6yWjSTq+vnz5y9fvvz07bvnbWtl0Cx/rypUwLQHa0E2l5q6ObCIiXuEUTiKY1sw2kspTrxXXR1CDqKwiJlpDgZiLDzy4XkhIlMNMWjukUwZ5B6InBKCKavhdru//Pxf3358/7Tv+7ZtLRdLYRy69DcJwfpxM3jvfdaL/LiMrsalFCbzOjrH0JJRwiTOE3v6cb93zXRgWe9ETpep5c6n2FfBUgKJcLQ/5/E8NtjHuY6NfIOUX58Qofd+/fSn77758O6y75dt37e9tVif9nc4Ruv360vs8hmt9RbnxBiXqoVHgRAxi4eBavwP1X4c6tMsFDYwtpeIiKjH/yQQmLmUnoVJwFcnavegvOpTygDqckzsYK/3zz/+x8dvPjxd9su+7/vupbDa6t8eHqYQRnuVLw1V14SKAvHjUfxASycYYD7j7w22pP2eO2np7M+SR1YyG818vkZknRckmFmYAsUs9DjczxVBzRcFma/7sv76c31+//H90+Vp3/eLS2Hbtn1vW5XTERWrt+BfEsKrVz9GQxbMTLUf3Ui1SE2+TcBBehiMioDIjISgvR89FlAOGWTdIJYdEDSMLhwjz4NQHDCHN4i+Rs4Vz+4Xu/qqhlRZ6x3S79eftu2y75fLvu/b3va27/u279vWfHHRubzNX/0PqjAzZRIa5kCAmluDUYdKdxRvEFAc1AYYQ8TvCujHvVsScqsuROMdisGbUWyqQnafpS7AV/4NbxiEIjyGdsw3ZW/hgN6ZTY9b8TUD+7bt27Zv+77vW/US8Wj2XBvA3nS0UowlzMzQn0q37mEfJmwGYlUmKkWYfUMbE4yUhNiO+5EYwW993QbkBUf/PnK8gvNw3LxGHleYK+Q5FyeSWT806hXj8DH2vNZsA8T0JiJ127eQwqVt29b2ffOqGdFC7r41B4Lpgbk9J4TQ1bSPkrp7As6ikyfQEGEYs5Eex+GxgROreJ4tOTbuWwiVkzcmPEKaEUkmJ5t7FgDth8XhQEuvr7c1EDFVNjHj4yosddv2/bL5Wq89u0iatxjMY0fWPtBq2u93+BaEvGFTc5qJBVZiXw1HPFN/AiLVmAmFAOuZNGBQADSb1pyzDb7OdFQhp8cKYjW35o2zTqNHxI6jo4zmRnbql5lND/8gcpTGxnS8shRvyd32fWvbHj0UbdtaLV9LLiqs93vsEyL2dL+r9m4wErECH8EgIQsvZqbgwhbelKA9MUI+zawYkC86c1UAM4uyaXq/0azJY65pVuKHCwB6PzqX3Kw6fokZ6EQE25CdViAQm/ZX5lK3fW/Za7Ztbdu3fWu1lbGIfgpBe8fI3YkFxqbH4UIopXqN3RdB52W5fwezFIFpVzWi3CyUT8xYRm+qf6gQM+ukC1In+cGRDD3w3Oo4Ony5mONCS69AsE6xV6CSbyhIJobteGURbzzcsnl9j8bU01xHXasDfnCVmS+WNgUJl9gEF3djNJ4ze3MBZ9fZWmhJu1prrnm4LZs6qc4YVeks7aRDWI6BgPV+aJ5GPkfn4atorTvqNqu1EA/H5B9sOK6+4npzYbhqtG3btr1FA7ZzjNEw70HLOFgPI4OySpe8NCICjSukXO5IsVfsYdMCcvPK4pV9BW2eb+UaNqd74jEsUvDC1tHNR+iGsrqnETGCQl0IICo8Yo1nckYG4k7J/O6+BLHtbXeH0V4Oo8q1bXtwOQEU0nPFf9kc82diLpXE92AJFSniB5v42bhLEHYsQbPtZkrBJwVHNKVhoZxmkTKAR4ZDH8uQw3MoYOJTyAZCGbv/aLkU//O4shcC3Ens+7bte/v3H1+OKrVtu5mfGZXZQz5EIgJbPAHJ7aB1cx7alyYqE1K9lwZEmKo/9ukvvcbMgjzE+asdiWPPracq/Ti6Cs2N+5FvMLEULyKLT9jDmh+TnG1muW12Quh+Y/GThz1+btu//denW+UReZLbxGLY019FbVlKa7W0wkw2Vsg6v5wiR7Tj+f4gxKSC6/ZoSjkjpfgbL19x+db70XX6mejlCuMhX4c89ixoIW8bD7y9PIN4wmxGB7H4ItF92/7vn74cNbJiZJi3GeWzWpBCiCCX1mlEHlXNrSH3gQzPZUzpFCiKaxGGWMIelniQf0beGDHWz/cSXtKt7KDwvF7FuQ71VSTsu2uK1JqdNEP6ywQJqN9fRUpr//6nTzE1H947Z6+mTUyVTTohoBzEyfPIOS2qyKGszFlFyqVkUwqjMJ/lOQdYc4XoWCIJIM84y1mA8LacC0ekFBNf2EIEskLWVW2c7+1F0+EYEpoxPALwlX/6+eWIzZzpTpz9zv2JSddhJD1SfHFQLlsno0ggmdjigc61pUjUuB5yNuV7OqlqaEXUaAhwXsPSKwYPHjftO6lQvFcYxp3IhNGPQ23ufvRqcs44Lj560AZmlRcT9a5SqA0MPFHM1NdhoM4x9d67EbEMiw+dXtxr4twYZxiSWRrmRsogPGj82MpNwdrGs/CBcTMnQKV4qx0DZkoCgt4OBZfSWomTdkotpVThqRE85MFEdSLvyIzWhxUCwslniTBZLHow2NGDO8sFXFPGuSxqxDQSI557/jEZ86wzD3aSgu1VUBxQRbNcnztPeXRUBOvkFYJ+9MMgxY+Zj0GOrLHH6UtTFxws5TOOqc1UG8oDB6JnJ83B96tGgo7jONRIKNaIzXP2zEyQLU9zdmNkaTSa50cLw1C1qLVY17GJpawBxZUEEnw2A9kPJVJFRMi0W2xI9d2e8UctpdZcZIs4aC81YQDh8eBGzIrPzYMhWUoVGuurTT1vSOwbv+W+xYN7RjRmG69bg1dqwDI/6dS/xrreMg/GixxJtXeFlBJeLqkGKAvXSOR9/jaem/fa+Ar1wHpjaW8dbmp2sYe5IpdAcgACh8qllpjwlRxuH64/QRUzw8wsaVIi34Ing0kZX/Ng72T0/Z2cgNZI1SbnGn1wRDBl0yyMDnrbmKXuTgQbgUjtIBJ2/xguwkNHLaWIrycN5ngJ2JIHGk655CkkuWI21Sf5FuIY5ErWMB3LiXKd3NH5izOLSFNkX0ajasjD3rJ0EyAhn5EKAzD2xQsAkTFz3WHAnSw/V6nfWbzY0dw6indjvty6VWJhsWwSWmLWwo4MtJCu33fam3/DgFx8TRz7UYKEmsRByizjC+j01otDQEauAGvzHL6HJmzzhaWcBwb7ihKKFl0ioqMPLAICqV9fKXE8W6211M+fr320+sM1No9gTxyYcG/GFval7IMecncatkpRaU0qfenqCATtCaBDEFeh6TOyYIPcvzh8xYzjM2hBybxjHEawGNpi8YkBl+jok0vtJI9kbhil1nZ9vWtNZ72ExtyBmH5oZKd+S6bwWogscFAyT8lkc62/nNqL/P4596pH+j/hC3lPj6paLrWloDdzI6aHaZgv5XChkTGzQVCEWLhsABHfdXSGpgYRQNT9EIpSaz/uvU7zXb6G+mYQB8+TA0eJ2CMZbJp7trxSMsVxC6nmvl1uepsR9UZaFEE1jnbk+S45pugxwrc7G8MwlmATmDRuUrhsRCxydI1YPtPJIG/VGyVMtcbnkY1jyMceQMr3dlOJx+zFkNysBlN/tXcvDY54tjqO94grsBGFlxpAhrms9qVsJUBShkcaomJmJiMenR+UxVRiYyYujaTU23GEMixlmCgKASCFGaqf1oRgOjJUDJyQYC47U10PVU3XnvXAi+n18hjSMlo/wilEWpAEzqBS5yrnOJQvFgCnIlE6D8538FMpDAxzIhMj6memVEhKqbdy+AEqI2OZr0jKpyYhECl53KpH+MmJsBcPfYIFZupNKnE5GIwD8vitcVhquhzL9wmtPvkJJAqJoQuKo9J8r+/AbMnMQPyHRubzxCNlDnYRgDB57229H900OpRP0xyDi6xnT4DFleZCl6VIPlKH4BdnYuwzjcbj7LFSg6LNPXojoPgjH09jctGU8TE9rNQSTP3gKCIb953G4G7jPf2QCIsHLEQMSbjZY8DgPNEyltnVcHOE0XxOmffPSwzal7yi5HmWALFHnQYuDKrYm+BrbB40yqmB5XAwGkp2hgijC4DimJfoDfNlnkwWsMLxNIEVURMxFikQEYjA13yCpcYDVIaRjRsPKcQ11AznMxwsMkjmc8n1vckI4RrDfwWrhOg8KJGrVBGCaFBSSZQsCBEnRaAcCIlQ41snswhMOqxiluMkRlWIiNnEpLBonotMQiAp1SGHMYslG3wWQxyuPZDiaJ3GRDphDvEsTUYeHpo7Usg8msGPwa3V+yE5NsWA3eGbgJbO75m7Ue5aYxCDuRQRX8lEmY3KYCXShJmELDpbTMRYRLgUTyb9YKUooTKT+IFej1+VEiOOjTbJOT682ltUxjOjWIrk0Xw43GEO4RkJxgozDkBKREySYD/uPhqZHtAd+ylW40CzgZaGKU3+Mv5m8ENcSxGN3DOQaJKEiwpMDFfzdjN3HA7I2YHZEscRITiq0kzef2lEUjKaRvOeUzmleKNW7OMjySxqeHIe6GBOICdEYY6FPuMVY+Sac150ag5NQotAZmKmfgK8B+CEtidTCNg+t+skGbxgdGCcFJys11g1wn6uUzDx0UMSTsNP53Yeh5gAkeL1FhCHe8oEeva10BKa4hQ9kVF0zPCUPhUPx3YtcY+IyIggIigWMH+RwZysYiKBGNWktyLFy1iZp2ZhoKVANEbe3zognMDbMi1cenQeukYwgQVSrABinHOhccsBfnh65sgsTIxpyiCXfbMwjb0Cc2gvLnxckQNYv0uheUYrR+6f6bgviQZZDY6ZPZRHHXWFjfNBSR5ylw1/Oe4fJ6XS0KWlnhplNUiEeBaIzXQrKiTZwpr/znP/EFWUEMSsYo/9pUMFBhALaBqSCDZyVZU5uV8MRBWEmNg9QYORNsR7M3nrYYmNIWAiYaiCY/NsEuyzWBNkrNfLBErELGRiYpSzs3GvGSwyC3CkRGRG2SKzZNOnjO9cz1uC7rqc4nTWII2ZCr8e1CzL+KMf3ibb7JNmMYavVONhpbDYMUnZ5pKyS58yop9AxARCzGLKluyKq3p0TMkkaYmYpEg46/BG853DBDCD2Io904n78DXFVsjwvBPwZJccGdeIiPFgkAl96CTx6G4k4sI5CwswwfpxKBUWhmIpRi7eF7Nvw5fcM4SjUW84OhjBWDBEAAg4z48OMsAgY1JgeNYlP17VHZG2esbNAEhScmENnLTG4Bizi2B9mu5BltoiSy21LK1YsON6mE909kP9iM8J8gHAZAYNQIxJQKyCMSzizWlxGCVn2z8ZiQipS1y9Ahh7vMnGosJhDJPnpvHMOINjPPglgoDYs7641br4FvBE9xTue1YlpLW2VVZKchdde3cHoUdXMPuhJYOqMp97noSTp52xjZajHY1gXfJEV3eGqkKSCjurFBbx3lLj5kKOZF+nUQAgMiLO7CYWk8ZL8ukABswFdPm7Y95iFh48FJZSiyMfH3wkK0FOaHduPPKrjDH5YLPfmeZJcpwcOQhGwgIPf8GaeYFiuE4apEjcKTNb0v3j4WY2vmAiDgLYYdByjhR7q7V5JKwYb80JKMIWJdom492Kn1mXY08ixNKI2lYMpt5LkjhSjM3EMg0J4kUmH0s82geCLx1MK0wna++/Ycy2EkhZc6NhwjFt6FNW+RIfLXD9kmgTw8J6BuaF5VlTmQ/lJv6hGZ7MgKW0WgupeM8Qi0DrJlQ3Vp2aE0FfYmnS9NzTBYzOXW8YNDEpTlRLnKN9kPoUjg6VoUGILSQb8unRcrkpgNFIwVzCeQohtJ6D43F3WSMnHtZHa0hFJMHEUrZ93/wU2HBsIJJCtTBKNT/2jsAhAjbTmYLDzSE7HEdYcGfpliKxS8OiId5g2iGFJLcaDnopk6LJzsX/D6qGmUWKlBiy9QINoq4Qp7MYnJdA9DaH1xgw9E0KSVRqa028A9sUZF1vh/m+UT/lqAgZjLTAxGLmP1KfxMgrSTWEzuLLZ/zJmXr7HAA7Dgi4BCrO1kX4makyItHiDzFSrVKqFCEQpLDBIDCMIlNkiwxmYkTxJRHy2Ad6FgNLrbX4YdUGaFdT9NvBUmM/klApoVwuAzaxgbuDh81qrcSH+UiPgBwW+FIVmBqxGEH7YSWnyHP3yiRi0yEiTTaxDclszVCWQkKmUBMQx3GvHMiZ1ywyrtZGxX3k6mGNUor4eQzaj6N3U/S7lpanabPvFxfnGFgAGMfBbyduStbJoAhwgqhzkun9dnT15+RTIGBhLum3Tq4/08dRVxkod0Y1Ii7EUJNuIAPHEljKvZjMPg0XzooEPDKALJ0mboZ2Pqwf/TiOww+r7NS4BLkEKUwicO/OxmyBbiQS7+EhneqgpH9Z2KKeAb3frjcF5ThoxHj2RjUHSxHBgtpYyJFJpGdsHTUe1yO13A477tf9Y01bTf+eLnJoiL+v9RsV+NyeKwIMzNVX7JMXQZWZAVMWYzaNp0UczbIx/S1xdmC0xvmR00yEftyv11v3SmLwBaSZlOU4qcOmDJ8DOcyYsRQBDZoxHcKiROxyCNlEEK9BcM2CMxIYD+KBiAn92kl778cc1kKBmSqBRPy0bxaCn91C3lPH7GcORw0DNFhC9xHwTnkmsn7cb9fr7SBpbnhdjQhknUVKycYlf8NErePpL8YQbYyBJTCarAQorMPaOYp59Ob41KXLaIZ+EIisC/Q4fH46jMrbwH32mx3WShIYzEpEZSa+MluUZp7FxOLmf7/fbq/X64G6gcnQjzF508WPJfbGl4RyPj2RvKDXCiMJD5WwbMVEepDoH5MAMZlaeBaZI1PA3D2ZHo2IiIysM+LpLN1DsC5ZLcu8V02znCyzCsiDKhkBMvgtAvr9fr1er9drh0ehOXUMAvWj1FKGr8qm8chFKOAcJpKIhXKmsMntCEDFH4bblFeHmbnSwMaIOa03uuEaNPDK7CxihvaciIlCtSgUahRH2BPBN04iw/BqyM7C9uN2vV6v19u9Uw1getxu9x5XA+1HrcGJBUFlRKWU2gf1N69XRErhGEHwVI8oCTkeoSO6zISJqumRCV4QPHB7dKXGIM4iRi/wwZdfaC7kc0rOj3ZSWD98Za9QTLYN1juo5ig0mbke3O6qYKAfdEDvt3v2mRCsi7D1XJQadSi93w/VuBn/n7CwmBbtRVwTmEvOSMa5ZHEGjddAiBRAtc6wOJwx8KLFCiBEaSnRMw8yNSVhRDCVsZmciaGqSgbtR82eFckVdMs/+UwJ2u+32+3uJ+hC9fDJm0PHAzA9CHosB23CXWk/lkL/gNOSXLsaJE4YCWzmHitX7bEQm4EqlKzzAKRE2YKBWSCgTCHc582SJRtLHzMp7kNM1WDQXrMfeZDKA3oRjQIxoP24Hz14QOhhxZVxqRyaErRkV2u4H5hqWsyI8aEP46iQ6C6KDS1ReR/VK/UH65IbvNZIJmcxaPHkNBKscS+D5MV8jcE7cSU515MCDFtKROKdOWnQQazkHuvF0PNc4RG4vLVrieRTHwZF7qT48EJY/oIBvv5/ycPUvJ7p2+0AAAAASUVORK5CYII=",
"text/plain": [
"260×260 Array{Gray{Float64},2} with eltype Gray{Float64}:\n",
" Gray{Float64}(0.0169412) … Gray{Float64}(0.0323137) \n",
" Gray{Float64}(0.0343529) Gray{Float64}(0.0655686) \n",
" Gray{Float64}(0.0520784) Gray{Float64}(0.0997647) \n",
" Gray{Float64}(0.0698039) Gray{Float64}(0.136) \n",
" Gray{Float64}(0.0872157) Gray{Float64}(0.171451) \n",
" Gray{Float64}(0.0872157) … Gray{Float64}(0.174745) \n",
" Gray{Float64}(0.0865882) Gray{Float64}(0.177412) \n",
" Gray{Float64}(0.0856471) Gray{Float64}(0.17349) \n",
" Gray{Float64}(0.0848627) Gray{Float64}(0.159843) \n",
" Gray{Float64}(0.0850196) Gray{Float64}(0.145725) \n",
" Gray{Float64}(0.0858039) … Gray{Float64}(0.131137) \n",
" Gray{Float64}(0.0861176) Gray{Float64}(0.116078) \n",
" Gray{Float64}(0.0864314) Gray{Float64}(0.111529) \n",
" ⋮ ⋱ \n",
" Gray{Float64}(0.0329412) Gray{Float64}(0.00831373)\n",
" Gray{Float64}(0.0332549) Gray{Float64}(0.00941176)\n",
" Gray{Float64}(0.0334118) … Gray{Float64}(0.00941176)\n",
" Gray{Float64}(0.0345098) Gray{Float64}(0.00941176)\n",
" Gray{Float64}(0.0357647) Gray{Float64}(0.00972549)\n",
" Gray{Float64}(0.0345098) Gray{Float64}(0.00956863)\n",
" Gray{Float64}(0.0326275) Gray{Float64}(0.00988235)\n",
" Gray{Float64}(0.032) … Gray{Float64}(0.0108235) \n",
" Gray{Float64}(0.025098) Gray{Float64}(0.00909804)\n",
" Gray{Float64}(0.0175686) Gray{Float64}(0.00690196)\n",
" Gray{Float64}(0.0112941) Gray{Float64}(0.00454902)\n",
" Gray{Float64}(0.00596078) Gray{Float64}(0.00235294)"
]
},
"execution_count": 12,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"kernalB = ones(5,5)/25\n",
"blurimage = matCov(image,kernalB)\n",
"Gray.(blurimage)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## LALFD 4.1 Fourier Transforms:Discrete and Continuous "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Load a discrete series f"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 13,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Intercept a signal of one period\n",
"f = ECGsample[500:699]\n",
"plot(f,legend = false)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Fourier Matrix F and DFT Matrix Ω"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"N×N Fourier matrix F contains powers of w = exp(2πi/N)"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"Fmat (generic function with 1 method)"
]
},
"execution_count": 14,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Function to build a Fourier matrix\n",
"function Fmat(N)\n",
" Fmat = ones(ComplexF64, N, N)\n",
" w = exp(2im*pi/N)\n",
" for i in 1:N\n",
" for j in 1:N\n",
" m = (i-1)*(j-1)\n",
" Fmat[i,j] = w^m\n",
" end\n",
" end\n",
" Fmat\n",
"end"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"N×N DFT matrix Ω contains powers of ω = exp(-2πi/N)"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"DFTmat (generic function with 1 method)"
]
},
"execution_count": 15,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Function to build a DFT matrix\n",
"function DFTmat(N)\n",
" Dmat = ones(ComplexF64, N, N)\n",
" w = exp(-2im*pi/N)\n",
" for i in 1:N\n",
" for j in 1:N\n",
" m = (i-1)*(j-1)\n",
" Dmat[i,j] = w^m\n",
" end\n",
" end\n",
" Dmat\n",
"end"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### DFT and inverse-DFT"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 16,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# generate the Fourier matrix and DFT matrix\n",
"Ω = DFTmat(200)\n",
"F = Fmat(200)\n",
"\n",
"# Ω times f produces discrete Fourier coefficients c\n",
"# That is also the discrete Fourier Transform\n",
"c = Ω * f\n",
"plot(c,legend = false)"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 17,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# F times c bring back f\n",
"# That is also inverse Fourier transform\n",
"a = F*c\n",
"plot(real(round.(a,digits = 5)),legend = false)"
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"200×200 Array{Complex{Float64},2}:\n",
" 200.0+0.0im 0.0+0.0im -0.0+0.0im … 0.0-0.0im 0.0-0.0im\n",
" 0.0-0.0im 200.0+0.0im -0.0+0.0im 0.0-0.0im 0.0-0.0im\n",
" -0.0-0.0im -0.0-0.0im 200.0+0.0im 0.0-0.0im 0.0-0.0im\n",
" -0.0-0.0im -0.0-0.0im -0.0-0.0im -0.0-0.0im -0.0-0.0im\n",
" -0.0-0.0im -0.0-0.0im -0.0-0.0im -0.0-0.0im -0.0-0.0im\n",
" -0.0-0.0im -0.0-0.0im -0.0-0.0im … -0.0-0.0im -0.0-0.0im\n",
" -0.0-0.0im -0.0-0.0im -0.0-0.0im -0.0-0.0im -0.0-0.0im\n",
" -0.0-0.0im -0.0-0.0im -0.0-0.0im -0.0-0.0im -0.0-0.0im\n",
" -0.0-0.0im -0.0-0.0im -0.0-0.0im -0.0-0.0im -0.0-0.0im\n",
" -0.0-0.0im -0.0-0.0im -0.0-0.0im -0.0-0.0im -0.0-0.0im\n",
" -0.0-0.0im -0.0-0.0im -0.0-0.0im … -0.0-0.0im -0.0-0.0im\n",
" -0.0-0.0im -0.0-0.0im -0.0-0.0im -0.0-0.0im -0.0-0.0im\n",
" -0.0-0.0im -0.0-0.0im -0.0-0.0im -0.0-0.0im -0.0-0.0im\n",
" ⋮ ⋱ \n",
" -0.0+0.0im -0.0+0.0im -0.0+0.0im -0.0+0.0im -0.0+0.0im\n",
" -0.0+0.0im -0.0+0.0im -0.0+0.0im -0.0+0.0im -0.0+0.0im\n",
" -0.0+0.0im -0.0+0.0im -0.0+0.0im … -0.0+0.0im -0.0+0.0im\n",
" -0.0+0.0im -0.0+0.0im -0.0+0.0im -0.0+0.0im -0.0+0.0im\n",
" -0.0+0.0im -0.0+0.0im -0.0+0.0im -0.0+0.0im -0.0+0.0im\n",
" -0.0+0.0im -0.0+0.0im -0.0+0.0im -0.0+0.0im -0.0+0.0im\n",
" -0.0+0.0im -0.0+0.0im -0.0+0.0im -0.0+0.0im -0.0+0.0im\n",
" -0.0+0.0im -0.0+0.0im -0.0+0.0im … -0.0+0.0im -0.0+0.0im\n",
" -0.0+0.0im -0.0+0.0im -0.0+0.0im -0.0+0.0im -0.0+0.0im\n",
" 0.0+0.0im -0.0+0.0im -0.0+0.0im -0.0+0.0im -0.0+0.0im\n",
" 0.0+0.0im 0.0+0.0im 0.0+0.0im 200.0+0.0im -0.0+0.0im\n",
" 0.0+0.0im 0.0+0.0im 0.0+0.0im -0.0-0.0im 200.0+0.0im"
]
},
"execution_count": 18,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# FΩ = NI\n",
"round.(F * Ω,digits = 5)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### The DFT Matrix Ω is a Permutation of the Fourier Matrix F"
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"Pmat (generic function with 1 method)"
]
},
"execution_count": 19,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Function to build a Permutation matrix \n",
"function Pmat(N)\n",
" P = zeros(N,N)\n",
" P[1,1] = 1\n",
" for i in 2:N\n",
" P[N+2-i,i] = 1\n",
" end\n",
" P\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"200×200 Array{Complex{Float64},2}:\n",
" 0.0+0.0im 0.0+0.0im 0.0+0.0im … 0.0+0.0im 0.0+0.0im 0.0+0.0im\n",
" 0.0+0.0im 0.0-0.0im 0.0-0.0im -0.0-0.0im -0.0-0.0im -0.0-0.0im\n",
" 0.0+0.0im 0.0-0.0im 0.0-0.0im -0.0-0.0im -0.0-0.0im -0.0-0.0im\n",
" 0.0+0.0im 0.0-0.0im 0.0-0.0im -0.0-0.0im -0.0-0.0im -0.0-0.0im\n",
" 0.0+0.0im 0.0-0.0im 0.0-0.0im -0.0-0.0im -0.0-0.0im -0.0-0.0im\n",
" 0.0+0.0im 0.0-0.0im 0.0-0.0im … -0.0-0.0im -0.0-0.0im -0.0-0.0im\n",
" 0.0+0.0im 0.0-0.0im 0.0-0.0im -0.0-0.0im -0.0-0.0im -0.0-0.0im\n",
" 0.0+0.0im 0.0-0.0im 0.0-0.0im -0.0-0.0im -0.0-0.0im -0.0-0.0im\n",
" 0.0+0.0im 0.0-0.0im 0.0-0.0im -0.0-0.0im -0.0-0.0im -0.0-0.0im\n",
" 0.0+0.0im 0.0-0.0im 0.0-0.0im -0.0-0.0im -0.0-0.0im -0.0-0.0im\n",
" 0.0+0.0im 0.0-0.0im 0.0-0.0im … -0.0-0.0im -0.0-0.0im -0.0-0.0im\n",
" 0.0+0.0im 0.0-0.0im 0.0-0.0im -0.0-0.0im -0.0-0.0im -0.0-0.0im\n",
" 0.0+0.0im 0.0-0.0im 0.0-0.0im -0.0-0.0im -0.0-0.0im -0.0-0.0im\n",
" ⋮ ⋱ \n",
" 0.0+0.0im 0.0+0.0im 0.0+0.0im -0.0+0.0im -0.0+0.0im -0.0+0.0im\n",
" 0.0+0.0im 0.0+0.0im 0.0+0.0im -0.0+0.0im -0.0+0.0im -0.0+0.0im\n",
" 0.0+0.0im 0.0+0.0im 0.0+0.0im … -0.0+0.0im -0.0+0.0im -0.0+0.0im\n",
" 0.0+0.0im 0.0+0.0im 0.0+0.0im -0.0+0.0im -0.0+0.0im -0.0+0.0im\n",
" 0.0+0.0im 0.0+0.0im 0.0+0.0im -0.0+0.0im -0.0+0.0im -0.0+0.0im\n",
" 0.0+0.0im 0.0+0.0im 0.0+0.0im -0.0+0.0im -0.0+0.0im -0.0+0.0im\n",
" 0.0+0.0im 0.0+0.0im 0.0+0.0im -0.0+0.0im -0.0+0.0im -0.0+0.0im\n",
" 0.0+0.0im 0.0+0.0im 0.0+0.0im … -0.0+0.0im -0.0+0.0im -0.0+0.0im\n",
" 0.0+0.0im 0.0+0.0im 0.0+0.0im -0.0+0.0im -0.0+0.0im -0.0+0.0im\n",
" 0.0+0.0im 0.0+0.0im 0.0+0.0im -0.0+0.0im -0.0+0.0im -0.0+0.0im\n",
" 0.0+0.0im 0.0-0.0im 0.0+0.0im -0.0+0.0im -0.0+0.0im -0.0-0.0im\n",
" 0.0+0.0im 0.0-0.0im 0.0-0.0im -0.0+0.0im -0.0-0.0im -0.0-0.0im"
]
},
"execution_count": 20,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"P = Pmat(200)\n",
"\n",
"# FP = Ω\n",
"round.(F*P-Ω,digits = 5)"
]
},
{
"cell_type": "code",
"execution_count": 21,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"200×200 Array{Complex{Float64},2}:\n",
" 0.0+0.0im 0.0+0.0im 0.0+0.0im … 0.0+0.0im 0.0+0.0im 0.0+0.0im\n",
" 0.0+0.0im 0.0+0.0im 0.0+0.0im -0.0+0.0im -0.0+0.0im -0.0+0.0im\n",
" 0.0+0.0im 0.0+0.0im 0.0+0.0im -0.0+0.0im -0.0+0.0im -0.0+0.0im\n",
" 0.0+0.0im 0.0+0.0im 0.0+0.0im -0.0+0.0im -0.0+0.0im -0.0+0.0im\n",
" 0.0+0.0im 0.0+0.0im 0.0+0.0im -0.0+0.0im -0.0+0.0im -0.0+0.0im\n",
" 0.0+0.0im 0.0+0.0im 0.0+0.0im … -0.0+0.0im -0.0+0.0im -0.0+0.0im\n",
" 0.0+0.0im 0.0+0.0im 0.0+0.0im -0.0+0.0im -0.0+0.0im -0.0+0.0im\n",
" 0.0+0.0im 0.0+0.0im 0.0+0.0im -0.0+0.0im -0.0+0.0im -0.0+0.0im\n",
" 0.0+0.0im 0.0+0.0im 0.0+0.0im -0.0+0.0im -0.0+0.0im -0.0+0.0im\n",
" 0.0+0.0im 0.0+0.0im 0.0+0.0im -0.0+0.0im -0.0+0.0im -0.0+0.0im\n",
" 0.0+0.0im 0.0+0.0im 0.0+0.0im … -0.0+0.0im -0.0+0.0im -0.0+0.0im\n",
" 0.0+0.0im 0.0+0.0im 0.0+0.0im -0.0+0.0im -0.0+0.0im -0.0+0.0im\n",
" 0.0+0.0im 0.0+0.0im 0.0+0.0im -0.0+0.0im -0.0+0.0im -0.0+0.0im\n",
" ⋮ ⋱ \n",
" 0.0+0.0im 0.0-0.0im 0.0-0.0im -0.0-0.0im -0.0-0.0im -0.0-0.0im\n",
" 0.0+0.0im 0.0-0.0im 0.0-0.0im -0.0-0.0im -0.0-0.0im -0.0-0.0im\n",
" 0.0+0.0im 0.0-0.0im 0.0-0.0im … -0.0-0.0im -0.0-0.0im -0.0-0.0im\n",
" 0.0+0.0im 0.0-0.0im 0.0-0.0im -0.0-0.0im -0.0-0.0im -0.0-0.0im\n",
" 0.0+0.0im 0.0-0.0im 0.0-0.0im -0.0-0.0im -0.0-0.0im -0.0-0.0im\n",
" 0.0+0.0im 0.0-0.0im 0.0-0.0im -0.0-0.0im -0.0-0.0im -0.0-0.0im\n",
" 0.0+0.0im 0.0-0.0im 0.0-0.0im -0.0-0.0im -0.0-0.0im -0.0-0.0im\n",
" 0.0+0.0im 0.0-0.0im 0.0-0.0im … -0.0-0.0im -0.0-0.0im -0.0-0.0im\n",
" 0.0+0.0im 0.0-0.0im 0.0-0.0im -0.0-0.0im -0.0-0.0im -0.0-0.0im\n",
" 0.0+0.0im 0.0-0.0im 0.0-0.0im -0.0-0.0im -0.0-0.0im -0.0-0.0im\n",
" 0.0+0.0im 0.0+0.0im 0.0-0.0im -0.0-0.0im -0.0-0.0im -0.0+0.0im\n",
" 0.0+0.0im 0.0+0.0im 0.0+0.0im -0.0-0.0im -0.0+0.0im -0.0+0.0im"
]
},
"execution_count": 21,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# PΩ = F\n",
"round.(Ω*P-F,digits = 5)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Fast Fourier Transform"
]
},
{
"cell_type": "code",
"execution_count": 22,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"Comat (generic function with 1 method)"
]
},
"execution_count": 22,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Combine matrix - C\n",
"using LinearAlgebra\n",
"function Comat(N)\n",
" n = floor(Int64,N/2)\n",
" hI = Matrix{Float64}(I, n, n)\n",
" \n",
" DDmat = zeros(ComplexF64,n,n)\n",
" w = exp(2im*pi/N)\n",
" for i in 1:n\n",
" DDmat[i,i] = w^(i-1)\n",
" end\n",
" Cm1 = hcat(hI,DDmat)\n",
" Cm2 = hcat(hI,-DDmat)\n",
" Cm = vcat(Cm1,Cm2)\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 23,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"HalfF (generic function with 1 method)"
]
},
"execution_count": 23,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Half-size transform matrix - H\n",
"function HalfF(N)\n",
" n = floor(Int64,N/2)\n",
" hF = Fmat(n)\n",
" zm = zeros(n,n)\n",
" HF1 = hcat(hF,zm)\n",
" HF2 = hcat(zm,hF)\n",
" HF = vcat(HF1,HF2)\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 24,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"PEOmat (generic function with 1 method)"
]
},
"execution_count": 24,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Even-Odd permutation matrix P - PP\n",
"function PEOmat(N)\n",
" P = zeros(N,N)\n",
" n = floor(Int64,N/2)\n",
" for i in 1:n\n",
" P[i,2*i-1] = 1\n",
" end\n",
" for i in n+1:N\n",
" P[i,2*(i-n)] = 1\n",
" end\n",
" P\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 25,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"1024×1024 Array{Complex{Float64},2}:\n",
" 0.0+0.0im 0.0+0.0im 0.0+0.0im … 0.0+0.0im 0.0+0.0im 0.0+0.0im\n",
" 0.0+0.0im 0.0+0.0im 0.0+0.0im -0.0+0.0im 0.0+0.0im 0.0-0.0im\n",
" 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0-0.0im 0.0+0.0im -0.0+0.0im\n",
" 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im -0.0+0.0im\n",
" 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0-0.0im 0.0+0.0im 0.0-0.0im\n",
" 0.0+0.0im 0.0+0.0im 0.0+0.0im … 0.0+0.0im 0.0+0.0im 0.0-0.0im\n",
" 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0-0.0im 0.0+0.0im 0.0+0.0im\n",
" 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im\n",
" 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im\n",
" 0.0+0.0im 0.0+0.0im 0.0+0.0im -0.0+0.0im 0.0+0.0im 0.0+0.0im\n",
" 0.0+0.0im 0.0+0.0im 0.0+0.0im … -0.0-0.0im 0.0+0.0im 0.0+0.0im\n",
" 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im\n",
" 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0-0.0im 0.0+0.0im 0.0-0.0im\n",
" ⋮ ⋱ \n",
" 0.0+0.0im 0.0-0.0im 0.0-0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im\n",
" 0.0+0.0im 0.0-0.0im 0.0-0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im\n",
" 0.0+0.0im 0.0-0.0im 0.0-0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im\n",
" 0.0+0.0im 0.0-0.0im 0.0-0.0im … 0.0+0.0im 0.0+0.0im 0.0+0.0im\n",
" 0.0+0.0im 0.0-0.0im 0.0-0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im\n",
" 0.0+0.0im 0.0-0.0im 0.0-0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im\n",
" 0.0+0.0im 0.0-0.0im 0.0-0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im\n",
" 0.0+0.0im 0.0-0.0im 0.0-0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im\n",
" 0.0+0.0im 0.0+0.0im 0.0-0.0im … 0.0+0.0im 0.0+0.0im 0.0+0.0im\n",
" 0.0+0.0im 0.0+0.0im 0.0-0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im\n",
" 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im\n",
" 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im"
]
},
"execution_count": 25,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# F = C*H*PP\n",
"C = Comat(1024)\n",
"H = HalfF(1024)\n",
"PP = PEOmat(1024)\n",
"FFT = C*H*PP\n",
"F1024 = Fmat(1024)\n",
"round.(FFT - F1024)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Full FFT by Recursion"
]
},
{
"cell_type": "code",
"execution_count": 26,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"fFFT (generic function with 1 method)"
]
},
"execution_count": 26,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Funtion for full FFT\n",
"function fFFT(N)\n",
" n = floor(Int64,N/2)\n",
" C = Comat(N)\n",
" PP = PEOmat(N)\n",
" \n",
" if N>4\n",
" hf = fFFT(n)\n",
" else\n",
" hf = Fmat(2)\n",
" end\n",
" \n",
" Z = zeros(n,n)\n",
" Hf1 = hcat(hf,Z)\n",
" Hf2 = hcat(Z,hf)\n",
" H = vcat(Hf1,Hf2)\n",
" \n",
" FFT = C*H*PP\n",
"end "
]
},
{
"cell_type": "code",
"execution_count": 27,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"1024×1024 Array{Complex{Float64},2}:\n",
" 0.0+0.0im 0.0+0.0im 0.0+0.0im … 0.0+0.0im 0.0+0.0im 0.0+0.0im\n",
" 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im\n",
" 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0-0.0im 0.0+0.0im 0.0+0.0im\n",
" 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0-0.0im 0.0-0.0im 0.0+0.0im\n",
" 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0-0.0im 0.0-0.0im 0.0+0.0im\n",
" 0.0+0.0im 0.0+0.0im 0.0+0.0im … 0.0-0.0im 0.0-0.0im 0.0-0.0im\n",
" 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0-0.0im 0.0-0.0im 0.0-0.0im\n",
" 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0-0.0im 0.0-0.0im 0.0-0.0im\n",
" 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0-0.0im 0.0-0.0im 0.0-0.0im\n",
" 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0-0.0im 0.0-0.0im 0.0-0.0im\n",
" 0.0+0.0im 0.0+0.0im 0.0+0.0im … 0.0-0.0im 0.0-0.0im 0.0-0.0im\n",
" 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0-0.0im 0.0-0.0im 0.0-0.0im\n",
" 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0-0.0im 0.0-0.0im 0.0-0.0im\n",
" ⋮ ⋱ \n",
" 0.0+0.0im 0.0-0.0im 0.0-0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im\n",
" 0.0+0.0im 0.0-0.0im 0.0-0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im\n",
" 0.0+0.0im 0.0-0.0im 0.0-0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im\n",
" 0.0+0.0im 0.0-0.0im 0.0-0.0im … 0.0+0.0im 0.0+0.0im 0.0+0.0im\n",
" 0.0+0.0im 0.0-0.0im 0.0-0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im\n",
" 0.0+0.0im 0.0-0.0im 0.0-0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im\n",
" 0.0+0.0im 0.0-0.0im 0.0-0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im\n",
" 0.0+0.0im 0.0-0.0im 0.0-0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im\n",
" 0.0+0.0im 0.0+0.0im 0.0-0.0im … 0.0+0.0im 0.0+0.0im 0.0+0.0im\n",
" 0.0+0.0im 0.0+0.0im 0.0-0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im\n",
" 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im\n",
" 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im"
]
},
"execution_count": 27,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# test\n",
"fFFT(1024)\n",
"round.(fFFT(1024) - F1024)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## LALFD 4.2 Shift Matrices and Circulant Matrices"
]
},
{
"cell_type": "code",
"execution_count": 28,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"Shiftm (generic function with 1 method)"
]
},
"execution_count": 28,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Function to build a shift matrix\n",
"function Shiftm(N)\n",
" mat = zeros(N,N)\n",
" for i in 1:N\n",
" if i == 1\n",
" mat[N,i] = 1\n",
" else\n",
" mat[i-1,i] = 1\n",
" end\n",
" end\n",
" mat\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 29,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"6-element Array{Float64,1}:\n",
" 3.0\n",
" 4.0\n",
" 5.0\n",
" 6.0\n",
" 1.0\n",
" 2.0"
]
},
"execution_count": 29,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"Ps = Shiftm(6)\n",
"x = [1,2,3,4,5,6]\n",
"Ps^2*x"
]
},
{
"cell_type": "code",
"execution_count": 30,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"Cirmat (generic function with 1 method)"
]
},
"execution_count": 30,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Function to build a Circulant matrix\n",
"function Cirmat(a)\n",
" n = length(a)\n",
" P = Shiftm(n)\n",
" Id = Matrix{Float64}(I, n, n)\n",
" mat = a[1]*Id\n",
" for i in 2:n\n",
" mat = mat+a[i]*P^(i-1)\n",
" end\n",
" mat\n",
"end "
]
},
{
"cell_type": "code",
"execution_count": 31,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"5×5 Array{Float64,2}:\n",
" 1.0 2.0 3.0 4.0 5.0\n",
" 5.0 1.0 2.0 3.0 4.0\n",
" 4.0 5.0 1.0 2.0 3.0\n",
" 3.0 4.0 5.0 1.0 2.0\n",
" 2.0 3.0 4.0 5.0 1.0"
]
},
"execution_count": 31,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"vecc = [1,2,3,4,5]\n",
"C = Cirmat(vecc)"
]
},
{
"cell_type": "code",
"execution_count": 32,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"5×5 Array{Float64,2}:\n",
" 67.0 94.0 81.0 78.0 85.0\n",
" 85.0 67.0 94.0 81.0 78.0\n",
" 78.0 85.0 67.0 94.0 81.0\n",
" 81.0 78.0 85.0 67.0 94.0\n",
" 94.0 81.0 78.0 85.0 67.0"
]
},
"execution_count": 32,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Product of two Circulant matrices is also circulant\n",
"vecd = [9,0,8,6,4]\n",
"D = Cirmat(vecd)\n",
"C*D\n",
"# C*D == D*C"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Cyclic Convolution"
]
},
{
"cell_type": "code",
"execution_count": 33,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"Cconv (generic function with 1 method)"
]
},
"execution_count": 33,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function Cconv(a,b)\n",
" A = Cirmat(a)\n",
" B = Cirmat(b)\n",
" m = A*B\n",
" m[1,:]\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 34,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"5-element Array{Float64,1}:\n",
" 67.0\n",
" 94.0\n",
" 81.0\n",
" 78.0\n",
" 85.0"
]
},
"execution_count": 34,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"Cconv(vecc,vecd)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Eigenvalues and Eigenvectors of P"
]
},
{
"cell_type": "code",
"execution_count": 35,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"4-element Array{Complex{Float64},1}:\n",
" -1.0000000000000004 + 0.0im \n",
" 8.326672684688674e-17 + 0.9999999999999996im\n",
" 8.326672684688674e-17 - 0.9999999999999996im\n",
" 0.9999999999999999 + 0.0im "
]
},
"execution_count": 35,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"eigvals(Shiftm(4))"
]
},
{
"cell_type": "code",
"execution_count": 36,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"4×4 Array{Complex{Float64},2}:\n",
" 1.0-0.0im 1.0+0.0im 1.0-0.0im 1.0-0.0im\n",
" 1.0-0.0im -0.0+1.0im -1.0-0.0im -0.0-1.0im\n",
" 1.0-0.0im -1.0-0.0im 1.0-0.0im -1.0+0.0im\n",
" 1.0-0.0im -0.0-1.0im -1.0-0.0im -0.0+1.0im"
]
},
"execution_count": 36,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# The eigenvector matrix for P is the Fourier matrix\n",
"eP = eigvecs(Shiftm(4))\n",
"\n",
"# Julia give one eigenvector situation. Eigenvector multiply a scalar is still eigenvector .\n",
"eigenP = hcat(eP[:,4]*(-2),eP[:,2]*2im,eP[:,1]*(-2),eP[:,3]*(-2im))\n",
"round.(eigenP,digits = 5)"
]
},
{
"cell_type": "code",
"execution_count": 37,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"4×4 Array{Complex{Float64},2}:\n",
" 1.0+0.0im 1.0+0.0im 1.0+0.0im 1.0+0.0im\n",
" 1.0+0.0im 0.0+1.0im -1.0+0.0im -0.0-1.0im\n",
" 1.0+0.0im -1.0+0.0im 1.0-0.0im -1.0+0.0im\n",
" 1.0+0.0im -0.0-1.0im -1.0+0.0im 0.0+1.0im"
]
},
"execution_count": 37,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"round.(Fmat(4),digits = 5)"
]
},
{
"cell_type": "code",
"execution_count": 38,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"4×4 Array{Float64,2}:\n",
" 0.0 1.0 0.0 0.0\n",
" 0.0 0.0 1.0 0.0\n",
" 0.0 0.0 0.0 1.0\n",
" 1.0 0.0 0.0 0.0"
]
},
"execution_count": 38,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"Shiftm(4)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Eigenvalues and Eigenvectors of a Circulant C"
]
},
{
"cell_type": "code",
"execution_count": 39,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"4×4 Array{Complex{Float64},2}:\n",
" -0.5+0.0im -0.5+0.0im -0.5-0.0im -0.5+0.0im\n",
" -0.5+0.0im -0.0+0.5im -0.0-0.5im 0.5+0.0im\n",
" -0.5+0.0im 0.5+0.0im 0.5-0.0im -0.5+0.0im\n",
" -0.5+0.0im 0.0-0.5im 0.0+0.5im 0.5+0.0im"
]
},
"execution_count": 39,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Eigenvectors of a circulant matrix\n",
"cc = [1,2,3,4]\n",
"Cir = Cirmat(cc)\n",
"round.(eigvecs(Cir),digits = 5)"
]
},
{
"cell_type": "code",
"execution_count": 40,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"4×4 Array{Complex{Float64},2}:\n",
" -0.5+0.0im -0.5+0.0im -0.5-0.0im -0.5+0.0im\n",
" -0.5+0.0im 0.0+0.5im 0.0-0.5im 0.5+0.0im\n",
" -0.5+0.0im 0.5-0.0im 0.5+0.0im -0.5+0.0im\n",
" -0.5+0.0im 0.0-0.5im 0.0+0.5im 0.5+0.0im"
]
},
"execution_count": 40,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Eigenvector of the permutation\n",
"# Julia give one eigenvector situation. Eigenvector multiply a scalar is still eigenvector .\n",
"eigenP2 = hcat(eP[:,4],eP[:,3]*im,eP[:,2]*-im,eP[:,1])\n",
"round.(eigenP2,digits = 5)\n",
"# Eigenvectors of a circulant matrix are the same as the eigenvectors of the permutation"
]
},
{
"cell_type": "code",
"execution_count": 41,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"4-element Array{Complex{Float64},1}:\n",
" 10.0 + 0.0im\n",
" -2.0 - 2.0im\n",
" -2.0 + 0.0im\n",
" -2.0 + 2.0im"
]
},
"execution_count": 41,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Multiply F times the vector c in the top row of C to find the eigenvalues\n",
"round.(Fmat(4)*Cir[1,:],digits = 5)"
]
},
{
"cell_type": "code",
"execution_count": 42,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"4-element Array{Complex{Float64},1}:\n",
" 10.0 + 0.0im\n",
" -2.0 + 2.0im\n",
" -2.0 - 2.0im\n",
" -2.0 + 0.0im"
]
},
"execution_count": 42,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"round.(eigvals(Cir),digits = 5)"
]
},
{
"cell_type": "code",
"execution_count": 43,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"1×4 Array{Complex{Float64},2}:\n",
" 10.0-0.0im -2.0+2.0im -2.0-0.0im -2.0-2.0im"
]
},
"execution_count": 43,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# The N eogenvalues of C are the components of Fc = inverse Fourier transform of c\n",
"round.(Cir[1,:]'*DFTmat(4),digits = 5)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### The Convolution Rule"
]
},
{
"cell_type": "code",
"execution_count": 44,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"5-element Array{Float64,1}:\n",
" 9.0\n",
" 0.0\n",
" 8.0\n",
" 6.0\n",
" 4.0"
]
},
"execution_count": 44,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Two circulant matrices C and D. Their top rows are the vectors c and d\n",
"c = C[1,:]\n",
"d = D[1,:]"
]
},
{
"cell_type": "code",
"execution_count": 45,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"5×5 Array{Float64,2}:\n",
" 67.0 94.0 81.0 78.0 85.0\n",
" 85.0 67.0 94.0 81.0 78.0\n",
" 78.0 85.0 67.0 94.0 81.0\n",
" 81.0 78.0 85.0 67.0 94.0\n",
" 94.0 81.0 78.0 85.0 67.0"
]
},
"execution_count": 45,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Top row of CD is the cyclic convolution of vectors c and d\n",
"C*D"
]
},
{
"cell_type": "code",
"execution_count": 46,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"5-element Array{Float64,1}:\n",
" 67.0\n",
" 94.0\n",
" 81.0\n",
" 78.0\n",
" 85.0"
]
},
"execution_count": 46,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"Cconv(c,d)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"F(c * d) = (Fc).×(Fd)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"LHS:Convolve c and d first, then transform by F"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"RHS:Transform by F first, then multiply Fc times Fd component by component"
]
},
{
"cell_type": "code",
"execution_count": 48,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"5-element Array{Complex{Float64},1}:\n",
" 405.0 + 0.0im \n",
" -6.3196601125010545 + 10.3228644035338im \n",
" -28.68033988749893 + 2.43689772174681im \n",
" -28.68033988749893 - 2.4368977217467886im\n",
" -6.3196601125010154 - 10.322864403533785im"
]
},
"execution_count": 48,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Convolution Rule\n",
"Fmat(5)*Cconv(c,d)"
]
},
{
"cell_type": "code",
"execution_count": 49,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"5-element Array{Complex{Float64},1}:\n",
" 405.0 + 0.0im \n",
" -6.319660112501046 + 10.3228644035338im \n",
" -28.68033988749893 + 2.4368977217467958im\n",
" -28.680339887498935 - 2.4368977217467833im\n",
" -6.319660112501056 - 10.322864403533783im"
]
},
"execution_count": 49,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"(Fmat(5)*c).*(Fmat(5)*d)"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Julia 1.0.3",
"language": "julia",
"name": "julia-1.0"
},
"language_info": {
"file_extension": ".jl",
"mimetype": "application/julia",
"name": "julia",
"version": "1.0.3"
}
},
"nbformat": 4,
"nbformat_minor": 2
}