{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {},
   "outputs": [],
   "source": [
    "# eda06_00 clear all variables and import vatious modules\n",
    "%reset -f\n",
    "import os\n",
    "from datetime import date\n",
    "from math import exp, pi, sin, sqrt, floor, ceil\n",
    "import numpy as np\n",
    "import scipy.linalg as la\n",
    "from matplotlib import pyplot as plt\n",
    "import matplotlib.cm as cm\n",
    "from matplotlib.colors import ListedColormap\n",
    "\n",
    "# eda_draw function makes a \"pictorial matrix equation\"\n",
    "# arguments are vectors, matrices and strings\n",
    "# which are plotted in the order that the appear\n",
    "# except that strings starting with 'title ' are plotted\n",
    "# under the subseqeunt matrix or vector\n",
    "# always returns a status of 1\n",
    "\n",
    "def eda_draw(*argv):\n",
    "    bw = np.zeros((256,4));\n",
    "    v = 0.9*(256 - np.linspace( 0, 255, 256 ))/255;\n",
    "    bw[:,0] = v;\n",
    "    bw[:,1] = v;\n",
    "    bw[:,2] = v;\n",
    "    bw[:,3] = np.ones(256);\n",
    "    bwcmap = ListedColormap(bw);\n",
    "    # size of plot\n",
    "    W = 16;\n",
    "    H = 4;\n",
    "    fig1 = plt.figure(1);\n",
    "    # figsize width and height in inches\n",
    "    fig1.set_size_inches(W,H);\n",
    "    ax1 = plt.subplot(1,1,1);\n",
    "    plt.axis([0, W, -H/2, H/2]);\n",
    "    plt.axis('off');\n",
    "    LM = W/6;    # matrix width and heoght\n",
    "    LV = W/40;   # vector width\n",
    "    FS = 0.12;    # character width\n",
    "    TO = 0.4;    # title vertical offset\n",
    "    SP = 0.2;    # space between objects\n",
    "    LS = 0.2;    # leading space\n",
    "    p = LS; # starting x-position\n",
    "    istitle=0; # flags presence of a title\n",
    "    for a in argv:\n",
    "        if isinstance(a,np.ndarray):\n",
    "            sh = np.shape(a);\n",
    "            if len(sh) == 1:  # conversion to nx1 array\n",
    "                n = sh[0];\n",
    "                m = 1;\n",
    "                ap = a;\n",
    "                a = np.zeros((n,1));\n",
    "                a[:,0] = ap;\n",
    "            else:\n",
    "                n = sh[0];\n",
    "                m = sh[1];\n",
    "            if m==1:\n",
    "                pold=p;\n",
    "                left=p;\n",
    "                right=p+LV;\n",
    "                bottom=-LM/2;\n",
    "                top=LM/2;\n",
    "                plt.imshow( a, cmap=bwcmap, vmin=np.min(a), vmax=np.max(a), extent=(left,right,bottom,top) );\n",
    "                p = p+LV;\n",
    "                pm = (p+pold)/2;\n",
    "                if istitle:\n",
    "                    plt.text(pm,-(LM/2)-TO,titlestr,horizontalalignment='center');\n",
    "                    istitle=0;\n",
    "                p = p+SP;\n",
    "            else:\n",
    "                pold=p;\n",
    "                left=p;\n",
    "                right=p+LM;\n",
    "                bottom=-LM/2;\n",
    "                top=LM/2;\n",
    "                plt.imshow( a, cmap=bwcmap, vmin=np.min(a), vmax=np.max(a), extent=(left,right,bottom,top) );\n",
    "                p = p+LM;\n",
    "                pm = (p+pold)/2;\n",
    "                if istitle:\n",
    "                    plt.text(pm,-(LM/2)-TO,titlestr,horizontalalignment='center');\n",
    "                    istitle=0;\n",
    "                p = p+SP;\n",
    "        elif isinstance(a,str):\n",
    "            ns = len(a);\n",
    "            istitle=0;\n",
    "            if( ns>=6 ):\n",
    "                if 'title ' in a[0:6]:\n",
    "                    istitle=1;\n",
    "                    titlestr=a[6:];\n",
    "            if( istitle != 1):\n",
    "                plt.text(p,0,a);\n",
    "                p = p + ns*FS + SP;\n",
    "    plt.show();\n",
    "    return 1;\n",
    "\n",
    "# function to make a numpy N-by-1 column vector\n",
    "# c=eda_cvec1(v) from entity v that is array-like,\n",
    "# including a number, a list of numbers, a tuple\n",
    "# of numbers, an N-by-0 np array and a N-by1 np array.\n",
    "# It promotes integers to floats, and integers\n",
    "# and floats to complex, by context. The function\n",
    "# also insures that, if v is an np array, that\n",
    "# c is a copy, as contrasted to a view, of it.\n",
    "# This version takes just one input argment.\n",
    "# whereas c=eda_cvec(v1,v2,...) concatenates\n",
    "# many argiments.\n",
    "def eda_cvec1(v):\n",
    "    if isinstance(v, int):\n",
    "        w = np.zeros((1,1),dtype=int);\n",
    "        w[0,0] = v;\n",
    "        return w;\n",
    "    elif isinstance(v, float):\n",
    "        w = np.zeros((1,1),dtype=float);\n",
    "        w[0,0] = v;\n",
    "        return w;\n",
    "    elif isinstance(v, complex):\n",
    "        w = np.zeros((1,1),dtype=complex);\n",
    "        w[0,0] = v;\n",
    "        return w;\n",
    "    elif isinstance(v, np.ndarray):\n",
    "        s = np.shape(v);\n",
    "        if len(s) == 1:\n",
    "            return np.copy(np.reshape(v,(s[0],1)));\n",
    "        else:\n",
    "            [r,c]=s;\n",
    "            if( c==1 ):\n",
    "                return(np.copy(v));\n",
    "            elif(r==1):\n",
    "                return(np.copy(v.T));\n",
    "            else:\n",
    "                raise TypeError(\"eda_cvec: %d by %d ndarray not allowed\" % (r, c));\n",
    "    elif isinstance(v, list):\n",
    "        r = len(v);\n",
    "        t = int;\n",
    "        for vi in v:\n",
    "            if isinstance(vi,int):\n",
    "                pass;\n",
    "            elif isinstance(vi,float):\n",
    "                t=float;\n",
    "            elif isinstance(vi,complex):\n",
    "                t=complex;\n",
    "                break;\n",
    "            else:\n",
    "                raise TypeError(\"eda_cvec: list contains unsupported type %s\" % type(vi));\n",
    "        w = np.zeros((r,1),dtype=t);\n",
    "        w[:,0] = v;\n",
    "        return w;\n",
    "    elif isinstance(v, tuple):\n",
    "        r = len(v);\n",
    "        t = int;\n",
    "        for vi in v:\n",
    "            if isinstance(vi,int):\n",
    "                pass;\n",
    "            elif isinstance(vi,float):\n",
    "                t=float;\n",
    "            elif isinstance(vi,complex):\n",
    "                t=complex;\n",
    "                break;\n",
    "            else:\n",
    "                raise TypeError(\"eda_cvec: tuple contains unsupported type %s\" % type(vi));\n",
    "        w = np.zeros((r,1),dtype=t);\n",
    "        w[:,0] = v;\n",
    "        return w;\n",
    "    else:\n",
    "        raise TypeError(\"eda_cvec: %s not supported\" % type(v));\n",
    "        \n",
    "# function to make a numpy N-by-1 column vector\n",
    "# c=eda_cvec(v1, v2, ...) from a list of several\n",
    "# array-like entities v1, v2, including a number\n",
    "# a list of numbers, a tuple of numbers, an N-by-0 np array\n",
    "# and a N-by-1 np array. The function\n",
    "# also insures that, if v is an np array, that\n",
    "# c is a copy, as contrasted to a view, of it\n",
    "# It promotes integers to floats, and integers\n",
    "# and floats to complex, by context.\n",
    "# This version concatenates many argments,\n",
    "# whereas c=eda_cvec1(v1) takes just one argiment.\n",
    "# I recommend always using eda_cvec(v1, v2, ...)\n",
    "def eda_cvec(*argv):\n",
    "    t = int;\n",
    "    Nt = 0;\n",
    "    for a in argv:\n",
    "        v = eda_cvec1(a);\n",
    "        N,M = np.shape(v);\n",
    "        Nt = Nt + N;\n",
    "        if (t==int) and isinstance(v[0,0],float):\n",
    "            t=float;\n",
    "        elif isinstance(v[0,0],complex):\n",
    "            t=complex;\n",
    "    w = np.zeros((Nt,1),dtype=t);\n",
    "    Nt = 0;\n",
    "    for a in argv:\n",
    "        v = eda_cvec1(a);\n",
    "        N,M = np.shape(v);\n",
    "        w[Nt:Nt+N,0] = v[0:N,0];\n",
    "        Nt = Nt + N;\n",
    "    return w;\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 24788 characters]\n",
      "text/plain": [
       "<Figure size 576x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# eda06_01\n",
    "# cosine wave\n",
    "\n",
    "# make the auxillary variable t\n",
    "N=100;\n",
    "Dt=1;\n",
    "tmin=0.0;\n",
    "tmax = tmin+Dt*(N-1);\n",
    "t = eda_cvec( np.linspace(tmin,tmax,N) );\n",
    "\n",
    "# a cosine of ampmitude C, time offset t0 and period T\n",
    "t0 = 10.0;\n",
    "T = 50;\n",
    "C= 2.0;\n",
    "d = C*np.cos( 2*pi*(t-t0)/T );\n",
    "\n",
    "# plot data\n",
    "fig1 = plt.figure(1,figsize=(8,4));\n",
    "ax1 = plt.subplot(1,1,1);\n",
    "plt.axis([tmin, tmax, 1.1*np.min(d), 1.1*np.max(d)])\n",
    "plt.plot(t,d,'k-');\n",
    "plt.xlabel('time t');\n",
    "plt.ylabel('d(t)');\n",
    "plt.title('Delayed cosine wave');\n",
    "plt.show();\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 35384 characters]\n",
      "text/plain": [
       "<Figure size 576x288 with 6 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# eda06_02\n",
    "# plot columns of G\n",
    "# (not normalized - see eda06_13 for properly normalized version)\n",
    "\n",
    "# square data kernel\n",
    "N=32; # must be an even number\n",
    "M=N;\n",
    "\n",
    "# make the auxillary variable t\n",
    "N=100;\n",
    "Dt = 1;\n",
    "tmin=0.0;\n",
    "tmax = tmin+Dt*(N-1);\n",
    "Dt = (tmax - tmin)/(N-1);\n",
    "t = eda_cvec( np.linspace(tmin,tmax,N) );\n",
    "\n",
    "# standard definition of nyquist frequency fmax\n",
    "# frequency increment Df, frequency vector f,\n",
    "# angular frequency vector w\n",
    "Nf=floor(N/2+1);\n",
    "fmax = 1/(2*Dt);\n",
    "Df = fmax/(N/2);\n",
    "f = eda_cvec( np.linspace(0,fmax,Nf) );\n",
    "Nw=Nf;\n",
    "wmax = 2*pi*fmax;\n",
    "Dw = wmax/(N/2);\n",
    "w = eda_cvec( np.linspace(0,wmax,Nw) );\n",
    "\n",
    "# set up data kernel G\n",
    "G = np.zeros((N,M));\n",
    "\n",
    "# zero frequency column\n",
    "G[:,0] = np.ones((N,1)).ravel();\n",
    "\n",
    "# interior M/2-1 columns.  Note indice arithmetic\n",
    "# to keep track of pairs of columns, one for cos(),\n",
    "# one for sin()\n",
    "Mo2 = floor(M/2);\n",
    "for i in range(1,Mo2):\n",
    "    j = 2*i-1;\n",
    "    k = j+1;\n",
    "    G[:,j]=np.cos(w[i,0]*t).ravel();\n",
    "    G[:,k]=np.sin(w[i,0]*t).ravel();\n",
    "\n",
    "# nyquist column\n",
    "G[:,M-1] = np.cos(w[Nw-1,0]*t).ravel();\n",
    "    \n",
    "# plot selected columns of dat kernel G\n",
    "fig1 = plt.figure(1,figsize=(8,4));\n",
    "plt.rcParams['xtick.bottom'] = False;\n",
    "plt.rcParams['xtick.labelbottom'] = False;\n",
    "plt.rcParams['xtick.top'] = True;\n",
    "plt.rcParams['xtick.labeltop'] = True;\n",
    "ip = 1;\n",
    "for i in [0,1,2,3,4,M-1]:\n",
    "    ax = plt.subplot(1,6,ip);\n",
    "    ax.invert_yaxis();\n",
    "    plt.axis([-2, 2, tmin, tmax ]);\n",
    "    plt.plot(G[:,i],t,'k-');\n",
    "    plt.plot(G[:,i],t,'k.');\n",
    "    mystr = \"col %d\" % (i+1);\n",
    "    plt.title(mystr);\n",
    "    if (i==0):\n",
    "        plt.ylabel('time t');\n",
    "    ip=ip+1;\n",
    "plt.show();\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "low frequency 0.12566370614359174\n",
      "high frequency 6.408849013323178\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 58400 characters]\n",
      "text/plain": [
       "<Figure size 1152x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# eda06_03\n",
    "# example of aliasing\n",
    "\n",
    "# make the auxillary variable t\n",
    "N=100;\n",
    "Dt=1;\n",
    "tmin=0.0;\n",
    "tmax = tmin+Dt*(N-1);\n",
    "Dt1 = (tmax - tmin)/(N-1);\n",
    "t = eda_cvec( np.linspace(tmin,tmax,N) );\n",
    "\n",
    "# standard definition of nyquist frequency fmax\n",
    "# frequency increment Df, frequency vector f,\n",
    "# angular frequency vector w\n",
    "Nf=floor(N/2+1);\n",
    "fmax = 1/(2*Dt1);\n",
    "Df = fmax/(N/2);\n",
    "f = eda_cvec( np.linspace(0,fmax,Nf) );\n",
    "Nw=Nf;\n",
    "wmax = 2*pi*fmax;\n",
    "Dw = wmax/(N/2);\n",
    "w = eda_cvec( np.linspace(0,wmax,Nw) );\n",
    "\n",
    "# low frequency oscillation\n",
    "w1= 2*Dw;\n",
    "print('low frequency', w1);\n",
    "d1 = np.cos( w1*t );\n",
    "\n",
    "# high frequency oscillation\n",
    "w2= (2+N)*Dw;\n",
    "print('high frequency', w2);\n",
    "d2 = np.cos( w2*t );\n",
    "\n",
    "# set up higher sample rate time\n",
    "N=500;\n",
    "Dt2=0.1;\n",
    "tmin=0.0;\n",
    "tmax = tmin+Dt2*(N-1);\n",
    "t3 = eda_cvec( np.linspace(tmin,tmax,N) );\n",
    "d3 = np.cos( w2*t3 );\n",
    "\n",
    "# plot signals\n",
    "fig1 = plt.figure(1,figsize=(16,4));\n",
    "plt.rcParams['xtick.bottom'] = False;\n",
    "plt.rcParams['xtick.labelbottom'] = False;\n",
    "plt.rcParams['xtick.top'] = True;\n",
    "plt.rcParams['xtick.labeltop'] = True;\n",
    "\n",
    "ax1 = plt.subplot(2,1,1);\n",
    "plt.axis([tmin, tmax, -2, 2 ]);\n",
    "plt.plot(t,d1,'k-');\n",
    "plt.plot(t,d2,'ko');\n",
    "plt.xlabel('time t');\n",
    "plt.ylabel('d(t)');\n",
    "\n",
    "ax1 = plt.subplot(2,1,2);\n",
    "plt.axis([tmin, tmax, -2, 2 ]);\n",
    "plt.plot(t,d2,'k-');\n",
    "plt.plot(t,d2,'ko');\n",
    "plt.plot(t3,d3,'k:');\n",
    "plt.xlabel('time t');\n",
    "plt.ylabel('d(t)');\n",
    "plt.show();\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "low frequency -0.25132741228718347\n",
      "high frequency 6.031857894892403\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 45308 characters]\n",
      "text/plain": [
       "<Figure size 576x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# eda06_03a\n",
    "# another example of aliasing\n",
    "\n",
    "# make the auxillary variable t\n",
    "N=50;\n",
    "Dt=1;\n",
    "tmin=0.0;\n",
    "tmax = tmin+Dt*(N-1);\n",
    "Dt1 = (tmax - tmin)/(N-1);\n",
    "t = eda_cvec( np.linspace(tmin,tmax,N) );\n",
    "\n",
    "# nyquist frequencies\n",
    "Nf=floor(N/2+1);\n",
    "fmax = 1/(2*Dt1);\n",
    "Df = fmax/(N/2);\n",
    "f = eda_cvec( np.linspace(0,fmax,Nf) );\n",
    "Nw=Nf;\n",
    "wmax = 2*pi*fmax;\n",
    "Dw = wmax/(N/2);\n",
    "w = eda_cvec( np.linspace(0,wmax,Nw) );\n",
    "\n",
    "# low frequency oscillation\n",
    "w1= -2*Dw;\n",
    "print('low frequency', w1);\n",
    "d1 = np.sin( w1*t );\n",
    "\n",
    "# high frequency oscillation\n",
    "w2= (N-2)*Dw;\n",
    "print('high frequency', w2);\n",
    "d2 = np.sin( w2*t );\n",
    "\n",
    "# set up higher sample rate time\n",
    "N=500;\n",
    "Dt2=0.1;\n",
    "tmin=0.0;\n",
    "tmax = tmin+Dt2*(N-1);\n",
    "t3 = eda_cvec( np.linspace(tmin,tmax,N) );\n",
    "d3 =  np.cos( w2*t3 );\n",
    "\n",
    "# plot signals\n",
    "fig1 = plt.figure(1,figsize=(8,4));\n",
    "plt.rcParams['xtick.bottom'] = False;\n",
    "plt.rcParams['xtick.labelbottom'] = False;\n",
    "plt.rcParams['xtick.top'] = True;\n",
    "plt.rcParams['xtick.labeltop'] = True;\n",
    "\n",
    "ax1 = plt.subplot(2,1,1);\n",
    "plt.axis([tmin, tmax, -2, 2 ]);\n",
    "plt.plot(t,d1,'k-');\n",
    "plt.plot(t,d2,'ko');\n",
    "plt.xlabel('time t');\n",
    "plt.ylabel('d(t)');\n",
    "\n",
    "ax1 = plt.subplot(2,1,2);\n",
    "plt.axis([tmin, tmax, -2, 2 ]);\n",
    "plt.plot(t,d2,'k-');\n",
    "plt.plot(t,d2,'ko');\n",
    "plt.plot(t3,d3,'k:');\n",
    "plt.xlabel('time t');\n",
    "plt.ylabel('d(t)');\n",
    "plt.show();\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 43072 characters]\n",
      "text/plain": [
       "<Figure size 720x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 37492 characters]\n",
      "text/plain": [
       "<Figure size 720x288 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Error [[2.87461543e-15]]\n"
     ]
    }
   ],
   "source": [
    "# eda06_04\n",
    "# amplitude spectral density of the Neuse River hydrograph\n",
    "# using least-squares\n",
    "# (not normalized - see eda06_13 for properly normalized version)\n",
    "\n",
    "# load data from file\n",
    "D = np.genfromtxt('neuse.txt', delimiter='\\t')\n",
    "[Nraw, K]=D.shape;\n",
    "traw = eda_cvec( D[:,0] );  # traw for \"t raw\"\n",
    "Dt = traw[1,0]-traw[0,0];  # sampling interval\n",
    "draw = eda_cvec( D[:,1] ); # draw for \"d raw\"\n",
    "\n",
    "# round off to even number of points\n",
    "No2 = floor(Nraw/2);\n",
    "N = 2*No2;\n",
    "d = eda_cvec( draw[0:N,0] );\n",
    "t = eda_cvec( traw[0:N,0] );\n",
    "\n",
    "# plot data\n",
    "fig1 = plt.figure(1,figsize=(10,4));\n",
    "plt.rcParams['xtick.bottom'] = True;\n",
    "plt.rcParams['xtick.labelbottom'] = True;\n",
    "plt.rcParams['xtick.top'] = False;\n",
    "plt.rcParams['xtick.labeltop'] = False;\n",
    "\n",
    "ax1 = plt.subplot(1,1,1);\n",
    "plt.axis([np.min(t), np.max(t), 0, 1.1*np.max(d) ]);\n",
    "plt.plot(t,d,'k-');\n",
    "plt.xlabel('time (days)');\n",
    "plt.ylabel('discharge (cfs)');\n",
    "plt.title('Neuse River Hydrograph');\n",
    "plt.show();\n",
    "\n",
    "# nyquist frequencies\n",
    "Nf=floor(N/2+1);\n",
    "fmax = 1/(2*Dt);\n",
    "Df = fmax/(N/2);\n",
    "f = np.zeros((Nf,1));\n",
    "f[:,0] = np.linspace(0,fmax,Nf);\n",
    "Nw=Nf;\n",
    "wmax = 2*pi*fmax;\n",
    "Dw = wmax/(N/2);\n",
    "w = np.zeros((Nw,1));\n",
    "w[:,0] = np.linspace(0,wmax,Nw);\n",
    "fpos=np.zeros((Nf,1));\n",
    "fpos[:,0] = Df * np.linspace(0,No2,Nf);\n",
    "wpos=np.zeros((Nf,1));\n",
    "wpos[:,0] = Dw * np.linspace(0,No2,Nf);\n",
    "\n",
    "# set up data kernel G\n",
    "M=N;\n",
    "Mo2=No2;\n",
    "G = np.zeros((N,M));\n",
    "\n",
    "# zero frequency column\n",
    "G[:,0] = np.ones((N,1)).ravel();\n",
    "\n",
    "# interior M/2-1 columns. Note index arithmetic\n",
    "# to handle pairs of columns, one for cos() and\n",
    "# one for sin()\n",
    "for i in range(1,Mo2):\n",
    "    j = 2*i-1;\n",
    "    k = j+1;\n",
    "    G[:,j]=np.cos(w[i,0]*t).ravel();\n",
    "    G[:,k]=np.sin(w[i,0]*t).ravel();\n",
    "\n",
    "# nyquist column\n",
    "G[:,M-1] = np.cos(w[Nw-1,0]*t).ravel();\n",
    "\n",
    "# solve least squares problem using\n",
    "# analytic formula for inv(GT*G)\n",
    "GTd = np.matmul( G.T, d );\n",
    "# since inv(GT G) is known analtically\n",
    "# and is a diaginal matrix, use a shortcut\n",
    "gtgi = 2*np.ones((M,1))/N;\n",
    "gtgi[0,0]=1/N;\n",
    "gtgi[M-1,0]=1/N;\n",
    "mest = np.multiply( gtgi , GTd );\n",
    "\n",
    "# compute spectral density\n",
    "s=np.zeros((Nw,1));\n",
    "\n",
    "# zero frequency\n",
    "s[0,0]=np.sqrt(mest[0,0]**2);\n",
    "\n",
    "# interior columns, where cos() and sin() are paired\n",
    "for i in range(1,Mo2):\n",
    "    j = 2*i-1;\n",
    "    k = j+1;\n",
    "    s[i,0] = sqrt(mest[j,0]**2 + mest[k,0]**2);\n",
    "\n",
    "# Nyquist frequency\n",
    "s[Nw-1,0]= sqrt(mest[M-1,0]**2);\n",
    "\n",
    "# plot spectrum\n",
    "fig2 = plt.figure(2,figsize=(10,4));\n",
    "plt.rcParams['xtick.bottom'] = True;\n",
    "plt.rcParams['xtick.labelbottom'] = True;\n",
    "plt.rcParams['xtick.top'] = False;\n",
    "plt.rcParams['xtick.labeltop'] = False;\n",
    "\n",
    "ax1 = plt.subplot(1,3,1);\n",
    "plt.axis([0, fmax, 0, 1.1*np.max(s[1:Nf,0]) ]);\n",
    "plt.plot(fpos[1:Nf,0],s[1:Nf,0],'k-');\n",
    "plt.xlabel('frequency f (cycles per day)');\n",
    "plt.ylabel('amplitude spectrum');\n",
    "plt.title('linear in frequency');\n",
    "\n",
    "ax1 = plt.subplot(1,3,2);\n",
    "Tmax = 1000;\n",
    "plt.axis([0, Tmax, 0, 1.1*np.max(s[1:Nf,0]) ]);\n",
    "plt.plot(np.reciprocal(fpos[1:Nf,0]),s[1:Nf,0],'k-');\n",
    "plt.xlabel('periot T (days)');\n",
    "plt.ylabel('amplitude spectrum');\n",
    "plt.title('linear in period');\n",
    "\n",
    "ax1 = plt.subplot(1,3,3);\n",
    "plt.semilogx(fpos[1:Nf,0], s[1:Nf,0], 'k-', basex=10);\n",
    "plt.axis([0.001, 1, 0, 1.1*np.max(s[1:Nf,0])])\n",
    "plt.xlabel('frequency f (cycles per day)');\n",
    "plt.ylabel('amplitude spectrum');\n",
    "plt.title('logarithmic in frequency');\n",
    "\n",
    "plt.show();\n",
    "\n",
    "# check results\n",
    "dpre =np.matmul( G, mest );\n",
    "e = d - dpre; # individual errors\n",
    "E = np.matmul( e.T, e ); # total error\n",
    "print('Error', E);\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 43072 characters]\n",
      "text/plain": [
       "<Figure size 720x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 32792 characters]\n",
      "text/plain": [
       "<Figure size 720x288 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Error of d->mest->dnew: [[1.92672359e-20]]\n"
     ]
    }
   ],
   "source": [
    "# eda06_05\n",
    "# amplitude spectral density of the Neuse River hydrograph\n",
    "# using fft command\n",
    "# (unnormalized spectral density - see eda06_14 for proper normalization)\n",
    "\n",
    "# load data from file\n",
    "D = np.genfromtxt('neuse.txt', delimiter='\\t')\n",
    "[Nraw, K]=D.shape;\n",
    "traw = eda_cvec( D[:,0] ); # traw for \"t raw\"\n",
    "Dt = traw[1,0]-traw[0,0]; # sampling interval\n",
    "draw = eda_cvec( D[:,1] ); # draw fro \"d raw\"\n",
    "\n",
    "# round off to even number of points\n",
    "No2 = floor(Nraw/2);\n",
    "N=2*No2;\n",
    "d = eda_cvec( draw[0:N,0] ); # truncate data vector\n",
    "t = eda_cvec( traw[0:N,0] ); # truncate time vector\n",
    "\n",
    "# plot data\n",
    "fig1 = plt.figure(1,figsize=(10,4));\n",
    "plt.rcParams['xtick.bottom'] = True;\n",
    "plt.rcParams['xtick.labelbottom'] = True;\n",
    "plt.rcParams['xtick.top'] = False;\n",
    "plt.rcParams['xtick.labeltop'] = False;\n",
    "\n",
    "ax1 = plt.subplot(1,1,1);\n",
    "plt.axis([np.min(t), np.max(t), 0, 1.1*np.max(d) ]);\n",
    "plt.plot(t,d,'k-');\n",
    "plt.xlabel('time (days)');\n",
    "plt.ylabel('discharge (cfs)');\n",
    "plt.title('Neuse River Hydrograph');\n",
    "plt.show();\n",
    "\n",
    "No2 = floor(N/2); # N/2 rounded off to an integer\n",
    "M=N; # number of model parameters equals number of data\n",
    "\n",
    "# set up time axis\n",
    "tmin = 0;\n",
    "tmax = Dt*(N-1);\n",
    "t= eda_cvec( np.linspace(tmin,tmax,N) );\n",
    "\n",
    "# generic frequency setuo\n",
    "fmax=1/(2.0*Dt);  # nyquist frequency\n",
    "Df=fmax/No2;      # frequency sampling\n",
    "Nf=No2+1;         # number of non-negative frequencies\n",
    "# f is full freqeuncy vector, [0, Df, 2Df ... fmax, -fmax+Df, ... -Df]\n",
    "f= eda_cvec(Df*np.concatenate((np.linspace(0,No2,Nf),np.linspace(-No2+1,-1,No2-1)), axis=0 ));\n",
    "Dw=2*pi*Df;       # angular frequency sampling\n",
    "w=2*pi*f;         # full angular fe=requency vector\n",
    "Nw=Nf;            # number of non-negative angular freqeuncies\n",
    "fpos = eda_cvec(Df * np.linspace(0,No2,Nf)); # non-negative frequencies\n",
    "wpos=2*pi*fpos;    # non-negative angular frequencies\n",
    "\n",
    "#  compute fourier coefficients\n",
    "# Don't forget the axis=0 (fft over elements of column).\n",
    "# The default is axis=1 (fft over elements of rows)\n",
    "mest = eda_cvec( np.fft.fft(d, axis=0) );\n",
    "\n",
    "# delete negative frequencies\n",
    "mestpos = eda_cvec( mest[0:Nf,0] );\n",
    "\n",
    "# ampltude spectrum (note is un-normalized)\n",
    "s = np.abs(mestpos);\n",
    "\n",
    "# plot spectrum\n",
    "fig2 = plt.figure(2,figsize=(10,4));\n",
    "plt.rcParams['xtick.bottom'] = True;\n",
    "plt.rcParams['xtick.labelbottom'] = True;\n",
    "plt.rcParams['xtick.top'] = False;\n",
    "plt.rcParams['xtick.labeltop'] = False;\n",
    "\n",
    "# linear in frequency plot\n",
    "ax1 = plt.subplot(1,3,1);\n",
    "plt.axis([0, fmax, 0, 1.1*np.max(s[1:Nf,0]) ]);\n",
    "plt.plot(fpos[1:Nf,0],s[1:Nf,0],'k-');\n",
    "plt.xlabel('frequency f (cycles per day)');\n",
    "plt.ylabel('amplitude spectrum');\n",
    "plt.title('linear in frequency');\n",
    "\n",
    "# linear in perior plot\n",
    "ax1 = plt.subplot(1,3,2);\n",
    "Tmax = 1000;\n",
    "plt.axis([0, Tmax, 0, 1.1*np.max(s[1:Nf,0]) ]);\n",
    "plt.plot(np.reciprocal(fpos[1:Nf,0]),s[1:Nf,0],'k-');\n",
    "plt.xlabel('periot T (days)');\n",
    "plt.ylabel('amplitude spectrum');\n",
    "plt.title('linear in period');\n",
    "\n",
    "# semilog plot\n",
    "ax1 = plt.subplot(1,3,3);\n",
    "plt.semilogx(fpos[1:Nf,0], s[1:Nf,0], 'k-', basex=10);\n",
    "plt.axis([0.001, 1, 0, 1.1*np.max(s[1:Nf,0])])\n",
    "plt.xlabel('frequency f (cycles per day)');\n",
    "plt.ylabel('amplitude spectrum');\n",
    "plt.title('logarithmic in frequency');\n",
    "\n",
    "plt.show();\n",
    "\n",
    "#  compute Error of d->mest->dnew\n",
    "# Don't forget the axis=0.\n",
    "dnew = eda_cvec( np.real(np.fft.ifft(mest,axis=0)));\n",
    "e = d - dnew; # individual errors\n",
    "E = np.matmul( e.T, e );  # total error\n",
    "print('Error of d->mest->dnew:', E );"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 4520 characters]\n",
      "text/plain": [
       "<Figure size 1152x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# eda06_06\n",
    "# spectra of Normal curves of varying widths\n",
    "\n",
    "# time series of length N=256 with sampling interval of Dt=1\n",
    "\n",
    "# generic time set up\n",
    "N=256;   # must be an even number\n",
    "Dt=1.0;\n",
    "tmin = 0;\n",
    "tmax = Dt*(N-1);\n",
    "t = eda_cvec( np.linspace(tmin,tmax,N) );\n",
    "\n",
    "# generic frequency set up\n",
    "No2 = floor(N/2);  # N divided by two as an integer\n",
    "fmax=1/(2.0*Dt);   # nyquist frequency\n",
    "Df=fmax/No2;       # frequency increment\n",
    "Nf=No2+1;          # number of non-negative frequencies\n",
    "# full vector of frequencies\n",
    "f=eda_cvec(Df*np.concatenate( ( np.linspace(0,No2,Nf), np.linspace(-No2+1,-1,No2-1) ), axis=0 ));\n",
    "Dw=2*pi*Df;        # angular frequency sampling\n",
    "w=2*pi*f;          # full vector of angular fequencies\n",
    "Nw=Nf;             # number of angular frequencies\n",
    "fpos=eda_cvec(Df*np.linspace(0,No2,Nf)); # non-negative frequencies\n",
    "wpos=2*pi*fpos;    # non-negative angular frequencies\n",
    "\n",
    "# stamdard deviations, chosen to make an aestheic plot\n",
    "Nsd=6;\n",
    "sd=eda_cvec( [1.25, 2.5, 5, 10, 20, 40] );\n",
    "\n",
    "# matrix D with each column a time series, one column for each sd\n",
    "D = np.zeros((N,Nsd));\n",
    "# matrix Ds with each column an amplitude spectral density\n",
    "Ds = np.zeros((Nf,Nsd));\n",
    "\n",
    "# populate D and Ds\n",
    "for j in range(Nsd):\n",
    "    t1 = tmax/2;     # mean\n",
    "    sd2 = sd[j]**2;  # variance\n",
    "    d = np.exp( -np.power(t-t1,2) / (2*sd2) );  # normal function\n",
    "    D[:,j] = d.ravel(); \n",
    "    \n",
    "    #  compute fourier coefficients\n",
    "    dfft = eda_cvec( np.fft.fft(d, axis=0) );\n",
    "\n",
    "    # delete negative frequencies\n",
    "    dfftpos = eda_cvec( dfft[0:Nf,0] );\n",
    "\n",
    "    # ampltude spectrum          \n",
    "    Ds[:,j]=np.abs(dfftpos).ravel();\n",
    "    \n",
    "eda_draw(D[:,0],D[:,1],D[:,2],D[:,3],D[:,4],D[:,5],'    ',Ds[:,0],Ds[:,1],Ds[:,2],Ds[:,3],Ds[:,4],Ds[:,5]);\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 12472 characters]\n",
      "text/plain": [
       "<Figure size 720x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# eda06_07 spike function\n",
    "\n",
    "# generic time set up\n",
    "N=256;  # length of time sereis; must be even\n",
    "Dt=1.0; # samppling interval\n",
    "tmin = 0;\n",
    "tmax = Dt*(N-1);\n",
    "t = eda_cvec( np.linspace(tmin,tmax,N) ); # time vector\n",
    "\n",
    "# spike time series of length N=256 with sampling interval of Dt=1\n",
    "d = np.zeros((N,1));\n",
    "d[0,0]=1.0;\n",
    "\n",
    "# generic frequency setup\n",
    "No2 = floor(N/2);\n",
    "fmax=1/(2.0*Dt); # nyquist frequency\n",
    "Df=fmax/No2; # frequency sampling\n",
    "Nf=No2+1; # number of non-negative frequencies\n",
    "# full (pos and neg) frequency vector\n",
    "f = eda_cvec( Df * np.concatenate( ( np.linspace(0,No2,Nf), np.linspace(-No2+1,-1,No2-1) ), axis=0 ));\n",
    "Dw=2*pi*Df; # angular frequency sampling\n",
    "w=2*pi*f; # full  (pos and neg) angular frequencies\n",
    "Nw=Nf; # number rof non-negative angular frequencies\n",
    "fpos = eda_cvec(Df * np.linspace(0,No2,Nf) ); # non-negative frequencies\n",
    "wpos= 2*pi*fpos; # non-negathive angular frequencies\n",
    "\n",
    "# take fourier transform\n",
    "# Note: Don't forget axis=0 option!\n",
    "dt = eda_cvec( np.fft.fft(d, axis=0) );\n",
    "\n",
    "# plot \n",
    "fig1 = plt.figure(1,figsize=(10,4));\n",
    "plt.rcParams['xtick.bottom'] = True;\n",
    "plt.rcParams['xtick.labelbottom'] = True;\n",
    "plt.rcParams['xtick.top'] = False;\n",
    "plt.rcParams['xtick.labeltop'] = False;\n",
    "\n",
    "ax1 = plt.subplot(2,1,1);\n",
    "plt.axis( [-tmax/10, tmax, -0.5, 1.5 ]);\n",
    "plt.plot(t,d,'k-');\n",
    "plt.plot(t,d,'k.');\n",
    "plt.xlabel('time t');\n",
    "plt.ylabel('d(t)');\n",
    "plt.title('time series');\n",
    "\n",
    "ax1 = plt.subplot(2,1,2);\n",
    "plt.axis( [-fmax, fmax, 0, 2] );\n",
    "plt.plot(f,np.real(dt),'k-');  # the transform should be purely real, so take real part\n",
    "plt.plot(f,np.real(dt),'k.');\n",
    "plt.xlabel('frequency f');\n",
    "plt.ylabel('d-tilde(f)');\n",
    "plt.title('time series');\n",
    "\n",
    "plt.show();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 18648 characters]\n",
      "text/plain": [
       "<Figure size 720x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# eda06_08 area under curve\n",
    "\n",
    "# time series of length N=256 with sampling interval of Dt=1\n",
    "\n",
    "# generic timeset up\n",
    "N=256; # length of time series, must be even\n",
    "Dt=1.0; # sampling interval\n",
    "tmin = 0; # minimum time\n",
    "tmax = Dt*(N-1); # maximum time\n",
    "t = eda_cvec( np.linspace(tmin,tmax,N) ); # time vector\n",
    "\n",
    "# generic frequency up\n",
    "No2 = floor(N/2);  # integer version of N/2\n",
    "fmax=1/(2.0*Dt);   # nyquist frequency \n",
    "Df=fmax/No2;       # frequency smaping\n",
    "Nf=No2+1;          # number on non-negative frequencies\n",
    "f=eda_cvec(Df*np.concatenate(( np.linspace(0,No2,Nf),np.linspace(-No2+1,-1,No2-1)),axis=0));\n",
    "Dw=2*pi*Df;        # angular frequeny increment\n",
    "w=2*pi*f;          # full (pos and neg) ana\n",
    "Nw=Nf;             # number of non-negative angular frequencies\n",
    "fpos=eda_cvec(Df * np.linspace(0,No2,Nf)); # non-negative frequencies\n",
    "wpos=2*pi*fpos;    # non-negative angular frequencies\n",
    "\n",
    "# exemplary function\n",
    "w0 = 2*pi*fmax/20;  # frequency of oscilation\n",
    "d= np.multiply(np.cos(w0*t),np.exp(-0.25*w0*t)); # product of cos(0 and delining exp()\n",
    "\n",
    "# Fourier transform; note normalization factor of Dt\n",
    "# Also, note axis=0 option, which must the there for this to work\n",
    "dt = eda_cvec( Dt*np.fft.fft(d, axis=0) );\n",
    "\n",
    "# area by Reimann summation\n",
    "area1 = Dt*np.sum(d);\n",
    "# area by fft\n",
    "area2 = np.real(dt[0,0]);\n",
    "\n",
    "# plot \n",
    "fig1 = plt.figure(1,figsize=(10,4));\n",
    "plt.rcParams['xtick.bottom'] = True;\n",
    "plt.rcParams['xtick.labelbottom'] = True;\n",
    "plt.rcParams['xtick.top'] = False;\n",
    "plt.rcParams['xtick.labeltop'] = False;\n",
    "\n",
    "ax1 = plt.subplot(1,1,1);\n",
    "plt.axis( [0, tmax, -1, 1] );\n",
    "plt.plot(t,d,'k-');\n",
    "plt.xlabel('time t');\n",
    "plt.ylabel('d(t)');\n",
    "titlestr = \"sum area %f     fft area %f\" % (area1, area2);\n",
    "plt.title(titlestr);\n",
    "\n",
    "plt.show();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 32708 characters]\n",
      "text/plain": [
       "<Figure size 720x504 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# eda06_09\n",
    "# time shift using fft example\n",
    "\n",
    "# generic time= set up\n",
    "N=256;\n",
    "Dt=1.0;\n",
    "tmin = 0.0;\n",
    "tmax = Dt*(N-1);\n",
    "t=eda_cvec(np.linspace(tmin,tmax,N));\n",
    "\n",
    "# generic time/frequency set up\n",
    "No2 = floor(N/2);\n",
    "fmax=1/(2.0*Dt);\n",
    "Df=fmax/No2;\n",
    "Nf=No2+1;\n",
    "f=eda_cvec(Df * np.concatenate( ( np.linspace(0,No2,Nf), np.linspace(-No2+1,-1,No2-1) ), axis=0 ));\n",
    "Dw=2*pi*Df;\n",
    "w=2*pi*f;\n",
    "Nw=Nf;\n",
    "fpos=eda_cvec(Df * np.linspace(0,No2,Nf));\n",
    "wpos=2*pi*fpos;\n",
    "\n",
    "# exemplary function\n",
    "d = np.zeros((N,1));\n",
    "w0 = 2*pi*fmax/20;\n",
    "d[:,0] = np.multiply(np.cos(w0*t),np.exp(-0.25*w0*t)).ravel();\n",
    "\n",
    "# times shift\n",
    "t0 = t[15,0];\n",
    "\n",
    "# take Fourier transform and mulpiply by exp( -w t0 )\n",
    "dt = np.zeros((N,1), dtype=np.cdouble);\n",
    "dt[:,0] = (np.multiply(  np.exp((complex(0,-1)*w*t0)), np.fft.fft(d, axis=0))).ravel();\n",
    "\n",
    "# inverse Fourier transform\n",
    "# result should be purely real, so take real part\n",
    "# to get rid of zero imaginary part\n",
    "ds = eda_cvec( np.real((np.fft.ifft(dt,axis=0))));\n",
    "           \n",
    "# plot \n",
    "fig1 = plt.figure(1,figsize=(10,7));\n",
    "plt.rcParams['xtick.bottom'] = True;\n",
    "plt.rcParams['xtick.labelbottom'] = True;\n",
    "plt.rcParams['xtick.top'] = False;\n",
    "plt.rcParams['xtick.labeltop'] = False;\n",
    "\n",
    "ax1 = plt.subplot(2,1,1);\n",
    "plt.axis( [0, tmax, -1, 1] );\n",
    "plt.plot(t,d,'k-');\n",
    "plt.ylabel('d(t)');\n",
    "plt.title('time series before shift');\n",
    "           \n",
    "ax1 = plt.subplot(2,1,2);\n",
    "plt.axis( [0, tmax, -1, 1] );\n",
    "plt.plot(t,ds,'k-');\n",
    "plt.xlabel('time t');\n",
    "plt.ylabel('d(t)');\n",
    "titlestr = \"time series after rshift of %.2f\" % (t0);\n",
    "plt.title(titlestr);\n",
    "\n",
    "plt.show();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 44536 characters]\n",
      "text/plain": [
       "<Figure size 720x576 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# eda06_10\n",
    "# derivative using fft example\n",
    "\n",
    "# time series of length N=256 with sampling interval of Dt=1\n",
    "N=256;\n",
    "Dt=1.0;\n",
    "tmin = 0;\n",
    "tmax = Dt*(N-1);\n",
    "t=eda_cvec(np.linspace(tmin,tmax,N));\n",
    "    \n",
    "# generic time/frequency set up\n",
    "No2 = floor(N/2);\n",
    "fmax=1/(2.0*Dt);\n",
    "Df=fmax/No2;\n",
    "Nf=No2+1;\n",
    "f=eda_cvec(Df * np.concatenate((np.linspace(0,No2,Nf),np.linspace(-No2+1,-1,No2-1)),axis=0));\n",
    "Dw=2*pi*Df;\n",
    "w=2*pi*f;\n",
    "Nw=Nf;\n",
    "fpos=eda_cvec(Df * np.linspace(0,No2,Nf));\n",
    "wpos=2*pi*wpos;\n",
    "\n",
    "# exemplary function\n",
    "d = np.zeros((N,1));\n",
    "t0 = tmax/2;\n",
    "sd = 50;\n",
    "d[:,0] = np.exp(-np.power((t-t0),2)/(2*sd**2)).ravel();\n",
    "\n",
    "# take Fourier transform and multiply by iw (except for zero freuency)\n",
    "dt = np.zeros((N,1), dtype=np.cdouble);\n",
    "dt[:,0] = (np.multiply(  complex(0,1)*w , np.fft.fft(d, axis=0))).ravel();\n",
    "\n",
    "# inverse Fourier transform\n",
    "dddt = np.zeros((N,1));\n",
    "dddt[:,0] = np.real((np.fft.ifft(dt,axis=0))).ravel();\n",
    "\n",
    "# derivative by finite differences\n",
    "dddt2 = np.zeros((N,1));\n",
    "dddt2[0:N-1,0]=(d[1:N,0]-d[0:N-1,0])/Dt;\n",
    "           \n",
    "# plot \n",
    "fig1 = plt.figure(1,figsize=(10,8));\n",
    "plt.rcParams['xtick.bottom'] = True;\n",
    "plt.rcParams['xtick.labelbottom'] = True;\n",
    "plt.rcParams['xtick.top'] = False;\n",
    "plt.rcParams['xtick.labeltop'] = False;\n",
    "\n",
    "ax1 = plt.subplot(3,1,1);\n",
    "plt.axis( [0, tmax, -1.1*np.max(d), 1.1*np.max(d)] );\n",
    "plt.plot(t,d,'k-');\n",
    "plt.ylabel('d(t)');\n",
    "           \n",
    "ax1 = plt.subplot(3,1,2);\n",
    "plt.axis( [0, tmax, -1.1*np.max(dddt), 1.1*np.max(dddt)] );\n",
    "plt.plot(t,dddt,'k-');\n",
    "plt.xlabel('time t');\n",
    "plt.ylabel('dd/dt by fft');\n",
    "\n",
    "ax1 = plt.subplot(3,1,3);\n",
    "plt.axis( [0, tmax, -1.1*np.max(dddt2), 1.1*np.max(dddt2)] );\n",
    "plt.plot(t,dddt2,'k-');\n",
    "plt.xlabel('time t');\n",
    "plt.ylabel('dd/dt by sum');\n",
    "plt.show();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 39856 characters]\n",
      "text/plain": [
       "<Figure size 720x576 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# eda06_11\n",
    "# integral using fft example\n",
    "\n",
    "# generic time/y set up\n",
    "N=256;\n",
    "Dt=1.0;\n",
    "tmin = 0;\n",
    "tmax = Dt*(N-1);\n",
    "t=eda_cvec(np.linspace(tmin,tmax,N));\n",
    "\n",
    "# generic frequency set up\n",
    "No2 = floor(N/2);\n",
    "fmax=1/(2.0*Dt);\n",
    "Df=fmax/No2;\n",
    "Nf=No2+1;\n",
    "f=eda_cvec(Df * np.concatenate( ( np.linspace(0,No2,Nf), np.linspace(-No2+1,-1,No2-1) ), axis=0 ));\n",
    "Dw=2*pi*Df;\n",
    "w=2*pi*f;\n",
    "Nw=Nf;\n",
    "fpos=eda_cvec(Df * np.linspace(0,No2,Nf));\n",
    "wpos=2*pi*fpos;\n",
    "\n",
    "# exemplary function\n",
    "sd = 50;\n",
    "t1 = tmax/4;\n",
    "t2 = 3*tmax/4;\n",
    "sd = 25;\n",
    "d = (np.exp(-np.power((t-t1),2)/(2*sd**2)) - np.exp(-np.power((t-t2),2)/(2*sd**2))  );\n",
    "\n",
    "# take Fourier transform and divide iw (except set zero frequency to zero)\n",
    "# note that data type of fourier trasnform is complex-double\n",
    "dt = np.zeros((N,1), dtype=np.cdouble);\n",
    "wr = np.zeros((N,1), dtype=np.cdouble);\n",
    "wr[1:N,0]=np.reciprocal(complex(0,1)*w[1:N,0]);\n",
    "wr[0,0]=0;\n",
    "dt[:,0] = (np.multiply(  wr, np.fft.fft(d, axis=0))).ravel();  # don't forget that axis=0\n",
    "\n",
    "# inverse Fourier transform, don't forget the axis=0\n",
    "# since result should be purely real, delete imaginary part\n",
    "integral1 = eda_cvec( np.real((np.fft.ifft(dt,axis=0)))); \n",
    "\n",
    "# integeral by Reimann sum\n",
    "integral2 = eda_cvec((Dt*np.cumsum(d)));\n",
    "\n",
    "# adjust so that integrals have the same integration constant (zero)\n",
    "integral1 = integral1 - integral1[0,0];\n",
    "integral2 = integral2 - integral2[0,0];\n",
    "           \n",
    "# plot \n",
    "fig1 = plt.figure(1,figsize=(10,8));\n",
    "plt.rcParams['xtick.bottom'] = True;\n",
    "plt.rcParams['xtick.labelbottom'] = True;\n",
    "plt.rcParams['xtick.top'] = False;\n",
    "plt.rcParams['xtick.labeltop'] = False;\n",
    "\n",
    "ax1 = plt.subplot(3,1,1);\n",
    "plt.axis( [0, tmax, -1.1*np.max(d), 1.1*np.max(d)] );\n",
    "plt.plot(t,d,'k-');\n",
    "plt.ylabel('d(t)');\n",
    "           \n",
    "ax1 = plt.subplot(3,1,2);\n",
    "plt.axis( [0, tmax, -1.1*np.max(integral1), 1.1*np.max(integral1)] );\n",
    "plt.plot(t,integral1,'k-');\n",
    "plt.xlabel('time t');\n",
    "plt.ylabel('integral by fft');\n",
    "\n",
    "ax1 = plt.subplot(3,1,3);\n",
    "plt.axis( [0, tmax, -1.1*np.max(integral2), 1.1*np.max(integral2)] );\n",
    "plt.plot(t,integral2,'k-');\n",
    "plt.xlabel('time t');\n",
    "plt.ylabel('dd/dt by diff');\n",
    "plt.show();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 78076 characters]\n",
      "text/plain": [
       "<Figure size 720x576 with 8 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "\n",
      "\n",
      "\n",
      "\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 105648 characters]\n",
      "text/plain": [
       "<Figure size 720x576 with 8 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# eda06_12\n",
    "# makes exemplary time series d(t),\n",
    "# amplitude spectral density s(f) pairs for\n",
    "# Crib Sheet 6.2\n",
    "\n",
    "# time series of length N=256 with sampling interval of Dt=1\n",
    "N=128;\n",
    "Dt=0.5;\n",
    "tmin = 0;\n",
    "tmax = Dt*(N-1);\n",
    "t=eda_cvec(np.linspace(tmin,tmax,N));\n",
    "\n",
    "# generic frequency set up\n",
    "No2 = floor(N/2);\n",
    "fmax=1/(2.0*Dt);\n",
    "Df=fmax/No2;\n",
    "Nf=No2+1;\n",
    "f=eda_cvec(Df * np.concatenate( ( np.linspace(0,No2,Nf), np.linspace(-No2+1,-1,No2-1) ), axis=0 ));\n",
    "Dw=2*pi*Df;\n",
    "w=2*pi*f;\n",
    "Nw=Nf;\n",
    "fpos=eda_cvec(Df * np.linspace(0,No2,Nf));\n",
    "wpos=2*pi*fpos;\n",
    "\n",
    "# duration of time series\n",
    "T = tmax-tmin;\n",
    "\n",
    "# plot 1\n",
    "fig1 = plt.figure(1,figsize=(10,8));\n",
    "plt.rcParams['xtick.bottom'] = True;\n",
    "plt.rcParams['xtick.labelbottom'] = True;\n",
    "plt.rcParams['xtick.top'] = False;\n",
    "plt.rcParams['xtick.labeltop'] = False;\n",
    "p=1;\n",
    "\n",
    "# narrow function\n",
    "sigmad = 1.0;\n",
    "d = (np.exp(-np.power(t-tmax/2,2)/(2*(sigmad**2)))+np.random.normal(0,0.1,(N,1)));\n",
    "          \n",
    "dfft = eda_cvec( Dt*np.fft.fft(d, axis=0) );   # fourier transform, note normalization of Dt\n",
    "dfftpos = eda_cvec( dfft[0:Nf,0] ) ; # positive fequencies, only\n",
    "ds=sqrt(2/T)*np.abs(dfftpos); # properly normalized amplitude spectral density\n",
    "\n",
    "ax1 = plt.subplot(4,2,p);\n",
    "plt.axis( [0, T, -1.1, 1.1] );\n",
    "plt.plot(t,d,'k-');\n",
    "plt.ylabel('d(t)');\n",
    "p=p+1;\n",
    "           \n",
    "ax1 = plt.subplot(4,2,p);\n",
    "plt.axis( [0, fmax, 0, 1.1*np.max(ds)] );\n",
    "plt.plot(fpos,ds,'k-');\n",
    "plt.ylabel('s(f)');\n",
    "p=p+1;\n",
    "\n",
    "# zero area wavelet\n",
    "sigmad1 = 2.5;\n",
    "sigmad2 = 1.0;\n",
    "d = np.multiply( np.exp(-np.power(t-T/2,2)/(2*(sigmad1**2))) , np.cos((t-T/2)/sigmad2));\n",
    "d = d + np.random.normal(0,0.1,(N,1));\n",
    "          \n",
    "dfft = eda_cvec( Dt*np.fft.fft(d, axis=0) );\n",
    "dfftpos = eda_cvec( dfft[0:Nf,0] );\n",
    "ds=sqrt(2/T)*np.abs(dfftpos);\n",
    "\n",
    "ax1 = plt.subplot(4,2,p);\n",
    "plt.axis( [0, T, -1.1, 1.1] );\n",
    "plt.plot(t,d,'k-');\n",
    "plt.ylabel('d(t)');\n",
    "p=p+1;\n",
    "           \n",
    "ax1 = plt.subplot(4,2,p);\n",
    "plt.axis( [0, fmax, 0, 1.1*np.max(ds)] );\n",
    "plt.plot(fpos,ds,'k-');\n",
    "plt.ylabel('s(f)');\n",
    "p=p+1;\n",
    "\n",
    "# wide function\n",
    "sigmad = 7.0;\n",
    "d = np.exp(-np.power(t-tmax/2,2)/(2*(sigmad**2))) + np.random.normal(0,0.1,(N,1));\n",
    "\n",
    "# fourier transform\n",
    "dfft = eda_cvec( Dt*np.fft.fft(d, axis=0) );  # note Dt normalization\n",
    "dfftpos = eda_cvec( dfft[0:Nf,0] );  # keep non-negative frequencies\n",
    "ds = sqrt(2/T)*np.abs(dfftpos); # amplitude spectral density\n",
    "\n",
    "ax1 = plt.subplot(4,2,p);\n",
    "plt.axis( [0, T, -1.1, 1.1] );\n",
    "plt.plot(t,d,'k-');\n",
    "plt.ylabel('d(t)');\n",
    "p=p+1;\n",
    "           \n",
    "ax1 = plt.subplot(4,2,p);\n",
    "plt.axis( [0, fmax, 0, 1.1*np.max(ds)] );\n",
    "plt.plot(fpos,ds,'k-');\n",
    "plt.ylabel('s(f)');\n",
    "p=p+1;\n",
    "\n",
    "# triangle\n",
    "No4 = floor(N/4);\n",
    "d = np.zeros((N,1));\n",
    "d[No4:2*No4,0] = np.linspace(0,1,No4);\n",
    "d[2*No4:3*No4,0] = np.linspace(1,0,No4);\n",
    "d[:,0] = (d + np.random.normal(0,0.01,(N,1))).ravel();\n",
    "        \n",
    "dfft = eda_cvec( Dt*np.fft.fft(d, axis=0) ); # fourier transform with Dt normalization\n",
    "dfftpos = eda_cvec( dfft[0:Nf,0] );  # keep positive frequencies\n",
    "ds = sqrt(2/T)*np.abs(dfftpos);  # amplitude spectral density\n",
    "\n",
    "ax1 = plt.subplot(4,2,p);\n",
    "plt.axis( [0, T, -1.1, 1.1] );\n",
    "plt.plot(t,d,'k-');\n",
    "plt.xlabel('time t');\n",
    "plt.ylabel('d(t)');\n",
    "p=p+1;\n",
    "           \n",
    "ax1 = plt.subplot(4,2,p);\n",
    "plt.axis( [0, fmax, 0, 1.1*np.max(ds)] );\n",
    "plt.plot(fpos,ds,'k-');\n",
    "plt.xlabel('frequency f');\n",
    "plt.ylabel('s(f)');\n",
    "p=p+1;\n",
    "\n",
    "plt.show();\n",
    "\n",
    "print(\"\\n\\n\\n\\n\");\n",
    "\n",
    "# plot 2\n",
    "fig2 = plt.figure(2,figsize=(10,8));\n",
    "plt.rcParams['xtick.bottom'] = True;\n",
    "plt.rcParams['xtick.labelbottom'] = True;\n",
    "plt.rcParams['xtick.top'] = False;\n",
    "plt.rcParams['xtick.labeltop'] = False;\n",
    "p=1;\n",
    "\n",
    "# approximately sinusoidal, with time varying frequency w0\n",
    "w0 = eda_cvec( (0.25*2*pi*(0.95+0.1*np.linspace(0,N-1,N)/(N-1))) );\n",
    "d = np.sin(np.multiply(w0,t)) + np.random.normal(0,0.2,(N,1));\n",
    "\n",
    "dfft = eda_cvec( Dt*np.fft.fft(d, axis=0) );\n",
    "dfftpos = eda_cvec( dfft[0:Nf,0] );\n",
    "ds = sqrt(2/T)*np.abs(dfftpos);\n",
    "\n",
    "ax1 = plt.subplot(4,2,p);\n",
    "plt.axis( [0, T, -1.1, 1.1] );\n",
    "plt.plot(t,d,'k-');\n",
    "plt.ylabel('d(t)');\n",
    "p=p+1;\n",
    "           \n",
    "ax1 = plt.subplot(4,2,p);\n",
    "plt.axis( [0, fmax, 0, 1.1*np.max(ds)] );\n",
    "plt.plot(fpos,ds,'k-');\n",
    "plt.ylabel('s(f)');\n",
    "p=p+1;\n",
    "\n",
    "# beating, two sinusoids with frequencies in 2:1 ratio\n",
    "w1 = 0.25*2*pi;\n",
    "w2 = w1/2;\n",
    "d = np.sin(w1*t) + 0.5*np.sin(w2*t) + np.random.normal(0,0.2,(N,1));\n",
    "\n",
    "dfft = eda_cvec( Dt*np.fft.fft(d, axis=0) ); # fourier trasnform with Dt normailzation\n",
    "dfftpos = eda_cvec( dfft[0:Nf,0] ); # non-negative fequencies only\n",
    "ds =sqrt(2/T)*np.abs(dfftpos).ravel(); # amplitude spectral density\n",
    "\n",
    "ax1 = plt.subplot(4,2,p);\n",
    "plt.axis( [0, T, -1.1, 1.1] );\n",
    "plt.plot(t,d,'k-');\n",
    "plt.ylabel('d(t)');\n",
    "p=p+1;\n",
    "           \n",
    "ax1 = plt.subplot(4,2,p);\n",
    "plt.axis( [0, fmax, 0, 1.1*np.max(ds)] );\n",
    "plt.plot(fpos,ds,'k-');\n",
    "plt.ylabel('s(f)');\n",
    "p=p+1;\n",
    "\n",
    "# chirp, frequency that increases with time\n",
    "w0 = eda_cvec(0.25*2*pi*(2*np.linspace(0,N-1,N)/(N-1)));\n",
    "d = np.sin(np.multiply(w0,t)) + np.random.normal(0,0.05,(N,1));\n",
    "\n",
    "dfft = eda_cvec( Dt*np.fft.fft(d, axis=0) ); # fourier transform with Dt normaization\n",
    "dfftpos = eda_cvec( dfft[0:Nf,0] );  # non-negative frequencies only\n",
    "ds = sqrt(2/T)*np.abs(dfftpos);\n",
    "\n",
    "ax1 = plt.subplot(4,2,p);\n",
    "plt.axis( [0, T, -1.1, 1.1] );\n",
    "plt.plot(t,d,'k-');\n",
    "plt.ylabel('d(t)');\n",
    "p=p+1;\n",
    "           \n",
    "ax1 = plt.subplot(4,2,p);\n",
    "plt.axis( [0, fmax, 0, 1.1*np.max(ds)] );\n",
    "plt.plot(fpos,ds,'k-');\n",
    "plt.ylabel('s(f)');\n",
    "p=p+1;\n",
    "\n",
    "# sawtooth\n",
    "d = ((np.mod(t,N/10)/(N/10)-0.5)+np.random.normal(0.0,0.005,(N,1)));\n",
    "\n",
    "dfft = eda_cvec( Dt*np.fft.fft(d, axis=0) ); # fourier transform with Dt normalization\n",
    "dfftpos = eda_cvec( dfft[0:Nf,0] ); # non-negative frequencies\n",
    "ds = sqrt(2/T)*np.abs(dfftpos); # amplitude spectral ratio\n",
    "\n",
    "ax1 = plt.subplot(4,2,p);\n",
    "plt.axis( [0, T, -1.1, 1.1] );\n",
    "plt.plot(t,d,'k-');\n",
    "plt.xlabel('time t');\n",
    "plt.ylabel('d(t)');\n",
    "p=p+1;\n",
    "           \n",
    "ax1 = plt.subplot(4,2,p);\n",
    "plt.axis( [0, fmax, 0, 1.1*np.max(ds)] );\n",
    "plt.plot(fpos,ds,'k-');\n",
    "plt.xlabel('frequency f');\n",
    "plt.ylabel('s(f)');\n",
    "p=p+1;\n",
    "\n",
    "plt.show();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "check agreement between sines/cosine and fft at frequency 10\n",
      "A2 B2                      -0.5604 -0.3851\n",
      "(2/N)real(C) -(2/N)imag(C) -0.5604 -0.3851\n",
      "variance estimated from time series            105.9732\n",
      "variance estimated from sine/cosines           105.9732\n",
      "variance estimated from fft coefficients       105.9732\n",
      "variance estimated from fourier transform      105.9732\n",
      "variance estimated from frequency integral     105.9732\n",
      "variance estimated from power spectral density 105.9732\n"
     ]
    }
   ],
   "source": [
    "# eda06_13\n",
    "# comparison of variances\n",
    "# between time series and its power spectral density\n",
    "\n",
    "# random time series with zero mean and variance sigmad^2\n",
    "N=256;   # length of time series, must be even\n",
    "Dt=0.01; # sampling interval\n",
    "tmin = 0;\n",
    "tmax = Dt*(N-1);\n",
    "t=eda_cvec(Dt*np.linspace(0,N-1,N));\n",
    "sigmad1=10; # standard deviation\n",
    "sigmad1sq=sigmad1**2; # variance\n",
    "d = np.random.normal(0,sigmad1,(N,1)); # random timeseroes\n",
    "# this code will also work for a timeseries with non-zero mean,\n",
    "# but then the \"variance\" is really \"power\"\n",
    "d = d-np.mean(d);  # this statement can be commented out for \"power\"\n",
    "\n",
    "# variance estimated from time series\n",
    "sigmad2sq = np.matmul(d.T,d)/N;\n",
    "\n",
    "# standard frequency setup for cos() and sin() representation of data\n",
    "Nf=floor(N/2+1); # number o fnon-negative frequencies\n",
    "fmax = 1/(2*Dt); # nyquist frequency\n",
    "Df = fmax/(N/2); # frequency interval\n",
    "f = eda_cvec( np.linspace(0,fmax,Nf) ); # positive frequencies\n",
    "Nw=Nf; # number of non-negative angular frequencies\n",
    "wmax = 2*pi*fmax; # maximum angular frequency\n",
    "Dw = 2*pi*Df; # angular frequency interval\n",
    "w = 2*pi*f; # non-negative angular frequencies\n",
    "\n",
    "# set up data kernel G for fourier coefficients\n",
    "M=N;\n",
    "G = np.zeros((N,M));\n",
    "\n",
    "# zero frequency column\n",
    "G[:,0] = np.ones((N,1)).ravel();\n",
    "\n",
    "# interior M/2-1 columns\n",
    "# note index arithmetic due to columns occuring in pairs\n",
    "Mo2 = floor(M/2);\n",
    "for i in range(1,Mo2):\n",
    "    j = 2*i-1;\n",
    "    k = j+1;\n",
    "    G[:,j]=np.cos(w[i,0]*t).ravel();\n",
    "    G[:,k]=np.sin(w[i,0]*t).ravel();\n",
    "\n",
    "# nyquist column\n",
    "G[:,M-1] = np.cos(w[Nw-1,0]*t).ravel();\n",
    "\n",
    "# solve least squares problem using\n",
    "# analytic formula for inv(GT*G)\n",
    "GTd = np.matmul( G.T, d );\n",
    "# since inv( GT G ) is known analytically, use shortcut\n",
    "gtgi = 2*np.ones((M,1))/N;\n",
    "gtgi[0,0]=1/N;\n",
    "gtgi[M-1,0]=1/N;\n",
    "mest = np.multiply( gtgi , GTd );\n",
    "\n",
    "e = d - np.matmul(G,mest); # individual errors\n",
    "E = np.matmul( e.T , e ); # overall erros\n",
    "# print('error', E);  # uncomment to check that error is zero\n",
    "\n",
    "# generic frequency setup for fourier transform\n",
    "No2 = floor(N/2);  # integer version of N/2\n",
    "fmax=1/(2.0*Dt);   # nyquist frequency\n",
    "Df=fmax/No2;       # frequency sampling\n",
    "Nf=No2+1;          # number of non-negative frequencies\n",
    "# full (pos, neg) frequencies\n",
    "f=eda_cvec(Df * np.concatenate( ( np.linspace(0,No2,Nf), np.linspace(-No2+1,-1,No2-1) ), axis=0 ));\n",
    "Dw=2*pi*Df;     # angular freuqncy samppling\n",
    "w=2*pi*f;       # angular frequencies\n",
    "Nw=Nf;          # number of non-negative angular frequencies\n",
    "fpos=eda_cvec(Df * np.linspace(0,No2,Nf)); # non-negative frequencies\n",
    "wpos=2*pi*fpos; # non-negative angular frequencies\n",
    "\n",
    "# discrete fourier transform, nno normalization\n",
    "C = eda_cvec( np.fft.fft(d, axis=0) );\n",
    "\n",
    "# inverse Fourier transform; result should be purely real, so take real part\n",
    "drecon = eda_cvec( np.real((np.fft.ifft(C,axis=0))) );\n",
    "e = d - drecon; # individual errors\n",
    "E = np.matmul( e.T, e ); # total error\n",
    "# print('error', E); # uncomment to check that error is zero\n",
    "\n",
    "# comparison of coefficients\n",
    "jf=10;\n",
    "kf=2*jf-1;\n",
    "printstr = \"check agreement between sines/cosine and fft at frequency %d\" % (jf);\n",
    "print(printstr);\n",
    "printstr = \"A2 B2                      %.4f %.4f\" % (mest[kf,0], mest[kf+1,0]);\n",
    "print(printstr);\n",
    "printstr = \"(2/N)real(C) -(2/N)imag(C) %.4f %.4f\" % ((2/N)*np.real(C[jf,0]), -(2/N)*np.imag(C[jf,0]));\n",
    "print(printstr);\n",
    "\n",
    "printstr = \"variance estimated from time series            %.4f\" % (sigmad2sq);\n",
    "print(printstr);\n",
    "\n",
    "sigmad3sq = 0.5*np.matmul(mest.T,mest);\n",
    "# normally one doesn't bother to correct the first and last frequency, but ...\n",
    "sigmad3sq = sigmad3sq  + 0.5*(mest[0,0]**2 + mest[N-1,0]**2)\n",
    "printstr = \"variance estimated from sine/cosines           %.4f\" % (sigmad3sq);\n",
    "print(printstr);\n",
    "\n",
    "# variance estimated from fft coefficients\n",
    "C2 = np.zeros((Nf,1))\n",
    "C2[:,0] = np.real(np.multiply(C[0:Nf,0].conj(),C[0:Nf,0])).ravel();\n",
    "sigmad4sq = (2/(N**2)) * np.sum(C2);\n",
    "# normally one doesn't bother to correct the first and last frequency, but ...\n",
    "sigmad4sq = sigmad4sq - 0.5*(2/(N**2))*C2[0,0] - 0.5*(2/(N**2))*C2[Nf-1,0];\n",
    "printstr = \"variance estimated from fft coefficients       %.4f\" % (sigmad4sq);\n",
    "print(printstr);\n",
    "\n",
    "# approximation to integral transform\n",
    "dbar=Dt*C;\n",
    "\n",
    "# variance estimated from fourier transform\n",
    "dbar2 = eda_cvec( np.real(np.multiply(dbar[0:Nf,0].conj(),dbar[0:Nf,0])));\n",
    "sigmad4sq = (2/((Dt**2)*(N**2))) * np.sum( dbar2 );\n",
    "# normally one doesn't bother to correct the first and last frequency, but ...\n",
    "sigmad4sq = sigmad4sq - 0.5*(2/((Dt**2)*(N**2)))*dbar2[0,0] - 0.5*(2/((Dt**2)*(N**2)))*dbar2[Nf-1,0];\n",
    "printstr = \"variance estimated from fourier transform      %.4f\" % (sigmad4sq);\n",
    "print(printstr);\n",
    "\n",
    "# integral from 0 to fny of |dbar(f)|2\n",
    "theintegral = Df * np.sum( dbar2 );\n",
    "sigmad5sq = (2/(Df*(Dt**2)*(N**2))) * theintegral;\n",
    "# normally one doesn't bother to correct the first and last frequency, but ...\n",
    "sigmad5sq = sigmad5sq - 0.5*Df*(2/(Df*(Dt**2)*(N**2)))*dbar2[0,0] - 0.5*Df*(2/(Df*(Dt**2)*(N**2)))*dbar2[Nf-1,0];\n",
    "printstr = \"variance estimated from frequency integral     %.4f\" % (sigmad5sq);\n",
    "print(printstr);\n",
    "\n",
    "# definition of power spectral density\n",
    "T = tmax-tmin+Dt;  # note the correction of the length of the window by Dt\n",
    "thepsd = (2/T) * dbar2;\n",
    "\n",
    "# integral of power spectral density\n",
    "sigmad6sq = Df * np.sum(thepsd);\n",
    "# normally one doesn't bother to correct the first and last frequency, but ...\n",
    "sigmad6sq  = sigmad6sq  - 0.5*Df*thepsd[0,0] - 0.5*Df*thepsd[Nf-1,0];\n",
    "printstr = \"variance estimated from power spectral density %.4f\" % (sigmad6sq);\n",
    "print(printstr);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 43124 characters]\n",
      "text/plain": [
       "<Figure size 720x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 26844 characters]\n",
      "text/plain": [
       "<Figure size 720x288 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "variance estimated from time series            8.396180e+06\n",
      "variance estimated from power spectral density 8.396180e+06\n"
     ]
    }
   ],
   "source": [
    "# eda06_14\n",
    "# properly normalized amplitude spectral density of the Neuse River hydrograph\n",
    "# using fft command\n",
    "# (unnormalized spectral density - see eda06_13 for proper normalization)\n",
    "\n",
    "# load data from file\n",
    "D = np.genfromtxt('neuse.txt', delimiter='\\t')\n",
    "[Nraw, K]=D.shape;\n",
    "traw = eda_cvec( D[:,0] ); # time\n",
    "Dt = traw[1,0]-traw[0,0];  # sampling interval\n",
    "draw = eda_cvec( D[:,1] ); # data\n",
    "\n",
    "# round off to even number of points\n",
    "No2 = floor(Nraw/2);\n",
    "N=2*No2;\n",
    "d = eda_cvec( draw[0:N,0] );\n",
    "t = eda_cvec( traw[0:N,0] );\n",
    "\n",
    "# generic time set up\n",
    "tmin = 0;\n",
    "tmax = Dt*(N-1);\n",
    "t=eda_cvec(np.linspace(tmin,tmax,N));\n",
    "\n",
    "# generic frequency set up\n",
    "No2 = floor(N/2);  # integer version of N/2\n",
    "M=N;               # number of model parameters equal number of data\n",
    "fmax=1/(2.0*Dt);   # nyquist frequency\n",
    "Df=fmax/No2;       # frequency sampling\n",
    "Nf=No2+1;          # number of non-negative frequenices\n",
    "# full (pos, neg) frequency vector\n",
    "f=eda_cvec(Df * np.concatenate( ( np.linspace(0,No2,Nf), np.linspace(-No2+1,-1,No2-1) ), axis=0 ));\n",
    "Dw=2*pi*Df;        # angular frequency sampling\n",
    "w=2*pi*f;          # full (pos, neg) anguar frequency vector\n",
    "Nw=Nf;             # number of non-negative angular frequencies\n",
    "fpos=eda_cvec(Df * np.linspace(0,No2,Nf)); # non-negative frequencies\n",
    "wpos=2*pi*fpos;    # non-negative angular frequencies\n",
    "\n",
    "#  compute fourier coefficients. Note normalization of Dt\n",
    "dbar = eda_cvec(Dt*np.fft.fft(d, axis=0));\n",
    "dbarpos = eda_cvec(dbar[0:Nf,0]); # non-negative frequencies\n",
    "\n",
    "# compute power spectral density. note the normalization of (2/T).\n",
    "T = N*Dt;  # time series length\n",
    "s2 = (2/T) * np.real(np.multiply(dbarpos.conj(), dbarpos)); # properly normalized power spectral density\n",
    "\n",
    "# plot data\n",
    "fig1 = plt.figure(1,figsize=(10,4));\n",
    "plt.rcParams['xtick.bottom'] = True;\n",
    "plt.rcParams['xtick.labelbottom'] = True;\n",
    "plt.rcParams['xtick.top'] = False;\n",
    "plt.rcParams['xtick.labeltop'] = False;\n",
    "\n",
    "ax1 = plt.subplot(1,1,1);\n",
    "plt.axis([np.min(t), np.max(t), 0, 1.1*np.max(d) ]);\n",
    "plt.plot(t,d,'k-');\n",
    "plt.xlabel('time (days)');\n",
    "plt.ylabel('discharge (cfs)');\n",
    "plt.title('Neuse River Hydrograph');\n",
    "plt.show();\n",
    "\n",
    "# plot power spectral density\n",
    "fig2 = plt.figure(2,figsize=(10,4));\n",
    "plt.rcParams['xtick.bottom'] = True;\n",
    "plt.rcParams['xtick.labelbottom'] = True;\n",
    "plt.rcParams['xtick.top'] = False;\n",
    "plt.rcParams['xtick.labeltop'] = False;\n",
    "\n",
    "ax1 = plt.subplot(1,3,1);\n",
    "plt.axis([0, fmax, 0, 1.1*np.max(s2[1:Nf,0]) ]);\n",
    "plt.plot(fpos[1:Nf,0],s2[1:Nf,0],'k-');\n",
    "plt.xlabel('frequency f (cycles per day)');\n",
    "plt.ylabel('power spectral density');\n",
    "plt.title('linear in frequency');\n",
    "\n",
    "ax1 = plt.subplot(1,3,2);\n",
    "Tmax = 1000;\n",
    "plt.axis([0, Tmax, 0, 1.1*np.max(s2[1:Nf,0]) ]);\n",
    "plt.plot(np.reciprocal(fpos[1:Nf,0]),s2[1:Nf,0],'k-');\n",
    "plt.xlabel('periot T (days)');\n",
    "plt.ylabel('power spectral density');\n",
    "plt.title('linear in period');\n",
    "\n",
    "ax1 = plt.subplot(1,3,3);\n",
    "plt.axis([0.001, 1, 0, 1.1*np.max(s2[1:Nf,0])]);\n",
    "plt.semilogx(fpos[1:Nf,0], s2[1:Nf,0], 'k-', basex=10);\n",
    "plt.xlabel('frequency f (cycles per day)');\n",
    "plt.ylabel('power spectral density');\n",
    "plt.title('logarithmic in frequency');\n",
    "\n",
    "plt.show();\n",
    "\n",
    "# variance calculation\n",
    "v1 = np.std(d)**2;\n",
    "printstr = \"variance estimated from time series            %.6e\" % (v1);\n",
    "print(printstr);\n",
    "v2 = Df * np.sum(s2[1:Nf,0]);\n",
    "printstr = \"variance estimated from power spectral density %.6e\" % (v2);\n",
    "print(printstr);\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 69,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 30440 characters]\n",
      "text/plain": [
       "<Figure size 720x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# eda06_15\n",
    "# make a random noise dataset\n",
    "\n",
    "# Note that he noise is saved to noise2.txt, not noise.txt\n",
    "# to avoid overwriting the standard distribution.  So change\n",
    "# the name mif you want to overwrite  noise.txt\n",
    "\n",
    "# generic time setup\n",
    "N=1024;          # length of time series, must be even\n",
    "Dt=1.0;          # sampling interval\n",
    "tmin = 0.0;      # start time\n",
    "tmax = Dt*(N-1); # end time\n",
    "t = eda_cvec(np.linspace(tmin,tmax,N)); # time vector\n",
    "\n",
    "# make vector d of random noise\n",
    "dbar=0.0;   # mean\n",
    "sigmad=1.0; # standard deviation\n",
    "d = np.random.normal( dbar, sigmad, (N,1) );\n",
    "\n",
    "# plot data\n",
    "fig1 = plt.figure(1,figsize=(10,4));\n",
    "plt.rcParams['xtick.bottom'] = True;\n",
    "plt.rcParams['xtick.labelbottom'] = True;\n",
    "plt.rcParams['xtick.top'] = False;\n",
    "plt.rcParams['xtick.labeltop'] = False;\n",
    "\n",
    "ax1 = plt.subplot(1,1,1);\n",
    "r = dbar+3*sigmad;  # this range encloses 99% of data\n",
    "plt.axis([tmin, tmax, -r, r ]);\n",
    "plt.plot(t,d,'k.');\n",
    "plt.plot([tmin,tmax],[0,0],'k-');\n",
    "plt.xlabel('time');\n",
    "plt.ylabel('d(t)');\n",
    "titlestr = \"Normally-distributed random timeseries witm mean %.2f and sigma %.2f\" % (dbar, sigmad);\n",
    "plt.title(titlestr);\n",
    "plt.show();\n",
    "\n",
    "# save matrix to file of tab-separated values\n",
    "# the file name has been changed from \"noise.txt\" to \"noise2.txt\"\n",
    "# so as not to overwrite the existing file in the distributuon\n",
    "Dm = np.concatenate( (t,d),axis=1);\n",
    "np.savetxt(\"noise2.txt\",Dm, delimiter=\"\\t\");\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
