edapy_07

In [1]:

    # edapy_07_00 clear all variables and import vatious modules
    %reset -f
    import os
    from datetime import date
    from math import exp, pi, sin, sqrt, floor, ceil
    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( N==0 ):
                continue; # skip vector of zero length
            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;e*h],
        # where G is a topplitz matrix built from the filter g
        # Note that "H" must be sparse and called "edaHsparse" and the
        # thefilter "g" bust be called edafilterg and must be a column vector
        # and that the parameter "e" 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): # convert to column-vector
            vv = np.zeros((s[0],1));
            vv[:,0] = v;
        else:
            vv=v;
        N, i = np.shape(edafilterg);
        temp1 = np.zeros((N+M-1,1));
        temp1[:,0] = np.convolve(edafilterg.ravel(),vv.ravel());
        a = np.zeros((N,1));
        a[:,0] = temp1[0:N,0];
        b = eda_e*edaHsparse * v;
        temp2 = np.zeros((N+N-1,1));
        temp2[:,0] = np.convolve(
            (np.flipud(edafilterg)).ravel(),a.ravel());
        a2 = temp2[N-1:N+M-1,0];
        b2 = eda_e*edaHsparse.transpose()*b;
        # a2+b2 = FT F v = GT G v + (e^**2) 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);

In [124]:

    # edapy_07_01: least squares estimation of filter input

    # example using impulse response of a heat-producing layer
    # with solution by damped least squares
    # Note: symbol q for theta, temperture
    # This is the non-sparse implementation.

    # generic time setup, except make two time arrays.
    # one in seconds for calulations, one in days for plotting
    M=1024;    # length of time series
    N=M;       # number of model parameters equal number of data
    Dtdays = 0.5; # sampling interval in days
    tdays = eda_cvec(np.linspace(0,Dtdays*(N-1),N)); # time vector (days)
    Dtseconds = Dtdays*3600*24; # samping interval in seconds
    tseconds = Dtseconds*tdays;  # time vector (seconds)

    # materal constants
    rho = 1500; # density, kg/m3
    cp = 800; # heat capacity, J / kg-K
    kappa = 1.25e-6; # thermal diffusivity, m2/s
    z = 1.0; # depth of observation, m

    # impulse response, which is a complicated analytic formula
    # time vector without zero
    tsecondswo0 = eda_cvec( tseconds[1:M,0:1] );
    P1 = (1.0/(rho*cp));
    P2 = (Dtseconds/sqrt(2*pi));
    P3 = np.reciprocal( np.sqrt(2*kappa*tsecondswo0));
    P4 = np.exp( -0.5*(z**2)*np.reciprocal(2*kappa*tsecondswo0) );
    g = eda_cvec( 0.0, P1 * P2 * np.multiply(P3 ,P4) );

    # true model parameters, heat production
    t1 = Dtdays*N/4;     # time of heat pulse 2
    t2 = Dtdays*7*N/16;  # time of heat pulse 2
    sd = Dtdays*N/256;   # width of pulses
    Ah=10;               # amplitude of pulses
    P5 = np.exp(-np.power(tdays-t1,2) / (2*sd**2) ); # gaussian pulse 1
    P6 = np.exp(-np.power(tdays-t2,2) / (2*sd**2) ); # gaussian pulse 2
    htrue = Ah*P5 + 0.5*Ah*P6;  # sum pulses

    # predict true temperature by convolving g and htrue
    # note that np.convolve() expects Nx0 arrays, so the
    # Nx1 columnn-vectors must be "raveled" before call
    # to np.convolve, and then turned back into a column-vector
    # by a call to eda_cvec()
    qtrue = eda_cvec(np.convolve(g.ravel(),htrue.ravel()));
    qtrue = np.copy(qtrue[0:N,0:1]);# keep only first N values

    # plot impulse response
    fig1 = plt.figure(1,figsize=(8,4));
    ax1 = plt.subplot(3,1,1);
    plt.axis([tdays[0,0], tdays[M-1,0], 0, 1.1*np.max(g)] );
    plt.plot(tdays[0:M],g,'k-');
    plt.xlabel('time t, days');
    plt.ylabel('g(t)');

    # 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:M,0:1],htrue,'k-');
    plt.xlabel('time t, days');
    plt.ylabel('htrue(t)');

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

    # Toeplitz matrix version of impulse response
    Gcol = eda_cvec( g, np.zeros((N-M,1)));
    Grow = eda_cvec( g[0,0], np.zeros((M-1,1)));
    G = la.toeplitz( Gcol, Grow.T );

    # temperature = G times heat_production
    qtrue2 = np.matmul(G,htrue);

    # create simulated observations, qobs = qtrue+noise
    sigmad=0.001*np.amax(qtrue);  # noise is 0.1% amplitude of signal
    qobs=qtrue+np.random.normal(0,sigmad,(N,1));

    # plot simulated observations
    fig2 = plt.figure(2,figsize=(8,4));
    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,0:1],qobs,'k-');
    plt.xlabel('time t, days');
    plt.ylabel('qobs(t)');

    # solve for h given qobs and G using damped least squares
    GTG = np.matmul(G.T,G);       # G-transpose times G
    GTGmax = np.max(np.abs(GTG)); # maximum value of G
    e2=1e-4*GTGmax;               # damping is a fraction of GTG
    hest1 = la.solve(GTG+e2*np.identity(M), np.matmul(G.T,qobs) );

    # 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,0:1],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,0:1],hest1,'k-');
    plt.xlabel('time t, days');
    plt.ylabel('hest1(t)');
    plt.show();

