\n", "The wavelet transform of a time signal $y(t)$,\n", " $$\n", "Y(s, \\tau) = \\int_{-\\infty}^{+\\infty} dt\\,\n", "\\psi_{s,\\tau}^*(t) y(t),\n", " $$\n", "uses a basis $\\psi_{s,\\tau}(t)$ that is localized in time and contains a small range of frequencies.\n", "We start with the mother wavelet\n", "\\begin{align}\n", "\\Psi (t) = \\sin(8t) e^{-t^2/2},\n", "\\end{align}\n", "and form wavelet basis functions \n", "$$\n", "\\psi_{s,\\tau}(t) = \\frac{1}{\\sqrt{s}}\n", "\\Psi\\left(\\frac{t-\\tau}{s}\\right) = \\frac{1}{\\sqrt{s}} \\sin\\left[\n", "\\frac{8(t-\\tau)}{s}\\right]\\, e^{-{(t-\\tau)^2} /{2s^2}}.\n", " $$\n", " The transform is applied to:\n", " $$\n", "y(t) = \\begin{cases}\n", " \\sin 2\\pi t, & \\hbox{for } \\ 0 \\leq t \\leq 2,\\\\[3pt]\n", " 5\\sin 2\\pi t + 10\\sin 4\\pi t , & \\hbox{for } \\ 2 \\leq t \\leq 8,\\\\[3pt]\n", " 2.5\\sin 2\\pi t + 6\\sin 4\\pi t + 10\\sin 6\\pi t, &\\hbox{for }\\ 8 \\leq t \\leq 12.\n", "\\end{cases}\n", " $$\n", " \n", " Based on a program written by Written by Zlatko Dimcovic." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true, "deletable": true, "editable": true, "scrolled": true }, "outputs": [], "source": [ "\"\"\" From \"COMPUTATIONAL PHYSICS\" & \"COMPUTER PROBLEMS in PHYSICS\"\n", " by RH Landau, MJ Paez, and CC Bordeianu (deceased).\n", " Copyright R Landau, Oregon State Unv, MJ Paez, Univ Antioquia, \n", " C Bordeianu, Univ Bucharest, 2021. \n", " Please respect copyright & acknowledge our work.\"\"\"\n", "\n", "# ContWavelet.ipynb Continuous Wavelet Transform. Written by Zlatko Dimcovic\n", "\n", "from vpython import *\n", "import numpy as np\n", "\n", "N = 240\n", "transfgr = canvas(x=0,y=0,width=600,height=200, title='Transform, not normalized')\n", "transf = curve(x=list(range(0,90)),display=transfgr,color=color.cyan)\n", "wavlgr = canvas(x=0,y=200,width=600,height=200,\n", " title='Morlet Wavelet at different scales, up to s=12.0')\n", "wavelet = curve(x=list(range(0,N)),display=wavlgr,color=color.yellow)\n", "invtrgr = canvas(x=0,y=400,width=600,height=200,\n", " title='Inverse TF, not normalized')\n", "invtr = curve(x=list(range(0,N)),display=invtrgr,color=color.green)\n", "wvlabel = label(pos=vec(0, -50,0), text='s=', box=0,display=wavlgr)\n", "iT = 0.0\n", "fT = 12.0\n", "W = fT - iT # i,f times\n", "h = W/N; # Steps\n", "noPtsSig = N\n", "noS = 30 \n", "noTau = 90 # of pts\n", "iTau = 0.\n", "iS = 0.1\n", "tau = iTau\n", "s = iS\n", "\n", "''' Need *very* small s steps for high-frequency, but only if s is small\n", "Thus increment s by multiplying by number close enough to 1 '''\n", "\n", "dTau = W/noTau\n", "dS = (W/iS)**(1./noS)\n", "sig = np.zeros((noPtsSig),float) # Signal\n", "Y = np.zeros((noS,noTau),float) # Transform\n", "maxY = 0.001;\n", "\n", "def signal(noPtsSig, y): # The signal\n", " t = 0.0\n", " hs = W / noPtsSig\n", " t1 = W/6.\n", " t2 = 4.*W/6.\n", " for i in range(0,noPtsSig): \n", " if t >= iT and t <= t1: y[i] = sin(2*pi*t)\n", " elif t >= t1 and t <= t2:\ty[i] = 5.*sin(2*pi*t) +10.*sin(4*pi*t);\n", " elif t >= t2 and t <= fT: y[i] = 2.5*sin(2*pi*t)+6.*sin(4*pi*t)+10.*sin(6*pi*t)\n", " else:\n", " print(\"In signal(...) : t out of range.\")\n", " sys.exit(1)\n", " t += hs \n", "signal(noPtsSig, sig)\n", "\n", "def morlet(t, s, tau): # Mother\n", " T = (t-tau)/s\n", " return sin(8*T) * exp(-T*T/2.)\n", " \n", "def transform(s, tau, sig, j): \n", " integral = 0.\n", " t = 0.0 # \"initial time\" = class variable\n", " if j%2==0: wvlabel.text = 's=%5.3f' %s\n", " wavelet.clear()\n", " for i in range(0,len(sig)):\n", " t +=h\n", " rate(1500)\n", " yy=morlet(t,s,tau)\n", " wavelet.append(pos=vec(110.0/3*t-240,40.0*yy,0))\n", " integral += sig[i]*yy*h\n", " output=integral/sqrt(s) \n", " wavelet.visible=True \n", " return output\n", "\n", "def invTransform(t,Y): # given the transform (from previous steps)\n", " s = iS # computes original signal\n", " tau = iTau\n", " recSig_t = 0 \n", " for i in range (0,noS):\n", " s *= dS \n", " tau = iTau\n", " for j in range (0,noTau):\n", " tau += dTau\n", " recSig_t += dTau*dS *(s**(-1.5))* Y[i,j] * morlet(t,s,tau)\n", " return recSig_t\n", " \n", "print(\"working, finding transform\")\n", "for i in range( 0, noS):\n", " s *= dS # Scaling s\n", " tau = 0.0\n", " transf.clear()\n", " for j in range(0,noTau):\n", " rate(1500)\n", " tau += dTau # Translation\n", " Y[i,j] = transform(s, tau, sig,i)\n", " transf.append(pos=vec(40./3.0*tau-80, 4.0*Y[i,j] ,0))\n", " transf.visible = True \n", "print(\"transform found\")\n", "print(\"finding inverse transform\")\t\t # Inverse TF\n", "recSigData = \"recSig.dat\" \n", "recSig = np.zeros(len(sig)) # Same resolution\n", "t = 0.0;\n", "kco = 0\n", "j = 0\n", "\n", "for rs in range(0, len(recSig)): \t\t\t # with inverse transform\n", " recSig[rs] = invTransform(t, Y) \t\t\t # find the original signal\n", " t += h \t\t # not normalized\n", " invtr.append(pos =vec( rs*2.0-220, 1.5*recSig[rs],0))\n", "print(\"nDone\") " ] } ], "metadata": { "kernelspec": { "display_name": "Python 2", "language": "python", "name": "python2" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 2 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython2", "version": "2.7.13" } }, "nbformat": 4, "nbformat_minor": 2 }