edapy_09

In [1]:

    # edapy_09_00 clear all variables and import vatious modules
    %reset -f
    import os
    from datetime import date
    from math import exp, pi, sin, cos, tan, sqrt, floor, ceil, log
    import numpy as np
    import scipy.sparse.linalg as las
    from scipy import sparse
    import scipy.linalg as la
    import scipy.signal as sg
    from matplotlib import pyplot as plt
    import matplotlib.cm as cm
    from matplotlib.colors import ListedColormap

    # eda_draw function makes a "pictorial matrix equation"
    # arguments are vectors, matrices and strings
    # which are plotted in the order that the appear
    # except that strings starting with 'title ' are plotted
    # under the subseqeunt matrix or vector
    # always returns a status of 1

    def eda_draw(*argv):
        bw = np.zeros((256,4));
        v = 0.9*(256 - np.linspace( 0, 255, 256 ))/255;
        bw[:,0] = v;
        bw[:,1] = v;
        bw[:,2] = v;
        bw[:,3] = np.ones(256);
        bwcmap = ListedColormap(bw);
        # size of plot
        W = 16;
        H = 4;
        fig1 = plt.figure(1);
        # figsize width and height in inches
        fig1.set_size_inches(W,H);
        ax1 = plt.subplot(1,1,1);
        plt.axis([0, W, -H/2, H/2]);
        plt.axis('off');
        LM = W/6;    # matrix width and heoght
        LV = W/40;   # vector width
        FS = 0.12;    # character width
        TO = 0.4;    # title vertical offset
        SP = 0.2;    # space between objects
        LS = 0.2;    # leading space
        p = LS; # starting x-position
        istitle=0; # flags presence of a title
        for a in argv:
            if isinstance(a,np.ndarray):
                sh = np.shape(a);
                if len(sh) == 1:  # conversion to nx1 array
                    n = sh[0];
                    m = 1;
                    ap = a;
                    a = np.zeros((n,1));
                    a[:,0] = ap;
                else:
                    n = sh[0];
                    m = sh[1];
                if m==1:
                    pold=p;
                    left=p;
                    right=p+LV;
                    bottom=-LM/2;
                    top=LM/2;
                    plt.imshow( a, cmap=bwcmap, vmin=np.min(a), vmax=np.max(a), extent=(left,right,bottom,top) );
                    p = p+LV;
                    pm = (p+pold)/2;
                    if istitle:
                        plt.text(pm,-(LM/2)-TO,titlestr,horizontalalignment='center');
                        istitle=0;
                    p = p+SP;
                else:
                    pold=p;
                    left=p;
                    right=p+LM;
                    bottom=-LM/2;
                    top=LM/2;
                    plt.imshow( a, cmap=bwcmap, vmin=np.min(a), vmax=np.max(a), extent=(left,right,bottom,top) );
                    p = p+LM;
                    pm = (p+pold)/2;
                    if istitle:
                        plt.text(pm,-(LM/2)-TO,titlestr,horizontalalignment='center');
                        istitle=0;
                    p = p+SP;
            elif isinstance(a,str):
                ns = len(a);
                istitle=0;
                if( ns>=6 ):
                    if 'title ' in a[0:6]:
                        istitle=1;
                        titlestr=a[6:];
                if( istitle != 1):
                    plt.text(p,0,a);
                    p = p + ns*FS + SP;
        plt.show();
        return 1;

    # function to make a numpy N-by-1 column vector
    # c=eda_cvec1(v) from entity v that is array-like,
    # including a number, a list of numbers, a tuple
    # of numbers, an N-by-0 np array and a N-by1 np array.
    # It promotes integers to floats, and integers
    # and floats to complex, by context. The function
    # also insures that, if v is an np array, that
    # c is a copy, as contrasted to a view, of it.
    # This version takes just one input argment.
    # whereas c=eda_cvec(v1,v2,...) concatenates
    # many argiments.
    def eda_cvec1(v):
        if isinstance(v, int):
            w = np.zeros((1,1),dtype=int);
            w[0,0] = v;
            return w;
        elif isinstance(v, float):
            w = np.zeros((1,1),dtype=float);
            w[0,0] = v;
            return w;
        elif isinstance(v, complex):
            w = np.zeros((1,1),dtype=complex);
            w[0,0] = v;
            return w;
        elif isinstance(v, np.ndarray):
            s = np.shape(v);
            if len(s) == 1:
                return np.copy(np.reshape(v,(s[0],1)));
            else:
                [r,c]=s;
                if( c==1 ):
                    return(np.copy(v));
                elif(r==1):
                    return(np.copy(v.T));
                else:
                    raise TypeError("eda_cvec: %d by %d ndarray not allowed" % (r, c));
        elif isinstance(v, list):
            r = len(v);
            t = int;
            for vi in v:
                if isinstance(vi,int):
                    pass;
                elif isinstance(vi,float):
                    t=float;
                elif isinstance(vi,complex):
                    t=complex;
                    break;
                else:
                    raise TypeError("eda_cvec: list contains unsupported type %s" % type(vi));
            w = np.zeros((r,1),dtype=t);
            w[:,0] = v;
            return w;
        elif isinstance(v, tuple):
            r = len(v);
            t = int;
            for vi in v:
                if isinstance(vi,int):
                    pass;
                elif isinstance(vi,float):
                    t=float;
                elif isinstance(vi,complex):
                    t=complex;
                    break;
                else:
                    raise TypeError("eda_cvec: tuple contains unsupported type %s" % type(vi));
            w = np.zeros((r,1),dtype=t);
            w[:,0] = v;
            return w;
        else:
            raise TypeError("eda_cvec: %s not supported" % type(v));

    # function to make a numpy N-by-1 column vector
    # c=eda_cvec(v1, v2, ...) from a list of several
    # array-like entities v1, v2, including a number
    # a list of numbers, a tuple of numbers, an N-by-0 np array
    # and a N-by-1 np array. The function
    # also insures that, if v is an np array, that
    # c is a copy, as contrasted to a view, of it
    # It promotes integers to floats, and integers
    # and floats to complex, by context.
    # This version concatenates many argments,
    # whereas c=eda_cvec1(v1) takes just one argiment.
    # I recommend always using eda_cvec(v1, v2, ...)
    def eda_cvec(*argv):
        t = int;
        Nt = 0;
        for a in argv:
            v = eda_cvec1(a);
            N,M = np.shape(v);
            Nt = Nt + N;
            if (t==int) and isinstance(v[0,0],float):
                t=float;
            elif isinstance(v[0,0],complex):
                t=complex;
        w = np.zeros((Nt,1),dtype=t);
        Nt = 0;
        for a in argv:
            v = eda_cvec1(a);
            N,M = np.shape(v);
            w[Nt:Nt+N,0] = v[0:N,0];
            Nt = Nt + N;
        return w;

    def FTFmul(v):
        # this function is used by the bicongugate gradient solver to solve the geneneralized least
        # squares problem Fm=f.  Note that "F" must be called "edaFsparse".
        global edaFsparse;
        s = np.shape(v);
        if(len(s)==1):
            vv = np.zeros((s[0],1));
            vv[:,0] = v;
        else:
            vv=v;
        temp = edaFsparse*vv;
        return edaFsparse.transpose()*temp;

    def GLSFilterMul(v):
        # this function is used by the bicongugate gradient solver to solve the
        # geneneralized least squares problem Fm=f with F=[G;eH] and f=[d,eh],
        # where G is a topplitz matrix built from the filter g
        # Note that "H" must be sparese and called "edaHsparse" and that
        # the filter "g" bust be called edafilterg and must be a column vector
        # the damping e = sigma_d/sigma_h must be called "eda_e"
        global eda_e;
        global edaHsparse;
        global edafilterg;
        # rearrange as column-vector
        s = np.shape(v);
        M = s[0];
        if(len(s)==1):
            vv = np.zeros((s[0],1));
            vv[:,0] = v;
        else:
            vv=v;
        N, i = np.shape(edafilterg);
        # implements (GT G + HT H) v
        temp1 = np.zeros((N+M-1,1));
        temp1[:,0] = np.convolve(edafilterg.ravel(),vv.ravel()); # G v is of length N
        a = np.zeros((N,1));
        a[:,0] = temp1[0:N,0];
        b = eda_e*edaHsparse * v; # H v is of length K
        temp2 = np.zeros((N+N-1,1));
        temp2[:,0] = np.convolve((np.flipud(edafilterg)).ravel(),a.ravel()); # GT (G v) is of length M
        a2 = temp2[N-1:N+M-1,0];
        b2 = eda_e*edaHsparse.transpose()*b; # HT (H v) is of length M
        # FT F v = GT G v + HT H v
        return (a2+b2);

    def GLSFilterMulRHS(g,dobs,e,H,h,M):
        # companion function to GLSFilterMul()
        # creates RHS of generalized least squares equattion
        # set up right hand side, F'f = GT qobs + e HT h
        # note that H must be sparse
        N,i = np.shape(dobs);
        temp = eda_cvec( np.convolve(np.flipud(g).ravel(),dobs.ravel()) );
        FTfa = eda_cvec( temp[N-1:N+M-1,0] );
        FTfb = (e**2)*H.T*h;
        FTf=FTfa+FTfb;
        return (FTf);

    def GLSautoMul(v):
        # this function is used by the bicongugate gradient solver to solve the
        # geneneralized least squares problem Fm=f with F=[G;H] and f=[d,h],
        # where GTG is built from the autocorrelation of the filter g
        global eda_e;
        global edaHsparse;
        global edaautoc;
        # rearrange as column-vector
        s = np.shape(v);
        M = s[0];
        if(len(s)==1):
            vv = np.zeros((s[0],1));
            vv[:,0] = v;
        else:
            vv=v;
        # FT F v = GT G v + (e^2) HT H v
        GTGv=np.zeros((M,1));
        for i in range(M):
            u = np.concatenate(
                (np.flipud(edaautoc[0:i+1,0]), edaautoc[1:M-i,0]),
                axis=0 );
            GTGv[i,0] =  np.matmul( u.T, vv );
        Hv = eda_e*edaHsparse*vv;
        HTHv = eda_e*edaHsparse.transpose()*Hv;
        return (GTGv + HTHv);

    def GLSautoMulRHS(g,d,e,H,h,M):
    # companion function to def GLSautoMul(), makes RHS
        # calculate FT d = GT d + (e^2) HT h
        # when G is a toeplitz matrix
        # built from a filter g, using only the
        # cross-correlation of g and d
        # M is the number of unnowns
        NA,i = np.shape(d);
        NB,i = np.shape(g);
        Nm = NB-1; # index of zero-lag element
        tmp = eda_cvec(np.correlate(d.ravel(),g.ravel(),mode='full'));
        c = tmp[Nm:Nm+M,0:1];    # c = GT d
        HTh =  edaHsparse.T * h; # HT h
        FTf = c + (e**2)*HTh;
        return (FTf);

    def chebyshevfilt(d, Dt, flow, fhigh):
        # chebyshevfilt
        # chebyshev IIR bandpass filter
        # d - input array of data
        # Dt - sampling interval
        # flow - low pass frequency, Hz
        # fhigh - high pass frequency, Hz
        # dout - output array of data
        # u - the numerator filter
        # v - the denominator filter
        # these filters can be used again using dout=filter(u,v,din);

        # make sure input timeseries is a column vector
        s = np.shape(d);
        N = s[0];
        if(N==1):
            dd = np.zeros((N,1));
            dd[:,0] = d;
        else:
            dd=d;

        # sampling rate
        rate=1/Dt;

        # ripple parameter, set to ten percent
        ripple=0.1;

        # normalise frequency
        fl=2.0*flow/rate;
        fh=2.0*fhigh/rate;

        # center frequency
        cf = 4 * tan( (fl*pi/2) ) * tan( (fh*pi/2) );

        # bandwidth
        bw = 2 * ( tan( (fh*pi/2) ) - tan( (fl*pi/2) ) );

        # ripple parameter factor
        rpf = sqrt((sqrt((1.0+1.0/(ripple*ripple))) + 1.0/ripple));
        a = 0.5*(rpf-1.0/rpf);
        b = 0.5*(rpf+1.0/rpf);

        u=np.zeros((5,1));
        v=np.zeros((5,1));
        theta = 3*pi/4;
        sr = a * cos(theta);
        si = b * sin(theta);
        es = sqrt(sr*sr+si*si);
        tmp= 16.0 - 16.0*bw*sr + 8.0*cf + 4.0*es*es*bw*bw - 4.0*bw*cf*sr + cf*cf;
        v[0,0] = 1.0;
        v[1,0] = 4.0*(-16.0 + 8.0*bw*sr - 2.0*bw*cf*sr + cf*cf)/tmp;
        v[2,0] = (96.0 - 16.0*cf - 8.0*es*es*bw*bw + 6.0*cf*cf)/tmp;
        v[3,0] = (-64.0 - 32.0*bw*sr + 8.0*bw*cf*sr + 4.0*cf*cf)/tmp;
        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;
        tmp = 4.0*es*es*bw*bw/tmp;
        u[0,0] = tmp;
        u[1,0] = 0.0;
        u[2,0] = -2.0*tmp;
        u[3,0] = 0.0;
        u[4,0] = tmp;

        dout = sg.lfilter(u.ravel(),v.ravel(),dd.ravel());
        return (dout,u,v);

    def eda_timedelay( uA, uB, Dt ):
    # time delay between two similar timeseries, by
    # cross-correlation followed by parabilic refinment.
    # delay is positive when pulse on B is later than pulse on A
    # input:
    # uA, uB, two signals to be compared, as N-by-1 np.arrays
    #         should have zero mean
    #         if of unequal length, the shorter is zero padded
    # Dt, sampling of both signals
    # output
    # tBmA_est, estimated time delay, positive when dB is delated with respect to uA
    # Cmax_est, cross-correlation at estimated delay

        # cross correlate
        c = eda_cvec( np.correlate(uA.ravel(), uB.ravel(), mode='full') );
        Nc,i = np.shape(c);

        # rough estimate of lag
        NA, i = np.shape(uA);  # length of u
        NB, i = np.shape(uB);  # length of v
        cmax = np.max(c);
        icmax = np.argmax(c);
        tBmArough = -Dt * (icmax-NB+1);

        # parabolic refinement
        # (c-cmax) = m(1) + m(2) dt + m(3) dt^2
        # d(c-cmax)/ddt = 0 = m(2) + 2 m(3) dt
        # dt = -0.5*m(2)/m(3)
        # c = cmax + m(1) + m(2) dt + m(3) dt^2
        Dt2 = Dt**2;
        g1 =     eda_cvec([1.0, 1.0, 1.0]);
        g2 =  Dt*eda_cvec([-1.0, 0.0, 1.0]);
        g3 = Dt2*eda_cvec([1.0, 0.0, 1.0]);
        G = np.concatenate( (g1,g2,g3), axis=1);
        d = eda_cvec( [ c[icmax-1,0]-cmax, c[icmax,0]-cmax, c[icmax+1,0]-cmax ] );
        m = la.solve( np.matmul(G.transpose(),G), np.matmul(G.transpose(),d) );
        dt = -0.5*m[1,0]/m[2,0];
        Cmax_est = cmax + m[0,0] + m[1,0]*dt + m[2,0]*Dt2;
        tBmA_est = tBmArough - dt;

        return tBmA_est, Cmax_est;

In [13]:

    # edapy_09_01
    # estimate covariance and etc of rorck data types

    # read data file
    S = np.genfromtxt('../Data/rocks.txt', delimiter='\t')
    [N, M]=S.shape;
    # 0: SiO2
    # 1: TiO2
    # 2: Al203
    # 3: FeO-total
    # 4: MgO
    # 5: CaO
    # 6: Na2O
    # 7: K2O

    # note use of a tupple to create array of variable names
    name = ( 'SiO2', 'TiO2', 'Al203', 'FeO', 'MgO', 'CaO', 'Na2O', 'K2O' );

    dbar = eda_cvec( np.mean(S,axis=0) );    # mean
    C = np.cov(S,rowvar=False);              # covariance
    sigma = eda_cvec( np.sqrt(np.diag(C)) ); # stddev
    R = np.corrcoef(S,rowvar=False);         # correlation coefficient

    # plot correlation coefficent matrix
    print('correlation coefficients:')
    eda_draw(np.abs(R));

    # display results; also find most correlated data
    print('mean and standard deviations:')
    for i in range(M):
        print( "    %s: %f %f" % (name[i], dbar[i,0], sigma[i,0]) );

    print(' ');
    icor=1;
    jcor=2;
    r=0;
    print('covariance and correlation coefficient:')
    for i in range(M):
        for j in range(i+1,M):
            print( "    %s vs %s: %f %f" % (name[i], name[j], C[i,j], R[i,j])) ;
            if( abs(R[i,j]) > r ):
                r = abs(R[i,j]);
                icor=i;
                jcor=j;

    fig2 = plt.figure(2,figsize=(6,6));
    ax1 = plt.subplot(1,1,1);
    plt.plot(S[:,icor],S[:,jcor],'ko');
    plt.xlabel(name[icor]);
    plt.ylabel(name[jcor]);
    plt.title('most correlated variables');
    plt.show();

    correlation coefficients:

    mean and standard deviations:
        SiO2: 50.388700 1.295814
        TiO2: 1.364637 0.377338
        Al203: 15.255543 1.601403
        FeO: 9.860720 1.322854
        MgO: 7.876269 3.209106
        CaO: 11.641090 1.332980
        Na2O: 2.471665 0.469354
        K2O: 0.190218 0.185841

    covariance and correlation coefficient:
        SiO2 vs TiO2: 0.176776 0.361535
        SiO2 vs Al203: 0.103385 0.049821
        SiO2 vs FeO: 0.427913 0.249633
        SiO2 vs MgO: -1.869543 -0.449582
        SiO2 vs CaO: 0.098530 0.057043
        SiO2 vs Na2O: 0.203303 0.334273
        SiO2 vs K2O: 0.021383 0.088796
        TiO2 vs Al203: 0.013893 0.022992
        TiO2 vs FeO: 0.251822 0.504489
        TiO2 vs MgO: -0.568284 -0.469300
        TiO2 vs CaO: -0.112748 -0.224158
        TiO2 vs Na2O: 0.129429 0.730804
        TiO2 vs K2O: 0.021179 0.302024
        Al203 vs FeO: -0.819501 -0.386845
        Al203 vs MgO: -3.517052 -0.684375
        Al203 vs CaO: 1.407959 0.659577
        Al203 vs Na2O: 0.291717 0.388115
        Al203 vs K2O: 0.024097 0.080970
        FeO vs MgO: -0.959504 -0.226022
        FeO vs CaO: -0.352446 -0.199874
        FeO vs Na2O: 0.051783 0.083402
        FeO vs K2O: -0.019208 -0.078130
        MgO vs CaO: -2.664767 -0.622947
        MgO vs Na2O: -0.791771 -0.525672
        MgO vs K2O: -0.104875 -0.175851
        CaO vs Na2O: -0.056522 -0.090343
        CaO vs K2O: -0.037734 -0.152324
        Na2O vs K2O: 0.019353 0.221870

In [15]:

    # edapy_09_02
    # load the neuse.txt data set
    # make scatter plots of lagged samples
    # in order to investigate correlation

    # red data file
    S = np.genfromtxt('../Data/neuse.txt', delimiter='\t')
    [N, M]=D.shape;

    # break out time and data vectors
    t = eda_cvec( S[:,0] );
    d = eda_cvec( S[:,1] );

    fig1 = plt.figure(1,figsize=(16,4));

    lag=1;
    ax1 = plt.subplot(1,3,1);
    plt.axis([0, 25000, 0, 25000]);
    plt.plot(d[0:N-lag,0], d[lag:N,0],'ko');
    plt.xlabel('discharge');
    plt.ylabel("discharge lagged by %d days" % (lag) );

    lag=3;
    ax2 = plt.subplot(1,3,2);
    plt.axis([0, 25000, 0, 25000]);
    plt.plot(d[0:N-lag,0], d[lag:N,0],'ko');
    plt.xlabel('discharge');
    plt.ylabel("discharge lagged by %d days" % (lag) );

    lag=30;
    ax3 = plt.subplot(1,3,3);
    plt.axis([0, 25000, 0, 25000]);
    plt.plot(d[0:N-lag,0], d[lag:N,0],'ko');
    plt.xlabel('discharge');
    plt.ylabel("discharge lagged by %d days" % (lag) );

    plt.show();

In [41]:

    # edapy_09_03, autocorrelation of Neuse hydrograph

    # read data file
    S = np.genfromtxt('../Data/neuse.txt', delimiter='\t')
    [N, M]=S.shape;

    # break out time and data vectors
    t = eda_cvec( S[:,0] );
    d = eda_cvec( S[:,1] );
    Dt=t[1,0]-t[0,0]; # sampling interval

    # demean the data
    d=d-np.mean(d);

    # compute the autocorrelation
    a = eda_cvec(np.correlate(d.ravel(),d.ravel(),mode='full'));
    # time axis for autocorrelation
    tau = eda_cvec( np.linspace(-N+1,N-1,2*N-1) );
    # number of points correlated for each lag
    Nk = eda_cvec(np.concatenate(
        (np.linspace(1,N,N), np.linspace(N-1,1,N-1))
         ));
    # correction for length of timeseries
    a = N*np.divide( a, Nk );

    fig1 = plt.figure(1,figsize=(16,6));

    lag=1;
    ax1 = plt.subplot(2,1,1);
    plt.axis([ -30, 30, 0, 1.1*np.max(a) ]);
    plt.plot( tau, a, 'k-');
    plt.xlabel('lag (days)');
    plt.ylabel('autocorrelation');

    ax2 = plt.subplot(2,1,2);
    plt.axis([ -365.25*10, 365.25*10, -1.1*np.max(a), 1.1*np.max(a) ]);
    plt.plot( tau, a,'k-');
    plt.xlabel('lag (days)');
    plt.ylabel('autocorrelation');

    plt.show();

In [46]:

    # edapy_09_04, align two timeseries using cross-correlation
    # this example uses synthetic data

    # time series paramaters
    N=1000; # number of data
    Dt=0.1; # sampling interval
    t = eda_cvec( Dt*np.linspace(0,N-1,N) );

    # synthetic time series, u and v
    t1 = Dt*4*N/8 + Dt/3.0;
    t2 = Dt*5*N/8;
    tlagtrue = t2-t1;
    w0 = 2*pi*0.5;
    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)));
    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)));

    fig1 = plt.figure(1,figsize=(16,6));
    ax1 = plt.subplot(2,1,1);
    plt.axis([t[0,0], t[N-1,0], -1.1, 1.1]);
    plt.plot( t, u+0.1, 'k-');
    plt.plot( t, v-0.1, 'k-');
    plt.xlabel('time');
    plt.ylabel('u(t) and v(t)');

    # cross correlation of u and v
    c = eda_cvec( np.correlate(u.ravel(),v.ravel(),mode='full'));
    NA, i = np.shape(u);  # length of u
    NB, i = np.shape(v);  # length of v
    Nc,i = np.shape(c);
    tau = eda_cvec(Dt*np.linspace(-NB+1,NA-1,Nc))

    # lag to nearest sample
    cmax = np.max(c);
    icmax = np.argmax(c);
    tlag = -Dt * (icmax-NB+1);

    # plot lagged timeseries
    ax2 = plt.subplot(2,1,2);
    plt.axis([t[0,0], t[N-1,0], -1.1, 1.1]);
    plt.plot( t, u+0.1, 'k-');
    plt.plot( t-tlag, v-0.1, 'k-');
    plt.xlabel('time');
    plt.ylabel('u(t) and v(t)');

    #  plot cross-correlation
    fig2 = plt.figure(2,figsize=(16,6));
    ax2 = plt.subplot(1,1,1);
    plt.axis( [-20, 20, -1.1*cmax, 1.1*cmax] );
    plt.plot( tau, c, 'k-');
    plt.plot( [-tlag, -tlag], [-cmax, cmax], 'k:');
    plt.xlabel('time')
    plt.ylabel('cross-correlation');
    plt.show();

    print( "lag: true: %.4f estimated: %.4f" % (tlagtrue, tlag) );
    tlag, cmax = eda_timedelay(u, v, Dt);
    print( "lag: true: %.4f refined:   %.4f" % (tlagtrue, tlag) );

    fig2 = plt.figure(2,figsize=(16,6));
    ax1 = plt.subplot(2,1,1);
    plt.axis([t[0,0], t[N-1,0], -1.1, 1.1]);
    plt.plot( t, u+0.1, 'k-');
    plt.plot( t, v-0.1, 'k-');
    plt.xlabel('time');
    plt.ylabel('u(t) and v(t)');

    # plot lagged timeseries
    ax2 = plt.subplot(2,1,2);
    plt.axis([t[0,0], t[N-1,0], -1.1, 1.1]);
    plt.plot( t, u+0.1, 'k-');
    plt.plot( t-tlag, v-0.1, 'k-');
    plt.xlabel('time');
    plt.ylabel('u(t) and v(t)');
    plt.show();

    lag: true: 12.4667 estimated: 12.5000
    lag: true: 12.4667 refined:   12.4571

In [19]:

    # edapy_09_05
    # using the cross correlation to align two timeseries
    # West Point ozone data

    # load data
    S = np.genfromtxt('../Data/ozone_nohead.txt', delimiter='\t')
    [N, M]=S.shape;

    # break out time, ozone, radiation and temperature
    t = eda_cvec( S[:,0] );  # time in days after 1:00 AM on Aug 1, 1993
    Dt=1/24; # sampling interval
    tmin = t[0,0]; # start time
    tmax = t[N-1,0]; # end time
    ozone = eda_cvec( S[:,1] );     # ozone in ppm
    radiation = eda_cvec(S[:,2]);   # solar radiation in W/m2
    temp = eda_cvec(S[:,3]);        # temperature in deg C

    # plot original timeseries
    print('unlagged data');
    fig1 = plt.figure(1,figsize=(10,6));

    # solar radiation
    ax1 = plt.subplot(3,1,1);
    plt.axis([tmin, tmax, 0, np.max(radiation)]);
    plt.plot( t, radiation, 'k-');
    plt.xlabel('time (days)');
    plt.ylabel('solar (W/m2)');

    # ozone
    ax2 = plt.subplot(3,1,2);
    plt.axis([tmin, tmax, 0, np.max(ozone)]);
    plt.plot( t, ozone, 'k-');
    plt.xlabel('time (days)');
    plt.ylabel('ozone (ppm)');

    # both
    ax3 = plt.subplot(3,1,3);
    plt.axis([tmin, tmax, 0, 1.1]);
    plt.plot( t, radiation/np.max(radiation), 'k-', color=(0.7,0.7,0.7));
    plt.plot( t, ozone/np.max(ozone), 'k-');
    plt.xlabel('time (days)');
    plt.ylabel('both');

    plt.show();

    # cross correlation
    Nc = 2*N-1;
    c = eda_cvec( np.correlate(radiation.ravel(), ozone.ravel(), mode='full') );

    # time asxis for autocorrelation
    tau = eda_cvec( Dt*np.linspace(-N+1,N-1,2*N-1) );

    # find best lag
    Ncenter = (Nc+1)/2;
    cmax = np.max(c);
    icmax = np.argmax(c);
    tlag = -Dt * (icmax-N+1);

    # unlagged plot
    print('lagged data');
    fig2 = plt.figure(2,figsize=(10,6));

    # solar radiation
    ax1 = plt.subplot(3,1,1);
    plt.axis([tmin, tmax, 0, np.max(radiation)]);
    plt.plot( t, radiation, 'k-');
    plt.xlabel('time (days)');
    plt.ylabel('solar (W/m2)');

    # ozone
    ax2 = plt.subplot(3,1,2);
    plt.axis([tmin, tmax, 0, np.max(ozone)]);
    plt.plot( t-tlag, ozone, 'k-');
    plt.xlabel('time (days)');
    plt.ylabel('ozone (ppm)');

    # both
    ax3 = plt.subplot(3,1,3);
    plt.axis([tmin, tmax, 0, 1.1]);
    plt.plot( t, radiation/np.max(radiation), 'k-', color=(0.7,0.7,0.7));
    plt.plot( t-tlag, ozone/np.max(ozone), 'k-');
    plt.xlabel('time (days)');
    plt.ylabel('both');

    plt.show();

    # lagged plot, zommed
    days=5;
    print('lagged data, zoomed');
    fig2 = plt.figure(2,figsize=(10,6));

    # solar radiation
    ax1 = plt.subplot(3,1,1);
    plt.axis([t[0,0], t[days*24-1,0], 0, np.max(radiation)]);
    plt.plot( t, radiation, 'k-');
    plt.xlabel('time (days)');
    plt.ylabel('solar (W/m2)');

    # ozone
    ax2 = plt.subplot(3,1,2);
    plt.axis([t[0,0], t[days*24-1,0], 0, np.max(ozone)]);
    plt.plot( t-tlag, ozone, 'k-');
    plt.xlabel('time (days)');
    plt.ylabel('ozone (ppm)');

    # both
    ax3 = plt.subplot(3,1,3);
    plt.axis([t[0,0], t[days*24-1,0], 0, 1.1]);
    plt.plot( t, radiation/np.max(radiation), 'k-', color=(0.7,0.7,0.7));
    plt.plot( t-tlag, ozone/np.max(ozone), 'k-');
    plt.xlabel('time (days)');
    plt.ylabel('both');

    plt.show();

    plt.show();

    #plot cross-correlation
    fig3 = plt.figure(3,figsize=(8,6));
    ax2 = plt.subplot(1,1,1);
    plt.axis( [-5, 5, 0, 1.1*cmax] );
    plt.plot( tau, c, 'k-');
    plt.plot( [-tlag, -tlag], [-cmax, cmax], 'k:');
    plt.xlabel('time (days)')
    plt.ylabel('cross-correlation');
    plt.show();

    print( "estimated lag: %.2f hours" % (24*tlag) );

    unlagged data

    lagged data

    lagged data, zoomed

    estimated lag: 3.00 hours

In [2]:

    # edapy_09_06, filter estimation using autocorrelation
    # example of a filter problem
    # using impulse repsonse of a heat-producing layer
    # with solution via generalized least squares
    # using bicongugate gradient method
    # with the prior information of smoothness
    # Note: symbol q for theta, temperture

    # time setup for g (length M)
    M=1024;
    N=M;
    Dtdays = 0.5;
    tdays = eda_cvec( Dtdays*np.linspace(0,Dtdays*(N-1),N) );
    Dtseconds = Dtdays*3600*24;
    tseconds = Dtseconds*tdays;

    # materal constants
    rho = 1500; # kg/m3
    cp = 800; # J / kg-K
    kappa = 1.25e-6; # m2/s
    z = 1; # m

    # impulse response; handle time zero separately from the rest
    g = np.zeros((M,1));
    g[0,0] = 0.0;
    P1 = (1/(rho*cp));
    P2 = (Dtseconds/sqrt(2*pi));
    P3 = np.reciprocal( np.sqrt(2*kappa*tseconds[1:M]));
    P4 = np.exp( -0.5*(z**2)*np.reciprocal(2*kappa*tseconds[1:M]));
    g[1:M,0] = (P1 * P2 * np.multiply(P3 ,P4)).ravel();
    edafilterg = g; # global defines in GLSFilterMul()

    # true heat production
    t1 = Dtdays*N/4;
    t2 = Dtdays*7*N/16;
    sd = Dtdays*N/256;
    Ah=10;
    P5 = np.exp(-np.power(tdays-t1,2) / (2*sd**2) );
    P6 = np.exp(-np.power(tdays-t2,2) / (2*sd**2) );
    htrue = Ah*P5 + 0.5*Ah*P6;

    # predict true temperature by convolving g and htrue; truncate to length N
    tmp = np.convolve(g.ravel(), htrue.ravel());
    qtrue=eda_cvec( tmp[0:N] );

    # create simulated observations, qobs = qtrue+noise
    sigmad=0.01*np.amax(qtrue);   # std dev is 1% of maximum value of qtrue
    qobs=qtrue+np.random.normal(0,sigmad,(N,1)); # observed = signal + noie

    # solve for h given qobs and g using generalized least squares
    # with a prior smoothness constraint using biconjugate gradient

    # autocorrelation, a, of g. Just keep positive values
    # it needs to be of at least length M, since g is of
    # length N, and ordinarily N>=M.
    tmp = np.correlate(g.ravel(),g.ravel(),mode='full');
    edaautoc = eda_cvec(tmp[M-1:2*M-1]);

    # set up H and h and f
    eda_e=1.0*np.amax(g);  # damping parameter e=sigma_d/sigma_h
    edaHsparse = sparse.diags(
        [1.0, -2.0, 1.0], [0, 1, 2], shape=(M-2, M));
    h = np.zeros((M-2,1)); # prior information: 2nd derivative is zero

    # set up RHS
    FTf = GLSautoMulRHS(edafilterg,qobs,eda_e,edaHsparse,h,M);

    # define linear operator needed for conjugate gradienet solver
    LO=las.LinearOperator(shape=(M,M),matvec=GLSautoMul,rmatvec=GLSautoMul);

    # solve for estimated model parameters using congugate gradient algrorithm
    # note solver returns tupple first member of which is solution
    q=las.cg(LO,FTf,tol=1e-6, maxiter=(3*(N+M)+100));
    mest = eda_cvec( q[0] );

    hest3 = mest;

    # figure 1
    fig1 = plt.figure(1,figsize=(8,4));

    # plot observations
    ax1 = plt.subplot(3,1,1);
    plt.axis([tdays[0,0], tdays[N-1,0], 0, 1.1*np.max(qobs)] );
    plt.plot(tdays[0:N],qobs,'k-');
    plt.xlabel('time t, days');
    plt.ylabel('qobs(t)');
    plt.title('Autocorrelation Method');

    # plot true heat production
    ax1 = plt.subplot(3,1,2);
    plt.axis([tdays[0,0], tdays[N-1,0], -1.1*np.max(htrue), 1.1*np.max(htrue)] );
    plt.plot(tdays[0:N],htrue,'k-');
    plt.xlabel('time t, days');
    plt.ylabel('htrue(t)');

    # plot estimated heat production
    ax1 = plt.subplot(3,1,3);
    plt.axis([tdays[0,0], tdays[N-1,0], -1.1*np.max(htrue), 1.1*np.max(htrue)] );
    plt.plot(tdays[0:N],hest3,'k-');
    plt.xlabel('time t, days');
    plt.ylabel('hest(t)');
    plt.show();

In [21]:

    # edapy_09_07, smooth neuse.txt

    # load data
    D = np.genfromtxt('../Data/neuse.txt', delimiter='\t')
    [N, M]=D.shape;

    # break out time and data into separate vectors
    t = eda_cvec( D[:,0] );         # time in days after 1:00 AM on Aug 1, 1993
    Dt=t[1,0]-t[0,0];               # sampling interval
    tmin = t[0,0];                  # start time
    tmax = t[N-1,0];                # end time
    d = eda_cvec( D[:,1] );         # discharge in cfs

    # construct 3-point smoothing filter
    L1 = 3; # width of filter; must be odd
    s1 = eda_cvec( [1, 2, 1] );
    s1 = s1/np.sum(s1);  # normalize to sum=1
    Lh1 = floor(L1/2); # (half width)-1

    # filter and time shift
    temp = np.convolve(s1.ravel(), d.ravel(), mode='full');
    d1 = eda_cvec( temp[Lh1:N+Lh1] );

    # construct 19-point smoothing filter
    L2 = 19; # width of fileter, must be odd
    s2 = eda_cvec( [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 9, 8, 7, 6, 5, 4, 3, 2, 1] );
    s2 = s2/np.sum(s1); # normalize
    Lh2 = floor(L2/2); # (half width) - 1

    # filter and time shift
    temp = np.convolve(s2.ravel(), d.ravel(), mode='full');
    d2 = eda_cvec( temp[Lh2:N+Lh2] );

    # plot only this time range (for clarity)
    N1=0;
    N2=500;

    fig1 = plt.figure(1,figsize=(10,9));

    # original
    ax1 = plt.subplot(3,1,1);
    plt.axis([t[N1,0], t[N2,0], 0, np.max(d)]);
    plt.plot( t, d, 'k-');
    plt.xlabel('time (days)');
    plt.ylabel('discharge (cfs)');

    # filtered with 3-pt filter
    ax2 = plt.subplot(3,1,2);
    plt.axis([t[N1,0], t[N2,0], 0, np.max(d1)]);
    plt.plot( t, d1, 'k-');
    plt.xlabel('time (days)');
    plt.ylabel("%d point filter" % L1);

    # filtered with 19-pt filter
    ax3 = plt.subplot(3,1,3);
    plt.axis([t[N1,0], t[N2,0], 0, np.max(d2)]);
    plt.plot( t, d2, 'k-');
    plt.xlabel('time (days)');
    plt.ylabel("%d point filter" % L2);

    plt.show();

In [22]:

    # edapy_09_08
    # amplitude spectrum of uniform smoothing filter

    # standard set up of time
    N=256;             # length of time series, must be even
    Dt=0.01;           # sampling interval
    tmin = 0;          # start time
    tmax = Dt*(N-1);   # end time
    t=eda_cvec(np.linspace(tmin,tmax,N)); # time vector

    # standard set up of frequency
    No2 = floor(N/2);  # integer version of N/2
    fmax=1/(2.0*Dt);   # Nyquist frequency
    Df=fmax/No2;       # frequency sampling
    Nf=No2+1;          # number of non-negative frequencies
    # full (pos and neg) frequency vector
    f=eda_cvec(Df * np.concatenate( ( np.linspace(0,No2,Nf), np.linspace(-No2+1,-1,No2-1) ), axis=0 ));
    Dw=2*pi*Df;        # angular frequency sampling
    w=2*pi*f;          # full (pos and neg) angular frequency vector
    Nw=Nf;             # number of non-negative angular frequencies
    fpos=eda_cvec(Df * np.linspace(0,No2,Nf) );  # vector non-negative frequencies
    wpos=2*pi*fpos;    # vector non-negative angular frequencies

    # evaluate fourier transform of filter of length L
    T1=3*Dt;                       # width in time domain
    s1 = np.sinc(wpos*T1/(2*pi));  # sinc function in frequency
    sa1 = np.abs(s1);              # amplitude spectral density

    # evaluate fourier transform of filter of length L
    T2=21*Dt;                      # width in time domain
    s2 = np.sinc(wpos*T2/(2*pi));  # sinc function in frequency
    sa2 = np.abs(s2);              # amplitude spectral density

    # plot
    fig1 = plt.figure(1,figsize=(10,9));

    ax1 = plt.subplot(2,1,1);
    plt.axis( [0,fmax, 0, 1.1]);
    plt.plot( fpos, sa1, 'k-');
    plt.xlabel('frequency (Hz)');
    plt.ylabel('|s| for L=3');

    ax2 = plt.subplot(2,1,2);
    plt.axis( [0,fmax, 0, 1.1]);
    plt.plot( fpos, sa2, 'k-');
    plt.xlabel('frequency (Hz)');
    plt.ylabel('|s| for L=21');

    plt.show();

In [23]:

    # edapy_09_09
    # amplitude spectrum of Normal smoothing filter

    # standard set up of time
    N=256;             # length of time series, must be even
    Dt=0.01;           # sampling interval
    tmin = 0;          # start time
    tmax = Dt*(N-1);   # end time
    t=eda_cvec(np.linspace(tmin,tmax,N)); # time vector

    # standard set up of frequency
    No2 = floor(N/2);  # integer version of N/2
    fmax=1/(2.0*Dt);   # Nyquist frequency
    Df=fmax/No2;       # frequency sampling
    Nf=No2+1;          # number of non-negative frequencies
    # full (pos and neg) frequency vector
    f=eda_cvec(Df * np.concatenate( ( np.linspace(0,No2,Nf), np.linspace(-No2+1,-1,No2-1) ), axis=0 ));
    Dw=2*pi*Df;        # angular frequency sampling
    w=2*pi*f;          # full (pos and neg) angular frequency vector
    Nw=Nf;             # number of non-negative angular frequencies
    fpos=eda_cvec(Df * np.linspace(0,No2,Nf) );  # vector non-negative frequencies
    wpos=2*pi*fpos;    # vector non-negative angular frequencies

    # evaluate fourier transform of filter, g, of length L
    L1=3;

    # compute variance of filter (as if it were distribution)
    g1 = np.ones((L1,1));                    # vector of 1's
    g1 = g1/np.sum(g1);                      # normalize
    x = eda_cvec( np.linspace(0,L1-1,L1) );  # vector of x's
    mg1 = np.sum( np.multiply( x, g1 ) ) ;   # sum of ( x times g1 )
    vg1 = np.sum( np.multiply( np.power(x-mg1,2),  g1) ); # variance = sum of ( (x-mean)**2 times g1)
    sg1 = Dt*sqrt(vg1);                      # Dt * sqrt of variance

    # fourier transform of gausian, using analytic formula
    sigma1=1/sg1;
    s1 = np.exp(-np.power(wpos,2)/(2*(sigma1**2)));
    sa1 = np.abs(s1);  # amplitude spreal density
    print( "length %d samples (%.3f seconds) sigma-t samples %.3f (%.3f seconds)" % (L1, Dt*L1, sg1/Dt, sg1) );

    # plot
    fig1 = plt.figure(1,figsize=(10,9));

    ax1 = plt.subplot(2,1,1);
    plt.axis( [0,fmax, 0, 1.1]);
    plt.plot( fpos, sa1, 'k-');
    plt.xlabel('frequency (Hz)');
    plt.ylabel('|s| for L=3');

    # evaluate fourier transform of filter, g, of length L
    L2=21;

    # compute variance of filter (as if it were distribution)
    g2 = np.ones((L2,1));
    g2 = g2/np.sum(g2);
    x = eda_cvec( np.linspace(0,L2-1,L2) );
    mg2 = np.sum( np.multiply( x, g2 ) );  # sum of ( x times g2 )
    vg2 = np.sum( np.multiply( np.power(x-mg2,2),  g2) ); # variance = sum of ( (x-mean)**2 times g2)
    sg2 = Dt*sqrt(vg2); # standard deviation

    # fourier transform of gausian, using analytic formula
    sigma2=1/sg2;  # variance in frequency domain
    s2 = np.exp(-np.power(wpos,2)/(2*(sigma2**2)));  # fourier transform
    sa2 = np.abs(s2); # amplitude spectral density
    print( "length %d samples (%.3f seconds) sigma-t samples %.3f (%.3f seconds)" % (L2, Dt*L2, sg2/Dt, sg2) );

    ax2 = plt.subplot(2,1,2);
    plt.axis( [0,fmax, 0, 1.1]);
    plt.plot( fpos, sa2, 'k-');
    plt.xlabel('frequency (Hz)');
    plt.ylabel('|s| for L=21');

    plt.show();

    length 3 samples (0.030 seconds) sigma-t samples 0.816 (0.008 seconds)
    length 21 samples (0.210 seconds) sigma-t samples 6.055 (0.061 seconds)

In [24]:

    # edapy_09_10
    # complex z-plane analysis of filters, case 1, highpass

    # standard set up of time
    N=32;              # length of time series, must be even
    Dt=1.0 ;           # sampling interval
    tmin = 0;          # start time
    tmax = Dt*(N-1);   # end time
    t=eda_cvec(np.linspace(tmin,tmax,N)); # time vector

    # standard set up of frequency
    No2 = floor(N/2);  # integer version of N/2
    fmax=1/(2.0*Dt);   # Nyquist frequency
    Df=fmax/No2;       # frequency sampling
    Nf=No2+1;          # number of non-negative frequencies
    # full (pos and neg) frequency vector
    f=eda_cvec(Df * np.concatenate( ( np.linspace(0,No2,Nf), np.linspace(-No2+1,-1,No2-1) ), axis=0 ));
    Dw=2*pi*Df;        # angular frequency sampling
    w=2*pi*f;          # full (pos and neg) angular frequency vector
    Nw=Nf;             # number of non-negative angular frequencies
    fpos=eda_cvec(Df * np.linspace(0,No2,Nf) );  # vector non-negative frequencies
    wpos=2*pi*fpos;    # vector non-negative angular frequencies

    # a monomial filter, and its one root
    F = np.zeros((2,1));
    F[:,0] = [1, -1.1];
    zp = -F[0,0]/F[1,0];  # 0=F1+F2*z  so z=-F1/F2
    print("zero at %.3f" % (zp) );

    # complex z-plane
    Nz=128;   # z-plane image is Nz by Nz
    zmin=-2.0;  # min z
    zmax=2.0;   # max z
    dz=(zmax-zmin)/(Nz-1); # sampling in z
    zr=eda_cvec(np.linspace(zmin,zmax,Nz)); # real-z vector
    zi=eda_cvec(np.linspace(zmin,zmax,Nz)); # imag-z vertor

    # evaluate filter polynomial everywhere in complex z-plane
    S = np.zeros((Nz,Nz));
    for n in range(Nz):
        for m in range(Nz):
            S[m,n] = abs( F[0,0] + F[1,0]*(zr[n,0]+complex(0,1)*zi[m,0]) )**2;

    # plot amplitude of polynomial
    fig1 = plt.figure(1,figsize=(12,6));

    # black and white colormap
    bw = np.zeros((256,4));
    v = 0.9*(256 - np.linspace( 0, 255, 256 ))/255;
    bw[:,0] = v;
    bw[:,1] = v;
    bw[:,2] = v;
    bw[:,3] = np.ones(256);
    bwcmap = ListedColormap(bw);

    ax1 = plt.subplot(1,2,1);
    plt.axis( [zmin,zmax,zmin,zmax]);

    plt.imshow( S, cmap=bwcmap, vmin=np.min(S), vmax=np.max(S), extent=(zmin,zmax,zmin,zmax) );
    plt.xlabel('real z');
    plt.ylabel('imag z');
    plt.title('|S|^2');

    # plot Fourier transform points on the unit circle
    plt.plot(0,0,'wo');
    z=np.exp(-complex(0,1)*w*Dt);
    plt.plot(np.real(z),np.imag(z),'k+');

    # plot zero of filter
    plt.plot(np.real(zp),np.imag(zp),'w*');

    # plot frequency response of filter

    Sz = np.abs( F[0,0]+F[1,0]*z )**2;

    ax2 = plt.subplot(1,2,2);
    plt.axis( [0, fmax, 0, max(Sz)] );

    plt.xlabel('frequency, Hz');
    plt.ylabel('|S|^2');
    plt.title('|S|^2');
    plt.plot(fpos,Sz[0:Nf,0],'k-');
    plt.plot(fpos,Sz[0:Nf,0],'ko');

    plt.show();

    zero at 0.909

In [25]:

    # edapy_09_11
    # complex z-plane analysis of filters, case 1, lowpass

    # standard set up of time
    N=32;              # length of time series, must be even
    Dt=1.0 ;           # sampling interval
    tmin = 0;          # start time
    tmax = Dt*(N-1);   # end time
    t=eda_cvec(np.linspace(tmin,tmax,N)); # time vector

    # standard set up of frequency
    No2 = floor(N/2);  # integer version of N/2
    fmax=1/(2.0*Dt);   # Nyquist frequency
    Df=fmax/No2;       # frequency sampling
    Nf=No2+1;          # number of non-negative frequencies
    # full (pos and neg) frequency vector
    f=eda_cvec(Df * np.concatenate( ( np.linspace(0,No2,Nf), np.linspace(-No2+1,-1,No2-1) ), axis=0 ));
    Dw=2*pi*Df;        # angular frequency sampling
    w=2*pi*f;          # full (pos and neg) angular frequency vector
    Nw=Nf;             # number of non-negative angular frequencies
    fpos=eda_cvec(Df * np.linspace(0,No2,Nf) );  # vector non-negative frequencies
    wpos=2*pi*fpos;    # vector non-negative angular frequencies

    # a monomial filter, and its one root
    F = eda_cvec( [1, 1.1] );
    zp = -F[0,0]/F[1,0];  # 0=F1+F2*z  so z=-F1/F2
    print("zero at %.3f" % (zp) );

    # complex z-plane
    Nz=128;   # z-plane image is Nz by Nz
    zmin=-2.0;  # min z
    zmax=2.0;   # max z
    dz=(zmax-zmin)/(Nz-1); # sampling in z
    zr=eda_cvec(np.linspace(zmin,zmax,Nz)); # real-z vector
    zi=eda_cvec(np.linspace(zmin,zmax,Nz)); # imag-z vertor

    # evaluate filter polynomial everywhere in complex z-plane
    S = np.zeros((Nz,Nz));
    for n in range(Nz):
        for m in range(Nz):
            S[m,n] = abs( F[0,0] + F[1,0]*(zr[n,0]+complex(0,1)*zi[m,0]) )**2;

    # plot amplitude of polynomial
    fig1 = plt.figure(1,figsize=(12,6));

    # black and white colormap
    bw = np.zeros((256,4));
    v = 0.9*(256 - np.linspace( 0, 255, 256 ))/255;
    bw[:,0] = v;
    bw[:,1] = v;
    bw[:,2] = v;
    bw[:,3] = np.ones(256);
    bwcmap = ListedColormap(bw);

    ax1 = plt.subplot(1,2,1);
    plt.axis( [zmin,zmax,zmin,zmax]);

    plt.imshow( S, cmap=bwcmap, vmin=np.min(S), vmax=np.max(S), extent=(zmin,zmax,zmin,zmax) );
    plt.xlabel('real z');
    plt.ylabel('imag z');
    plt.title('|S|^2');

    # plot Fourier transform points on the unit circle
    plt.plot(0,0,'wo');
    z=np.exp(-complex(0,1)*w*Dt);
    plt.plot(np.real(z),np.imag(z),'k+');

    # plot zero of filter
    plt.plot(np.real(zp),np.imag(zp),'w*');

    # plot frequency response of filter

    Sz = np.abs( F[0,0]+F[1,0]*z )**2;

    ax2 = plt.subplot(1,2,2);
    plt.axis( [0, fmax, 0, max(Sz)] );

    plt.xlabel('frequency, Hz');
    plt.ylabel('|S|^2');
    plt.title('|S|^2');
    plt.plot(fpos,Sz[0:Nf,0],'k-');
    plt.plot(fpos,Sz[0:Nf,0],'ko');

    plt.show();

    zero at -0.909

In [26]:

    # edapy_09_12
    # complex z-plane analysis of filters, case 3, narrow band

    # standard set up of time
    N=32;              # length of time series, must be even
    Dt=1.0 ;           # sampling interval
    tmin = 0;          # start time
    tmax = Dt*(N-1);   # end time
    t=eda_cvec(np.linspace(tmin,tmax,N)); # time vector

    # standard set up of frequency
    No2 = floor(N/2);  # integer version of N/2
    fmax=1/(2.0*Dt);   # Nyquist frequency
    Df=fmax/No2;       # frequency sampling
    Nf=No2+1;          # number of non-negative frequencies
    # full (pos and neg) frequency vector
    f=eda_cvec(Df * np.concatenate( ( np.linspace(0,No2,Nf), np.linspace(-No2+1,-1,No2-1) ), axis=0 ));
    Dw=2*pi*Df;        # angular frequency sampling
    w=2*pi*f;          # full (pos and neg) angular frequency vector
    Nw=Nf;             # number of non-negative angular frequencies
    fpos=eda_cvec(Df * np.linspace(0,No2,Nf) );  # vector non-negative frequencies
    wpos=2*pi*fpos;    # vector non-negative angular frequencies

    # a monomial filter, and its one root
    F1 = eda_cvec( [1, 0.6*(1-1.1*complex(0.0,1.0))] );
    zp1 = -F1[0,0]/F1[1,0];  # 0=F1+F2*z  so z=-F1/F2
    F2 = np.conj(F1);
    zp2 = -F2[0,0]/F2[1,0];
    print("poles at (%.3f %.3f) and (%.3f %.3f)" % (np.real(zp1),np.imag(zp1),np.real(zp2),np.imag(zp2)) );

    # complex z-plane
    Nz=128;     # z-plane image is Nz by Nz
    zmin=-2.0;  # min z
    zmax=2.0;   # max z
    dz=(zmax-zmin)/(Nz-1); # sampling in z
    zr=eda_cvec(np.linspace(zmin,zmax,Nz)); # real-z vector
    zi=eda_cvec(np.linspace(zmin,zmax,Nz)); # imag-z vertor

    # evaluate filter polynomial everywhere in complex z-plane
    S = np.zeros((Nz,Nz));
    for n in range(Nz):
        for m in range(Nz):
            v1 = 1/(F1[0,0] + F1[1,0]*(zr[n,0]+complex(0.0,1.0)*zi[m,0]));
            v2 = 1/(F2[0,0] + F2[1,0]*(zr[n,0]+complex(0.0,1.0)*zi[m,0]));
            S[m,n] = log( abs( v1*v2 )**2);  # take the log for plot-scaling  purposes

    # plot amplitude of polynomial
    fig1 = plt.figure(1,figsize=(12,6));

    # black and white colormap
    bw = np.zeros((256,4));
    v = 0.9*(256 - np.linspace( 0, 255, 256 ))/255;
    bw[:,0] = v;
    bw[:,1] = v;
    bw[:,2] = v;
    bw[:,3] = np.ones(256);
    bwcmap = ListedColormap(bw);

    ax1 = plt.subplot(1,2,1);
    plt.axis( [zmin,zmax,zmin,zmax]);

    plt.imshow( S, cmap=bwcmap, vmin=np.min(S), vmax=np.max(S), extent=(zmin,zmax,zmin,zmax) );
    plt.xlabel('real z');
    plt.ylabel('imag z');
    plt.title('|S|^2');

    # plot Fourier transform points on the unit circle
    plt.plot(0,0,'wo');
    z=np.exp(-complex(0,1)*w*Dt);
    plt.plot(np.real(z),np.imag(z),'k+');

    # plot poles of filter
    plt.plot(np.real(zp1),np.imag(zp1),'w*');
    plt.plot(np.real(zp2),np.imag(zp2),'w*');

    # plot frequency response of filter
    Sz = np.abs( np.reciprocal( np.multiply(F1[0,0]+F1[1,0]*z,F2[0,0]+F2[1,0]*z)) )**2;

    ax2 = plt.subplot(1,2,2);
    plt.axis( [0, fmax, 0, max(Sz)] );

    plt.xlabel('frequency, Hz');
    plt.ylabel('|S|^2');
    plt.title('|S|^2');
    plt.plot(fpos,Sz[0:Nf,0],'k-');
    plt.plot(fpos,Sz[0:Nf,0],'ko');

    plt.show();

    poles at (-0.754 -0.830) and (-0.754 0.830)

In [27]:

    # edapy_09_13
    # complex z-plane analysis of filters, case 4, notch

    # standard set up of time
    N=32;              # length of time series, must be even
    Dt=1.0 ;           # sampling interval
    tmin = 0;          # start time
    tmax = Dt*(N-1);   # end time
    t=eda_cvec(np.linspace(tmin,tmax,N)); # time vector

    # standard set up of frequency
    No2 = floor(N/2);  # integer version of N/2
    fmax=1/(2.0*Dt);   # Nyquist frequency
    Df=fmax/No2;       # frequency sampling
    Nf=No2+1;          # number of non-negative frequencies
    # full (pos and neg) frequency vector
    f=eda_cvec(Df * np.concatenate( ( np.linspace(0,No2,Nf), np.linspace(-No2+1,-1,No2-1) ), axis=0 ));
    Dw=2*pi*Df;        # angular frequency sampling
    w=2*pi*f;          # full (pos and neg) angular frequency vector
    Nw=Nf;             # number of non-negative angular frequencies
    fpos=eda_cvec(Df * np.linspace(0,No2,Nf) );  # vector non-negative frequencies
    wpos=2*pi*fpos;    # vector non-negative angular frequencies

    # a monomial filter, and its one root
    F1 = eda_cvec( [1, 0.9*complex(0.0,1.0)] );
    zp1 = -F1[0,0]/F1[1,0];  # 0=F1+F2*z  so z=-F1/F2
    F2 = np.conj(F1);
    zp2 = -F2[0,0]/F2[1,0];
    print("zeros at (%.3f %.3f) and (%.3f %.3f)" % (np.real(zp1),np.imag(zp1),np.real(zp2),np.imag(zp2)) );

    F3 = eda_cvec( [1, 0.8*complex(0.0,1.0)] );
    zp3 = -F3[0,0]/F3[1,0];  # 0=F1+F2*z  so z=-F1/F2
    F4 = np.conj(F3);
    zp4 = -F4[0,0]/F4[1,0];
    print("poles at (%.3f %.3f) and (%.3f %.3f)" % (np.real(zp3),np.imag(zp3),np.real(zp4),np.imag(zp4)) );

    # complex z-plane
    Nz=128;     # z-plane image is Nz by Nz
    zmin=-2.0;  # min z
    zmax=2.0;   # max z
    dz=(zmax-zmin)/(Nz-1); # sampling in z
    zr=eda_cvec(np.linspace(zmin,zmax,Nz)); # real-z vector
    zi=eda_cvec(np.linspace(zmin,zmax,Nz)); # imag-z vertor

    # evaluate filter polynomial everywhere in complex z-plane
    S = np.zeros((Nz,Nz));
    for n in range(Nz):
        for m in range(Nz):
            v1 = (F1[0,0] + F1[1,0]*(zr[n,0]+complex(0.0,1.0)*zi[m,0]));
            v2 = (F2[0,0] + F2[1,0]*(zr[n,0]+complex(0.0,1.0)*zi[m,0]));
            v3 = (F3[0,0] + F3[1,0]*(zr[n,0]+complex(0.0,1.0)*zi[m,0]));
            v4 = (F4[0,0] + F4[1,0]*(zr[n,0]+complex(0.0,1.0)*zi[m,0]));
            S[m,n] = log( abs( (v1*v2)/(v3*v4) )**2);  # take the log for plot-scaling  purposes

    # plot amplitude of polynomial
    fig1 = plt.figure(1,figsize=(12,6));

    # black and white colormap
    bw = np.zeros((256,4));
    v = 0.9*(256 - np.linspace( 0, 255, 256 ))/255;
    bw[:,0] = v;
    bw[:,1] = v;
    bw[:,2] = v;
    bw[:,3] = np.ones(256);
    bwcmap = ListedColormap(bw);

    ax1 = plt.subplot(1,2,1);
    plt.axis( [zmin,zmax,zmin,zmax]);

    plt.imshow( S, cmap=bwcmap, vmin=np.min(S), vmax=np.max(S), extent=(zmin,zmax,zmin,zmax) );
    plt.xlabel('real z');
    plt.ylabel('imag z');
    plt.title('|S|^2');

    # plot Fourier transform points on the unit circle
    plt.plot(0,0,'wo');
    z=np.exp(-complex(0,1)*w*Dt);
    plt.plot(np.real(z),np.imag(z),'k+');

    # plot poles of filter
    plt.plot(np.real(zp1),np.imag(zp1),'w*');
    plt.plot(np.real(zp2),np.imag(zp2),'w*');
    plt.plot(np.real(zp3),np.imag(zp3),'w+');
    plt.plot(np.real(zp4),np.imag(zp4),'w+');

    # plot frequency response of filter
    partA = np.multiply(  F1[0,0]+F1[1,0]*z,  F2[0,0]+F2[1,0]*z);
    partB = np.multiply(  F3[0,0]+F3[1,0]*z,  F4[0,0]+F4[1,0]*z);
    Sz = np.abs(np.divide(partA, partB))**2;

    ax2 = plt.subplot(1,2,2);
    plt.axis( [0, fmax, 0, max(Sz)] );

    plt.xlabel('frequency, Hz');
    plt.ylabel('|S|^2');
    plt.title('|S|^2');
    plt.plot(fpos,Sz[0:Nf,0],'k-');
    plt.plot(fpos,Sz[0:Nf,0],'ko');

    plt.show();

    zeros at (-0.000 1.111) and (-0.000 -1.111)
    poles at (-0.000 1.250) and (-0.000 -1.250)

In [28]:

    # edapy_09_14, Chebychev bandpass filter

    # limits of bandpass filter
    flow=5;
    fhigh=10;

    # standard set up of time
    N=1024;            # length of time series, must be even
    Dt=0.01 ;          # sampling interval
    tmin = 0;          # start time
    tmax = Dt*(N-1);   # end time
    t=eda_cvec(np.linspace(tmin,tmax,N)); # time vector

    # standard set up of frequency
    No2 = floor(N/2);  # integer version of N/2
    fmax=1/(2.0*Dt);   # Nyquist frequency
    Df=fmax/No2;       # frequency sampling
    Nf=No2+1;          # number of non-negative frequencies
    # full (pos and neg) frequency vector
    f=eda_cvec(Df * np.concatenate( ( np.linspace(0,No2,Nf), np.linspace(-No2+1,-1,No2-1) ), axis=0 ));
    Dw=2*pi*Df;        # angular frequency sampling
    w=2*pi*f;          # full (pos and neg) angular frequency vector
    Nw=Nf;             # number of non-negative angular frequencies
    fpos=eda_cvec(Df * np.linspace(0,No2,Nf) );  # vector non-negative frequencies
    wpos=2*pi*fpos;    # vector non-negative angular frequencies

    # input timeseries is a spike
    din=np.zeros((N,1));
    din[100,0]=1;

    # filter
    [dout, u, v] = chebyshevfilt( din, Dt, flow, fhigh);

    #  spectrum of input
    dint = eda_cvec( np.fft.fft(din, axis=0) );  # fourier transform
    dintpos = dint[0:Nf,0]; # keep only non-negative frequencies
    dinspec = np.abs(dintpos); # amplitude spectral density

    #  spectrum of output
    doutt = eda_cvec( np.fft.fft(dout, axis=0) ); # fourier transform
    douttpos = eda_cvec( doutt[0:Nf,0] ); # keep only non-negative frequencies
    doutspec = np.abs(douttpos); # amplitude spectral density

    # plot unfiltered and filtered timeseries
    fig1 = plt.figure(1,figsize=(12,6));
    ax1 = plt.subplot(2,1,1);
    plt.axis( [0.75, 1.25, -1, 1 ]);
    plt.plot(t, din, 'k-' );
    plt.xlabel('time (s)');
    plt.ylabel('din(t)');
    plt.title('timeseries');
    ax2 = plt.subplot(2,1,2);
    plt.axis( [0.75, 1.25, -1, 1 ]);
    plt.plot(t, dout, 'k-' );
    plt.xlabel('time (s)');
    plt.ylabel('din(t)');
    plt.show();

    # plot spectra
    fig2 = plt.figure(2,figsize=(12,6));
    ax1 = plt.subplot(1,2,1);
    specmax=np.max(dinspec);
    plt.axis( [0, fmax, 0, 1.1*specmax ]);
    plt.plot(fpos, dinspec, 'k-' );
    plt.xlabel('frequency (Hx)');
    plt.ylabel('din(f)');
    plt.title('input spectra');
    ax2 = plt.subplot(1,2,2);
    plt.axis( [0, fmax, 0, 1.1*specmax ]);
    plt.plot(fpos, doutspec, 'k-' );
    plt.xlabel('frequency (Hz)');
    plt.ylabel('dout(f)');
    plt.title('output spectra');
    plt.plot( [flow, flow], [0, 0.2*specmax], 'k:' );
    plt.plot( [fhigh, fhigh], [0, 0.2*specmax], 'k:' );
    plt.show();

    # roots of filters
    [Nu, i] = np.shape(u);           # length of u
    Nru = Nu-1;                      # number or roots
    pol = (np.flipud(u)).ravel();    # polynomial coefficents
    ru = eda_cvec( np.roots(pol) );  # vector of roots

    [Nv,i] = np.shape(v);           # length of u
    Nrv = Nv-1;                     # number or roots
    pol = (np.flipud(v)).ravel();   # polynomial coefficents
    rv = eda_cvec( np.roots(pol) ); # vector of roots

    # complex z-plane
    Nz=128;     # z-plane image is Nz by Nz
    zmin=-2.0;  # min z
    zmax=2.0;   # max z
    dz=(zmax-zmin)/(Nz-1); # sampling in z
    zr=eda_cvec(np.linspace(zmin,zmax,Nz)); # real-z vector
    zi=eda_cvec(np.linspace(zmin,zmax,Nz)); # imag-z vertor

    # evaluate filter polynomials everywhere in complex z-plane
    S = np.zeros((Nz,Nz));
    for n in range(Nz):
        for m in range(Nz):
            z = zr[n,0]+complex(0.0,1.0)*zi[m,0];
            uu=0.0;
            zpow=1;
            for k in range(Nu):
                uu = uu + u[k,0]*zpow;
                zpow=zpow*z;
            vv=0.0;
            zpow=1;
            for k in range(Nv):
                vv = vv + v[k,0]*zpow;
                zpow=zpow*z;
            SS = uu/vv;
            S[m,n]=np.abs(SS);

    # plot amplitude of polynomial
    fig3 = plt.figure(3,figsize=(12,6));

    # black and white colormap
    bw = np.zeros((256,4));
    v = 0.9*(256 - np.linspace( 0, 255, 256 ))/255;
    bw[:,0] = v;
    bw[:,1] = v;
    bw[:,2] = v;
    bw[:,3] = np.ones(256);
    bwcmap = ListedColormap(bw);

    ax1 = plt.subplot(1,2,1);
    plt.axis( [zmin,zmax,zmin,zmax]);

    plt.imshow( S, cmap=bwcmap, vmin=np.min(S), vmax=np.max(S), extent=(zmin,zmax,zmin,zmax) );
    plt.xlabel('real z');
    plt.ylabel('imag z');
    plt.title('|S|^2');

    # plot Fourier transform points on the unit circle
    plt.plot(0,0,'wo');
    z=np.exp(-complex(0,1)*w*Dt);
    plt.plot(np.real(z),np.imag(z),'k-');

    # plot origin
    plt.plot(0,0,'ko');
    # plot unit circle
    z=np.exp(-complex(0.0,1.0)*w*Dt);
    plt.plot(np.real(z),np.imag(z),'k-');

    # plot zeros and poles
    for k in range(Nru):
        plt.plot(np.real(ru[k,0]),np.imag(ru[k,0]),'k*');

    for k in range(Nrv):
        plt.plot(np.real(rv[k,0]),np.imag(rv[k,0]),'k+');

    plt.show();

In [31]:

    # edapy_09_15
    # read Reynolds Channel Water Quality dataset and plot it
    # columns are days, precip, air_temp, water_temp, salinity, turbidity, chlorophyl

    # load data from file
    S = np.genfromtxt('../Data/reynolds_interpolated.txt', delimiter='\t')
    [N, M]=S.shape;

    # pull out time vector
    t = eda_cvec( D[:,0] )  # time in days
    Dt=t[1,0]-t[0,0];       # sampling interval
    tmin = t[0,0];          # start time
    tmax = t[N-1,0];        # end time

    # pull out data vectors
    p =  eda_cvec(S[:,1]);       # precipitation, inches
    at = eda_cvec(S[:,2]);       # air temperature, deg C
    wt = eda_cvec(S[:,3]);       # water temperature, deg C
    s =  eda_cvec(S[:,4]);       # salinity, practical salinity units
    tb = eda_cvec(S[:,5]);       # turbidity, (no units given)
    c =  eda_cvec(S[:,6]);       # chlorophyl, micrograms per liter

    fig1 = plt.figure(1,figsize=(12,6));

    # plotting limits
    itmin=0;
    itmax=N;
    tmin=t[itmin];
    tmax=t[itmax-1];

    ax1 = plt.subplot(6,1,1);
    plt.axis( [tmin, tmax, np.min(p[itmin:itmax,0]), np.max(p[itmin:itmax,0])]);
    plt.plot(t, p, 'k-' );
    plt.ylabel('precipitation');
    plt.title('unfiltered data');

    ax2 = plt.subplot(6,1,2);
    plt.axis( [tmin, tmax, np.min(at[itmin:itmax,0]), np.max(at[itmin:itmax,0])]);
    plt.plot(t, at, 'k-' );
    plt.ylabel('T-air');

    ax3 = plt.subplot(6,1,3);
    plt.axis( [tmin, tmax, np.min(wt[itmin:itmax,0]), np.max(wt[itmin:itmax,0])]);
    plt.plot(t, wt, 'k-' );
    plt.ylabel('T-water');

    ax4 = plt.subplot(6,1,4);
    plt.axis( [tmin, tmax, np.min(s[itmin:itmax,0]), np.max(s[itmin:itmax,0])]);
    plt.plot(t, s, 'k-' );
    plt.ylabel('salinity');

    ax5 = plt.subplot(6,1,5);
    plt.axis( [tmin, tmax, np.min(tb[itmin:itmax,0]), np.max(tb[itmin:itmax,0])]);
    plt.plot(t, tb, 'k-' );
    plt.ylabel('turbidity');

    ax6 = plt.subplot(6,1,6);
    plt.axis( [tmin, tmax, np.min(c[itmin:itmax,0]), np.max(c[itmin:itmax,0])]);
    plt.plot(t, c, 'k-' );
    plt.ylabel('chlorophyl');

In [32]:

    # edapy_09_16
    # read Reynolds Channel Water Quality dataset and plot it
    # columns are days, precip, air_temp, water_temp, salinity, turbidity, chlorophyl

    # load data from file
    S = np.genfromtxt('../Data/reynolds_interpolated.txt', delimiter='\t')
    [N, M]=S.shape;

    # pull out time vector
    t = eda_cvec( D[:,0] ); # time in days
    Dt=t[1,0]-t[0,0];       # sampling interval
    tmin = t[0,0];          # start time
    tmax = t[N-1,0];        # end time

    # low-ish bandpass filter parameters
    fcenter = 1/365.25;
    flow = 0.75*fcenter;
    fhigh =1.25*fcenter;

    # bandpass filter
    d = np.zeros((N,1));
    Df = np.zeros((N,M));
    Df[:,0] = S[:,0];
    for i in range(1,M):
        d[:,0] = D[:,i];
        [z, u, v] = chebyshevfilt(d-np.mean(d), Dt, flow, fhigh);
        Df[:,i] = z.ravel();

    fig1 = plt.figure(1,figsize=(12,6));

    # pull out data vectors
    p =  eda_cvec(Df[:,1]);       # precipitation, inches
    at = eda_cvec(Df[:,2]);       # air temperature, deg C
    wt = eda_cvec(Df[:,3]);       # water temperature, deg C
    s =  eda_cvec(Df[:,4]);       # salinity, practical salinity units
    tb = eda_cvec(Df[:,5]);       # turbidity, (no units given)
    c =  eda_cvec(Df[:,6]);       # chlorophyl, micrograms per liter

    # plotting limits
    itmin=0;
    itmax=N;
    tmin=t[itmin];
    tmax=t[itmax-1];

    ax1 = plt.subplot(6,1,1);
    plt.axis( [tmin, tmax, np.min(p[itmin:itmax,0]), np.max(p[itmin:itmax,0])]);
    plt.plot(t, p, 'k-' );
    plt.ylabel('precipitation');
    plt.title("data filtered bwetween periods of %.1f and %.1f days" % (1/fhigh,1/flow) );

    ax2 = plt.subplot(6,1,2);
    plt.axis( [tmin, tmax, np.min(at[itmin:itmax,0]), np.max(at[itmin:itmax,0])]);
    plt.plot(t, at, 'k-' );
    plt.ylabel('T-air');

    ax3 = plt.subplot(6,1,3);
    plt.axis( [tmin, tmax, np.min(wt[itmin:itmax,0]), np.max(wt[itmin:itmax,0])]);
    plt.plot(t, wt, 'k-' );
    plt.ylabel('T-water');

    ax4 = plt.subplot(6,1,4);
    plt.axis( [tmin, tmax, np.min(s[itmin:itmax,0]), np.max(s[itmin:itmax,0])]);
    plt.plot(t, s, 'k-' );
    plt.ylabel('salinity');

    ax5 = plt.subplot(6,1,5);
    plt.axis( [tmin, tmax, np.min(tb[itmin:itmax,0]), np.max(tb[itmin:itmax,0])]);
    plt.plot(t, tb, 'k-' );
    plt.ylabel('turbidity');

    ax6 = plt.subplot(6,1,6);
    plt.axis( [tmin, tmax, np.min(c[itmin:itmax,0]), np.max(c[itmin:itmax,0])]);
    plt.plot(t, c, 'k-' );
    plt.ylabel('chlorophyl');

    plt.show();

    # high-ish bandpass filter parameters
    fcenter = 1.0/5.0;
    flow = 0.75*fcenter;
    fhigh =1.25*fcenter;

    # bandpass filter
    d = np.zeros((N,1));
    Df = np.zeros((N,M));
    Df[:,0] = S[:,0];
    for i in range(1,M):
        d[:,0] = D[:,i];
        [z, u, v] = chebyshevfilt(d-np.mean(d), Dt, flow, fhigh);
        Df[:,i] = z.ravel();

    fig2 = plt.figure(2,figsize=(12,6));

    # pull out data vectors
    p =  eda_cvec(Df[:,1]);       # precipitation, inches
    at = eda_cvec(Df[:,2]);       # air temperature, deg C
    wt = eda_cvec(Df[:,3]);       # water temperature, deg C
    s =  eda_cvec(Df[:,4]);       # salinity, practical salinity units
    tb = eda_cvec(Df[:,5]);       # turbidity, (no units given)
    c =  eda_cvec(Df[:,6]);       # chlorophyl, micrograms per liter

    # plotting limits
    itmin=0;
    itmax=N;
    tmin=t[itmin];
    tmax=t[itmax-1];

    ax1 = plt.subplot(6,1,1);
    plt.axis( [tmin, tmax, np.min(p[itmin:itmax,0]), np.max(p[itmin:itmax,0])]);
    plt.plot(t, p, 'k-' );
    plt.ylabel('precipitation');
    plt.title("data filtered bwetween periods of %.1f and %.1f days" % (1/fhigh,1/flow) );

    ax2 = plt.subplot(6,1,2);
    plt.axis( [tmin, tmax, np.min(at[itmin:itmax,0]), np.max(at[itmin:itmax,0])]);
    plt.plot(t, at, 'k-' );
    plt.ylabel('T-air');

    ax3 = plt.subplot(6,1,3);
    plt.axis( [tmin, tmax, np.min(wt[itmin:itmax,0]), np.max(wt[itmin:itmax,0])]);
    plt.plot(t, wt, 'k-' );
    plt.ylabel('T-water');

    ax4 = plt.subplot(6,1,4);
    plt.axis( [tmin, tmax, np.min(s[itmin:itmax,0]), np.max(s[itmin:itmax,0])]);
    plt.plot(t, s, 'k-' );
    plt.ylabel('salinity');

    ax5 = plt.subplot(6,1,5);
    plt.axis( [tmin, tmax, np.min(tb[itmin:itmax,0]), np.max(tb[itmin:itmax,0])]);
    plt.plot(t, tb, 'k-' );
    plt.ylabel('turbidity');

    ax6 = plt.subplot(6,1,6);
    plt.axis( [tmin, tmax, np.min(c[itmin:itmax,0]), np.max(c[itmin:itmax,0])]);
    plt.plot(t, c, 'k-' );
    plt.ylabel('chlorophyl');

    plt.show();

In [33]:

    # edapy_09_17
    # cohernce calculation (single frequency band only)
    # as a check of the formula for C

    # make a time vector
    N=256;  # length of time series
    Dt=1.0; # sampling interval
    t = eda_cvec( Dt*np.linspace(0,N-1,N) ); # time vector

    # random time series
    rts1=np.random.normal( 0.0, 1.0, (256, 1) );
    rts2=np.random.normal( 0.0, 1.0, (256, 1) );
    rts3=np.random.normal( 0.0, 1.0, (256, 1) );

    # create partially correlated time series, by adding
    # together uncorrelated ones
    uorig = rts1 + rts2;
    vorig = rts1 + rts3;

    # standard set up of frequency
    No2 = floor(N/2);  # integer version of N/2
    fmax=1/(2.0*Dt);   # Nyquist frequency
    Df=fmax/No2;       # frequency sampling
    Nf=No2+1;          # number of non-negative frequencies
    # full (pos and neg) frequency vector
    f=eda_cvec(Df * np.concatenate( ( np.linspace(0,No2,Nf), np.linspace(-No2+1,-1,No2-1) ), axis=0 ));
    Dw=2*pi*Df;        # angular frequency sampling
    w=2*pi*f;          # full (pos and neg) angular frequency vector
    Nw=Nf;             # number of non-negative angular frequencies
    fpos=eda_cvec(Df * np.linspace(0,No2,Nf) );  # vector non-negative frequencies
    wpos=2*pi*fpos;    # vector non-negative angular frequencies

    # frequency and bandwith for calculation
    f0 = fmax/2;            # center of frequency [0,nyquist] range
    lowside = 0.75;         # fraction below nyquist
    highside =1.25;         # ractiona bve nyquist
    flow = f0*lowside;      # low side of frequency band
    fhigh= f0*highside;     # high side of frequency band
    ilow = floor(flow/Df);    # corresponding index
    ihigh = floor(fhigh/Df);  # corresponding index

    # box car filter matching the parameters above
    bp = np.zeros((N,1));
    bp[ilow:ihigh,0]=1.0;          # positive frequencies
    bp[N-ihigh+1:N-ilow+1,0]=1.0;  # negative frequencies

    # bandpass filter the data, by tahing inverse fourier transform of product fourier transforms
    ufilt = eda_cvec( np.real( np.fft.ifft( np.multiply( bp, np.fft.fft(uorig,axis=0) ), axis=0 ) ) );
    vfilt = eda_cvec( np.real( np.fft.ifft( np.multiply( bp, np.fft.fft(vorig,axis=0) ), axis=0 ) ) );

    # calculation of power and zero-lag cross-correlation
    Pu = np.sum( np.power(ufilt,2) );
    Pv = np.sum( np.power(vfilt,2) );
    x = np.correlate( ufilt.ravel(), vfilt.ravel(), mode='full' );
    x0 = eda_cvec( x[N-1] );

    # plot time series
    fig1 = plt.figure(1,figsize=(12,6));
    ax1 = plt.subplot(1,1,1);
    plt.axis( [t[0,0], t[N-1,0], -10, 10 ] );
    plt.plot(t, ufilt+3, 'k-' );
    plt.plot(t, vfilt-3, 'k-' );
    plt.xlabel('time (s)');
    plt.ylabel('u(t) and v(t)');
    plt.title("pass-band: %.2f to %.2f Hz" % (flow, fhigh));
    plt.show();

    # display results
    print('Calculated directly from the band-passed time series:');
    print("    Power in u %.4f" % (Pu) );
    print("    Power in v %.4f" % (Pv) );
    print("    Zero-lag cross correlation %.4f" % (x0) );
    print("    Curly C squared %.4f" % ((x0**2)/(Pu*Pv)) );
    print(' ');

    # Fourier transform
    ut = eda_cvec( np.fft.fft(uorig, axis=0) ); # fourier trasnform
    u =  eda_cvec( ut[0:Nf,0] ); # delete negative frequencies

    vt = eda_cvec( np.fft.fft(vorig,axis=0) ); # fourier trasnform
    v =  eda_cvec(vt[0:Nf,0]); # delete negative frequencies

    # spectral densities
    usv = np.multiply( np.conj(u), v ); # cross spectral density of u and v
    usu = np.real( np.multiply( np.conj(u), u ) ); # power spectral density of u
    vsv = np.real( np.multiply( np.conj(v), v ) ); # power spectral density of v

    # band-averaged quantities
    usva = np.sum( (usv[ilow:ihigh,0]) ) / (Nf-1);
    usvar = np.sum( np.real( (usv[ilow:ihigh,0] ) ) )/ (Nf-1);
    usua = np.sum( usu[ilow:ihigh,0] )/ (Nf-1);
    vsva = np.sum( vsv[ilow:ihigh,0] )/ (Nf-1);
    curlyc2 = np.multiply( np.conj(usvar), usvar) / (usua * vsva);
    coherence = np.real(np.multiply(np.conj(usva),usva)) / (usua * vsva);

    # display quantities
    print('Calculated from the cross spectral density');
    print("    Pover in u %.4f" % (usua));
    print("    Pover in v %.4f" % (vsva));
    print("    Zero-lag cross correlation %.4f" % (usvar));
    print("    Curly C squared %.4f" % (curlyc2));
    print("    Coherence %.4f" % (coherence));

    Calculated directly from the band-passed time series:
        Power in u 124.7036
        Power in v 128.6350
        Zero-lag cross correlation 46.3628
        Curly C squared 0.1340

    Calculated from the cross spectral density
        Pover in u 124.7036
        Pover in v 128.6350
        Zero-lag cross correlation 46.3628
        Curly C squared 0.1340
        Coherence 0.1582

In [52]:

    # edapy_09_18, coherence of Reynolds Channel data set
    # columns are days, precip, air_temp, water_temp, salinity, turbidity, chlorophyl
    # compute coherence and plot it

    # set DOCOHEREENCE to 1 to add results from the scipy.signal.coherence() function
    # (which are plotted in red).  Set to 0 to omit them.
    DOCOHEREENCE = 1;

    # this is an example of tupple used here to hold the names of the variables
    names = ('time', 'precipitation', 'air-temp', 'water-temp', 'salinity', 'turbidity', 'chlorophyll');

    # load data
    Sraw = np.genfromtxt('../Data/reynolds_interpolated.txt', delimiter='\t')
    [N, M] = Sraw.shape;
    N = 2*floor(N/2); # truncate to enev length
    S = np.zeros((N,M));
    S[:,:] = Sraw[0:N,:];

    # time in days
    t = np.zeros((N,1));
    t[:,0] = S[:,0];
    Dt=t[1,0]-t[0,0];
    tmin = t[0,0];
    tmax = t[N-1,0];

    # standard frequency setup
    Nf = floor(N/2+1);
    fny = 1/(2*Dt);
    Df = fny/(N/2);
    f = np.zeros((Nf,1)); # positive frequencies
    f[:,0] = np.linspace(0.0,fny,Nf);

    # bandwidth factors
    lowside = 0.75;
    highside =1.25;

    ifig=1;
    fig1 = plt.figure(1,figsize=(16,16));

    # loop over all pairs of columns, computing coherence for each
    for icol in range(1,M):
        for jcol in range(icol+1,M-1):
            # compute cross spectral density and power spectral density
            # no need to normalize, since all normalizations cancel
            ut = eda_cvec( np.fft.fft( S[:,icol] ) );
            vt = eda_cvec( np.fft.fft( S[:,jcol] ) );
            u = np.zeros((Nf,1), dtype=np.complex );
            v = np.zeros((Nf,1), dtype=np.complex );
            u = eda_cvec( ut[0:Nf,0] ); # delete negative frequencies
            v = eda_cvec( vt[0:Nf,0] );
            usv = np.multiply( np.conj(u) , v ); # cross spectral density of u and v
            usu = np.real( np.multiply( np.conj(u) , u ) ); # power spectral density of u
            vsv = np.real( np.multiply( np.conj(v) , v ) ); # power spectral density of v
            # average over band
            usva = np.zeros((Nf,1), dtype=np.complex );
            usua = np.zeros((Nf,1));
            vsva = np.zeros((Nf,1));
            c = np.zeros((Nf,1));
            for i in range(Nf):
                fi = i*Df; # center frequency
                flow = lowside*fi;
                fhigh = highside*fi;
                ilow = floor(flow/Df);
                ihigh = floor(fhigh/Df);
                if( ihigh==ilow):
                    ihigh=ilow+1;
                if (ilow < 0):
                    ilow=0;
                elif (ilow > Nf):
                    ilow=Nf-1;
                if (ihigh < 0):
                    ihigh=0;
                elif (ihigh > Nf):
                    ihigh=Nf;
                tmp1 = np.mean( usv[ilow:ihigh,0] );
                tmp2 = np.mean( usu[ilow:ihigh,0] );
                tmp3 = np.mean( vsv[ilow:ihigh,0] );
                usva[i,0] = tmp1;
                usua[i,0] = tmp2;
                vsva[i,0] = tmp3;

                top = np.real( np.conj(tmp1) * tmp1 );
                bottom = (tmp2 * tmp3);
                if( bottom > 0.0 ):
                    c[i,0] =  top / bottom;
                else:
                    c[i,0] = 0.0;
            # plot coherence
            ax1 = plt.subplot(4,3,ifig);
            plt.axis( [0, 1/5, 0, 1.1] );
            plt.plot(f, c, 'k-' );
            plt.ylabel("coherent of %s vs %s" % (names[icol],names[jcol]));
            plt.xlabel('frequency, cycles per day');
            if( DOCOHEREENCE == 1):
                # result using mscohere() (plotted in red)
                # the black and red traced will not agree
                # perfectly, because the bandwidth parameters
                # cannot be made to match, expecially at low frequencies.
                # figure our various lengths needed by np.coherence()
                N2 = 2;
                while (N2 < N ): # find smallest power of 2 greater than N
                    N2 = 2*N2;
                N3 = floor(N2/16);
                N4 = 3;
                # f2 is frequency axis, C is coherence
                [f2, C2] = sg.coherence( S[:,icol], S[:,jcol],
                        nfft=N2, window='hamming', nperseg=N3, noverlap=N3-N4 );
                f2max = np.max(f2);
                plt.plot( (fny/f2max)*f2, C2, 'r:');
            ifig=ifig+1;

    plt.show();

In [37]:

    # edapy_09_19, boxcar window
    # demonsrate effect of windowing on power spectral density
    # boxcar window

    # standard time/frequency definitions
    N=512;  # length of time serires, must be even
    Dt=1.0; # sampling interval
    t = eda_cvec( Dt*np.linspace(0,N-1,N) ); # time vector
    tmin = t[0,0]; # start time
    tmax=t[N-1,0]; # end time

    # standard set up of frequency
    No2 = floor(N/2);  # integer version of N/2
    fmax=1/(2.0*Dt);   # Nyquist frequency
    Df=fmax/No2;       # frequency sampling
    Nf=No2+1;          # number of non-negative frequencies
    # full (pos and neg) frequency vector
    f=eda_cvec(Df * np.concatenate( ( np.linspace(0,No2,Nf), np.linspace(-No2+1,-1,No2-1) ), axis=0 ));
    Dw=2*pi*Df;        # angular frequency sampling
    w=2*pi*f;          # full (pos and neg) angular frequency vector
    Nw=Nf;             # number of non-negative angular frequencies
    fpos=eda_cvec(Df * np.linspace(0,No2,Nf) );  # vector non-negative frequencies
    wpos=2*pi*fpos;    # vector non-negative angular frequencies

    # boxcar window function, W
    W = np.zeros((N,1));
    No8 = floor(N/8);
    W[0:No8,0]=1.0;
    Wtmax = tmax/8.0;

    # timeseriee
    w0 = 40.0*Dw;
    d = np.sin( w0 * t );

    # windowed timeseries
    Wd = np.multiply( W, d );

    # fourier transform
    dfft = eda_cvec( np.fft.fft(d, axis=0) );
    Wfft = eda_cvec( np.fft.fft(W, axis=0) );
    Wdfft = eda_cvec( np.fft.fft(Wd, axis=0) );

    # keep non-negative frequencies, only
    dfftpos = eda_cvec(dfft[0:Nf,0]);
    Wfftpos = eda_cvec(Wfft[0:Nf,0]);
    Wdfftpos = eda_cvec(Wdfft[0:Nf,0]);

    # power spectral density
    norm = (Dt**2)/tmax;
    ds=norm*np.power(np.abs(dfftpos),2);
    Ws=norm*np.power(np.abs(Wfftpos),2);
    Wds=norm*np.power(np.abs(Wdfftpos),2);

    # time domain plot
    fig1 = plt.figure(1,figsize=(16,16));

    ax1 = plt.subplot(3,2,1);
    plt.axis( [0, tmax, -1.1, 1.1] );
    plt.plot(t, d, 'k-' );
    plt.xlabel('time t (s)');
    plt.ylabel('d(t)');

    ax3 = plt.subplot(3,2,3);
    plt.axis( [0, tmax, -1.1, 1.1] );
    plt.plot(t, W, 'k-' );
    plt.xlabel('time t (s)');
    plt.ylabel('W(t)');

    ax5 = plt.subplot(3,2,5);
    plt.axis( [0, tmax, -1.1, 1.1] );
    plt.plot(t, Wd, 'k-' );
    plt.xlabel('time t (s)');
    plt.ylabel('W(t)*d(t)');

    # frequency domain plot
    ax2 = plt.subplot(3,2,2);
    plt.axis( [0, fmax, 0, 1.1*np.max(np.sqrt(ds))] );
    plt.plot(fpos, np.sqrt(ds), 'k-' );
    plt.xlabel('frequency f (Hz)');
    plt.ylabel('ASD of d');

    ax4 = plt.subplot(3,2,4);
    plt.axis( [0, fmax, 0, 1.1*np.max(np.sqrt(Ws))] );
    plt.plot(fpos, np.sqrt(Ws), 'k-' );
    plt.xlabel('frequency f (Hz)');
    plt.ylabel('ASD of W');

    ax5 = plt.subplot(3,2,6);
    plt.axis( [0, fmax, 0, 1.1*np.max(np.sqrt(Wds))] );
    plt.plot( fpos, np.sqrt(Wds), 'k-' );
    plt.xlabel('frequency f (Hz))');
    plt.ylabel('ASD of Wd');

    plt.show();

In [36]:

    # edapy_09_20, Hamming window
    # demonsrate effect of windowing on power spectral density

    # standard time/frequency definitions
    N=512;
    Dt=1.0;
    t = eda_cvec( Dt*np.linspace(0,N-1,N) ); # time vector
    tmin = t[0,0]; # start time
    tmax=t[N-1,0]; # end time

    # standard set up of frequency
    No2 = floor(N/2);  # integer version of N/2
    fmax=1/(2.0*Dt);   # Nyquist frequency
    Df=fmax/No2;       # frequency sampling
    Nf=No2+1;          # number of non-negative frequencies
    # full (pos and neg) frequency vector
    f=eda_cvec(Df * np.concatenate( ( np.linspace(0,No2,Nf), np.linspace(-No2+1,-1,No2-1) ), axis=0 ));
    Dw=2*pi*Df;        # angular frequency sampling
    w=2*pi*f;          # full (pos and neg) angular frequency vector
    Nw=Nf;             # number of non-negative angular frequencies
    fpos=eda_cvec(Df * np.linspace(0,No2,Nf) );  # vector non-negative frequencies
    wpos=2*pi*fpos;    # vector non-negative angular frequencies

    # Hamming window function, W
    W = np.zeros((N,1));
    Nw = floor(N/8);
    W[0:Nw,0]=0.54-0.46*np.cos(2*pi*np.linspace(0,Nw-1,Nw)/(Nw-1));

    # timeseriee
    w0 = 40*Dw;
    d = np.sin( w0 * t );

    # windowed timeseries
    Wd = np.multiply( W, d );

    # fourier transform
    dfft = np.zeros((N,1), dtype=np.complex );
    Wfft = np.zeros((N,1), dtype=np.complex );
    Wdfft = np.zeros((N,1), dtype=np.complex );
    dfft[:,0] = np.fft.fft(d.ravel());
    Wfft[:,0] = np.fft.fft(W.ravel());
    Wdfft[:,0] = np.fft.fft(Wd.ravel());

    # power spectral density
    norm = (Dt**2)/tmax;
    ds = np.zeros((Nf,1));
    Ws = np.zeros((Nf,1));
    Wds = np.zeros((Nf,1));
    ds[:,0]=norm*np.power(abs(dfft[0:Nf,0]),2);
    Ws[:,0]=norm*np.power(abs(Wfft[0:Nf,0]),2);
    Wds[:,0]=norm*np.power(abs(Wdfft[0:Nf,0]),2);

    # time domain plot
    fig1 = plt.figure(1,figsize=(16,16));

    ax1 = plt.subplot(3,2,1);
    plt.axis( [0, tmax, -1.1, 1.1] );
    plt.plot(t, d, 'k-' );
    plt.xlabel('time t (s)');
    plt.ylabel('d(t)');

    ax3 = plt.subplot(3,2,3);
    plt.axis( [0, tmax, -1.1, 1.1] );
    plt.plot(t, W, 'k-' );
    plt.xlabel('time t (s)');
    plt.ylabel('W(t)');

    ax5 = plt.subplot(3,2,5);
    plt.axis( [0, tmax, -1.1, 1.1] );
    plt.plot(t, Wd, 'k-' );
    plt.xlabel('time t (s)');
    plt.ylabel('W(t)*d(t)');

    # frequency domain plot
    ax2 = plt.subplot(3,2,2);
    plt.axis( [0, fmax, 0, 1.1*np.max(np.sqrt(ds))] );
    plt.plot(fpos, np.sqrt(ds), 'k-' );
    plt.xlabel('frequency f (Hz)');
    plt.ylabel('ASD of d');

    ax4 = plt.subplot(3,2,4);
    plt.axis( [0, fmax, 0, 1.1*np.max(np.sqrt(Ws))] );
    plt.plot(fpos, np.sqrt(Ws), 'k-' );
    plt.xlabel('frequency f (Hz)');
    plt.ylabel('ASD of W');

    ax5 = plt.subplot(3,2,6);
    plt.axis( [0, fmax, 0, 1.1*np.max(np.sqrt(Wds))] );
    plt.plot( fpos, np.sqrt(Wds), 'k-' );
    plt.xlabel('frequency f (Hz))');
    plt.ylabel('ASD of Wd');

    plt.show();