In [125]:

    # edapy_07_02: estimation of filter input, with smoothing

    # using impulse repsonse of a heat-producing layer
    # with solution via generalized least squares
    # with the prior information of smoothness
    # Note: symbol q for theta, temperture.
    # This is the non-sparse implementation.

    # generic time setup, except make two time arrays.
    # one in seconds for calulations, one in days for plotting
    M=1024;    # length of time series
    N=M;       # number of model parameters equal number of data
    Dtdays = 0.5; # sampling interval in days
    tdays = eda_cvec(np.linspace(0,Dtdays*(N-1),N)); # time vector (days)
    Dtseconds = Dtdays*3600*24; # samping interval in seconds
    tseconds = Dtseconds*tdays;  # time vector (seconds)

    # materal constants
    rho = 1500.0; # density, kg/m3
    cp = 800.0;   # heat capacity, J / kg-K
    kappa = 1.25e-6; # thermal diffusivity, m2/s
    z = 1.0; # depth of observation, m

    # impulse response, which is a analysic formula that's
    # fairly complicated, so break it into three pieces,
    # P1, P2 and P3
    # impulse respone is g, but time=0 will be handled separately
    tsecondswo0 = eda_cvec( tseconds[1:M,0] ); # time vector without zero
    P1 = (1.0/(rho*cp));
    P2 = (Dtseconds/sqrt(2*pi));
    P3 = np.reciprocal( np.sqrt(2*kappa*tsecondswo0));
    P4 = np.exp( -0.5*(z**2)*np.reciprocal(2*kappa*tsecondswo0) );
    g = eda_cvec( 0.0, P1 * P2 * np.multiply(P3 ,P4) );

    #true model parameters, heat production
    t1 = Dtdays*N/4;     # time of heat pulse 1
    t2 = Dtdays*7*N/16;  # time of heat pulse 2
    sd = Dtdays*N/256;   # width of pulses
    Ah = 10.0;           # amplitude of pulses
    P5 = np.exp(-np.power(tdays-t1,2) / (2*sd**2) ); # gaussian pulse 1
    P6 = np.exp(-np.power(tdays-t2,2) / (2*sd**2) ); # gaussian pulse 2
    htrue = Ah*P5 + 0.5*Ah*P6;  # sum pulses

    # Toeplitz matrix version of impulse response
    Gcol = eda_cvec( g, np.zeros((N-M,1)));
    Grow = eda_cvec( g[0,0], np.zeros((M-1,1)));
    G = la.toeplitz( Gcol, Grow.T );

    # temperature = G times heat_production
    qtrue2 = np.matmul(G,htrue);

    # create simulated observations, qobs = qtrue+noise
    # noise leve is 1% of maximum temperature
    sigmad=0.01*np.max(np.abs(qtrue));
    qobs=qtrue+np.random.normal(0,sigmad,(N,1));

    # solve for h given qobs and g using damped least squared
    GTG = np.matmul(G.T,G);  # G_transpose time G
    GTGmax = np.amax(GTG);   # maximum value of GTG
    e2=1e-4*GTGmax;  # damping is fraction of GTG
    # damped least squares solution
    hest2 = la.solve(GTG+e2*np.identity(M), np.matmul(G.T,qobs) );

    # generalized least squares with F = [G; H]  and  f = [qobs, 0]
    # second derivative smoothing in interior of interval; nothing on ends
    # H is a Toeplitz matrix with diaginals 1, -2 and 1.  Note
    # that he factor of 1/(Dt^2) is omitted because the r.h.s. is zero.
    # The damping cnsnstant e = sigma_h/sigma_d determines the
    # relative weighting of Gm=d and Hm=h. It is being tuned
    # here by trail and error to achieve the desired degree of smoothness
    e=10*np.max(np.abs(G));# e is fraction of maximum value of G
    Hrow = eda_cvec( 1.0, -2.0,  1.0, np.zeros((M-3,1)));
    Hcol = eda_cvec( 1.0, np.zeros((M-3,1)));
    H = la.toeplitz( Hcol, Hrow.T ); # construct H
    h = np.zeros((M-2,1)); # since prior info is that 2nd derivative is zero, h=0
    F = np.concatenate( (G,e*H), axis=0 );  # combine G and H into F
    f = np.concatenate( (qobs,e*h), axis=0 ); # combine qobs and h into f

    # deneralized least squares solution
    FTF = np.matmul(F.T,F); # F_transpose times F
    FTf = np.matmul(F.T,f); # F_transpose times f
    hest3 = la.solve(FTF, FTf);  # generalized least squares solution

    # 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,0:1],qobs,'k-');
    plt.xlabel('time t, days');
    plt.ylabel('htrue(t)');
    plt.title('Damped Least Squares');

    # 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,0:1],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,0:1],hest2,'k-');
    plt.xlabel('time t, days');
    plt.ylabel('hest(t)');
    plt.show();

    # figure 2
    fig2 = plt.figure(2,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,0:1],qobs,'k-');
    plt.xlabel('time t, days');
    plt.ylabel('htrue(t)');
    plt.title('Generalized Least Squares with 2nd Derivative Smoothing');

    # 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,0:1],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,0:1],hest3,'k-');
    plt.xlabel('time t, days');
    plt.ylabel('hest1(t)');
    plt.show();

In [2]:

    # edapy_07_03: sparse, least sq estimation of filter input

    # using impulse repsonse of a heat-producing layer
    # with solution via generalized least squares
    # with the prior information of smoothness
    # Note: symbol q for theta, temperture
    # This is the sparse implementation.

    # generic time setup, except make two time arrays.
    # one in seconds for calulations, one in days for plotting
    M=1024;    # length of time series
    N=M;       # number of model parameters equal number of data
    Dtdays = 0.5; # sampling interval in days
    tdays = eda_cvec(np.linspace(0,Dtdays*(N-1),N)); # time vector (days)
    Dtseconds = Dtdays*3600*24; # samping interval in seconds
    tseconds = Dtseconds*tdays;  # time vector (seconds)

    # materal constants
    rho = 1500.0;    # density, kg/m3
    cp = 800.0;      # heat capacity, J / kg-K
    kappa = 1.25e-6; # thermal diffusivity, m2/s
    z = 1.0;         # depth of observation, m

    # impulse response, which is a analysic formula that's
    # fairly complicated, so break it into three pieces,
    # P1, P2 and P3
    # impulse respone is g, but time=0 will be handled separately
    tsecondswo0 = eda_cvec( tseconds[1:M,0] ); # time vector without zero
    P1 = (1.0/(rho*cp));
    P2 = (Dtseconds/sqrt(2*pi));
    P3 = np.reciprocal( np.sqrt(2*kappa*tsecondswo0));
    P4 = np.exp( -0.5*(z**2)*np.reciprocal(2*kappa*tsecondswo0) );
    g = eda_cvec( 0.0, P1 * P2 * np.multiply(P3 ,P4) );
    edafilterg = g;  # edafilterg is a global variable defined in GLSFilterMul()

    # true heat production, sum of two pulses
    t1 = Dtdays*N/4;      # time of first pulse
    t2 = Dtdays*7*N/16;   # time of second pulse
    sd = Dtdays*N/256;    # variance of pulse
    Ah=10;                # amplitude of pulse
    # sum of two gaussian pulses
    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
    # note that np.convolve() expects Nx0 arrays, so the
    # Nx1 columnn-vectors must be "raveled" before call
    # to np.convolve, and then turned back into a column-vector
    # by a call to eda_cvec()
    qtrue = eda_cvec(np.convolve(g.ravel(), htrue.ravel()));
    qtrue = eda_cvec(qtrue[0:N]); # keep only first N values

    # create noisy simulated observations, qobs = qtrue+noise
    sigmad=0.01*np.max(np.abs(qtrue));  # fraction of signal
    qobs=qtrue+np.random.normal(0,sigmad,(N,1))

    # for comparison purposes, solve for h given qobs and G using damped least squares

    # Toeplitz matrix version of impulse response
    Gcol = eda_cvec( g, np.zeros((N-M,1)));
    Grow = eda_cvec( g[0,0], np.zeros((M-1,1)));
    G = la.toeplitz( Gcol, Grow.T );

    # predict true temperature from true heat prodcution
    qtrue2 = np.matmul(G,htrue);
    # solve with damped least squares
    GTG = np.matmul(G.T,G); # G_transpose times G
    GTGmax = np.max(np.abs(GTG));  # maximum value of GTG
    e2=1e-4*GTGmax; # damping is fraction of maximum value of GTG
    # damped least squares ssolution
    hest2 = la.solve(GTG+e2*np.identity(M), np.matmul(G.T,qobs) );

    # generalized leasr squares with F = [G; H]  and  f = [qobs, 0]
    # second derivative smoothing in interior of interval; nothing on ends
    eda_e=10.0*np.amax(G); # note damoing is a global variable

    # Implement prior information of smoothness (2nd derivative = 0).
    # Note that a factor of 1/(Dt**2) is omitted because the r.h.s. is zero.
    # H is Toeplitz, so use sparse.spdiags() to create a sparse version of it.
    # The [0,1,2] means the non-zero diagonals are the main diagonal
    # (named 0) and two immediateley the right of it (the 1 and 2), The
    # (D,-2.0*D,D) are the corresponding values of the diagonals.
    # The sparse matrix is K by M.
    K=M-2;
    D = np.ones((1,M));
    edaHsparse = sparse.spdiags(
        np.concatenate((D,-2.0*D,D),axis=0), [0, 1, 2], K, M);
    h = np.zeros((M-2,1));

    # set up right hand side, F'f = GT qobs + e HT h
    FTf = GLSFilterMulRHS(edafilterg,qobs,eda_e,edaHsparse,h,M)

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

    # solve for estimated model parameters using congugate gradient algrorithm
    hest3 = np.zeros((M,1));
    q=las.cg(LO,FTf,tol=1e-12, maxiter=(3*(N+M)+100));
    # note that q is a tupple, the first element of which is the solution
    hest3 = eda_cvec(q[0]);

    # 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('htrue(t)');
    plt.title('Damped Least Squares');

    # 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],hest2,'k-');
    plt.xlabel('time t, days');
    plt.ylabel('hest(t)');
    plt.show();

    # figure 2
    fig2 = plt.figure(2,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('htrue(t)');
    plt.title('Generalized Least Squares with 2nd Derivative Smoothing');

    # 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('hest1(t)');
    plt.show();

In [6]:

    # edapy_07_04: filter estimation problem

    # example of estimating a filter r that predicts river
    # discharge d from precipitation data g; that is, g*r=d

    # in this example, the filter length is chosen to be
    # M<<N, to embody the idea that the response of the river
    # (meaning the length M of the filter r) is short

    M=25; # length of filter r

    # load merged discharge and precipitation data
    # column 1: time in inter count of days starting 1/1/2002
    # column 2: Hudson River discharge in cfs
    # column 3: Albany NY precipitation in inches
    D = np.genfromtxt('../Data/precip_discharge_data.txt', delimiter='\t')
    [Nlong, K]=D.shape;
    tlong = eda_cvec(D[:,0]); # time in days
    dlong = eda_cvec(D[:,1]/35.3146); # discharge in cubic meters per second
    glong = eda_cvec(D[:,2]*25.4) # precipitation, im millimeters pper day;
    Dt = tlong[1,0]-tlong[0,0]; # sampling interval

    # in this example, we use only a 101 day portion of the
    # dataset, mainly so that when we plot it, storm events
    # are clearly visible

    # select small segment of time, discharge and precip data
    N=101;
    istart = 900-1;
    t = eda_cvec( tlong[istart:istart+N,0] );
    d = eda_cvec( dlong[istart:istart+N,0] );
    g = eda_cvec( glong[istart:istart+N,0] );
    tmin = t[0,0];    # start time
    tmax = t[N-1,0];  # end time
    edafilterg = g; # edafilterg is a global variable defined in GLSFilterMul()

    # figure 1: discharge and precipitation data
    fig1 = plt.figure(1,figsize=(8,4));

    # plot observations
    ax1 = plt.subplot(2,1,1);
    plt.axis([tmin, tmax, -20, 60 ] );
    plt.plot(t,g,'k-');
    plt.xlabel('time t, days');
    plt.ylabel('precipitation (mm/day)');

    ax1 = plt.subplot(2,1,2);
    plt.axis([tmin, tmax, -1000, 2000 ] );
    plt.plot(t,d,'k-');
    plt.xlabel('time t, days');
    plt.ylabel('discharge (m3/s)');
    plt.show();

    # Given the filter equation g*r=d, solve for r given d and g
    # using genralized least squares

    # Note that the convolution relationship g*r=d is equivalent
    # to the matrix equation G r = d where G is the (g*) operator,
    # Since we are using filterfun(), we do not need to construct
    # G.  We do, however, need to implement the prior information
    # equation Hr=h.

    # This code can toggle between the prior information of r
    # being small, its first derivative being small and its second
    # derivative being small, depending on whether DERIVATIVE is 0,
    # 1 or 2. (However, the results turn out not to be especially
    # sensitive to the choice).
    DERIVATIVE=1;
    if( DERIVATIVE==2 ):
        # The prior information equation H*r = h is implemented
        # for the smoothness condition of the second derivative
        # being close to zero
        eda_e=2*np.max(np.abs(g)); # the coefficient was manually tuned
                                   # to produce a reasonably smooth filter
        K=M-2;
        L=N+K;
        # second derivative operator; no end conditions
        D = np.ones((1,M));
        edaHsparse = sparse.spdiags( np.concatenate((D,-2.0*D,D),axis=0), [0, 1, 2], K, M);
        h = np.zeros((K,1));
    elif (DERIVATIVE==1):
        # The prior information equation H*r = h is implemented
        # for the smoothness condition of the first derivative
        # being close to zero
        eda_e=4.0*np.max(np.abs(g)); # the coefficient was manually tuned
                                     # % to produce a reasonably smooth filter
        K=M-1;
        L=N+K;
        D = np.ones((1,M));
        edaHsparse = sparse.spdiags( np.concatenate((D,-D),axis=0), [0, 1], K, M);
        h = np.zeros((K,1));
    else: # DERIVATIVE == 0
        # The prior information equation H*r = h is implemented
        # for the smallness condition of the filter is close to zero
        eda_e=0.1*np.max(np.abs(g)); # the coefficient was manually tuned
                                     # to produce a reasonably small filter
        K=M;
        L=N+K;
        D = np.ones((1,M));
        edaHsparse = sparse.spdiags( D, [0], K, M);
        h = np.zeros((K,1));

    # set up right hand side, F'f = GT qobs + e HT h
    FTf = GLSFilterMulRHS(edafilterg,d,eda_e,edaHsparse,h,M);

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

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

    # filter time axis
    tr = eda_cvec( Dt*np.linspace(0,M-1,M) );

    # predicted discharge, convolve and truncate to length N
    dpre = eda_cvec( np.convolve(g.ravel(),r.ravel()) );
    dpre = dpre[0:N,0:1];
    e = d-dpre;
    E = np.matmul(e.T,e);
    print("total error E: %e" % (E));

    # figure 2: filter and predicted discharge
    fig2 = plt.figure(2,figsize=(8,4));

    # plot observations
    ax1 = plt.subplot(2,1,1);
    plt.axis([0, M, -5, 10 ] );
    plt.plot(tr,r,'k-');
    plt.plot([tr[0,0], tr[M-1,0]],[0,0],'k:');
    plt.xlabel('time t, days');
    plt.ylabel('filter r(t)');

    ax1 = plt.subplot(2,1,2);
    plt.axis([tmin, tmax, -1000, 2000 ] );
    plt.plot(t,d,'k:');
    plt.plot(t,dpre,'k-');
    plt.xlabel('time t, days');
    plt.ylabel('discharge (m3/s)');
    plt.show();

    total error E: 2.316860e+06

In [128]:

    # edapy_07_05: prediction error filter (pef)

    # this examples uses the for Neuse Riverr Hydrograph
    # this code uses generalized least squares to implement
    # the prior information that pef(1)=(-1).  Two solutions
    # are calculated, one via standard matrix algebra,
    # (FTF)\(FTf) and the other via the biconjugate gradient method

    # read in data
    D = np.genfromtxt('../Data/neuse.txt', delimiter='\t')
    [N, K]=D.shape;
    t = eda_cvec( D[:,0] ); # time in days
    d = eda_cvec( D[:,1] ); # discharge in cubic meters per second
    Dt = t[1,0]-t[0,0];     # sampling interval
    tmin = t[0,0];          # start time
    tmax = t[N-1,0];        # end time

    # sizes of various arrays
    M=100; # length of pef is user tunable
    K=1;   # one piece of prior information
    L=N+K; # total number of equations

    # for a prediction error filter, pef, the data is the filter, g
    g = d;
    edafilterg = d; # edafilterg is a global variable defined in GLSFilterMul()

    # prior information that first element of the pef is (-1)
    eda_e=(10^1)*np.amax(g); # large weight
    H=np.zeros((K,M));   # the pror information equation, Hm=h, such that m(1)=(-1);
    H[0,0]=1.0;
    h=eda_cvec( (-1.0) );

    # Data Kernel (needed only for standard matrix algebra calculation
    # data kernel is Toeoplitz, so use standard setup
    # Toeplitz matrix version of impulse response
    Gcol = g;
    Grow = eda_cvec( g[0,0], np.zeros((M-1,1)));
    G = la.toeplitz( Gcol, Grow.T );
    # for pef, Gm=rhs=0

    # F needed only for standard matrix algebra calculation
    F = np.concatenate( (G, eda_e*H), axis=0 );
    rhs = np.zeros((N,1));
    f = np.concatenate( (rhs, eda_e*h), axis=0 );

    # standard GLS solution
    FTF = np.matmul( F.T, F );  # F_transpose times F
    FTf = np.matmul( F.T, f );  # F_transpose times f
    pef1 = la.solve(FTF, FTf ); # least squares solution

    # now set up for biconjugate gradient method

    # Hm=h
    eda_e=(10^1)*np.amax(g); # large weight
    # edaHsparse is global variable defined in GLSFilterMul()
    edaHsparse = sparse.diags([1.0], [0], shape=(1, M));
    h=np.zeros((1,1));
    h[0,0]=(-1.0);

    # FT f can be done analytucally:
    # FTf = GT 0 + e HT e h = 0 + [e 0 0 ... 0]T [-e] = [-e^2 0 0 ... 0]T
    FTf = np.zeros((M,1));
    FTf[0,0] = -(eda_e**2);

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

    # solve for estimated model parameters using congugate gradient algrorithm
    r = np.zeros((M,1));
    # conjugate gradient solver returns soln vector as first element of a tupple
    q=las.cg(LO,FTf,tol=1e-12, maxiter=(3*(N+M)+100));
    pef2 = eda_cvec(q[0]);

    # compute prediction error by colvolving g with pef, then truncate to length N
    temp = np.convolve(pef2.ravel(),g.ravel());
    perror = eda_cvec( temp[0:N] );

    # short time axis
    tpef = eda_cvec( t[0:M,0] );

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

    # plot prediction error filters
    ax1 = plt.subplot(2,1,1);
    plt.axis([ 0, Dt*(M-1), np.min(pef1), np.max(pef1)] );
    plt.plot(tpef,pef1,'k-');
    plt.xlabel('time t, days');
    plt.ylabel('pef(t) (std)');

    ax1 = plt.subplot(2,1,2);
    plt.axis([ 0, Dt*(M-1), np.min(pef1), np.max(pef1)] );
    plt.plot(tpef,pef2,'k-');
    plt.xlabel('time t, days');
    plt.ylabel('pef(t) (std)');
    plt.show();

    # plot hydrograph data and prediction error for first year
    fig2 = plt.figure(2,figsize=(8,4));

    Np = 365;
    ax1 = plt.subplot(2,1,1);
    plt.axis([tmin, tmin+Dt*Np, -5000, 15000]);
    plt.plot(t[0:Np,0],d[0:Np,0],'k-');
    plt.xlabel('time t, days');
    plt.ylabel('d(t)');

    ax1 = plt.subplot(2,1,2);
    plt.axis([tmin, tmin+Dt*Np, -5000, 15000]);
    plt.plot(t[0:Np,0],perror[0:Np,0],'k-');
    plt.xlabel('time t, days');
    plt.ylabel('e(t)');
    plt.show();

In [4]:

    # edapy_07_06 Script supporting CribSheet 7.2
    # Least squares filter estimation of filter m
    # obeying g*m=d

    # standard set-up of time axis, t
    N=101;    # number of data
    Dt = 1.0; # sampling interval
    t = eda_cvec(np.linspace(0,Dt*(N-1),N));
    tmin = t[0,0];   # start time
    tmax = t[N-1,0]; # end time

    # make some synthetic timeseries g
    # draw from an exponential distribution so its all positive
    g = np.random.exponential(1,(N,1));
    # edafilterg is a global variable defined in GLSFilterMul()
    edafilterg = g;

    # exemplary true filter is decaying exponential
    M=10;     # length of filter
    c=5.0*Dt; # decay constant
    mtrue = eda_cvec(0.3*np.exp(-t[0:M,0]/c));

    # the true data is g covolved with m; then truncate to length N
    dtrue = eda_cvec(np.convolve(g.ravel(),mtrue.ravel()));
    dtrue = eda_cvec(dtrue[0:N,0:1]);

    # the observed data is the true data plus noise
    # the noise is Normally distributed with zero mean and variance sigmad^2
    sigmad = 0.1;  # standard deviation of noise
    noise = np.random.normal(0,sigmad,(N,1)); # random noise
    d = dtrue + noise; # data = signal + noise

    # now invert for the filter using the observed data and the time series g

    # prior infoirmation of smoothness, with no end conditions
    eda_e = 0.001;         # 1 damping factor = sigma_h / sigma_d
    # make Toeplitz matrix of second derivatives
    K = M-2;
    D = np.ones((1,M))/(Dt**2);
    DIAGS = np.concatenate((D,-2.0*D,D),axis=0);
    edaHsparse = sparse.spdiags(DIAGS, [0, 1, 2], K, M);
    h = np.zeros((M-2,1)); # prior info that second derivative is zero,

    # set up right hand side, F'f = GT qobs + e HT h
    FTf = GLSFilterMulRHS(edafilterg,d,eda_e,edaHsparse,h,M);

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

    # solve for estimated model parameters using congugate gradient algrorithm
    # returned value 1 is a tuple whose first element is the solution vector
    q=las.cg(LO,FTf,tol=1e-12, maxiter=(3*(N+M)+100));
    mest = eda_cvec(q[0]);

    # predicted data by colvolving mest and g, then truncat to length N
    dpre = eda_cvec(np.convolve(g.ravel(),mest.ravel()));
    dpre = eda_cvec(dpre[0:N]);

    # compute the error
    e = d-dpre;             # individual errors
    E = np.matmul(e.T, e);  # total error

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

    ax1 = plt.subplot(2,1,1);
    plt.axis([tmin, tmax, -5, 5] );
    plt.plot(t,g,'k-');
    plt.xlabel('t');
    plt.ylabel('g(t)');

    ax1 = plt.subplot(2,1,2);
    plt.axis([tmin, tmax, -5, 5] );
    plt.plot(t,d,'k-',lw=1);
    plt.plot(t,dpre,'k:',lw=3);
    plt.xlabel('t');
    plt.ylabel('d(t)');
    plt.show();

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

    ax1 = plt.subplot(1,1,1);
    plt.axis([t[0,0], t[M-1,0], 0, 0.5] );
    plt.plot( t[0:M], mtrue.ravel(),'k-');
    plt.plot( t[0:M], mest.ravel(), 'k-', color = (0.8,0.8,0.8));
    plt.xlabel('t');
    plt.ylabel('m(t)');
    plt.show();

In [130]:

    # edapy_07_07: inverse filter via z-transforms

    # use z-transforms to:
    #    check if a filter is minumum phase
    #    represent it as a cascade of length-2 filters
    #    compute the inverse filter

    # define a filter
    g = eda_cvec( [6, 5.5, 5.3, 5, 4, 3, 2] );
    N,K = np.shape(g);
    gN = g[N-1,0];

    # find roots of g
    pol = (np.flipud(g)).ravel(); # pol is a list of coefficients
    r = np.roots(pol);  # r is a list of roots

    # check size of roots
    badroots=0;
    for j in range(N-1):
        if( abs(r[j]) <= 1 ):
            badroots=badroots+1;

    printstr = "%d roots not outside the unit circle" % (badroots);
    print(printstr);

    # rebuild filter from its roots
    gp = np.zeros((1,1), dtype=np.complex);
    gp[0,0]=gN;
    for j in range(N-1):
        temp = np.convolve(gp.ravel(),[-r[j],1]);
        K = len(temp);
        gp = np.zeros((K,1), dtype=np.complex);
        gp[:,0] = temp;

    # there should be no imiginary part
    gp = np.real(gp);

    # check that it worked by computing error
    E1 = np.matmul( (g-gp).T , (g-gp) );
    printstr = "error of filter rebuilt from its roots: %f" % (E1);
    print(printstr);

    # construct inverse filter, gi, of length Ni
    Ni = 50;
    gi = np.zeros((Ni,1),dtype=np.complex);
    gi[0,0]=1/gN;
    # filter cascade, one filter per loop
    for j in range(N-1):
        # construct inverse filter of a length-2 filter
        temp = np.zeros((Ni,1),dtype=np.complex);
        temp[0,0] = 1;
        for k  in range(1,Ni):
            temp[k,0]=temp[k-1,0]/r[j];
        temp = -temp/r[j];
        temp2 = eda_cvec(np.convolve(gi.ravel(),temp.ravel()));
        gi[0:Ni,0:1]=temp2[0:Ni,0:1];
    # delete imaginary part (which should be zero)
    gi = np.real(gi);

    # test quality of inverse filter; truncate to length Ni
    gig=eda_cvec(np.convolve(g.ravel(),gi.ravel()));
    gig=eda_cvec( gig[0:Ni,0:1] );

    # impulse function
    d = np.zeros((Ni,1));
    d[0,0]=1.0;
    E2 = np.matmul( (gig-d).T, (gig-d) );  # error
    printstr = "error of inverse filter: %f" % (E2);
    print(printstr);

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

    # plot filter
    ax1 = plt.subplot(3,1,1);
    plt.axis([0, Ni, np.min(g), np.max(g)] );
    gx = np.concatenate( (g, np.zeros((Ni-N,1))), axis=0);
    plt.plot( np.linspace(0,Ni-1,Ni), gx.ravel(),'k-');
    plt.xlabel('element j');
    plt.ylabel('g(j)');

    ax1 = plt.subplot(3,1,2);
    plt.axis([0, Ni, np.min(gi), np.max(gi)] );
    plt.plot( np.linspace(0,Ni-1,Ni), gi.ravel(),'k-');
    plt.xlabel('element j');
    plt.ylabel('ginv(j)');

    ax1 = plt.subplot(3,1,3);
    plt.axis([0, Ni, -1.2, 1.2] );
    plt.plot( np.linspace(0,Ni-1,Ni), gig.ravel(),'k-');
    plt.xlabel('element j');
    plt.ylabel('[g*ginv](j)');

    0 roots not outside the unit circle
    error of filter rebuilt from its roots: 0.000000
    error of inverse filter: 0.000000

In [131]:

    # edapy_07_08, recursive smoothing filter, by for-loops

    # example of recursive smoothing filter
    # this version used for-loops, for demonstration purposes
    # the method in edapy_07_09 is prefered

    # standard set-up of time axis, t
    N=101;   # number of data
    Dt=1.0;  # sampling interval
    t = eda_cvec( Dt*np.linspace(0,N-1,N) );  # time vector

    # create two synthetic datasets
    # h1, a unit spike
    # h2, random noise with zero mean and unit variance
    h1=np.zeros((N,1)); # h1 is a time series with a spike in its middle
    No2 = floor(N/2);   # middle at N/2
    h1[No2,0]=1.0;      # spike has unit amplitude
    h2=np.random.normal(0.0,1.0,(N,1));  # h2 is a time series of random noise

    # apply recursive filter to h1 to yield q1
    q1=np.zeros((N,1));
    q1[0,0]=0.5*h1[0,0];
    for j in range(2,N):
        q1[j,0]=0.5*(h1[j,0]+q1[j-1,0]);

    # apply recursive filter to h2 to yield q2
    q2=np.zeros((N,1));
    q2[0,0]=0.5*h2[0,0];
    for j in range(2,N):
        q2[j,0]=0.5*(h2[j,0]+q2[j-1,0]);

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

    ax1 = plt.subplot(2,1,1);
    plt.axis( [t[0,0], t[N-1,0], -1.1, 1.1] );
    plt.plot( t, h1, 'k.');
    plt.plot( t, q1, 'k-', color=(0.7,0.7,0.7));
    plt.xlabel('time t');
    plt.ylabel('h1(t) and q1(t)');

    ax1 = plt.subplot(2,1,2);
    plt.axis([t[0,0], t[N-1,0], -5, 5]);
    plt.plot( t, h2, 'k.');
    plt.plot( t, q2, 'k-', color=(0.7,0.7,0.7));
    plt.xlabel('time t');
    plt.ylabel('h2(t) and q2(t)');
    plt.show();

    displaystr = "sum of elements of the unit spike: %f" % (np.sum(h1));
    print(displaystr);
    displaystr = "sum of elements of the filtered spike: %f" % (np.sum(q1));
    print(displaystr);

    sum of elements of the unit spike: 1.000000
    sum of elements of the filtered spike: 1.000000

In [132]:

    # edapy_07_09: recursive smoothing filter, by sg.lfilter()

    # example of recursive smoothing filter
    # this version used scipy.signal.lfilter()

    # standard set-up of time axis, t
    N=101;   # number of data
    Dt=1.0;  # sampling interval
    t = eda_cvec(np.linspace(0,Dt*(N-1),N)); # time vector

    # create two synthetic datasets
    # h1, a unit spike
    # h2, random noise with zero mean and unit variance
    h1=np.zeros((N,1)); # h1 is a time series with a spike in its middle
    No2 = floor(N/2);   # middle position
    h1[No2,0]=1.0;      # unit spike
    h2=np.random.normal(0.0,1.0,(N,1)); # h2 is a random time series

    Nu = 2;                      # length of filter u
    Nv = 2;                      # length of filter v
    u = eda_cvec( [0.5, 0.0] );  # filter u
    v = eda_cvec([1.0, -0.5] );  # filter v
    # filter to h1 to q1
    q1 = eda_cvec(sg.lfilter(u.ravel(),v.ravel(),h1.ravel()));
    q1 = eda_cvec( q1[0:N,0:1] ); # truncate to length N
    # apply filter to h2 to q2
    q2 = eda_cvec(sg.lfilter(u.ravel(),v.ravel(),h2.ravel()));
    q2 = eda_cvec( q2[0:N,0:1] ); # truncate to length N

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

    ax1 = plt.subplot(2,1,1);
    plt.axis( [t[0,0], t[N-1,0], -1.1, 1.1] );
    plt.plot( t, h1, 'k.');
    plt.plot( t, q1, 'k-', color=(0.7,0.7,0.7));
    plt.xlabel('time t');
    plt.ylabel('h1(t) and q1(t)');

    ax1 = plt.subplot(2,1,2);
    plt.axis([t[0,0], t[N-1,0], -5, 5]);
    plt.plot( t, h2, 'k.');
    plt.plot( t, q2, 'k-', color=(0.7,0.7,0.7));
    plt.xlabel('time t');
    plt.ylabel('h2(t) and q2(t)');
    plt.show();

    displaystr = "sum of elements of the unit spike: %f" % (np.sum(h1));
    print(displaystr);
    displaystr = "sum of elements of the filtered spike: %f" % (np.sum(q1));
    print(displaystr);

    sum of elements of the unit spike: 1.000000
    sum of elements of the filtered spike: 1.000000

In [ ]: