{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# eda09_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, cos, tan, sqrt, floor, ceil, log\n",
    "import numpy as np\n",
    "import scipy.sparse.linalg as las\n",
    "from scipy import sparse\n",
    "import scipy.linalg as la\n",
    "import scipy.signal as sg\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",
    "\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",
    "\n",
    "\n",
    "def FTFmul(v):\n",
    "    # this function is used by the bicongugate gradient solver to solve the geneneralized least\n",
    "    # squares problem Fm=f.  Note that \"F\" must be called \"edaFsparse\".\n",
    "    global edaFsparse;\n",
    "    s = np.shape(v);\n",
    "    if(len(s)==1):\n",
    "        vv = np.zeros((s[0],1));\n",
    "        vv[:,0] = v;\n",
    "    else:\n",
    "        vv=v;\n",
    "    temp = edaFsparse*vv;\n",
    "    return edaFsparse.transpose()*temp;\n",
    "\n",
    "\n",
    "def GLSFilterMul(v):\n",
    "    # this function is used by the bicongugate gradient solver to solve the\n",
    "    # geneneralized least squares problem Fm=f with F=[G;H] and f=[d,h],\n",
    "    # where G is a topplitz matrix built from the filter g\n",
    "    # Note that \"H\" must be sparese and called \"edaHsparse\" and that\n",
    "    # the filter \"g\" bust be called edafilterg and must be a column vector\n",
    "    global edaHsparse;\n",
    "    global edafilterg;\n",
    "    # rearrange as column-vector\n",
    "    s = np.shape(v);\n",
    "    M = s[0];\n",
    "    if(len(s)==1):\n",
    "        vv = np.zeros((s[0],1));\n",
    "        vv[:,0] = v;\n",
    "    else:\n",
    "        vv=v;\n",
    "    N, i = np.shape(edafilterg);\n",
    "    # implements (GT G + HT H) v\n",
    "    temp1 = np.zeros((N+M-1,1));\n",
    "    temp1[:,0] = np.convolve(edafilterg.ravel(),vv.ravel()); # G v is of length N\n",
    "    a = np.zeros((N,1));\n",
    "    a[:,0] = temp1[0:N,0];\n",
    "    b = edaHsparse * v; # H v is of length K\n",
    "    temp2 = np.zeros((N+N-1,1));\n",
    "    temp2[:,0] = np.convolve((np.flipud(edafilterg)).ravel(),a.ravel()); # GT (G v) is of length M\n",
    "    a2 = temp2[N-1:N+M-1,0];\n",
    "    b2 = edaHsparse.transpose()*b; # HT (H v) is of length M\n",
    "    # FT F v = GT G v + HT H v\n",
    "    return (a2+b2);\n",
    "\n",
    "\n",
    "def GLSautoMul(v):\n",
    "    # this function is used by the bicongugate gradient solver to solve the\n",
    "    # geneneralized least squares problem Fm=f with F=[G;H] and f=[d,h],\n",
    "    # where GTG is built from the autocorrelation of the filter g\n",
    "    global edaHsparse;\n",
    "    global edaautoc;\n",
    "    # rearrange as column-vector\n",
    "    s = np.shape(v);\n",
    "    M = s[0];\n",
    "    if(len(s)==1):\n",
    "        vv = np.zeros((s[0],1));\n",
    "        vv[:,0] = v;\n",
    "    else:\n",
    "        vv=v;\n",
    "    # FT F v = GT G v + HT H v\n",
    "    GTGv=np.zeros((M,1));\n",
    "    for i in range(M):\n",
    "        u = np.concatenate( (np.flipud(edaautoc[0:i+1,0]), edaautoc[1:M-i,0]), axis=0 );\n",
    "        GTGv[i,0] =  np.matmul( u.T, vv );\n",
    "    Hv = edaHsparse*vv;\n",
    "    HTHv = edaHsparse.transpose()*Hv;\n",
    "    return (GTGv + HTHv);\n",
    "\n",
    "\n",
    "def chebyshevfilt(d, Dt, flow, fhigh):\n",
    "    #chebyshevfilt\n",
    "    # chebyshev IIR bandpass filter\n",
    "    # d - input array of data\n",
    "    # Dt - sampling interval\n",
    "    # flow - low pass frequency, Hz\n",
    "    # fhigh - high pass frequency, Hz\n",
    "    # dout - output array of data\n",
    "    # u - the numerator filter\n",
    "    # v - the denominator filter\n",
    "    # these filters can be used again using dout=filter(u,v,din);\n",
    "\n",
    "    # make sure input timeseries is a column vector\n",
    "    s = np.shape(d);\n",
    "    N = s[0];\n",
    "    if(N==1):\n",
    "        dd = np.zeros((N,1));\n",
    "        dd[:,0] = d;\n",
    "    else:\n",
    "        dd=d;\n",
    "        \n",
    "    # sampling rate\n",
    "    rate=1/Dt;\n",
    "\n",
    "    # ripple parameter, set to ten percent\n",
    "    ripple=0.1;  \n",
    "\n",
    "    # normalise frequency\n",
    "    fl=2.0*flow/rate;\n",
    "    fh=2.0*fhigh/rate;\n",
    "\n",
    "    # center frequency \n",
    "    cf = 4 * tan( (fl*pi/2) ) * tan( (fh*pi/2) );\n",
    "\n",
    "    # bandwidth\n",
    "    bw = 2 * ( tan( (fh*pi/2) ) - tan( (fl*pi/2) ) );\n",
    "\n",
    "    # ripple parameter factor\n",
    "    rpf = sqrt((sqrt((1.0+1.0/(ripple*ripple))) + 1.0/ripple));\n",
    "    a = 0.5*(rpf-1.0/rpf);\n",
    "    b = 0.5*(rpf+1.0/rpf);\n",
    "\n",
    "    u=np.zeros((5,1));\n",
    "    v=np.zeros((5,1));\n",
    "    theta = 3*pi/4;\n",
    "    sr = a * cos(theta);\n",
    "    si = b * sin(theta);\n",
    "    es = sqrt(sr*sr+si*si);\n",
    "    tmp= 16.0 - 16.0*bw*sr + 8.0*cf + 4.0*es*es*bw*bw - 4.0*bw*cf*sr + cf*cf;\n",
    "    v[0,0] = 1.0;\n",
    "    v[1,0] = 4.0*(-16.0 + 8.0*bw*sr - 2.0*bw*cf*sr + cf*cf)/tmp;\n",
    "    v[2,0] = (96.0 - 16.0*cf - 8.0*es*es*bw*bw + 6.0*cf*cf)/tmp;\n",
    "    v[3,0] = (-64.0 - 32.0*bw*sr + 8.0*bw*cf*sr + 4.0*cf*cf)/tmp;\n",
    "    v[4,0] = (16.0 + 16.0*bw*sr + 8.0*cf + 4.0*es*es*bw*bw + 4.0*bw*cf*sr + cf*cf)/tmp;\n",
    "    tmp = 4.0*es*es*bw*bw/tmp;\n",
    "    u[0,0] = tmp;\n",
    "    u[1,0] = 0.0;\n",
    "    u[2,0] = -2.0*tmp;\n",
    "    u[3,0] = 0.0;\n",
    "    u[4,0] = tmp;\n",
    "\n",
    "    dout = sg.lfilter(u.ravel(),v.ravel(),dd.ravel());\n",
    "    return (dout,u,v);\n",
    "\n",
    "\n",
    "def eda_timedelay( uA, uB, Dt ):\n",
    "# time delay between two similar timeseries, by\n",
    "# cross-correlation followed by parabilic refinment.\n",
    "# delay is positive when pulse on B is later than pulse on A\n",
    "# input:\n",
    "# uA, uB, two signals to be compared, as N-by-1 np.arrays\n",
    "#         should have zero mean\n",
    "#         if of unequal length, the shorter is zero padded\n",
    "# Dt, sampling of both signals\n",
    "# output\n",
    "# tBmA_est, estimated time delay, positive when dB is delated with respect to uA\n",
    "# Cmax_est, cross-correlation at estimated delay\n",
    "\n",
    "    # lengths\n",
    "    sh = np.shape(uA);\n",
    "    NA = sh[0];\n",
    "    sh = np.shape(uB);\n",
    "    NB = sh[0];\n",
    "    N = np.max( eda_cvec([NA,NB]) ) ;\n",
    "\n",
    "    # cross correlate\n",
    "    c = eda_cvec( np.correlate(uA.ravel(), uB.ravel(), mode='full') );\n",
    "\n",
    "    # rough estimate of lag\n",
    "    Nc = 2*N-1;\n",
    "    Ncenter = (Nc+1)/2;\n",
    "    cmax = np.max(c);\n",
    "    icmax = np.argmax(c);\n",
    "    tBmArough = -Dt * (icmax-N+1);\n",
    "    \n",
    "    # parabolic refinement\n",
    "    # (c-cmax) = m(1) + m(2) dt + m(3) dt^2\n",
    "    # d(c-cmax)/ddt = 0 = m(2) + 2 m(3) dt\n",
    "    # dt = -0.5*m(2)/m(3)\n",
    "    # c = cmax + m(1) + m(2) dt + m(3) dt^2\n",
    "    Dt2 = Dt**2;\n",
    "    g1 =     eda_cvec([1.0, 1.0, 1.0]);\n",
    "    g2 =  Dt*eda_cvec([-1.0, 0.0, 1.0]);\n",
    "    g3 = Dt2*eda_cvec([1.0, 0.0, 1.0]);\n",
    "    G = np.concatenate( (g1,g2,g3), axis=1);\n",
    "    d = eda_cvec( [ c[icmax-1,0]-cmax, c[icmax,0]-cmax, c[icmax+1,0]-cmax ] );\n",
    "    m = la.solve( np.matmul(G.transpose(),G), np.matmul(G.transpose(),d) );\n",
    "    dt = -0.5*m[1,0]/m[2,0];\n",
    "    Cmax_est = cmax + m[0,0] + m[1,0]*dt + m[2,0]*Dt2;\n",
    "    tBmA_est = tBmArough - dt;\n",
    "\n",
    "    return tBmA_est, Cmax_est;\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "correlation coefficients:\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA3QAAADnCAYAAABWihTLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAG9UlEQVR4nO3ZIW4VWBiG4XMnTW+6gibUl+JRGCzFgavoAnAsBVdfg0eUpKKqFtEgei2GBEVYwZ01zORPZt70eRbw5eS4N/9mv98vAAAAev76rx8AAADAvyPoAAAAogQdAABAlKADAACIEnQAAABRgg4AACBK0AEAAEQJOgAAgChBBwAAECXoAAAAogQdAABAlKADAACIEnQAAABRgg4AACBK0AEAAEQJOgAAgChBBwAAECXoAAAAogQdAABAlKADAACIEnQAAABRgg4AACBK0AEAAEQJOgAAgChBBwAAECXoAAAAogQdAABAlKADAACIEnQAAABRgg4AACBK0AEAAEQJOgAAgChBBwAAECXoAAAAogQdAABAlKADAACIEnQAAABRgg4AACBK0AEAAEQJOgAAgChBBwAAECXoAAAAog6mhjabzX5i5+rqamJmnZ+fj+w8PDz8r3bWWuvy8nJk5+7ubmTn1atXIztnZ2ebkSEAAHgiXOgAAACiBB0AAECUoAMAAIgSdAAAAFGCDgAAIErQAQAARAk6AACAKEEHAAAQJegAAACiBB0AAECUoAMAAIgSdAAAAFGCDgAAIErQAQAARAk6AACAKEEHAAAQJegAAACiDqaGrq6uRnY+fPgwsvPz58+Rnfv7+5Gdo6OjkZ211trtdiM719fXIzvb7XZk5+zsbGQHAACeChc6AACAKEEHAAAQJegAAACiBB0AAECUoAMAAIgSdAAAAFGCDgAAIErQAQAARAk6AACAKEEHAAAQJegAAACiBB0AAECUoAMAAIgSdAAAAFGCDgAAIErQAQAARAk6AACAqM1+vx8Z+vHjx8jQ4eHhxMx69uzZyM6nT59Gdm5vb0d21lpru92O7Lx+/Xpk582bNyM7z58/34wMAQDAE+FCBwAAECXoAAAAogQdAABAlKADAACIEnQAAABRgg4AACBK0AEAAEQJOgAAgChBBwAAECXoAAAAogQdAABAlKADAACIEnQAAABRgg4AACBK0AEAAEQJOgAAgChBBwAAELXZ7/cjQ1++fBkZur+/n5hZJycnIzsfP34c2bm5uRnZWWut7XY7sjP1R5vNZmTn9PR0ZggAAJ4IFzoAAIAoQQcAABAl6AAAAKIEHQAAQJSgAwAAiBJ0AAAAUYIOAAAgStABAABECToAAIAoQQcAABAl6AAAAKIEHQAAQJSgAwAAiBJ0AAAAUYIOAAAgStABAABECToAAICog6mhh4eHkZ2jo6ORndvb25Gdm5ubkZ23b9+O7Ky11sXFxcjO8fHxyM779+9Hdk5PT0d2AADgqXChAwAAiBJ0AAAAUYIOAAAgStABAABECToAAIAoQQcAABAl6AAAAKIEHQAAQJSgAwAAiBJ0AAAAUYIOAAAgStABAABECToAAIAoQQcAABAl6AAAAKIEHQAAQJSgAwAAiBJ0AAAAUQdTQ5eXlyM7u91uZOf79+8jO9vtdmTn4uJiZGettT5//jyy8/j4OLLz58+fkR0AAOCfcaEDAACIEnQAAABRgg4AACBK0AEAAEQJOgAAgChBBwAAECXoAAAAogQdAABAlKADAACIEnQAAABRgg4AACBK0AEAAEQJOgAAgChBBwAAECXoAAAAogQdAABAlKADAACIOpgauru7G9m5vr4e2Xn37t3IzsnJycjO8fHxyM5aaz0+Po7svHjxYmTn69evIzsAAMA/40IHAAAQJegAAACiBB0AAECUoAMAAIgSdAAAAFGCDgAAIErQAQAARAk6AACAKEEHAAAQJegAAACiBB0AAECUoAMAAIgSdAAAAFGCDgAAIErQAQAARAk6AACAKEEHAAAQtdnv9yNDu91uZOjbt28TM+vly5cjO5vNZmTn169fIztrrXV4eDiy8/v375Gd8/PzkZ39fj/z2QAA8ES40AEAAEQJOgAAgChBBwAAECXoAAAAogQdAABAlKADAACIEnQAAABRgg4AACBK0AEAAEQJOgAAgChBBwAAECXoAAAAogQdAABAlKADAACIEnQAAABRgg4AACBK0AEAAERt9vv9f/0GAAAA/gUXOgAAgChBBwAAECXoAAAAogQdAABAlKADAACIEnQAAABRgg4AACBK0AEAAEQJOgAAgChBBwAAECXoAAAAogQdAABAlKADAACIEnQAAABRgg4AACBK0AEAAEQJOgAAgChBBwAAECXoAAAAogQdAABAlKADAACIEnQAAABRgg4AACBK0AEAAEQJOgAAgChBBwAAECXoAAAAogQdAABAlKADAACIEnQAAABRgg4AACBK0AEAAEQJOgAAgChBBwAAECXoAAAAogQdAABAlKADAACIEnQAAABRgg4AACBK0AEAAEQJOgAAgChBBwAAEPU3xIhyoFWo2DsAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 1152x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "mean and standard deviations:\n",
      "    SiO2: 50.388700 7.209022\n",
      "    TiO2: 1.364637 1.208305\n",
      "    Al203: 15.255543 3.949684\n",
      "    FeO: 9.860720 3.112876\n",
      "    MgO: 7.876269 2.853069\n",
      "    CaO: 11.641090 3.511410\n",
      "    Na2O: 2.471665 1.469694\n",
      "    K2O: 0.190218 0.387298\n",
      " \n",
      "covariance and correlation coefficient:\n",
      "    SiO2 vs TiO2: 0.176776 0.361535\n",
      "    SiO2 vs Al203: 0.103385 0.049821\n",
      "    SiO2 vs FeO: 0.427913 0.249633\n",
      "    SiO2 vs MgO: -1.869543 -0.449582\n",
      "    SiO2 vs CaO: 0.098530 0.057043\n",
      "    SiO2 vs Na2O: 0.203303 0.334273\n",
      "    SiO2 vs K2O: 0.021383 0.088796\n",
      "    TiO2 vs Al203: 0.013893 0.022992\n",
      "    TiO2 vs FeO: 0.251822 0.504489\n",
      "    TiO2 vs MgO: -0.568284 -0.469300\n",
      "    TiO2 vs CaO: -0.112748 -0.224158\n",
      "    TiO2 vs Na2O: 0.129429 0.730804\n",
      "    TiO2 vs K2O: 0.021179 0.302024\n",
      "    Al203 vs FeO: -0.819501 -0.386845\n",
      "    Al203 vs MgO: -3.517052 -0.684375\n",
      "    Al203 vs CaO: 1.407959 0.659577\n",
      "    Al203 vs Na2O: 0.291717 0.388115\n",
      "    Al203 vs K2O: 0.024097 0.080970\n",
      "    FeO vs MgO: -0.959504 -0.226022\n",
      "    FeO vs CaO: -0.352446 -0.199874\n",
      "    FeO vs Na2O: 0.051783 0.083402\n",
      "    FeO vs K2O: -0.019208 -0.078130\n",
      "    MgO vs CaO: -2.664767 -0.622947\n",
      "    MgO vs Na2O: -0.791771 -0.525672\n",
      "    MgO vs K2O: -0.104875 -0.175851\n",
      "    CaO vs Na2O: -0.056522 -0.090343\n",
      "    CaO vs K2O: -0.037734 -0.152324\n",
      "    Na2O vs K2O: 0.019353 0.221870\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 16464 characters]\n",
      "text/plain": [
       "<Figure size 432x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# eda09_01\n",
    "# estimate covariance and etc of rorck data types\n",
    "\n",
    "# read data file\n",
    "D = np.genfromtxt('rocks.txt', delimiter='\\t')\n",
    "[N, M]=D.shape;\n",
    "# 0: SiO2\n",
    "# 1: TiO2\n",
    "# 2: Al203\n",
    "# 3: FeO-total\n",
    "# 4: MgO\n",
    "# 5: CaO\n",
    "# 6: Na2O\n",
    "# 7: K2O\n",
    "\n",
    "# note use of a tupple to create array of variable names\n",
    "name = ( 'SiO2', 'TiO2', 'Al203', 'FeO', 'MgO', 'CaO', 'Na2O', 'K2O' );\n",
    "\n",
    "dbar = eda_cvec( np.mean(D,axis=0) ); # mean\n",
    "C = np.cov(D,rowvar=False); # covariance\n",
    "sigma = eda_cvec( np.sqrt(np.diag(D)) ); # stddev\n",
    "R = np.corrcoef(D,rowvar=False); # correlation coefficient\n",
    "\n",
    "# plot correlation coefficent matrix\n",
    "print('correlation coefficients:')\n",
    "eda_draw(np.abs(R));\n",
    "\n",
    "# display results; also find most correlated data\n",
    "print('mean and standard deviations:')\n",
    "for i in range(M):\n",
    "    print( \"    %s: %f %f\" % (name[i], dbar[i,0], sigma[i,0]) );\n",
    "\n",
    "print(' ');\n",
    "icor=1;\n",
    "jcor=2;\n",
    "r=0;\n",
    "print('covariance and correlation coefficient:')\n",
    "for i in range(M):\n",
    "    for j in range(i+1,M):\n",
    "        print( \"    %s vs %s: %f %f\" % (name[i], name[j], C[i,j], R[i,j])) ;\n",
    "        if( abs(R[i,j]) > r ):\n",
    "            r = abs(R[i,j]);\n",
    "            icor=i;\n",
    "            jcor=j;\n",
    "        \n",
    "fig2 = plt.figure(2,figsize=(6,6));\n",
    "ax1 = plt.subplot(1,1,1);\n",
    "plt.plot(D[:,icor],D[:,jcor],'ko');\n",
    "plt.xlabel(name[icor]);\n",
    "plt.ylabel(name[jcor]);\n",
    "plt.title('most correlated variables');\n",
    "plt.show();\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 32728 characters]\n",
      "text/plain": [
       "<Figure size 1152x288 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# eda09_02\n",
    "# load the neuse.txt data set\n",
    "# make scatter plots of lagged samples\n",
    "# in order to investigate correlation\n",
    "\n",
    "# red data file\n",
    "D = np.genfromtxt('neuse.txt', delimiter='\\t')\n",
    "[N, M]=D.shape;\n",
    "\n",
    "# break out time and data vectors\n",
    "t = eda_cvec( D[:,0] );\n",
    "d = eda_cvec( D[:,1] ); \n",
    "\n",
    "fig1 = plt.figure(1,figsize=(16,4));\n",
    "\n",
    "lag=1;\n",
    "ax1 = plt.subplot(1,3,1);\n",
    "plt.axis([0, 25000, 0, 25000]);\n",
    "plt.plot(d[0:N-lag,0], d[lag:N,0],'ko');\n",
    "plt.xlabel('discharge');\n",
    "plt.ylabel(\"discharge lagged by %d days\" % (lag) );\n",
    "\n",
    "lag=3;\n",
    "ax2 = plt.subplot(1,3,2);\n",
    "plt.axis([0, 25000, 0, 25000]);\n",
    "plt.plot(d[0:N-lag,0], d[lag:N,0],'ko');\n",
    "plt.xlabel('discharge');\n",
    "plt.ylabel(\"discharge lagged by %d days\" % (lag) );\n",
    "\n",
    "lag=30;\n",
    "ax3 = plt.subplot(1,3,3);\n",
    "plt.axis([0, 25000, 0, 25000]);\n",
    "plt.plot(d[0:N-lag,0], d[lag:N,0],'ko');\n",
    "plt.xlabel('discharge');\n",
    "plt.ylabel(\"discharge lagged by %d days\" % (lag) );\n",
    "\n",
    "plt.show();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 54964 characters]\n",
      "text/plain": [
       "<Figure size 1152x432 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# eda09_03\n",
    "# load the neuse.txt data set\n",
    "# compute and plot autocorrelation\n",
    "\n",
    "# read data file\n",
    "D = np.genfromtxt('neuse.txt', delimiter='\\t')\n",
    "[N, M]=D.shape;\n",
    "\n",
    "# breal out time and data vectors\n",
    "t = eda_cvec( D[:,0] );\n",
    "d = eda_cvec( D[:,1] );\n",
    "Dt=t[1,0]-t[0,0]; # sampling interval\n",
    "\n",
    "# demean the data\n",
    "d=d-np.mean(d);\n",
    "\n",
    "# compute the autocorrelation\n",
    "a = np.correlate(d.ravel(), d.ravel(), mode='full');\n",
    "\n",
    "# time axis for autocorrelation\n",
    "tau = eda_cvec( np.linspace(-N+1,N-1,2*N-1) );\n",
    "\n",
    "# number pf points correlated for each lag\n",
    "Nk = np.concatenate( ( np.linspace(1,N,N), np.linspace(N-1,1,N-1) ) );\n",
    "\n",
    "# correction for length of timeseries\n",
    "a = N*np.divide( a, Nk );\n",
    "\n",
    "fig1 = plt.figure(1,figsize=(16,6));\n",
    "\n",
    "lag=1;\n",
    "ax1 = plt.subplot(2,1,1);\n",
    "plt.axis([ -30, 30, 0, 1.1*np.max(a) ]);\n",
    "plt.plot( tau, a, 'k-');\n",
    "plt.xlabel('lag (days)');\n",
    "plt.ylabel('autocorrelation');\n",
    "\n",
    "ax2 = plt.subplot(2,1,2);\n",
    "plt.axis([ -365.25*10, 365.25*10, -1.1*np.max(a), 1.1*np.max(a) ]);\n",
    "plt.plot( tau, a,'k-');\n",
    "plt.xlabel('lag (days)');\n",
    "plt.ylabel('autocorrelation');\n",
    "\n",
    "plt.show();\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 79252 characters]\n",
      "text/plain": [
       "<Figure size 1152x432 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 34760 characters]\n",
      "text/plain": [
       "<Figure size 1152x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "lag: true: 12.49 estimated: 12.50\n"
     ]
    }
   ],
   "source": [
    "# eda09_04\n",
    "# using the cross correlation to align two timeseries\n",
    "# this example uses synthetic data\n",
    "\n",
    "# time series paramaters\n",
    "N=1000; # number of data\n",
    "Dt=0.1; # sampling interval\n",
    "t = eda_cvec( Dt*np.linspace(0,N-1,N) );\n",
    "\n",
    "# synthetic time series, u and v\n",
    "t1 = Dt*4*N/8 + Dt/10.0;\n",
    "t2 = Dt*5*N/8;\n",
    "tlagtrue = t2-t1;\n",
    "w0 = 2*pi*0.5;\n",
    "u = (np.multiply( np.cos( w0 * (t-t1) ), np.exp( -np.power((t-t1),2)/2 )) + np.random.normal(0.0,0.02,(N,1))).ravel();\n",
    "v = (np.multiply( np.cos( w0 * (t-t2) ), np.exp( -np.power((t-t2),2)/2 )) + np.random.normal(0.0,0.02,(N,1))).ravel();\n",
    "\n",
    "fig1 = plt.figure(1,figsize=(16,6));\n",
    "ax1 = plt.subplot(2,1,1);\n",
    "plt.axis([t[0,0], t[N-1,0], -1.1, 1.1]);\n",
    "plt.plot( t, u+0.1, 'k-');\n",
    "plt.plot( t, v-0.1, 'k-');\n",
    "plt.xlabel('time');\n",
    "plt.ylabel('u(t) and v(t)');\n",
    "\n",
    "# cross correlation\n",
    "c = eda_cvec( np.correlate(u.ravel(), v.ravel(), mode='full') );\n",
    "\n",
    "# time asxis for autocorrelation\n",
    "tau = np.zeros((2*N-1,1));\n",
    "tau[:,0] = Dt*np.linspace(-N+1,N-1,2*N-1);\n",
    "\n",
    "# find best lag\n",
    "Nc = 2*N-1;\n",
    "Ncenter = (Nc+1)/2;\n",
    "cmax = np.max(c);\n",
    "icmax = np.argmax(c);\n",
    "tlag = -Dt * (icmax-N+1);\n",
    "\n",
    "# plot lagged timeseries\n",
    "ax2 = plt.subplot(2,1,2);\n",
    "plt.axis([t[0,0], t[N-1,0], -1.1, 1.1]);\n",
    "plt.plot( t, u+0.1, 'k-');\n",
    "plt.plot( t-tlag, v-0.1, 'k-');\n",
    "plt.xlabel('time');\n",
    "plt.ylabel('u(t) and v(t)');\n",
    "\n",
    "#  plot cross-correlation\n",
    "fig2 = plt.figure(2,figsize=(16,6));\n",
    "ax2 = plt.subplot(1,1,1);\n",
    "plt.axis( [-20, 20, -1.1*cmax, 1.1*cmax] );\n",
    "plt.plot( tau, c, 'k-');\n",
    "plt.plot( [-tlag, -tlag], [-cmax, cmax], 'k:');\n",
    "plt.xlabel('time')\n",
    "plt.ylabel('cross-correlation');\n",
    "plt.show();\n",
    "\n",
    "print( \"lag: true: %.2f estimated: %.2f\" % (tlagtrue, tlag) );"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "lag: true: 12.4900 estimated: 12.4889\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 80488 characters]\n",
      "text/plain": [
       "<Figure size 1152x432 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# eda09_04a\n",
    "# using the cross correlation to align two timeseries\n",
    "# this example uses eda_timedelay() function\n",
    "\n",
    "# time series paramaters\n",
    "N=1000; # number of data\n",
    "Dt=0.1; # sampling interval\n",
    "t = eda_cvec( Dt*np.linspace(0,N-1,N) );\n",
    "\n",
    "# synthetic time series, u and v\n",
    "t1 = Dt*4*N/8 + Dt/10.0;\n",
    "t2 = Dt*5*N/8;\n",
    "tlagtrue = t2-t1;\n",
    "w0 = 2*pi*0.5;\n",
    "u = (np.multiply( np.cos( w0 * (t-t1) ), np.exp( -np.power((t-t1),2)/2 )) + np.random.normal(0.0,0.02,(N,1))).ravel();\n",
    "v = (np.multiply( np.cos( w0 * (t-t2) ), np.exp( -np.power((t-t2),2)/2 )) + np.random.normal(0.0,0.02,(N,1))).ravel();\n",
    "\n",
    "fig1 = plt.figure(1,figsize=(16,6));\n",
    "ax1 = plt.subplot(2,1,1);\n",
    "plt.axis([t[0,0], t[N-1,0], -1.1, 1.1]);\n",
    "plt.plot( t, u+0.1, 'k-');\n",
    "plt.plot( t, v-0.1, 'k-');\n",
    "plt.xlabel('time');\n",
    "plt.ylabel('u(t) and v(t)');\n",
    "\n",
    "tlag, cmax = eda_timedelay(u, v, Dt)\n",
    "\n",
    "# plot lagged timeseries\n",
    "ax2 = plt.subplot(2,1,2);\n",
    "plt.axis([t[0,0], t[N-1,0], -1.1, 1.1]);\n",
    "plt.plot( t, u+0.1, 'k-');\n",
    "plt.plot( t-tlag, v-0.1, 'k-');\n",
    "plt.xlabel('time');\n",
    "plt.ylabel('u(t) and v(t)');\n",
    "\n",
    "print( \"lag: true: %.4f estimated: %.4f\" % (tlagtrue, tlag) );"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "unlagged data\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 94416 characters]\n",
      "text/plain": [
       "<Figure size 720x432 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "lagged data\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 93968 characters]\n",
      "text/plain": [
       "<Figure size 720x432 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "lagged data, zoomed\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 69720 characters]\n",
      "text/plain": [
       "<Figure size 720x432 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 41024 characters]\n",
      "text/plain": [
       "<Figure size 576x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "estimated lag: 3.00 hours\n"
     ]
    }
   ],
   "source": [
    "# eda09_05\n",
    "# using the cross correlation to align two timeseries\n",
    "# West Point ozone data\n",
    "\n",
    "# load data\n",
    "D = np.genfromtxt('ozone_nohead.txt', delimiter='\\t')\n",
    "[N, M]=D.shape;\n",
    "\n",
    "# break out time, ozone, radiation and temperature\n",
    "t = eda_cvec( D[:,0] );  # time in days after 1:00 AM on Aug 1, 1993\n",
    "Dt=1/24; # sampling interval\n",
    "tmin = t[0,0]; # start time\n",
    "tmax = t[N-1,0]; # end time\n",
    "ozone = eda_cvec( D[:,1] );     # ozone in ppm\n",
    "radiation = eda_cvec(D[:,2]);   # solar radiation in W/m2\n",
    "temp = eda_cvec(D[:,3]);        # temperature in deg C\n",
    "\n",
    "# plot original timeseries\n",
    "print('unlagged data');\n",
    "fig1 = plt.figure(1,figsize=(10,6));\n",
    "\n",
    "# solar radiation\n",
    "ax1 = plt.subplot(3,1,1);\n",
    "plt.axis([tmin, tmax, 0, np.max(radiation)]);\n",
    "plt.plot( t, radiation, 'k-');\n",
    "plt.xlabel('time (days)');\n",
    "plt.ylabel('solar (W/m2)');\n",
    "\n",
    "# ozone\n",
    "ax2 = plt.subplot(3,1,2);\n",
    "plt.axis([tmin, tmax, 0, np.max(ozone)]);\n",
    "plt.plot( t, ozone, 'k-');\n",
    "plt.xlabel('time (days)');\n",
    "plt.ylabel('ozone (ppm)');\n",
    "\n",
    "# both\n",
    "ax3 = plt.subplot(3,1,3);\n",
    "plt.axis([tmin, tmax, 0, 1.1]);\n",
    "plt.plot( t, radiation/np.max(radiation), 'k-', color=(0.7,0.7,0.7));\n",
    "plt.plot( t, ozone/np.max(ozone), 'k-');\n",
    "plt.xlabel('time (days)');\n",
    "plt.ylabel('both');\n",
    "\n",
    "plt.show();\n",
    "\n",
    "# cross correlation\n",
    "Nc = 2*N-1;\n",
    "c = eda_cvec( np.correlate(radiation.ravel(), ozone.ravel(), mode='full') );\n",
    "\n",
    "# time asxis for autocorrelation\n",
    "tau = eda_cvec( Dt*np.linspace(-N+1,N-1,2*N-1) );\n",
    "\n",
    "# find best lag\n",
    "Ncenter = (Nc+1)/2;\n",
    "cmax = np.max(c);\n",
    "icmax = np.argmax(c);\n",
    "tlag = -Dt * (icmax-N+1);\n",
    "\n",
    "# unlagged plot\n",
    "print('lagged data');\n",
    "fig2 = plt.figure(2,figsize=(10,6));\n",
    "\n",
    "# solar radiation\n",
    "ax1 = plt.subplot(3,1,1);\n",
    "plt.axis([tmin, tmax, 0, np.max(radiation)]);\n",
    "plt.plot( t, radiation, 'k-');\n",
    "plt.xlabel('time (days)');\n",
    "plt.ylabel('solar (W/m2)');\n",
    "\n",
    "# ozone\n",
    "ax2 = plt.subplot(3,1,2);\n",
    "plt.axis([tmin, tmax, 0, np.max(ozone)]);\n",
    "plt.plot( t-tlag, ozone, 'k-');\n",
    "plt.xlabel('time (days)');\n",
    "plt.ylabel('ozone (ppm)');\n",
    "\n",
    "# both\n",
    "ax3 = plt.subplot(3,1,3);\n",
    "plt.axis([tmin, tmax, 0, 1.1]);\n",
    "plt.plot( t, radiation/np.max(radiation), 'k-', color=(0.7,0.7,0.7));\n",
    "plt.plot( t-tlag, ozone/np.max(ozone), 'k-');\n",
    "plt.xlabel('time (days)');\n",
    "plt.ylabel('both');\n",
    "\n",
    "plt.show();\n",
    "\n",
    "# lagged plot, zommed\n",
    "days=5;\n",
    "print('lagged data, zoomed');\n",
    "fig2 = plt.figure(2,figsize=(10,6));\n",
    "\n",
    "# solar radiation\n",
    "ax1 = plt.subplot(3,1,1);\n",
    "plt.axis([t[0,0], t[days*24-1,0], 0, np.max(radiation)]);\n",
    "plt.plot( t, radiation, 'k-');\n",
    "plt.xlabel('time (days)');\n",
    "plt.ylabel('solar (W/m2)');\n",
    "\n",
    "# ozone\n",
    "ax2 = plt.subplot(3,1,2);\n",
    "plt.axis([t[0,0], t[days*24-1,0], 0, np.max(ozone)]);\n",
    "plt.plot( t-tlag, ozone, 'k-');\n",
    "plt.xlabel('time (days)');\n",
    "plt.ylabel('ozone (ppm)');\n",
    "\n",
    "# both\n",
    "ax3 = plt.subplot(3,1,3);\n",
    "plt.axis([t[0,0], t[days*24-1,0], 0, 1.1]);\n",
    "plt.plot( t, radiation/np.max(radiation), 'k-', color=(0.7,0.7,0.7));\n",
    "plt.plot( t-tlag, ozone/np.max(ozone), 'k-');\n",
    "plt.xlabel('time (days)');\n",
    "plt.ylabel('both');\n",
    "\n",
    "plt.show();\n",
    "\n",
    "plt.show();\n",
    "\n",
    "#plot cross-correlation\n",
    "fig3 = plt.figure(3,figsize=(8,6));\n",
    "ax2 = plt.subplot(1,1,1);\n",
    "plt.axis( [-5, 5, 0, 1.1*cmax] );\n",
    "plt.plot( tau, c, 'k-');\n",
    "plt.plot( [-tlag, -tlag], [-cmax, cmax], 'k:');\n",
    "plt.xlabel('time (days)')\n",
    "plt.ylabel('cross-correlation');\n",
    "plt.show();\n",
    "\n",
    "print( \"estimated lag: %.2f hours\" % (24*tlag) );\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 29296 characters]\n",
      "text/plain": [
       "<Figure size 576x288 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# eda09_06\n",
    "# example of a filter problem\n",
    "# using impulse repsonse of a heat-producing layer\n",
    "# with solution via generalized least squares\n",
    "# using bicongugate gradient method\n",
    "# with the prior information of smoothness\n",
    "# Note: symbol q for theta, temperture\n",
    "\n",
    "# time setup for g (length M)\n",
    "M=1024;\n",
    "N=M;\n",
    "Dtdays = 0.5;\n",
    "tdays = eda_cvec( Dtdays*np.linspace(0,Dtdays*(N-1),N) );\n",
    "Dtseconds = Dtdays*3600*24;\n",
    "tseconds = Dtseconds*tdays;\n",
    "\n",
    "# materal constants\n",
    "rho = 1500; # kg/m3\n",
    "cp = 800; # J / kg-K\n",
    "kappa = 1.25e-6; # m2/s\n",
    "z = 1; # m\n",
    "\n",
    "# impulse response; handle time zero separately from the rest\n",
    "g = np.zeros((M,1));\n",
    "g[0,0] = 0.0;\n",
    "P1 = (1/(rho*cp));\n",
    "P2 = (Dtseconds/sqrt(2*pi));\n",
    "P3 = np.reciprocal( np.sqrt(2*kappa*tseconds[1:M]));\n",
    "P4 = np.exp( -0.5*(z**2)*np.reciprocal(2*kappa*tseconds[1:M]));\n",
    "g[1:M,0] = (P1 * P2 * np.multiply(P3 ,P4)).ravel();\n",
    "edafilterg = g; # global defines in GLSFilterMul()\n",
    "\n",
    "# true heat production\n",
    "t1 = Dtdays*N/4;\n",
    "t2 = Dtdays*7*N/16;\n",
    "sd = Dtdays*N/256;\n",
    "Ah=10;\n",
    "P5 = np.exp(-np.power(tdays-t1,2) / (2*sd**2) );\n",
    "P6 = np.exp(-np.power(tdays-t2,2) / (2*sd**2) );\n",
    "htrue = Ah*P5 + 0.5*Ah*P6;\n",
    "\n",
    "# predict true temperature by convolving g and htrue; truncate to length N\n",
    "tmp = np.convolve(g.ravel(), htrue.ravel());\n",
    "qtrue=eda_cvec( tmp[0:N] );\n",
    "\n",
    "# create simulated observations, qobs = qtrue+noise\n",
    "sigmad=0.01*np.amax(qtrue);   # std dev is 1% of maximum value of qtrue\n",
    "qobs=qtrue+np.random.normal(0,sigmad,(N,1)); # observed = signal + noie\n",
    "\n",
    "# solve for h given qobs and g using generalized least squares\n",
    "# with a prior smoothness constraint using biconjugate gradient\n",
    "\n",
    "# autocorrelation, a, of g. Just keep positive values\n",
    "# it needs to be of at least length M, since g is of\n",
    "# length N, and ordinarily N>=M.\n",
    "tmp = np.correlate(g.ravel(),g.ravel(),mode='full');\n",
    "edaautoc = eda_cvec(tmp[M-1:2*M-1]);\n",
    "\n",
    "# cross-correlation, c, of qobs2 and g. Just keep positive values\n",
    "# also, it needs to be of length M\n",
    "tmp = np.correlate(qobs.ravel(), g.ravel(),mode='full');\n",
    "Ntmp = len(tmp)\n",
    "Ntmp2 = floor(Ntmp/2);\n",
    "Nc = Ntmp-Ntmp2;\n",
    "tmp2 = eda_cvec( tmp[Ntmp2: Ntmp] );\n",
    "c = eda_cvec( tmp2[0:M,0] );\n",
    "\n",
    "# set up H and h and f\n",
    "e=1.0*np.amax(g);  # e is a damping parameter that plays role of 1/sqrt(variance)\n",
    "edaHsparse = sparse.diags([e, -2*e, e], [0, 1, 2], shape=(M-2, M));\n",
    "h = np.zeros((M-2,1)); # prior information is that 2nd derivative is zero\n",
    "\n",
    "# set up F'f = GT qobs + HT h\n",
    "# GT qobs is c=qobs2*g\n",
    "HTh =  edaHsparse.transpose() * h;\n",
    "FTf = c + HTh;\n",
    "\n",
    "# define linear operator needed for conjugate gradienet solver\n",
    "LO=las.LinearOperator(shape=(M,M),matvec=GLSautoMul,rmatvec=GLSautoMul);\n",
    "\n",
    "# solve for estimated model parameters using congugate gradient algrorithm\n",
    "# note solver returns tupple first member of which is solution\n",
    "q=las.cg(LO,FTf,tol=1e-6, maxiter=(3*(N+M)+100));\n",
    "hest3 = eda_cvec( q[0] );\n",
    "\n",
    "# figure 1\n",
    "fig1 = plt.figure(1,figsize=(8,4));\n",
    "\n",
    "# plot observations\n",
    "ax1 = plt.subplot(3,1,1);\n",
    "plt.axis([tdays[0,0], tdays[N-1,0], 0, 1.1*np.max(qobs)] );\n",
    "plt.plot(tdays[0:N],qobs,'k-');\n",
    "plt.xlabel('time t, days');\n",
    "plt.ylabel('qobs(t)');\n",
    "plt.title('Autocorrelation Method');\n",
    "\n",
    "# plot true heat production\n",
    "ax1 = plt.subplot(3,1,2);\n",
    "plt.axis([tdays[0,0], tdays[N-1,0], -1.1*np.max(htrue), 1.1*np.max(htrue)] );\n",
    "plt.plot(tdays[0:N],htrue,'k-');\n",
    "plt.xlabel('time t, days');\n",
    "plt.ylabel('htrue(t)');\n",
    "\n",
    "# plot estimated heat production\n",
    "ax1 = plt.subplot(3,1,3);\n",
    "plt.axis([tdays[0,0], tdays[N-1,0], -1.1*np.max(htrue), 1.1*np.max(htrue)] );\n",
    "plt.plot(tdays[0:N],hest3,'k-');\n",
    "plt.xlabel('time t, days');\n",
    "plt.ylabel('hest(t)');\n",
    "plt.show();\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 76124 characters]\n",
      "text/plain": [
       "<Figure size 720x648 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# eda09_07\n",
    "# smooth neuse.txt\n",
    "\n",
    "# load data\n",
    "D = np.genfromtxt('neuse.txt', delimiter='\\t')\n",
    "[N, M]=D.shape;\n",
    "\n",
    "# break out time and data into separate vectors\n",
    "t = eda_cvec( D[:,0] );         # time in days after 1:00 AM on Aug 1, 1993\n",
    "Dt=t[1,0]-t[0,0];               # sampling interval\n",
    "tmin = t[0,0];                  # start time\n",
    "tmax = t[N-1,0];                # end time\n",
    "d = eda_cvec( D[:,1] );         # discharge in cfs\n",
    "\n",
    "# construct 3-point smoothing filter\n",
    "L1 = 3; # width of filter; must be odd\n",
    "s1 = eda_cvec( [1, 2, 1] );\n",
    "s1 = s1/np.sum(s1);  # normalize to sum=1\n",
    "Lh1 = floor(L1/2); # (half width)-1\n",
    "\n",
    "# filter and time shift\n",
    "temp = np.convolve(s1.ravel(), d.ravel(), mode='full');\n",
    "d1 = eda_cvec( temp[Lh1:N+Lh1] );\n",
    "\n",
    "# construct 19-point smoothing filter\n",
    "L2 = 19; # width of fileter, must be odd\n",
    "s2 = eda_cvec( [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 9, 8, 7, 6, 5, 4, 3, 2, 1] );\n",
    "s2 = s2/np.sum(s1); # normalize\n",
    "Lh2 = floor(L2/2); # (half width) - 1\n",
    "\n",
    "# filter and time shift\n",
    "temp = np.convolve(s2.ravel(), d.ravel(), mode='full');\n",
    "d2 = eda_cvec( temp[Lh2:N+Lh2] );\n",
    "\n",
    "# plot only this time range (for clarity)\n",
    "N1=0;\n",
    "N2=500;\n",
    "\n",
    "fig1 = plt.figure(1,figsize=(10,9));\n",
    "\n",
    "# original\n",
    "ax1 = plt.subplot(3,1,1);\n",
    "plt.axis([t[N1,0], t[N2,0], 0, np.max(d)]);\n",
    "plt.plot( t, d, 'k-');\n",
    "plt.xlabel('time (days)');\n",
    "plt.ylabel('discharge (cfs)');\n",
    "\n",
    "# filtered with 3-pt filter\n",
    "ax2 = plt.subplot(3,1,2);\n",
    "plt.axis([t[N1,0], t[N2,0], 0, np.max(d1)]);\n",
    "plt.plot( t, d1, 'k-');\n",
    "plt.xlabel('time (days)');\n",
    "plt.ylabel(\"%d point filter\" % L1);\n",
    "\n",
    "# filtered with 19-pt filter\n",
    "ax3 = plt.subplot(3,1,3);\n",
    "plt.axis([t[N1,0], t[N2,0], 0, np.max(d2)]);\n",
    "plt.plot( t, d2, 'k-');\n",
    "plt.xlabel('time (days)');\n",
    "plt.ylabel(\"%d point filter\" % L2);\n",
    "\n",
    "plt.show();\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 37356 characters]\n",
      "text/plain": [
       "<Figure size 720x648 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# eda09_08\n",
    "# amplitude spectrum of uniform smoothing filter\n",
    "\n",
    "# standard set up of time\n",
    "N=256;             # length of time series, must be even\n",
    "Dt=0.01;           # sampling interval\n",
    "tmin = 0;          # start time\n",
    "tmax = Dt*(N-1);   # end time\n",
    "t=eda_cvec(np.linspace(tmin,tmax,N)); # time vector\n",
    "\n",
    "# standard set up of frequency\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 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 frequency vector\n",
    "Nw=Nf;             # number of non-negative angular frequencies\n",
    "fpos=eda_cvec(Df * np.linspace(0,No2,Nf) );  # vector non-negative frequencies\n",
    "wpos=2*pi*fpos;    # vector non-negative angular frequencies\n",
    "\n",
    "# evaluate fourier transform of filter of length L\n",
    "T1=3*Dt;                       # width in time domain\n",
    "s1 = np.sinc(wpos*T1/(2*pi));  # sinc function in frequency\n",
    "sa1 = np.abs(s1);              # amplitude spectral density\n",
    "\n",
    "# evaluate fourier transform of filter of length L\n",
    "T2=21*Dt;                      # width in time domain\n",
    "s2 = np.sinc(wpos*T2/(2*pi));  # sinc function in frequency\n",
    "sa2 = np.abs(s2);              # amplitude spectral density\n",
    "\n",
    "# plot\n",
    "fig1 = plt.figure(1,figsize=(10,9));\n",
    "\n",
    "ax1 = plt.subplot(2,1,1);\n",
    "plt.axis( [0,fmax, 0, 1.1]);\n",
    "plt.plot( fpos, sa1, 'k-');\n",
    "plt.xlabel('frequency (Hz)');\n",
    "plt.ylabel('|s| for L=3');\n",
    "\n",
    "ax2 = plt.subplot(2,1,2);\n",
    "plt.axis( [0,fmax, 0, 1.1]);\n",
    "plt.plot( fpos, sa2, 'k-');\n",
    "plt.xlabel('frequency (Hz)');\n",
    "plt.ylabel('|s| for L=21');\n",
    "\n",
    "plt.show();\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 69,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "length 3 samples (0.030 seconds) sigma-t samples 0.816 (0.008 seconds)\n",
      "length 21 samples (0.210 seconds) sigma-t samples 6.055 (0.061 seconds)\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 28900 characters]\n",
      "text/plain": [
       "<Figure size 720x648 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# eda09_09\n",
    "# amplitude spectrum of Normal smoothing filter\n",
    "\n",
    "# standard set up of time\n",
    "N=256;             # length of time series, must be even\n",
    "Dt=0.01;           # sampling interval\n",
    "tmin = 0;          # start time\n",
    "tmax = Dt*(N-1);   # end time\n",
    "t=eda_cvec(np.linspace(tmin,tmax,N)); # time vector\n",
    "\n",
    "# standard set up of frequency\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 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 frequency vector\n",
    "Nw=Nf;             # number of non-negative angular frequencies\n",
    "fpos=eda_cvec(Df * np.linspace(0,No2,Nf) );  # vector non-negative frequencies\n",
    "wpos=2*pi*fpos;    # vector non-negative angular frequencies\n",
    "\n",
    "# evaluate fourier transform of filter, g, of length L\n",
    "L1=3;\n",
    "\n",
    "# compute variance of filter (as if it were distribution)\n",
    "g1 = np.ones((L1,1));                    # vector of 1's\n",
    "g1 = g1/np.sum(g1);                      # normalize\n",
    "x = eda_cvec( np.linspace(0,L1-1,L1) );  # vector of x's\n",
    "mg1 = np.sum( np.multiply( x, g1 ) ) ;   # sum of ( x times g1 )\n",
    "vg1 = np.sum( np.multiply( np.power(x-mg1,2),  g1) ); # variance = sum of ( (x-mean)**2 times g1)\n",
    "sg1 = Dt*sqrt(vg1);                      # Dt * sqrt of variance\n",
    "\n",
    "# fourier transform of gausian, using analytic formula\n",
    "sigma1=1/sg1;\n",
    "s1 = np.exp(-np.power(wpos,2)/(2*(sigma1**2)));\n",
    "sa1 = np.abs(s1);  # amplitude spreal density\n",
    "print( \"length %d samples (%.3f seconds) sigma-t samples %.3f (%.3f seconds)\" % (L1, Dt*L1, sg1/Dt, sg1) );\n",
    "\n",
    "# plot\n",
    "fig1 = plt.figure(1,figsize=(10,9));\n",
    "\n",
    "ax1 = plt.subplot(2,1,1);\n",
    "plt.axis( [0,fmax, 0, 1.1]);\n",
    "plt.plot( fpos, sa1, 'k-');\n",
    "plt.xlabel('frequency (Hz)');\n",
    "plt.ylabel('|s| for L=3');\n",
    "\n",
    "# evaluate fourier transform of filter, g, of length L\n",
    "L2=21;\n",
    "\n",
    "# compute variance of filter (as if it were distribution)\n",
    "g2 = np.ones((L2,1));\n",
    "g2 = g2/np.sum(g2);\n",
    "x = eda_cvec( np.linspace(0,L2-1,L2) );\n",
    "mg2 = np.sum( np.multiply( x, g2 ) );  # sum of ( x times g2 )\n",
    "vg2 = np.sum( np.multiply( np.power(x-mg2,2),  g2) ); # variance = sum of ( (x-mean)**2 times g2)\n",
    "sg2 = Dt*sqrt(vg2); # standard deviation \n",
    "\n",
    "# fourier transform of gausian, using analytic formula\n",
    "sigma2=1/sg2;  # variance in frequency domain\n",
    "s2 = np.exp(-np.power(wpos,2)/(2*(sigma2**2)));  # fourier transform\n",
    "sa2 = np.abs(s2); # amplitude spectral density\n",
    "print( \"length %d samples (%.3f seconds) sigma-t samples %.3f (%.3f seconds)\" % (L2, Dt*L2, sg2/Dt, sg2) );\n",
    "\n",
    "ax2 = plt.subplot(2,1,2);\n",
    "plt.axis( [0,fmax, 0, 1.1]);\n",
    "plt.plot( fpos, sa2, 'k-');\n",
    "plt.xlabel('frequency (Hz)');\n",
    "plt.ylabel('|s| for L=21');\n",
    "\n",
    "plt.show();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 78,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "zero at 0.909\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 59904 characters]\n",
      "text/plain": [
       "<Figure size 864x432 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# eda09_10\n",
    "# complex z-plane analysis of filters, case 1, highpass\n",
    "\n",
    "# standard set up of time\n",
    "N=32;              # length of time series, must be even\n",
    "Dt=1.0 ;           # sampling interval\n",
    "tmin = 0;          # start time\n",
    "tmax = Dt*(N-1);   # end time\n",
    "t=eda_cvec(np.linspace(tmin,tmax,N)); # time vector\n",
    "\n",
    "# standard set up of frequency\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 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 frequency vector\n",
    "Nw=Nf;             # number of non-negative angular frequencies\n",
    "fpos=eda_cvec(Df * np.linspace(0,No2,Nf) );  # vector non-negative frequencies\n",
    "wpos=2*pi*fpos;    # vector non-negative angular frequencies\n",
    "\n",
    "# a monomial filter, and its one root\n",
    "F = np.zeros((2,1));\n",
    "F[:,0] = [1, -1.1];\n",
    "zp = -F[0,0]/F[1,0];  # 0=F1+F2*z  so z=-F1/F2\n",
    "print(\"zero at %.3f\" % (zp) );\n",
    "\n",
    "# complex z-plane\n",
    "Nz=128;   # z-plane image is Nz by Nz\n",
    "zmin=-2.0;  # min z\n",
    "zmax=2.0;   # max z\n",
    "dz=(zmax-zmin)/(Nz-1); # sampling in z\n",
    "zr=eda_cvec(np.linspace(zmin,zmax,Nz)); # real-z vector\n",
    "zi=eda_cvec(np.linspace(zmin,zmax,Nz)); # imag-z vertor\n",
    "\n",
    "# evaluate filter polynomial everywhere in complex z-plane\n",
    "S = np.zeros((Nz,Nz));\n",
    "for n in range(Nz):\n",
    "    for m in range(Nz):\n",
    "        S[m,n] = abs( F[0,0] + F[1,0]*(zr[n,0]+complex(0,1)*zi[m,0]) )**2;\n",
    "\n",
    "# plot amplitude of polynomial\n",
    "fig1 = plt.figure(1,figsize=(12,6));\n",
    "\n",
    "# black and white colormap\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",
    "\n",
    "ax1 = plt.subplot(1,2,1);\n",
    "plt.axis( [zmin,zmax,zmin,zmax]);\n",
    "\n",
    "plt.imshow( S, cmap=bwcmap, vmin=np.min(S), vmax=np.max(S), extent=(zmin,zmax,zmin,zmax) );\n",
    "plt.xlabel('real z');\n",
    "plt.ylabel('imag z');\n",
    "plt.title('|S|^2');\n",
    "\n",
    "# plot Fourier transform points on the unit circle\n",
    "plt.plot(0,0,'wo');\n",
    "z=np.exp(-complex(0,1)*w*Dt);\n",
    "plt.plot(np.real(z),np.imag(z),'k+');\n",
    "\n",
    "# plot zero of filter\n",
    "plt.plot(np.real(zp),np.imag(zp),'w*');\n",
    "\n",
    "# plot frequency response of filter\n",
    "\n",
    "Sz = np.abs( F[0,0]+F[1,0]*z )**2;\n",
    "\n",
    "ax2 = plt.subplot(1,2,2);\n",
    "plt.axis( [0, fmax, 0, max(Sz)] );\n",
    "\n",
    "plt.xlabel('frequency, Hz');\n",
    "plt.ylabel('|S|^2');\n",
    "plt.title('|S|^2');\n",
    "plt.plot(fpos,Sz[0:Nf,0],'k-');\n",
    "plt.plot(fpos,Sz[0:Nf,0],'ko');\n",
    "\n",
    "plt.show();\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 79,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "zero at -0.909\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 60328 characters]\n",
      "text/plain": [
       "<Figure size 864x432 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# eda09_11\n",
    "# complex z-plane analysis of filters, case 1, lowpass\n",
    "\n",
    "# standard set up of time\n",
    "N=32;              # length of time series, must be even\n",
    "Dt=1.0 ;           # sampling interval\n",
    "tmin = 0;          # start time\n",
    "tmax = Dt*(N-1);   # end time\n",
    "t=eda_cvec(np.linspace(tmin,tmax,N)); # time vector\n",
    "\n",
    "# standard set up of frequency\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 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 frequency vector\n",
    "Nw=Nf;             # number of non-negative angular frequencies\n",
    "fpos=eda_cvec(Df * np.linspace(0,No2,Nf) );  # vector non-negative frequencies\n",
    "wpos=2*pi*fpos;    # vector non-negative angular frequencies\n",
    "\n",
    "# a monomial filter, and its one root\n",
    "F = eda_cvec( [1, 1.1] );\n",
    "zp = -F[0,0]/F[1,0];  # 0=F1+F2*z  so z=-F1/F2\n",
    "print(\"zero at %.3f\" % (zp) );\n",
    "\n",
    "# complex z-plane\n",
    "Nz=128;   # z-plane image is Nz by Nz\n",
    "zmin=-2.0;  # min z\n",
    "zmax=2.0;   # max z\n",
    "dz=(zmax-zmin)/(Nz-1); # sampling in z\n",
    "zr=eda_cvec(np.linspace(zmin,zmax,Nz)); # real-z vector\n",
    "zi=eda_cvec(np.linspace(zmin,zmax,Nz)); # imag-z vertor\n",
    "\n",
    "# evaluate filter polynomial everywhere in complex z-plane\n",
    "S = np.zeros((Nz,Nz));\n",
    "for n in range(Nz):\n",
    "    for m in range(Nz):\n",
    "        S[m,n] = abs( F[0,0] + F[1,0]*(zr[n,0]+complex(0,1)*zi[m,0]) )**2;\n",
    "\n",
    "# plot amplitude of polynomial\n",
    "fig1 = plt.figure(1,figsize=(12,6));\n",
    "\n",
    "# black and white colormap\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",
    "\n",
    "ax1 = plt.subplot(1,2,1);\n",
    "plt.axis( [zmin,zmax,zmin,zmax]);\n",
    "\n",
    "plt.imshow( S, cmap=bwcmap, vmin=np.min(S), vmax=np.max(S), extent=(zmin,zmax,zmin,zmax) );\n",
    "plt.xlabel('real z');\n",
    "plt.ylabel('imag z');\n",
    "plt.title('|S|^2');\n",
    "\n",
    "# plot Fourier transform points on the unit circle\n",
    "plt.plot(0,0,'wo');\n",
    "z=np.exp(-complex(0,1)*w*Dt);\n",
    "plt.plot(np.real(z),np.imag(z),'k+');\n",
    "\n",
    "# plot zero of filter\n",
    "plt.plot(np.real(zp),np.imag(zp),'w*');\n",
    "\n",
    "# plot frequency response of filter\n",
    "\n",
    "Sz = np.abs( F[0,0]+F[1,0]*z )**2;\n",
    "\n",
    "ax2 = plt.subplot(1,2,2);\n",
    "plt.axis( [0, fmax, 0, max(Sz)] );\n",
    "\n",
    "plt.xlabel('frequency, Hz');\n",
    "plt.ylabel('|S|^2');\n",
    "plt.title('|S|^2');\n",
    "plt.plot(fpos,Sz[0:Nf,0],'k-');\n",
    "plt.plot(fpos,Sz[0:Nf,0],'ko');\n",
    "\n",
    "plt.show();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 82,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "poles at (-0.754 -0.830) and (-0.754 0.830)\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 54056 characters]\n",
      "text/plain": [
       "<Figure size 864x432 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# eda09_12\n",
    "# complex z-plane analysis of filters, case 3, narrow band\n",
    "\n",
    "# standard set up of time\n",
    "N=32;              # length of time series, must be even\n",
    "Dt=1.0 ;           # sampling interval\n",
    "tmin = 0;          # start time\n",
    "tmax = Dt*(N-1);   # end time\n",
    "t=eda_cvec(np.linspace(tmin,tmax,N)); # time vector\n",
    "\n",
    "# standard set up of frequency\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 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 frequency vector\n",
    "Nw=Nf;             # number of non-negative angular frequencies\n",
    "fpos=eda_cvec(Df * np.linspace(0,No2,Nf) );  # vector non-negative frequencies\n",
    "wpos=2*pi*fpos;    # vector non-negative angular frequencies\n",
    "\n",
    "# a monomial filter, and its one root\n",
    "F1 = eda_cvec( [1, 0.6*(1-1.1*complex(0.0,1.0))] );\n",
    "zp1 = -F1[0,0]/F1[1,0];  # 0=F1+F2*z  so z=-F1/F2\n",
    "F2 = np.conj(F1);\n",
    "zp2 = -F2[0,0]/F2[1,0];\n",
    "print(\"poles at (%.3f %.3f) and (%.3f %.3f)\" % (np.real(zp1),np.imag(zp1),np.real(zp2),np.imag(zp2)) );\n",
    "\n",
    "# complex z-plane\n",
    "Nz=128;     # z-plane image is Nz by Nz\n",
    "zmin=-2.0;  # min z\n",
    "zmax=2.0;   # max z\n",
    "dz=(zmax-zmin)/(Nz-1); # sampling in z\n",
    "zr=eda_cvec(np.linspace(zmin,zmax,Nz)); # real-z vector\n",
    "zi=eda_cvec(np.linspace(zmin,zmax,Nz)); # imag-z vertor\n",
    "\n",
    "# evaluate filter polynomial everywhere in complex z-plane\n",
    "S = np.zeros((Nz,Nz));\n",
    "for n in range(Nz):\n",
    "    for m in range(Nz):\n",
    "        v1 = 1/(F1[0,0] + F1[1,0]*(zr[n,0]+complex(0.0,1.0)*zi[m,0]));\n",
    "        v2 = 1/(F2[0,0] + F2[1,0]*(zr[n,0]+complex(0.0,1.0)*zi[m,0]));\n",
    "        S[m,n] = log( abs( v1*v2 )**2);  # take the log for plot-scaling  purposes\n",
    "\n",
    "# plot amplitude of polynomial\n",
    "fig1 = plt.figure(1,figsize=(12,6));\n",
    "\n",
    "# black and white colormap\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",
    "\n",
    "ax1 = plt.subplot(1,2,1);\n",
    "plt.axis( [zmin,zmax,zmin,zmax]);\n",
    "\n",
    "plt.imshow( S, cmap=bwcmap, vmin=np.min(S), vmax=np.max(S), extent=(zmin,zmax,zmin,zmax) );\n",
    "plt.xlabel('real z');\n",
    "plt.ylabel('imag z');\n",
    "plt.title('|S|^2');\n",
    "\n",
    "# plot Fourier transform points on the unit circle\n",
    "plt.plot(0,0,'wo');\n",
    "z=np.exp(-complex(0,1)*w*Dt);\n",
    "plt.plot(np.real(z),np.imag(z),'k+');\n",
    "\n",
    "# plot poles of filter\n",
    "plt.plot(np.real(zp1),np.imag(zp1),'w*');\n",
    "plt.plot(np.real(zp2),np.imag(zp2),'w*');\n",
    "\n",
    "# plot frequency response of filter\n",
    "Sz = np.abs( np.reciprocal( np.multiply(F1[0,0]+F1[1,0]*z,F2[0,0]+F2[1,0]*z)) )**2;\n",
    "\n",
    "ax2 = plt.subplot(1,2,2);\n",
    "plt.axis( [0, fmax, 0, max(Sz)] );\n",
    "\n",
    "plt.xlabel('frequency, Hz');\n",
    "plt.ylabel('|S|^2');\n",
    "plt.title('|S|^2');\n",
    "plt.plot(fpos,Sz[0:Nf,0],'k-');\n",
    "plt.plot(fpos,Sz[0:Nf,0],'ko');\n",
    "\n",
    "plt.show();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 88,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "zeros at (-0.000 1.111) and (-0.000 -1.111)\n",
      "poles at (-0.000 1.250) and (-0.000 -1.250)\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 40620 characters]\n",
      "text/plain": [
       "<Figure size 864x432 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# eda09_13\n",
    "# complex z-plane analysis of filters, case 4, notch\n",
    "\n",
    "# standard set up of time\n",
    "N=32;              # length of time series, must be even\n",
    "Dt=1.0 ;           # sampling interval\n",
    "tmin = 0;          # start time\n",
    "tmax = Dt*(N-1);   # end time\n",
    "t=eda_cvec(np.linspace(tmin,tmax,N)); # time vector\n",
    "\n",
    "# standard set up of frequency\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 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 frequency vector\n",
    "Nw=Nf;             # number of non-negative angular frequencies\n",
    "fpos=eda_cvec(Df * np.linspace(0,No2,Nf) );  # vector non-negative frequencies\n",
    "wpos=2*pi*fpos;    # vector non-negative angular frequencies\n",
    "\n",
    "# a monomial filter, and its one root\n",
    "F1 = eda_cvec( [1, 0.9*complex(0.0,1.0)] );\n",
    "zp1 = -F1[0,0]/F1[1,0];  # 0=F1+F2*z  so z=-F1/F2\n",
    "F2 = np.conj(F1);\n",
    "zp2 = -F2[0,0]/F2[1,0];\n",
    "print(\"zeros at (%.3f %.3f) and (%.3f %.3f)\" % (np.real(zp1),np.imag(zp1),np.real(zp2),np.imag(zp2)) );\n",
    "\n",
    "F3 = eda_cvec( [1, 0.8*complex(0.0,1.0)] );\n",
    "zp3 = -F3[0,0]/F3[1,0];  # 0=F1+F2*z  so z=-F1/F2\n",
    "F4 = np.conj(F3);\n",
    "zp4 = -F4[0,0]/F4[1,0];\n",
    "print(\"poles at (%.3f %.3f) and (%.3f %.3f)\" % (np.real(zp3),np.imag(zp3),np.real(zp4),np.imag(zp4)) );\n",
    "\n",
    "# complex z-plane\n",
    "Nz=128;     # z-plane image is Nz by Nz\n",
    "zmin=-2.0;  # min z\n",
    "zmax=2.0;   # max z\n",
    "dz=(zmax-zmin)/(Nz-1); # sampling in z\n",
    "zr=eda_cvec(np.linspace(zmin,zmax,Nz)); # real-z vector\n",
    "zi=eda_cvec(np.linspace(zmin,zmax,Nz)); # imag-z vertor\n",
    "\n",
    "# evaluate filter polynomial everywhere in complex z-plane\n",
    "S = np.zeros((Nz,Nz));\n",
    "for n in range(Nz):\n",
    "    for m in range(Nz):\n",
    "        v1 = (F1[0,0] + F1[1,0]*(zr[n,0]+complex(0.0,1.0)*zi[m,0]));\n",
    "        v2 = (F2[0,0] + F2[1,0]*(zr[n,0]+complex(0.0,1.0)*zi[m,0]));\n",
    "        v3 = (F3[0,0] + F3[1,0]*(zr[n,0]+complex(0.0,1.0)*zi[m,0]));\n",
    "        v4 = (F4[0,0] + F4[1,0]*(zr[n,0]+complex(0.0,1.0)*zi[m,0]));\n",
    "        S[m,n] = log( abs( (v1*v2)/(v3*v4) )**2);  # take the log for plot-scaling  purposes\n",
    "\n",
    "# plot amplitude of polynomial\n",
    "fig1 = plt.figure(1,figsize=(12,6));\n",
    "\n",
    "# black and white colormap\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",
    "\n",
    "ax1 = plt.subplot(1,2,1);\n",
    "plt.axis( [zmin,zmax,zmin,zmax]);\n",
    "\n",
    "plt.imshow( S, cmap=bwcmap, vmin=np.min(S), vmax=np.max(S), extent=(zmin,zmax,zmin,zmax) );\n",
    "plt.xlabel('real z');\n",
    "plt.ylabel('imag z');\n",
    "plt.title('|S|^2');\n",
    "\n",
    "# plot Fourier transform points on the unit circle\n",
    "plt.plot(0,0,'wo');\n",
    "z=np.exp(-complex(0,1)*w*Dt);\n",
    "plt.plot(np.real(z),np.imag(z),'k+');\n",
    "\n",
    "# plot poles of filter\n",
    "plt.plot(np.real(zp1),np.imag(zp1),'w*');\n",
    "plt.plot(np.real(zp2),np.imag(zp2),'w*');\n",
    "plt.plot(np.real(zp3),np.imag(zp3),'w+');\n",
    "plt.plot(np.real(zp4),np.imag(zp4),'w+');\n",
    "\n",
    "# plot frequency response of filter\n",
    "partA = np.multiply(  F1[0,0]+F1[1,0]*z,  F2[0,0]+F2[1,0]*z);\n",
    "partB = np.multiply(  F3[0,0]+F3[1,0]*z,  F4[0,0]+F4[1,0]*z);\n",
    "Sz = np.abs(np.divide(partA, partB))**2;\n",
    "\n",
    "ax2 = plt.subplot(1,2,2);\n",
    "plt.axis( [0, fmax, 0, max(Sz)] );\n",
    "\n",
    "plt.xlabel('frequency, Hz');\n",
    "plt.ylabel('|S|^2');\n",
    "plt.title('|S|^2');\n",
    "plt.plot(fpos,Sz[0:Nf,0],'k-');\n",
    "plt.plot(fpos,Sz[0:Nf,0],'ko');\n",
    "\n",
    "plt.show();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 90,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 17640 characters]\n",
      "text/plain": [
       "<Figure size 864x432 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 23812 characters]\n",
      "text/plain": [
       "<Figure size 864x432 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 21852 characters]\n",
      "text/plain": [
       "<Figure size 864x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# eda09_14\n",
    "# example of Chebychev bandpass filter\n",
    "\n",
    "# limits of bandpass filter\n",
    "flow=5;\n",
    "fhigh=10;\n",
    "\n",
    "# standard set up of time\n",
    "N=1024;            # length of time series, must be even\n",
    "Dt=0.01 ;          # sampling interval\n",
    "tmin = 0;          # start time\n",
    "tmax = Dt*(N-1);   # end time\n",
    "t=eda_cvec(np.linspace(tmin,tmax,N)); # time vector\n",
    "\n",
    "# standard set up of frequency\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 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 frequency vector\n",
    "Nw=Nf;             # number of non-negative angular frequencies\n",
    "fpos=eda_cvec(Df * np.linspace(0,No2,Nf) );  # vector non-negative frequencies\n",
    "wpos=2*pi*fpos;    # vector non-negative angular frequencies\n",
    "\n",
    "# input timeseries is a spike\n",
    "din=np.zeros((N,1));\n",
    "din[100,0]=1;\n",
    "\n",
    "# filter\n",
    "[dout, u, v] = chebyshevfilt( din, Dt, flow, fhigh);\n",
    "\n",
    "#  spectrum of input\n",
    "dint = eda_cvec( np.fft.fft(din, axis=0) );  # fourier transform\n",
    "dintpos = dint[0:Nf,0]; # keep only non-negative frequencies\n",
    "dinspec = np.abs(dintpos); # amplitude spectral density\n",
    "\n",
    "#  spectrum of output\n",
    "doutt = eda_cvec( np.fft.fft(dout, axis=0) ); # fourier transform\n",
    "douttpos = eda_cvec( doutt[0:Nf,0] ); # keep only non-negative frequencies\n",
    "doutspec = np.abs(douttpos); # amplitude spectral density\n",
    "\n",
    "# plot unfiltered and filtered timeseries\n",
    "fig1 = plt.figure(1,figsize=(12,6));\n",
    "ax1 = plt.subplot(2,1,1);\n",
    "plt.axis( [0.75, 1.25, -1, 1 ]);\n",
    "plt.plot(t, din, 'k-' );\n",
    "plt.xlabel('time (s)');\n",
    "plt.ylabel('din(t)');\n",
    "plt.title('timeseries');\n",
    "ax2 = plt.subplot(2,1,2);\n",
    "plt.axis( [0.75, 1.25, -1, 1 ]);\n",
    "plt.plot(t, dout, 'k-' );\n",
    "plt.xlabel('time (s)');\n",
    "plt.ylabel('din(t)');\n",
    "plt.show();\n",
    "\n",
    "# plot spectra\n",
    "fig2 = plt.figure(2,figsize=(12,6));\n",
    "ax1 = plt.subplot(1,2,1);\n",
    "specmax=np.max(dinspec);\n",
    "plt.axis( [0, fmax, 0, 1.1*specmax ]);\n",
    "plt.plot(fpos, dinspec, 'k-' );\n",
    "plt.xlabel('frequency (Hx)');\n",
    "plt.ylabel('din(f)');\n",
    "plt.title('input spectra');\n",
    "ax2 = plt.subplot(1,2,2);\n",
    "plt.axis( [0, fmax, 0, 1.1*specmax ]);\n",
    "plt.plot(fpos, doutspec, 'k-' );\n",
    "plt.xlabel('frequency (Hz)');\n",
    "plt.ylabel('dout(f)');\n",
    "plt.title('output spectra');\n",
    "plt.plot( [flow, flow], [0, 0.2*specmax], 'k:' );\n",
    "plt.plot( [fhigh, fhigh], [0, 0.2*specmax], 'k:' );\n",
    "plt.show();\n",
    "\n",
    "# roots of filters\n",
    "[Nu, i] = np.shape(u);           # length of u\n",
    "Nru = Nu-1;                      # number or roots\n",
    "pol = (np.flipud(u)).ravel();    # polynomial coefficents\n",
    "ru = eda_cvec( np.roots(pol) );  # vector of roots\n",
    "\n",
    "[Nv,i] = np.shape(v);           # length of u\n",
    "Nrv = Nv-1;                     # number or roots\n",
    "pol = (np.flipud(v)).ravel();   # polynomial coefficents\n",
    "rv = eda_cvec( np.roots(pol) ); # vector of roots\n",
    "\n",
    "# complex z-plane\n",
    "Nz=128;     # z-plane image is Nz by Nz\n",
    "zmin=-2.0;  # min z\n",
    "zmax=2.0;   # max z\n",
    "dz=(zmax-zmin)/(Nz-1); # sampling in z\n",
    "zr=eda_cvec(np.linspace(zmin,zmax,Nz)); # real-z vector\n",
    "zi=eda_cvec(np.linspace(zmin,zmax,Nz)); # imag-z vertor\n",
    "\n",
    "# evaluate filter polynomials everywhere in complex z-plane\n",
    "S = np.zeros((Nz,Nz));\n",
    "for n in range(Nz):\n",
    "    for m in range(Nz):\n",
    "        z = zr[n,0]+complex(0.0,1.0)*zi[m,0];\n",
    "        uu=0.0;\n",
    "        zpow=1;\n",
    "        for k in range(Nu):\n",
    "            uu = uu + u[k,0]*zpow;\n",
    "            zpow=zpow*z;\n",
    "        vv=0.0;\n",
    "        zpow=1;\n",
    "        for k in range(Nv):\n",
    "            vv = vv + v[k,0]*zpow;\n",
    "            zpow=zpow*z;\n",
    "        SS = uu/vv;\n",
    "        S[m,n]=np.abs(SS);\n",
    "\n",
    "# plot amplitude of polynomial\n",
    "fig3 = plt.figure(3,figsize=(12,6));\n",
    "\n",
    "# black and white colormap\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",
    "\n",
    "ax1 = plt.subplot(1,2,1);\n",
    "plt.axis( [zmin,zmax,zmin,zmax]);\n",
    "\n",
    "plt.imshow( S, cmap=bwcmap, vmin=np.min(S), vmax=np.max(S), extent=(zmin,zmax,zmin,zmax) );\n",
    "plt.xlabel('real z');\n",
    "plt.ylabel('imag z');\n",
    "plt.title('|S|^2');\n",
    "\n",
    "# plot Fourier transform points on the unit circle\n",
    "plt.plot(0,0,'wo');\n",
    "z=np.exp(-complex(0,1)*w*Dt);\n",
    "plt.plot(np.real(z),np.imag(z),'k-');\n",
    "\n",
    "# plot origin\n",
    "plt.plot(0,0,'ko');\n",
    "# plot unit circle\n",
    "z=np.exp(-complex(0.0,1.0)*w*Dt);\n",
    "plt.plot(np.real(z),np.imag(z),'k-');\n",
    "\n",
    "# plot zeros and poles\n",
    "for k in range(Nru):\n",
    "    plt.plot(np.real(ru[k,0]),np.imag(ru[k,0]),'k*');\n",
    "\n",
    "for k in range(Nrv):\n",
    "    plt.plot(np.real(rv[k,0]),np.imag(rv[k,0]),'k+');\n",
    "\n",
    "plt.show();\n",
    "        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 93,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 107164 characters]\n",
      "text/plain": [
       "<Figure size 864x432 with 6 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# eda09_15\n",
    "# read Reynolds Channel Water Quality dataset and plot it\n",
    "# columns are days, precip, air_temp, water_temp, salinity, turbidity, chlorophyl\n",
    "\n",
    "# load data from file\n",
    "D = np.genfromtxt('reynolds_interpolated.txt', delimiter='\\t')\n",
    "[N, M]=D.shape;\n",
    "\n",
    "# pull out time vector\n",
    "t = eda_cvec( D[:,0] )  # time in days \n",
    "Dt=t[1,0]-t[0,0];       # sampling interval\n",
    "tmin = t[0,0];          # start time\n",
    "tmax = t[N-1,0];        # end time\n",
    "\n",
    "# pull out data vectors\n",
    "p =  eda_cvec(D[:,1]);       # precipitation, inches\n",
    "at = eda_cvec(D[:,2]);       # air temperature, deg C\n",
    "wt = eda_cvec(D[:,3]);       # water temperature, deg C\n",
    "s =  eda_cvec(D[:,4]);       # salinity, practical salinity units\n",
    "tb = eda_cvec(D[:,5]);       # turbidity, (no units given)\n",
    "c =  eda_cvec(D[:,6]);       # chlorophyl, micrograms per liter\n",
    "\n",
    "fig1 = plt.figure(1,figsize=(12,6));\n",
    "\n",
    "# plotting limits\n",
    "itmin=0;\n",
    "itmax=N;\n",
    "tmin=t[itmin];\n",
    "tmax=t[itmax-1];\n",
    "\n",
    "ax1 = plt.subplot(6,1,1);\n",
    "plt.axis( [tmin, tmax, np.min(p[itmin:itmax,0]), np.max(p[itmin:itmax,0])]);\n",
    "plt.plot(t, p, 'k-' );\n",
    "plt.ylabel('precipitation');\n",
    "plt.title('unfiltered data');\n",
    "\n",
    "ax2 = plt.subplot(6,1,2);\n",
    "plt.axis( [tmin, tmax, np.min(at[itmin:itmax,0]), np.max(at[itmin:itmax,0])]);\n",
    "plt.plot(t, at, 'k-' );\n",
    "plt.ylabel('T-air');\n",
    "      \n",
    "ax3 = plt.subplot(6,1,3);\n",
    "plt.axis( [tmin, tmax, np.min(wt[itmin:itmax,0]), np.max(wt[itmin:itmax,0])]);\n",
    "plt.plot(t, wt, 'k-' );\n",
    "plt.ylabel('T-water');\n",
    "\n",
    "ax4 = plt.subplot(6,1,4);\n",
    "plt.axis( [tmin, tmax, np.min(s[itmin:itmax,0]), np.max(s[itmin:itmax,0])]);\n",
    "plt.plot(t, s, 'k-' );\n",
    "plt.ylabel('salinity');\n",
    "\n",
    "ax5 = plt.subplot(6,1,5);\n",
    "plt.axis( [tmin, tmax, np.min(tb[itmin:itmax,0]), np.max(tb[itmin:itmax,0])]);\n",
    "plt.plot(t, tb, 'k-' );\n",
    "plt.ylabel('turbidity');\n",
    "\n",
    "ax6 = plt.subplot(6,1,6);\n",
    "plt.axis( [tmin, tmax, np.min(c[itmin:itmax,0]), np.max(c[itmin:itmax,0])]);\n",
    "plt.plot(t, c, 'k-' );\n",
    "plt.ylabel('chlorophyl');\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 99,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 85400 characters]\n",
      "text/plain": [
       "<Figure size 864x432 with 6 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 128096 characters]\n",
      "text/plain": [
       "<Figure size 864x432 with 6 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# eda09_16\n",
    "# read Reynolds Channel Water Quality dataset and plot it\n",
    "# columns are days, precip, air_temp, water_temp, salinity, turbidity, chlorophyl\n",
    "\n",
    "# load data from file\n",
    "D = np.genfromtxt('reynolds_interpolated.txt', delimiter='\\t')\n",
    "[N, M]=D.shape;\n",
    "\n",
    "# pull out time vector\n",
    "t = eda_cvec( D[:,0] ); # time in days \n",
    "Dt=t[1,0]-t[0,0];       # sampling interval\n",
    "tmin = t[0,0];          # start time\n",
    "tmax = t[N-1,0];        # end time\n",
    "\n",
    "# low-ish bandpass filter parameters\n",
    "fcenter = 1/365.25;\n",
    "flow = 0.75*fcenter;\n",
    "fhigh =1.25*fcenter;\n",
    "\n",
    "# bandpass filter\n",
    "d = np.zeros((N,1));\n",
    "Df = np.zeros((N,M));\n",
    "Df[:,0] = D[:,0];\n",
    "for i in range(1,M):\n",
    "    d[:,0] = D[:,i];\n",
    "    [z, u, v] = chebyshevfilt(d-np.mean(d), Dt, flow, fhigh);\n",
    "    Df[:,i] = z.ravel();\n",
    "\n",
    "fig1 = plt.figure(1,figsize=(12,6));\n",
    "\n",
    "# pull out data vectors\n",
    "p =  eda_cvec(Df[:,1]);       # precipitation, inches\n",
    "at = eda_cvec(Df[:,2]);       # air temperature, deg C\n",
    "wt = eda_cvec(Df[:,3]);       # water temperature, deg C\n",
    "s =  eda_cvec(Df[:,4]);       # salinity, practical salinity units\n",
    "tb = eda_cvec(Df[:,5]);       # turbidity, (no units given)\n",
    "c =  eda_cvec(Df[:,6]);       # chlorophyl, micrograms per liter\n",
    "\n",
    "# plotting limits\n",
    "itmin=0;\n",
    "itmax=N;\n",
    "tmin=t[itmin];\n",
    "tmax=t[itmax-1];\n",
    "\n",
    "ax1 = plt.subplot(6,1,1);\n",
    "plt.axis( [tmin, tmax, np.min(p[itmin:itmax,0]), np.max(p[itmin:itmax,0])]);\n",
    "plt.plot(t, p, 'k-' );\n",
    "plt.ylabel('precipitation');\n",
    "plt.title(\"data filtered bwetween periods of %.1f and %.1f days\" % (1/fhigh,1/flow) );\n",
    "\n",
    "ax2 = plt.subplot(6,1,2);\n",
    "plt.axis( [tmin, tmax, np.min(at[itmin:itmax,0]), np.max(at[itmin:itmax,0])]);\n",
    "plt.plot(t, at, 'k-' );\n",
    "plt.ylabel('T-air');\n",
    "      \n",
    "ax3 = plt.subplot(6,1,3);\n",
    "plt.axis( [tmin, tmax, np.min(wt[itmin:itmax,0]), np.max(wt[itmin:itmax,0])]);\n",
    "plt.plot(t, wt, 'k-' );\n",
    "plt.ylabel('T-water');\n",
    "\n",
    "ax4 = plt.subplot(6,1,4);\n",
    "plt.axis( [tmin, tmax, np.min(s[itmin:itmax,0]), np.max(s[itmin:itmax,0])]);\n",
    "plt.plot(t, s, 'k-' );\n",
    "plt.ylabel('salinity');\n",
    "\n",
    "ax5 = plt.subplot(6,1,5);\n",
    "plt.axis( [tmin, tmax, np.min(tb[itmin:itmax,0]), np.max(tb[itmin:itmax,0])]);\n",
    "plt.plot(t, tb, 'k-' );\n",
    "plt.ylabel('turbidity');\n",
    "\n",
    "ax6 = plt.subplot(6,1,6);\n",
    "plt.axis( [tmin, tmax, np.min(c[itmin:itmax,0]), np.max(c[itmin:itmax,0])]);\n",
    "plt.plot(t, c, 'k-' );\n",
    "plt.ylabel('chlorophyl');\n",
    "\n",
    "plt.show();\n",
    "\n",
    "# high-ish bandpass filter parameters\n",
    "fcenter = 1.0/5.0;\n",
    "flow = 0.75*fcenter;\n",
    "fhigh =1.25*fcenter;\n",
    "\n",
    "# bandpass filter\n",
    "d = np.zeros((N,1));\n",
    "Df = np.zeros((N,M));\n",
    "Df[:,0] = D[:,0];\n",
    "for i in range(1,M):\n",
    "    d[:,0] = D[:,i];\n",
    "    [z, u, v] = chebyshevfilt(d-np.mean(d), Dt, flow, fhigh);\n",
    "    Df[:,i] = z.ravel();\n",
    "\n",
    "fig2 = plt.figure(2,figsize=(12,6));\n",
    "\n",
    "# pull out data vectors\n",
    "p =  eda_cvec(Df[:,1]);       # precipitation, inches\n",
    "at = eda_cvec(Df[:,2]);       # air temperature, deg C\n",
    "wt = eda_cvec(Df[:,3]);       # water temperature, deg C\n",
    "s =  eda_cvec(Df[:,4]);       # salinity, practical salinity units\n",
    "tb = eda_cvec(Df[:,5]);       # turbidity, (no units given)\n",
    "c =  eda_cvec(Df[:,6]);       # chlorophyl, micrograms per liter\n",
    "\n",
    "# plotting limits\n",
    "itmin=0;\n",
    "itmax=N;\n",
    "tmin=t[itmin];\n",
    "tmax=t[itmax-1];\n",
    "\n",
    "ax1 = plt.subplot(6,1,1);\n",
    "plt.axis( [tmin, tmax, np.min(p[itmin:itmax,0]), np.max(p[itmin:itmax,0])]);\n",
    "plt.plot(t, p, 'k-' );\n",
    "plt.ylabel('precipitation');\n",
    "plt.title(\"data filtered bwetween periods of %.1f and %.1f days\" % (1/fhigh,1/flow) );\n",
    "\n",
    "ax2 = plt.subplot(6,1,2);\n",
    "plt.axis( [tmin, tmax, np.min(at[itmin:itmax,0]), np.max(at[itmin:itmax,0])]);\n",
    "plt.plot(t, at, 'k-' );\n",
    "plt.ylabel('T-air');\n",
    "      \n",
    "ax3 = plt.subplot(6,1,3);\n",
    "plt.axis( [tmin, tmax, np.min(wt[itmin:itmax,0]), np.max(wt[itmin:itmax,0])]);\n",
    "plt.plot(t, wt, 'k-' );\n",
    "plt.ylabel('T-water');\n",
    "\n",
    "ax4 = plt.subplot(6,1,4);\n",
    "plt.axis( [tmin, tmax, np.min(s[itmin:itmax,0]), np.max(s[itmin:itmax,0])]);\n",
    "plt.plot(t, s, 'k-' );\n",
    "plt.ylabel('salinity');\n",
    "\n",
    "ax5 = plt.subplot(6,1,5);\n",
    "plt.axis( [tmin, tmax, np.min(tb[itmin:itmax,0]), np.max(tb[itmin:itmax,0])]);\n",
    "plt.plot(t, tb, 'k-' );\n",
    "plt.ylabel('turbidity');\n",
    "\n",
    "ax6 = plt.subplot(6,1,6);\n",
    "plt.axis( [tmin, tmax, np.min(c[itmin:itmax,0]), np.max(c[itmin:itmax,0])]);\n",
    "plt.plot(t, c, 'k-' );\n",
    "plt.ylabel('chlorophyl');\n",
    "\n",
    "plt.show();\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 104,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 76252 characters]\n",
      "text/plain": [
       "<Figure size 864x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Calculated directly from the band-passed time series:\n",
      "    Power in u 187.7300\n",
      "    Power in v 115.6218\n",
      "    Zero-lag cross correlation 77.9754\n",
      "    Curly C squared 0.2801\n",
      " \n",
      "Calculated from the cross spectral density\n",
      "    Pover in u 187.7300\n",
      "    Pover in v 115.6218\n",
      "    Zero-lag cross correlation 77.9754\n",
      "    Curly C squared 0.2801\n",
      "    Coherence 0.2861\n"
     ]
    }
   ],
   "source": [
    "# eda09_17\n",
    "# cohernce calculation (single frequency band only)\n",
    "# as a check of the formula for C\n",
    "\n",
    "# make a time vector\n",
    "N=256;  # length of time series\n",
    "Dt=1.0; # sampling interval\n",
    "t = eda_cvec( Dt*np.linspace(0,N-1,N) ); # time vector\n",
    "\n",
    "# random time series\n",
    "rts1=np.random.normal( 0.0, 1.0, (256, 1) );\n",
    "rts2=np.random.normal( 0.0, 1.0, (256, 1) );\n",
    "rts3=np.random.normal( 0.0, 1.0, (256, 1) );\n",
    "\n",
    "# create partially correlated time series, by adding\n",
    "# together uncorrelated ones\n",
    "uorig = rts1 + rts2;\n",
    "vorig = rts1 + rts3;\n",
    "\n",
    "# standard set up of frequency\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 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 frequency vector\n",
    "Nw=Nf;             # number of non-negative angular frequencies\n",
    "fpos=eda_cvec(Df * np.linspace(0,No2,Nf) );  # vector non-negative frequencies\n",
    "wpos=2*pi*fpos;    # vector non-negative angular frequencies\n",
    "\n",
    "# frequency and bandwith for calculation\n",
    "f0 = fmax/2;            # center of frequency [0,nyquist] range\n",
    "lowside = 0.75;         # fraction below nyquist\n",
    "highside =1.25;         # ractiona bve nyquist\n",
    "flow = f0*lowside;      # low side of frequency band\n",
    "fhigh= f0*highside;     # high side of frequency band\n",
    "ilow = floor(flow/Df);    # corresponding index\n",
    "ihigh = floor(fhigh/Df);  # corresponding index\n",
    "\n",
    "# box car filter matching the parameters above\n",
    "bp = np.zeros((N,1));\n",
    "bp[ilow:ihigh,0]=1.0;          # positive frequencies\n",
    "bp[N-ihigh+1:N-ilow+1,0]=1.0;  # negative frequencies\n",
    "\n",
    "# bandpass filter the data, by tahing inverse fourier transform of product fourier transforms\n",
    "ufilt = eda_cvec( np.real( np.fft.ifft( np.multiply( bp, np.fft.fft(uorig,axis=0) ), axis=0 ) ) );\n",
    "vfilt = eda_cvec( np.real( np.fft.ifft( np.multiply( bp, np.fft.fft(vorig,axis=0) ), axis=0 ) ) );\n",
    "\n",
    "# calculation of power and zero-lag cross-correlation\n",
    "Pu = np.sum( np.power(ufilt,2) );\n",
    "Pv = np.sum( np.power(vfilt,2) );\n",
    "x = np.correlate( ufilt.ravel(), vfilt.ravel(), mode='full' );\n",
    "x0 = eda_cvec( x[N-1] );\n",
    "\n",
    "# plot time series\n",
    "fig1 = plt.figure(1,figsize=(12,6));\n",
    "ax1 = plt.subplot(1,1,1);\n",
    "plt.axis( [t[0,0], t[N-1,0], -10, 10 ] );\n",
    "plt.plot(t, ufilt+3, 'k-' );\n",
    "plt.plot(t, vfilt-3, 'k-' );\n",
    "plt.xlabel('time (s)');\n",
    "plt.ylabel('u(t) and v(t)');\n",
    "plt.title(\"pass-band: %.2f to %.2f Hz\" % (flow, fhigh));\n",
    "plt.show();\n",
    "\n",
    "# display results\n",
    "print('Calculated directly from the band-passed time series:');\n",
    "print(\"    Power in u %.4f\" % (Pu) );\n",
    "print(\"    Power in v %.4f\" % (Pv) );\n",
    "print(\"    Zero-lag cross correlation %.4f\" % (x0) );\n",
    "print(\"    Curly C squared %.4f\" % ((x0**2)/(Pu*Pv)) );\n",
    "print(' ');\n",
    "\n",
    "# Fourier transform\n",
    "ut = eda_cvec( np.fft.fft(uorig, axis=0) ); # fourier trasnform\n",
    "u =  eda_cvec( ut[0:Nf,0] ); # delete negative frequencies\n",
    "\n",
    "vt = eda_cvec( np.fft.fft(vorig,axis=0) ); # fourier trasnform\n",
    "v =  eda_cvec(vt[0:Nf,0]); # delete negative frequencies\n",
    "\n",
    "# spectral densities\n",
    "usv = np.multiply( np.conj(u), v ); # cross spectral density of u and v\n",
    "usu = np.real( np.multiply( np.conj(u), u ) ); # power spectral density of u\n",
    "vsv = np.real( np.multiply( np.conj(v), v ) ); # power spectral density of v\n",
    "\n",
    "# band-averaged quantities\n",
    "usva = np.sum( (usv[ilow:ihigh,0]) ) / (Nf-1);\n",
    "usvar = np.sum( np.real( (usv[ilow:ihigh,0] ) ) )/ (Nf-1);\n",
    "usua = np.sum( usu[ilow:ihigh,0] )/ (Nf-1);\n",
    "vsva = np.sum( vsv[ilow:ihigh,0] )/ (Nf-1);\n",
    "curlyc2 = np.multiply( np.conj(usvar), usvar) / (usua * vsva);\n",
    "coherence = np.real(np.multiply(np.conj(usva),usva)) / (usua * vsva);\n",
    "\n",
    "# display quantities\n",
    "print('Calculated from the cross spectral density');\n",
    "print(\"    Pover in u %.4f\" % (usua));\n",
    "print(\"    Pover in v %.4f\" % (vsva));\n",
    "print(\"    Zero-lag cross correlation %.4f\" % (usvar));\n",
    "print(\"    Curly C squared %.4f\" % (curlyc2));\n",
    "print(\"    Coherence %.4f\" % (coherence));"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 138948 characters]\n",
      "text/plain": [
       "<Figure size 1152x1152 with 10 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# eda09_18\n",
    "# read Reynolds Channel Water Quality dataset\n",
    "# columns are days, precip, air_temp, water_temp, salinity, turbidity, chlorophyl\n",
    "# compute coherence and plot it\n",
    "\n",
    "# set DOCOHEREENCE to 1 to add results from the scipy.signal.coherence() function\n",
    "# (which are plotted in red).  Set to 0 to omit them.\n",
    "DOCOHEREENCE = 0;\n",
    "\n",
    "# this is an example of tupple used here to hold the names of the variables\n",
    "names = ('time', 'precipitation', 'air-temp', 'water-temp', 'salinity', 'turbidity', 'chlorophyll');\n",
    "\n",
    "# load data\n",
    "Draw = np.genfromtxt('reynolds_interpolated.txt', delimiter='\\t')\n",
    "[N, M] = Draw.shape;\n",
    "N = 2*floor(N/2); # truncate to enev length\n",
    "D = np.zeros((N,M));\n",
    "D[:,:] = Draw[0:N,:];\n",
    "\n",
    "# time in days \n",
    "t = np.zeros((N,1));\n",
    "t[:,0] = D[:,0];\n",
    "Dt=t[1,0]-t[0,0];\n",
    "tmin = t[0,0];\n",
    "tmax = t[N-1,0];\n",
    "\n",
    "# standard frequency setup\n",
    "Nf = floor(N/2+1);\n",
    "fny = 1/(2*Dt);\n",
    "Df = fny/(N/2);\n",
    "f = np.zeros((Nf,1)); # positive frequencies\n",
    "f[:,0] = np.linspace(0.0,fny,Nf);\n",
    "\n",
    "# bandwidth factors\n",
    "lowside = 0.75;\n",
    "highside =1.25;\n",
    "\n",
    "ifig=1;\n",
    "fig1 = plt.figure(1,figsize=(16,16));\n",
    "\n",
    "# loop over all pairs of columns, computing coherence for each\n",
    "for icol in range(1,M):\n",
    "    for jcol in range(icol+1,M-1):\n",
    "        # compute cross spectral density and power spectral density\n",
    "        # no need to normalize, since all normalizations cancel\n",
    "        ut = eda_cvec( np.fft.fft( D[:,icol] ) );\n",
    "        vt = eda_cvec( np.fft.fft( D[:,jcol] ) );\n",
    "        u = np.zeros((Nf,1), dtype=np.complex );\n",
    "        v = np.zeros((Nf,1), dtype=np.complex );\n",
    "        u = eda_cvec( ut[0:Nf,0] ); # delete negative frequencies\n",
    "        v = eda_cvec( vt[0:Nf,0] );\n",
    "        usv = np.multiply( np.conj(u) , v ); # cross spectral density of u and v\n",
    "        usu = np.real( np.multiply( np.conj(u) , u ) ); # power spectral density of u\n",
    "        vsv = np.real( np.multiply( np.conj(v) , v ) ); # power spectral density of v\n",
    "        # average over band\n",
    "        usva = np.zeros((Nf,1), dtype=np.complex );\n",
    "        usua = np.zeros((Nf,1));\n",
    "        vsva = np.zeros((Nf,1)); \n",
    "        c = np.zeros((Nf,1)); \n",
    "        for i in range(Nf):\n",
    "            fi = i*Df; # center frequency\n",
    "            flow = lowside*fi;\n",
    "            fhigh = highside*fi;\n",
    "            ilow = floor(flow/Df);\n",
    "            ihigh = floor(fhigh/Df);\n",
    "            if( ihigh==ilow):\n",
    "                ihigh=ilow+1;\n",
    "            if (ilow < 0):\n",
    "                ilow=0;\n",
    "            elif (ilow > Nf):\n",
    "                ilow=Nf-1;\n",
    "            if (ihigh < 0):\n",
    "                ihigh=0;\n",
    "            elif (ihigh > Nf):\n",
    "                ihigh=Nf;\n",
    "            tmp1 = np.mean( usv[ilow:ihigh,0] );\n",
    "            tmp2 = np.mean( usu[ilow:ihigh,0] );\n",
    "            tmp3 = np.mean( vsv[ilow:ihigh,0] );\n",
    "            usva[i,0] = tmp1;\n",
    "            usua[i,0] = tmp2;\n",
    "            vsva[i,0] = tmp3;\n",
    "            \n",
    "            top = np.real( np.conj(tmp1) * tmp1 );\n",
    "            bottom = (tmp2 * tmp3);\n",
    "            if( bottom > 0.0 ):\n",
    "                c[i,0] =  top / bottom;\n",
    "            else:\n",
    "                c[i,0] = 0.0;\n",
    "        # plot coherence\n",
    "        ax1 = plt.subplot(4,3,ifig);\n",
    "        plt.axis( [0, 1/5, 0, 1.1] );\n",
    "        plt.plot(f, c, 'k-' );\n",
    "        plt.ylabel(\"coherent of %s vs %s\" % (names[icol],names[jcol]));\n",
    "        plt.xlabel('frequency, cycles per day');\n",
    "        if( DOCOHEREENCE == 1):\n",
    "            # result using mscohere() (plotted in red)\n",
    "            # the black and red traced will not agree\n",
    "            # perfectly, because the bandwidth parameters\n",
    "            # cannot be made to match, expecially at low frequencies.\n",
    "            N2 = 2;\n",
    "            while (N2 < N ): # find smallest power of 2 greater than N\n",
    "                N2 = 2*N2;\n",
    "            N3 = floor(N2/16);\n",
    "            N4 = 3;\n",
    "            [f2, C2] = sg.coherence( D[:,icol], D[:,jcol], nfft=N2, window='hamming', nperseg=N3, noverlap=N3-N4 );\n",
    "            f2max = np.max(f2);\n",
    "            plt.plot( (fny/f2max)*f2, C2, 'r:');\n",
    "        ifig=ifig+1;\n",
    "\n",
    "plt.show();\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 109784 characters]\n",
      "text/plain": [
       "<Figure size 1152x1152 with 6 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# eda09_19\n",
    "# demonsrate effect of windowing on power spectral density\n",
    "# boxcar window\n",
    "\n",
    "# standard time/frequency definitions\n",
    "N=512;  # length of time serires, must be even\n",
    "Dt=1.0; # sampling interval\n",
    "t = eda_cvec( Dt*np.linspace(0,N-1,N) ); # time vector\n",
    "tmin = t[0,0]; # start time\n",
    "tmax=t[N-1,0]; # end time\n",
    "\n",
    "# standard set up of frequency\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 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 frequency vector\n",
    "Nw=Nf;             # number of non-negative angular frequencies\n",
    "fpos=eda_cvec(Df * np.linspace(0,No2,Nf) );  # vector non-negative frequencies\n",
    "wpos=2*pi*fpos;    # vector non-negative angular frequencies\n",
    "\n",
    "# boxcar window function, W\n",
    "W = np.zeros((N,1));\n",
    "No8 = floor(N/8);\n",
    "W[0:No8,0]=1.0;\n",
    "Wtmax = tmax/8.0;\n",
    "\n",
    "# timeseriee\n",
    "w0 = 40.0*Dw;\n",
    "d = np.sin( w0 * t );\n",
    "\n",
    "# windowed timeseries\n",
    "Wd = np.multiply( W, d );\n",
    "\n",
    "# fourier transform\n",
    "dfft = eda_cvec( np.fft.fft(d, axis=0) );\n",
    "Wfft = eda_cvec( np.fft.fft(W, axis=0) );\n",
    "Wdfft = eda_cvec( np.fft.fft(Wd, axis=0) );\n",
    "\n",
    "# keep non-negative frequencies, only\n",
    "dfftpos = eda_cvec(dfft[0:Nf,0]);\n",
    "Wfftpos = eda_cvec(Wfft[0:Nf,0]);\n",
    "Wdfftpos = eda_cvec(Wdfft[0:Nf,0]);\n",
    "\n",
    "# power spectral density\n",
    "norm = (Dt**2)/tmax;\n",
    "ds=norm*np.power(np.abs(dfftpos),2);\n",
    "Ws=norm*np.power(np.abs(Wfftpos),2);\n",
    "Wds=norm*np.power(np.abs(Wdfftpos),2);\n",
    "\n",
    "# time domain plot\n",
    "fig1 = plt.figure(1,figsize=(16,16));\n",
    "                       \n",
    "ax1 = plt.subplot(3,2,1);\n",
    "plt.axis( [0, tmax, -1.1, 1.1] );\n",
    "plt.plot(t, d, 'k-' );\n",
    "plt.xlabel('time t (s)');                 \n",
    "plt.ylabel('d(t)');\n",
    "                       \n",
    "ax3 = plt.subplot(3,2,3);\n",
    "plt.axis( [0, tmax, -1.1, 1.1] );\n",
    "plt.plot(t, W, 'k-' );\n",
    "plt.xlabel('time t (s)');                 \n",
    "plt.ylabel('W(t)');\n",
    "                       \n",
    "ax5 = plt.subplot(3,2,5);\n",
    "plt.axis( [0, tmax, -1.1, 1.1] );\n",
    "plt.plot(t, Wd, 'k-' );\n",
    "plt.xlabel('time t (s)');                 \n",
    "plt.ylabel('W(t)*d(t)');    \n",
    "            \n",
    "# frequency domain plot\n",
    "ax2 = plt.subplot(3,2,2);\n",
    "plt.axis( [0, fmax, 0, 1.1*np.max(np.sqrt(ds))] );\n",
    "plt.plot(fpos, np.sqrt(ds), 'k-' );\n",
    "plt.xlabel('frequency f (Hz)');                 \n",
    "plt.ylabel('ASD of d');\n",
    "                       \n",
    "ax4 = plt.subplot(3,2,4);\n",
    "plt.axis( [0, fmax, 0, 1.1*np.max(np.sqrt(Ws))] );\n",
    "plt.plot(fpos, np.sqrt(Ws), 'k-' );\n",
    "plt.xlabel('frequency f (Hz)');                 \n",
    "plt.ylabel('ASD of W');\n",
    "                       \n",
    "ax5 = plt.subplot(3,2,6);\n",
    "plt.axis( [0, fmax, 0, 1.1*np.max(np.sqrt(Wds))] );\n",
    "plt.plot( fpos, np.sqrt(Wds), 'k-' );\n",
    "plt.xlabel('frequency f (Hz))');                 \n",
    "plt.ylabel('ASD of Wd');                     \n",
    "                       \n",
    "plt.show();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 100752 characters]\n",
      "text/plain": [
       "<Figure size 1152x1152 with 6 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# eda09_20\n",
    "# demonsrate effect of windowing on power spectral density\n",
    "# Hamming window\n",
    "\n",
    "# standard time/frequency definitions\n",
    "N=512;\n",
    "Dt=1.0;\n",
    "t = eda_cvec( Dt*np.linspace(0,N-1,N) ); # time vector\n",
    "tmin = t[0,0]; # start time\n",
    "tmax=t[N-1,0]; # end time\n",
    "\n",
    "# standard set up of frequency\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 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 frequency vector\n",
    "Nw=Nf;             # number of non-negative angular frequencies\n",
    "fpos=eda_cvec(Df * np.linspace(0,No2,Nf) );  # vector non-negative frequencies\n",
    "wpos=2*pi*fpos;    # vector non-negative angular frequencies\n",
    "\n",
    "# Hamming window function, W\n",
    "W = np.zeros((N,1));\n",
    "Nw = floor(N/8);\n",
    "W[0:Nw,0]=0.54-0.46*np.cos(2*pi*np.linspace(0,Nw-1,Nw)/(Nw-1));\n",
    "\n",
    "# timeseriee\n",
    "w0 = 40*Dw;\n",
    "d = np.sin( w0 * t );\n",
    "\n",
    "# windowed timeseries\n",
    "Wd = np.multiply( W, d );\n",
    "\n",
    "# fourier transform\n",
    "dfft = np.zeros((N,1), dtype=np.complex );\n",
    "Wfft = np.zeros((N,1), dtype=np.complex );\n",
    "Wdfft = np.zeros((N,1), dtype=np.complex );\n",
    "dfft[:,0] = np.fft.fft(d.ravel());\n",
    "Wfft[:,0] = np.fft.fft(W.ravel());\n",
    "Wdfft[:,0] = np.fft.fft(Wd.ravel());\n",
    "\n",
    "# power spectral density\n",
    "norm = (Dt**2)/tmax;\n",
    "ds = np.zeros((Nf,1));\n",
    "Ws = np.zeros((Nf,1));\n",
    "Wds = np.zeros((Nf,1));\n",
    "ds[:,0]=norm*np.power(abs(dfft[0:Nf,0]),2);\n",
    "Ws[:,0]=norm*np.power(abs(Wfft[0:Nf,0]),2);\n",
    "Wds[:,0]=norm*np.power(abs(Wdfft[0:Nf,0]),2);\n",
    "\n",
    "# time domain plot\n",
    "fig1 = plt.figure(1,figsize=(16,16));\n",
    "                       \n",
    "ax1 = plt.subplot(3,2,1);\n",
    "plt.axis( [0, tmax, -1.1, 1.1] );\n",
    "plt.plot(t, d, 'k-' );\n",
    "plt.xlabel('time t (s)');                 \n",
    "plt.ylabel('d(t)');\n",
    "                       \n",
    "ax3 = plt.subplot(3,2,3);\n",
    "plt.axis( [0, tmax, -1.1, 1.1] );\n",
    "plt.plot(t, W, 'k-' );\n",
    "plt.xlabel('time t (s)');                 \n",
    "plt.ylabel('W(t)');\n",
    "                       \n",
    "ax5 = plt.subplot(3,2,5);\n",
    "plt.axis( [0, tmax, -1.1, 1.1] );\n",
    "plt.plot(t, Wd, 'k-' );\n",
    "plt.xlabel('time t (s)');                 \n",
    "plt.ylabel('W(t)*d(t)');    \n",
    "    \n",
    "# frequency domain plot\n",
    "ax2 = plt.subplot(3,2,2);\n",
    "plt.axis( [0, fmax, 0, 1.1*np.max(np.sqrt(ds))] );\n",
    "plt.plot(fpos, np.sqrt(ds), 'k-' );\n",
    "plt.xlabel('frequency f (Hz)');                 \n",
    "plt.ylabel('ASD of d');\n",
    "                       \n",
    "ax4 = plt.subplot(3,2,4);\n",
    "plt.axis( [0, fmax, 0, 1.1*np.max(np.sqrt(Ws))] );\n",
    "plt.plot(fpos, np.sqrt(Ws), 'k-' );\n",
    "plt.xlabel('frequency f (Hz)');                 \n",
    "plt.ylabel('ASD of W');\n",
    "                       \n",
    "ax5 = plt.subplot(3,2,6);\n",
    "plt.axis( [0, fmax, 0, 1.1*np.max(np.sqrt(Wds))] );\n",
    "plt.plot( fpos, np.sqrt(Wds), 'k-' );\n",
    "plt.xlabel('frequency f (Hz))');                 \n",
    "plt.ylabel('ASD of Wd');                     \n",
    "                       \n",
    "plt.show();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "6 largest eigenvalues:\n",
      "    59: 1.724369\n",
      "    60: 4.536603\n",
      "    61: 6.029180\n",
      "    62: 6.268016\n",
      "    63: 6.282832\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 66768 characters]\n",
      "text/plain": [
       "<Figure size 864x1008 with 4 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 46076 characters]\n",
      "text/plain": [
       "<Figure size 864x1008 with 4 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 64148 characters]\n",
      "text/plain": [
       "<Figure size 864x1008 with 5 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 92568 characters]\n",
      "text/plain": [
       "<Figure size 864x1008 with 6 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# eda09_21\n",
    "# MultiTaper analysis\n",
    "\n",
    "# standard time set up\n",
    "N=64;\n",
    "Dt=1.0;\n",
    "t = eda_cvec( Dt*np.linspace(0,N-1,N) ); # time vector\n",
    "tmin = t[0,0]; # start time\n",
    "tmax=t[N-1,0]; # end time\n",
    "\n",
    "# standard set up of frequency\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 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 frequency vector\n",
    "Nw=Nf;             # number of non-negative angular frequencies\n",
    "fpos=eda_cvec(Df * np.linspace(0,No2,Nf) );  # vector non-negative frequencies\n",
    "wpos=2*pi*fpos;    # vector non-negative angular frequencies\n",
    "\n",
    "# choose bandwidth to be +/- some multiple of Dw\n",
    "f0 = 2*Df;\n",
    "w0 = 2*pi*f0;\n",
    "\n",
    "# matrix, M, of sinc functions\n",
    "M = np.zeros((N,N));\n",
    "for m in range(N):\n",
    "    for n in range(N):\n",
    "        M[m,n] = 2*w0*np.sinc( w0*(m-n)*Dt / pi);\n",
    "\n",
    "# solve algebric eigenvalue problem\n",
    "# eig sorts in ascending order\n",
    "[lam, V] = la.eigh(M);\n",
    "print('6 largest eigenvalues:');\n",
    "for m in range(N-5,N):\n",
    "    print(\"    %d: %f\" % (m, lam[m]) );\n",
    "\n",
    "# four largest eigenvectors\n",
    "w1 = eda_cvec( V[:,N-1] );\n",
    "w2 = eda_cvec( V[:,N-2] );\n",
    "w3 = eda_cvec( V[:,N-3] );\n",
    "w4 = eda_cvec( V[:,N-4] );\n",
    "\n",
    "# power spectral density of window functions\n",
    "s1 = np.zeros((Nf,1));\n",
    "s2 = np.zeros((Nf,1));\n",
    "s3 = np.zeros((Nf,1));\n",
    "s4 = np.zeros((Nf,1));\n",
    "norm = (Dt**2)/tmax;\n",
    "temp1 = norm * np.power(np.abs(np.fft.fft(w1,axis=0)),2);\n",
    "temp2 = norm * np.power(np.abs(np.fft.fft(w2,axis=0)),2);\n",
    "temp3 = norm * np.power(np.abs(np.fft.fft(w3,axis=0)),2);\n",
    "temp4 = norm * np.power(np.abs(np.fft.fft(w4,axis=0)),2);\n",
    "s1[:,0] = temp1[0:Nf,0];\n",
    "s2[:,0] = temp2[0:Nf,0];\n",
    "s3[:,0] = temp3[0:Nf,0];\n",
    "s4[:,0] = temp4[0:Nf,0];\n",
    "\n",
    "# time domain plot\n",
    "fig1 = plt.figure(1,figsize=(12,14));\n",
    "winmin = -0.25;\n",
    "winmax = -winmin;\n",
    "\n",
    "ax1 = plt.subplot(4,1,1);\n",
    "plt.axis( [0, tmax, winmin, winmax] );\n",
    "plt.plot(t, w1, 'k-' );\n",
    "plt.xlabel('time t (s)');                 \n",
    "plt.ylabel('w1(t)');\n",
    "plt.title('tapers');\n",
    "\n",
    "ax2 = plt.subplot(4,1,2);\n",
    "plt.axis( [0, tmax, winmin, winmax] );\n",
    "plt.plot(t, w2, 'k-' );\n",
    "plt.xlabel('time t (s)');                 \n",
    "plt.ylabel('w1(t)');\n",
    "\n",
    "ax3 = plt.subplot(4,1,3);\n",
    "plt.axis( [0, tmax, winmin, winmax] );\n",
    "plt.plot(t, w3, 'k-' );\n",
    "plt.xlabel('time t (s)');                 \n",
    "plt.ylabel('w1(t)');\n",
    "\n",
    "ax4 = plt.subplot(4,1,4);\n",
    "plt.axis( [0, tmax, winmin, winmax] );\n",
    "plt.plot(t, w4, 'k-' );\n",
    "plt.xlabel('time t (s)');                 \n",
    "plt.ylabel('w1(t)');\n",
    "\n",
    "plt.show();\n",
    "\n",
    "# plot amplitude spectral density of window functions\n",
    "fig2 = plt.figure(2,figsize=(12,14));\n",
    "\n",
    "ax1 = plt.subplot(4,1,1);\n",
    "plt.axis( [0, fmax, 0, 1.1*np.max(np.sqrt(s1))] );\n",
    "plt.plot(fpos,s1, 'k-' );\n",
    "plt.plot( [f0, f0], [0, np.max(np.sqrt(s1))], 'k:' );\n",
    "plt.xlabel('frequency f (Hz)');                 \n",
    "plt.ylabel('w1(f)');\n",
    "plt.title('ASD of tapers');\n",
    "\n",
    "ax2 = plt.subplot(4,1,2);\n",
    "plt.axis( [0, fmax, 0, 1.1*np.max(np.sqrt(s2))] );\n",
    "plt.plot(fpos,s2, 'k-' );\n",
    "plt.plot( [f0, f0], [0, np.max(np.sqrt(s3))], 'k:' );\n",
    "plt.xlabel('frequency f (Hz)');                 \n",
    "plt.ylabel('w1(f)');\n",
    "\n",
    "ax3 = plt.subplot(4,1,3);\n",
    "plt.axis( [0, fmax, 0, 1.1*np.max(np.sqrt(s3))] );\n",
    "plt.plot(fpos,s3, 'k-' );\n",
    "plt.plot( [f0, f0], [0, np.max(np.sqrt(s3))], 'k:' );\n",
    "plt.xlabel('frequency f (Hz)');                 \n",
    "plt.ylabel('w1(f)');\n",
    "\n",
    "ax4 = plt.subplot(4,1,4);\n",
    "plt.axis( [0, fmax, 0, 1.1*np.max(np.sqrt(s4))] );\n",
    "plt.plot(fpos,s4, 'k-' );\n",
    "plt.plot( [f0, f0], [0, np.max(np.sqrt(s4))], 'k:' );\n",
    "plt.xlabel('frequency f (Hz)');                 \n",
    "plt.ylabel('w1(f)');\n",
    "\n",
    "plt.show();\n",
    "\n",
    "# new, longer time interval, with same sampling interval\n",
    "N2=512;                 # length of longer interval\n",
    "Dt2 = 1.0;              # sampling interval\n",
    "N2o2 = floor(N2/2);     # interger version of N2/2\n",
    "Nf2=N2o2+1;             # number of non-negative frequencies\n",
    "t2 = eda_cvec( Dt*np.linspace(0,N2-1,N2) ); # time vector\n",
    "tmin2 = t2[0,0];        # start time\n",
    "tmax2 = t2[N2-1,0];     # end time\n",
    "fmax2=1/(2.0*Dt2);      # nyquist freqency\n",
    "Df2=fmax2/N2o2;         # frequency sampling\n",
    "Nf2=N2o2+1;             # number of non-negative frequencies\n",
    "# full (pos and neg) vector of freqeuncies\n",
    "f2=eda_cvec(Df2 * np.concatenate( ( np.linspace(0,N2o2,Nf2), np.linspace(-N2o2+1,-1,N2o2-1) ), axis=0 ));\n",
    "Dw2=2*pi*Df2; # angular frequency sampling;\n",
    "# full (pos and neg) vector of angular freqeuncies\n",
    "w_2 = 2 * pi * f2;\n",
    "Nw2=Nf2; # number of angular frequencies\n",
    "fpos2=eda_cvec(Df2 * np.linspace(0,N2o2,Nf2)); # vector of non-negative frequencies\n",
    "wpos2=2*pi*fpos2; # vector of non-negative angular frequencies\n",
    "\n",
    "# create time series and window it\n",
    "ff= fmax/3.3;\n",
    "ww= 2*pi*ff;\n",
    "d=np.cos(ww*t2);\n",
    "\n",
    "# w0 is a uniform window function\n",
    "w0 = np.ones((N,1));\n",
    "\n",
    "# we is a bounch of zero that will be padded onto the\n",
    "# end of a short time series to make a long one\n",
    "we = np.zeros((N2-N,1));\n",
    "\n",
    "# zero-padded window functions\n",
    "w0we = eda_cvec( np.concatenate( (w0,we), axis=0 ));\n",
    "w1we = eda_cvec( np.concatenate( (w1,we), axis=0 ));\n",
    "w2we = eda_cvec( np.concatenate( (w2,we), axis=0 ));\n",
    "w3we = eda_cvec( np.concatenate( (w3,we), axis=0 ));\n",
    "w4we = eda_cvec( np.concatenate( (w4,we), axis=0 ));\n",
    "\n",
    "# windowed data\n",
    "dw0 = eda_cvec( np.multiply( w0we, d ) );\n",
    "dw1 = eda_cvec( np.multiply( w1we, d ) );\n",
    "dw2 = eda_cvec( np.multiply( w2we, d ) );\n",
    "dw3 = eda_cvec( np.multiply( w3we, d ) );\n",
    "dw4 = eda_cvec( np.multiply( w4we, d ) );\n",
    "\n",
    "# plot windowed time series\n",
    "fig2 = plt.figure(2,figsize=(12,14));\n",
    "\n",
    "ax1 = plt.subplot(5,1,1);\n",
    "plt.axis( [tmin2, tmax2, -1.1, 1.1] );\n",
    "plt.plot(t2, dw0, 'k-' );\n",
    "plt.xlabel('time t (s)');                 \n",
    "plt.ylabel('d(t)');\n",
    "plt.title('windowed time series');\n",
    "\n",
    "ax2 = plt.subplot(5,1,2);\n",
    "plt.axis( [tmin2, tmax2, winmin, winmax] );\n",
    "plt.plot(t2, dw1, 'k-' );\n",
    "plt.xlabel('time t (s)');                 \n",
    "plt.ylabel('d*w1')\n",
    "\n",
    "ax3 = plt.subplot(5,1,3);\n",
    "plt.axis( [tmin2, tmax2, winmin, winmax] );\n",
    "plt.plot(t2, dw2, 'k-' );\n",
    "plt.xlabel('time t (s)');                 \n",
    "plt.ylabel('d*w2')\n",
    "\n",
    "ax4 = plt.subplot(5,1,4);\n",
    "plt.axis( [tmin2, tmax2, winmin, winmax] );\n",
    "plt.plot(t2, dw3, 'k-' );\n",
    "plt.xlabel('time t (s)');                 \n",
    "plt.ylabel('d*w3')\n",
    "\n",
    "ax5 = plt.subplot(5,1,5);\n",
    "plt.axis( [tmin2, tmax2, winmin, winmax] );\n",
    "plt.plot(t2, dw4, 'k-' );\n",
    "plt.xlabel('time t (s)');                 \n",
    "plt.ylabel('d*w4')\n",
    "\n",
    "plt.show();\n",
    "\n",
    "# power spectral density of time series\n",
    "norm = (Dt2**2)/tmax2;\n",
    "temp = norm * np.power( np.abs( eda_cvec(np.fft.fft(d,axis=0))),2);  # power spectrum\n",
    "ds = eda_cvec( temp[0:Nf2,0] ); # delete negative frequencies\n",
    "\n",
    "# power spectral density of windowed time series\n",
    "norm = (Dt2**2)/tmax;\n",
    "temp0 = norm * np.power( np.abs( eda_cvec(np.fft.fft(dw0,axis=0))),2); # fourier transform and etc.\n",
    "temp1 = norm * np.power( np.abs( eda_cvec(np.fft.fft(dw1,axis=0))),2);\n",
    "temp2 = norm * np.power( np.abs( eda_cvec(np.fft.fft(dw2,axis=0))),2);\n",
    "temp3 = norm * np.power( np.abs( eda_cvec(np.fft.fft(dw3,axis=0))),2);\n",
    "temp4 = norm * np.power( np.abs( eda_cvec(np.fft.fft(dw4,axis=0))),2);\n",
    "dw0s = eda_cvec( temp0[0:Nf2,0] ); # keep only non-negative frequencies\n",
    "dw1s = eda_cvec( temp1[0:Nf2,0] );\n",
    "dw2s = eda_cvec( temp2[0:Nf2,0] );\n",
    "dw3s = eda_cvec( temp3[0:Nf2,0] );\n",
    "dw4s = eda_cvec( temp4[0:Nf2,0] );\n",
    "\n",
    "# average power spectral densities\n",
    "dws = (1.0/3.0)*( dw1s + dw2s + dw3s  );\n",
    "\n",
    "# plot amplitude spectral density of window functions\n",
    "fig4 = plt.figure(2,figsize=(12,14));\n",
    "\n",
    "ax1 = plt.subplot(6,1,1);\n",
    "plt.axis( [0, fmax2, 0, 1.1*np.max(np.sqrt(ds))] );\n",
    "plt.plot(fpos2,np.sqrt(ds), 'k-' );\n",
    "plt.plot( [ff, ff], [0, np.max(np.sqrt(ds))], 'k:' );\n",
    "plt.xlabel('frequency f (Hz)');                 \n",
    "plt.ylabel('ASD of d');\n",
    "plt.ylabel('ASD of taperd version of d(t)');\n",
    "\n",
    "ax2 = plt.subplot(6,1,2);\n",
    "plt.axis( [0, fmax2, 0, 1.1*np.max(np.sqrt(dw0s))] );\n",
    "plt.plot(fpos2,np.sqrt(dw0s), 'k-' );\n",
    "plt.plot( [ff, ff], [0, np.max(np.sqrt(dw0s))], 'k:' );\n",
    "plt.xlabel('frequency f (Hz)');                 \n",
    "plt.ylabel('ASD of d*w0');\n",
    "\n",
    "ax3 = plt.subplot(6,1,3);\n",
    "plt.axis( [0, fmax2, 0, 1.1*np.max(np.sqrt(dw1s))] );\n",
    "plt.plot(fpos2,np.sqrt(dw1s), 'k-' );\n",
    "plt.plot( [ff, ff], [0, np.max(np.sqrt(dw1s))], 'k:' );\n",
    "plt.xlabel('frequency f (Hz)');                 \n",
    "plt.ylabel('ASD of d*w1');\n",
    "\n",
    "ax4 = plt.subplot(6,1,4);\n",
    "plt.axis( [0, fmax2, 0, 1.1*np.max(np.sqrt(dw2s))] );\n",
    "plt.plot(fpos2,np.sqrt(dw2s), 'k-' );\n",
    "plt.plot( [ff, ff], [0, np.max(np.sqrt(dw2s))], 'k:' );\n",
    "plt.xlabel('frequency f (Hz)');                 \n",
    "plt.ylabel('ASD of d*w2');\n",
    "\n",
    "ax5 = plt.subplot(6,1,5);\n",
    "plt.axis( [0, fmax2, 0, 1.1*np.max(np.sqrt(dw3s))] );\n",
    "plt.plot(fpos2,np.sqrt(dw3s), 'k-' );\n",
    "plt.plot( [ff, ff], [0, np.max(np.sqrt(dw3s))], 'k:' );\n",
    "plt.xlabel('frequency f (Hz)');                 \n",
    "plt.ylabel('ASD of d*w3');\n",
    "\n",
    "ax6 = plt.subplot(6,1,6);\n",
    "plt.axis( [0, fmax2, 0, 1.1*np.max(np.sqrt(dws))] );\n",
    "plt.plot(fpos2,np.sqrt(dws), 'k-' );\n",
    "plt.plot( [ff, ff], [0, np.max(np.sqrt(dws))], 'k:' );\n",
    "plt.xlabel('frequency f (Hz)');                 \n",
    "plt.ylabel('mean ASD)');\n",
    "\n",
    "plt.show();\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
}