In [39]:

    # edapy_09_21, multiTaper analysis

    # standard time set up
    N=64;
    Dt=1.0;
    t = eda_cvec( Dt*np.linspace(0,N-1,N) ); # time vector
    tmin = t[0,0]; # start time
    tmax=t[N-1,0]; # end time

    # standard set up of frequency
    No2 = floor(N/2);  # integer version of N/2
    fmax=1/(2.0*Dt);   # Nyquist frequency
    Df=fmax/No2;       # frequency sampling
    Nf=No2+1;          # number of non-negative frequencies
    # full (pos and neg) frequency vector
    f=eda_cvec(Df * np.concatenate( ( np.linspace(0,No2,Nf), np.linspace(-No2+1,-1,No2-1) ), axis=0 ));
    Dw=2*pi*Df;        # angular frequency sampling
    w=2*pi*f;          # full (pos and neg) angular frequency vector
    Nw=Nf;             # number of non-negative angular frequencies
    fpos=eda_cvec(Df * np.linspace(0,No2,Nf) );  # vector non-negative frequencies
    wpos=2*pi*fpos;    # vector non-negative angular frequencies

    # choose bandwidth to be +/- some multiple of Dw
    f0 = 2*Df;
    w0 = 2*pi*f0;

    # matrix, M, of sinc functions
    M = np.zeros((N,N));
    for m in range(N):
        for n in range(N):
            M[m,n] = 2*w0*np.sinc( w0*(m-n)*Dt / pi);

    # solve algebric eigenvalue problem
    # eig sorts in ascending order
    [lam, V] = la.eigh(M);
    print('6 largest eigenvalues:');
    for m in range(N-5,N):
        print("    %d: %f" % (m, lam[m]) );

    # four largest eigenvectors
    w1 = eda_cvec( V[:,N-1] );
    w2 = eda_cvec( V[:,N-2] );
    w3 = eda_cvec( V[:,N-3] );
    w4 = eda_cvec( V[:,N-4] );

    # power spectral density of window functions
    s1 = np.zeros((Nf,1));
    s2 = np.zeros((Nf,1));
    s3 = np.zeros((Nf,1));
    s4 = np.zeros((Nf,1));
    norm = (Dt**2)/tmax;
    temp1 = norm * np.power(np.abs(np.fft.fft(w1,axis=0)),2);
    temp2 = norm * np.power(np.abs(np.fft.fft(w2,axis=0)),2);
    temp3 = norm * np.power(np.abs(np.fft.fft(w3,axis=0)),2);
    temp4 = norm * np.power(np.abs(np.fft.fft(w4,axis=0)),2);
    s1[:,0] = temp1[0:Nf,0];
    s2[:,0] = temp2[0:Nf,0];
    s3[:,0] = temp3[0:Nf,0];
    s4[:,0] = temp4[0:Nf,0];

    # time domain plot
    fig1 = plt.figure(1,figsize=(12,14));
    winmin = -0.25;
    winmax = -winmin;

    ax1 = plt.subplot(4,1,1);
    plt.axis( [0, tmax, winmin, winmax] );
    plt.plot(t, w1, 'k-' );
    plt.xlabel('time t (s)');
    plt.ylabel('w1(t)');
    plt.title('tapers');

    ax2 = plt.subplot(4,1,2);
    plt.axis( [0, tmax, winmin, winmax] );
    plt.plot(t, w2, 'k-' );
    plt.xlabel('time t (s)');
    plt.ylabel('w1(t)');

    ax3 = plt.subplot(4,1,3);
    plt.axis( [0, tmax, winmin, winmax] );
    plt.plot(t, w3, 'k-' );
    plt.xlabel('time t (s)');
    plt.ylabel('w1(t)');

    ax4 = plt.subplot(4,1,4);
    plt.axis( [0, tmax, winmin, winmax] );
    plt.plot(t, w4, 'k-' );
    plt.xlabel('time t (s)');
    plt.ylabel('w1(t)');

    plt.show();

    # plot amplitude spectral density of window functions
    fig2 = plt.figure(2,figsize=(12,14));

    ax1 = plt.subplot(4,1,1);
    plt.axis( [0, fmax, 0, 1.1*np.max(np.sqrt(s1))] );
    plt.plot(fpos,s1, 'k-' );
    plt.plot( [f0, f0], [0, np.max(np.sqrt(s1))], 'k:' );
    plt.xlabel('frequency f (Hz)');
    plt.ylabel('w1(f)');
    plt.title('ASD of tapers');

    ax2 = plt.subplot(4,1,2);
    plt.axis( [0, fmax, 0, 1.1*np.max(np.sqrt(s2))] );
    plt.plot(fpos,s2, 'k-' );
    plt.plot( [f0, f0], [0, np.max(np.sqrt(s3))], 'k:' );
    plt.xlabel('frequency f (Hz)');
    plt.ylabel('w1(f)');

    ax3 = plt.subplot(4,1,3);
    plt.axis( [0, fmax, 0, 1.1*np.max(np.sqrt(s3))] );
    plt.plot(fpos,s3, 'k-' );
    plt.plot( [f0, f0], [0, np.max(np.sqrt(s3))], 'k:' );
    plt.xlabel('frequency f (Hz)');
    plt.ylabel('w1(f)');

    ax4 = plt.subplot(4,1,4);
    plt.axis( [0, fmax, 0, 1.1*np.max(np.sqrt(s4))] );
    plt.plot(fpos,s4, 'k-' );
    plt.plot( [f0, f0], [0, np.max(np.sqrt(s4))], 'k:' );
    plt.xlabel('frequency f (Hz)');
    plt.ylabel('w1(f)');

    plt.show();

    # new, longer time interval, with same sampling interval
    N2=512;                 # length of longer interval
    Dt2 = 1.0;              # sampling interval
    N2o2 = floor(N2/2);     # interger version of N2/2
    Nf2=N2o2+1;             # number of non-negative frequencies
    t2 = eda_cvec( Dt*np.linspace(0,N2-1,N2) ); # time vector
    tmin2 = t2[0,0];        # start time
    tmax2 = t2[N2-1,0];     # end time
    fmax2=1/(2.0*Dt2);      # nyquist freqency
    Df2=fmax2/N2o2;         # frequency sampling
    Nf2=N2o2+1;             # number of non-negative frequencies
    # full (pos and neg) vector of freqeuncies
    f2=eda_cvec(Df2 * np.concatenate( ( np.linspace(0,N2o2,Nf2), np.linspace(-N2o2+1,-1,N2o2-1) ), axis=0 ));
    Dw2=2*pi*Df2; # angular frequency sampling;
    # full (pos and neg) vector of angular freqeuncies
    w_2 = 2 * pi * f2;
    Nw2=Nf2; # number of angular frequencies
    fpos2=eda_cvec(Df2 * np.linspace(0,N2o2,Nf2)); # vector of non-negative frequencies
    wpos2=2*pi*fpos2; # vector of non-negative angular frequencies

    # create time series and window it
    ff= fmax/3.3;
    ww= 2*pi*ff;
    d=np.cos(ww*t2);

    # w0 is a uniform window function
    w0 = np.ones((N,1));

    # we is a bounch of zero that will be padded onto the
    # end of a short time series to make a long one
    we = np.zeros((N2-N,1));

    # zero-padded window functions
    w0we = eda_cvec( np.concatenate( (w0,we), axis=0 ));
    w1we = eda_cvec( np.concatenate( (w1,we), axis=0 ));
    w2we = eda_cvec( np.concatenate( (w2,we), axis=0 ));
    w3we = eda_cvec( np.concatenate( (w3,we), axis=0 ));
    w4we = eda_cvec( np.concatenate( (w4,we), axis=0 ));

    # windowed data
    dw0 = eda_cvec( np.multiply( w0we, d ) );
    dw1 = eda_cvec( np.multiply( w1we, d ) );
    dw2 = eda_cvec( np.multiply( w2we, d ) );
    dw3 = eda_cvec( np.multiply( w3we, d ) );
    dw4 = eda_cvec( np.multiply( w4we, d ) );

    # plot windowed time series
    fig2 = plt.figure(2,figsize=(12,14));

    ax1 = plt.subplot(5,1,1);
    plt.axis( [tmin2, tmax2, -1.1, 1.1] );
    plt.plot(t2, dw0, 'k-' );
    plt.xlabel('time t (s)');
    plt.ylabel('d(t)');
    plt.title('windowed time series');

    ax2 = plt.subplot(5,1,2);
    plt.axis( [tmin2, tmax2, winmin, winmax] );
    plt.plot(t2, dw1, 'k-' );
    plt.xlabel('time t (s)');
    plt.ylabel('d*w1')

    ax3 = plt.subplot(5,1,3);
    plt.axis( [tmin2, tmax2, winmin, winmax] );
    plt.plot(t2, dw2, 'k-' );
    plt.xlabel('time t (s)');
    plt.ylabel('d*w2')

    ax4 = plt.subplot(5,1,4);
    plt.axis( [tmin2, tmax2, winmin, winmax] );
    plt.plot(t2, dw3, 'k-' );
    plt.xlabel('time t (s)');
    plt.ylabel('d*w3')

    ax5 = plt.subplot(5,1,5);
    plt.axis( [tmin2, tmax2, winmin, winmax] );
    plt.plot(t2, dw4, 'k-' );
    plt.xlabel('time t (s)');
    plt.ylabel('d*w4')

    plt.show();

    # power spectral density of time series
    norm = (Dt2**2)/tmax2;
    temp = norm * np.power( np.abs( eda_cvec(np.fft.fft(d,axis=0))),2);  # power spectrum
    ds = eda_cvec( temp[0:Nf2,0] ); # delete negative frequencies

    # power spectral density of windowed time series
    norm = (Dt2**2)/tmax;
    temp0 = norm * np.power( np.abs( eda_cvec(np.fft.fft(dw0,axis=0))),2); # fourier transform and etc.
    temp1 = norm * np.power( np.abs( eda_cvec(np.fft.fft(dw1,axis=0))),2);
    temp2 = norm * np.power( np.abs( eda_cvec(np.fft.fft(dw2,axis=0))),2);
    temp3 = norm * np.power( np.abs( eda_cvec(np.fft.fft(dw3,axis=0))),2);
    temp4 = norm * np.power( np.abs( eda_cvec(np.fft.fft(dw4,axis=0))),2);
    dw0s = eda_cvec( temp0[0:Nf2,0] ); # keep only non-negative frequencies
    dw1s = eda_cvec( temp1[0:Nf2,0] );
    dw2s = eda_cvec( temp2[0:Nf2,0] );
    dw3s = eda_cvec( temp3[0:Nf2,0] );
    dw4s = eda_cvec( temp4[0:Nf2,0] );

    # average power spectral densities
    dws = (1.0/3.0)*( dw1s + dw2s + dw3s  );

    # plot amplitude spectral density of window functions
    fig4 = plt.figure(2,figsize=(12,14));

    ax1 = plt.subplot(6,1,1);
    plt.axis( [0, fmax2, 0, 1.1*np.max(np.sqrt(ds))] );
    plt.plot(fpos2,np.sqrt(ds), 'k-' );
    plt.plot( [ff, ff], [0, np.max(np.sqrt(ds))], 'k:' );
    plt.xlabel('frequency f (Hz)');
    plt.ylabel('ASD of d');
    plt.ylabel('ASD of taperd version of d(t)');

    ax2 = plt.subplot(6,1,2);
    plt.axis( [0, fmax2, 0, 1.1*np.max(np.sqrt(dw0s))] );
    plt.plot(fpos2,np.sqrt(dw0s), 'k-' );
    plt.plot( [ff, ff], [0, np.max(np.sqrt(dw0s))], 'k:' );
    plt.xlabel('frequency f (Hz)');
    plt.ylabel('ASD of d*w0');

    ax3 = plt.subplot(6,1,3);
    plt.axis( [0, fmax2, 0, 1.1*np.max(np.sqrt(dw1s))] );
    plt.plot(fpos2,np.sqrt(dw1s), 'k-' );
    plt.plot( [ff, ff], [0, np.max(np.sqrt(dw1s))], 'k:' );
    plt.xlabel('frequency f (Hz)');
    plt.ylabel('ASD of d*w1');

    ax4 = plt.subplot(6,1,4);
    plt.axis( [0, fmax2, 0, 1.1*np.max(np.sqrt(dw2s))] );
    plt.plot(fpos2,np.sqrt(dw2s), 'k-' );
    plt.plot( [ff, ff], [0, np.max(np.sqrt(dw2s))], 'k:' );
    plt.xlabel('frequency f (Hz)');
    plt.ylabel('ASD of d*w2');

    ax5 = plt.subplot(6,1,5);
    plt.axis( [0, fmax2, 0, 1.1*np.max(np.sqrt(dw3s))] );
    plt.plot(fpos2,np.sqrt(dw3s), 'k-' );
    plt.plot( [ff, ff], [0, np.max(np.sqrt(dw3s))], 'k:' );
    plt.xlabel('frequency f (Hz)');
    plt.ylabel('ASD of d*w3');

    ax6 = plt.subplot(6,1,6);
    plt.axis( [0, fmax2, 0, 1.1*np.max(np.sqrt(dws))] );
    plt.plot(fpos2,np.sqrt(dws), 'k-' );
    plt.plot( [ff, ff], [0, np.max(np.sqrt(dws))], 'k:' );
    plt.xlabel('frequency f (Hz)');
    plt.ylabel('mean ASD)');

    plt.show();

    6 largest eigenvalues:
        59: 1.724369
        60: 4.536603
        61: 6.029180
        62: 6.268016
        63: 6.282832

In [ ]: